#pragma once
#define TOLERANCE 0.05

__device__ 
float oneNorm(const float A[3][3])
{
  float norm = 0.0;
  #pragma unroll 3
  for (int i=0; i<3; i++) 
  {
    float columnAbsSum = fabsf(A[0][i]) + fabsf(A[1][i]) + fabsf(A[2][i]);
    if (columnAbsSum > norm) 
      norm = columnAbsSum;
  }
  return norm;
}

__device__
float infNorm(const float A[3][3])
{
  float norm = 0.0;
  #pragma unroll 3
  for (int i=0; i<3; i++) 
  {
    float rowSum = fabsf(A[i][0]) + fabsf(A[i][1]) + fabsf(A[i][2]);
    if (rowSum > norm) 
      norm = rowSum;
  }
  return norm;
}

// cross product: c = a x b
__device__
void crossProduct(const float* a, const float* b, float* c)
{
  c[0] = a[1] * b[2] - a[2] * b[1];
  c[1] = a[2] * b[0] - a[0] * b[2];
  c[2] = a[0] * b[1] - a[1] * b[0];
}

// Input: M (3x3 mtx)
// Output: Q (3x3 rotation mtx), S (3x3 symmetric mtx)
__device__
void gpuComputePolarDecomposition(float Mk[3][3])
{
  //float Mk[3][3];
  float Ek[3][3];
  float det, M_oneNorm, M_infNorm, E_oneNorm;

  // Mk = M^T
  //for(int i=0; i<3; i++)
  //  for(int j=0; j<3; j++)
  //    Mk[i][j] = M[j][i];

  M_oneNorm = oneNorm(Mk); 
  M_infNorm = infNorm(Mk);

  do 
  {
    float MadjTk[3][3];
 
    // row 2 x row 3
    crossProduct(Mk[1], Mk[2], MadjTk[0]); 
    // row 3 x row 1
    crossProduct(Mk[2], Mk[0], MadjTk[1]);
    // row 1 x row 2
    crossProduct(Mk[0], Mk[1], MadjTk[2]);

    det = Mk[0][0] * MadjTk[0][0] + Mk[0][1] * MadjTk[0][1] + Mk[0][2] * MadjTk[0][2];
    if (det == 0.0) 
    {
      //printf("Warning (polarDecomposition) : zero determinant encountered.\n");
      break;
    }

    float MadjT_one = oneNorm(MadjTk); 
    float MadjT_inf = infNorm(MadjTk);

    float gamma = sqrtf(sqrtf((MadjT_one * MadjT_inf) / (M_oneNorm * M_infNorm)) / fabsf(det));
    float g1 = gamma * 0.5;
    float g2 = 0.5 / (gamma * det);

	#pragma unroll 3
    for(int i=0; i<3; i++)
		#pragma unroll 3
		for(int j=0; j<3; j++)
		{
			Ek[i][j] = Mk[i][j];
			Mk[i][j] = g1 * Mk[i][j] + g2 * MadjTk[i][j];
			Ek[i][j] -= Mk[i][j];
		}

    E_oneNorm = oneNorm(Ek);
    M_oneNorm = oneNorm(Mk);  
    M_infNorm = infNorm(Mk);
  }
  while ( E_oneNorm > M_oneNorm * TOLERANCE );

  // Q = Mk^T 
  //for(int i=0; i<3; i++)
  //  for(int j=0; j<3; j++)
  //    Q[i][j] = Mk[j][i];

  //return (det);
}

	.file	"tmpxft_002a8d59_00000000-6_GPUPolarDecompose.cudafe1.cpp"
	.text
#APP

#NO_APP
	.type	_ZL26__cudaUnregisterBinaryUtilv, @function
_ZL26__cudaUnregisterBinaryUtilv:
.LFB2014:
	.cfi_startproc
	endbr64
	movq	_ZL20__cudaFatCubinHandle(%rip), %rdi
	jmp	__cudaUnregisterFatBinary@PLT
	.cfi_endproc
.LFE2014:
	.size	_ZL26__cudaUnregisterBinaryUtilv, .-_ZL26__cudaUnregisterBinaryUtilv
	.globl	_Z7oneNormPA3_Kf
	.type	_Z7oneNormPA3_Kf, @function
_Z7oneNormPA3_Kf:
.LFB2008:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	popq	%rax
	.cfi_def_cfa_offset 8
	subq	$24, %rsp
	.cfi_def_cfa_offset 32
	movl	$1, 12(%rsp)
	movl	12(%rsp), %edi
	call	exit@PLT
	.cfi_endproc
.LFE2008:
	.size	_Z7oneNormPA3_Kf, .-_Z7oneNormPA3_Kf
	.globl	_Z7infNormPA3_Kf
	.type	_Z7infNormPA3_Kf, @function
_Z7infNormPA3_Kf:
.LFB2153:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	popq	%rax
	.cfi_def_cfa_offset 8
	pushq	%rax
	.cfi_def_cfa_offset 16
	call	_Z7oneNormPA3_Kf
	.cfi_endproc
.LFE2153:
	.size	_Z7infNormPA3_Kf, .-_Z7infNormPA3_Kf
	.globl	_Z12crossProductPKfS0_Pf
	.type	_Z12crossProductPKfS0_Pf, @function
_Z12crossProductPKfS0_Pf:
.LFB2010:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	popq	%rax
	.cfi_def_cfa_offset 8
	subq	$24, %rsp
	.cfi_def_cfa_offset 32
	movl	$1, 12(%rsp)
	movl	12(%rsp), %edi
	call	exit@PLT
	.cfi_endproc
.LFE2010:
	.size	_Z12crossProductPKfS0_Pf, .-_Z12crossProductPKfS0_Pf
	.globl	_Z28gpuComputePolarDecompositionPA3_f
	.type	_Z28gpuComputePolarDecompositionPA3_f, @function
_Z28gpuComputePolarDecompositionPA3_f:
.LFB2011:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	popq	%rax
	.cfi_def_cfa_offset 8
	subq	$24, %rsp
	.cfi_def_cfa_offset 32
	movl	$1, 12(%rsp)
	movl	12(%rsp), %edi
	call	exit@PLT
	.cfi_endproc
.LFE2011:
	.size	_Z28gpuComputePolarDecompositionPA3_f, .-_Z28gpuComputePolarDecompositionPA3_f
	.section	.text.startup,"ax",@progbits
	.type	_ZL24__sti____cudaRegisterAllv, @function
_ZL24__sti____cudaRegisterAllv:
.LFB2037:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	leaq	_ZL15__fatDeviceText(%rip), %rdi
	call	__cudaRegisterFatBinary@PLT
	movq	%rax, %rdi
	movq	%rax, _ZL20__cudaFatCubinHandle(%rip)
	call	__cudaRegisterFatBinaryEnd@PLT
	leaq	_ZL26__cudaUnregisterBinaryUtilv(%rip), %rdi
	popq	%rdx
	.cfi_def_cfa_offset 8
	jmp	atexit@PLT
	.cfi_endproc
.LFE2037:
	.size	_ZL24__sti____cudaRegisterAllv, .-_ZL24__sti____cudaRegisterAllv
	.section	.init_array,"aw"
	.align 8
	.quad	_ZL24__sti____cudaRegisterAllv
	.section	.nvFatBinSegment,"aw"
	.align 8
	.type	_ZL15__fatDeviceText, @object
	.size	_ZL15__fatDeviceText, 24
_ZL15__fatDeviceText:
	.long	1180844977
	.long	1
	.quad	fatbinData
	.quad	0
	.local	_ZL20__cudaFatCubinHandle
	.comm	_ZL20__cudaFatCubinHandle,8,8
	.ident	"GCC: (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0"
	.section	.note.GNU-stack,"",@progbits
	.section	.note.gnu.property,"a"
	.align 8
	.long	1f - 0f
	.long	4f - 1f
	.long	5
0:
	.string	"GNU"
1:
	.align 8
	.long	0xc0000002
	.long	3f - 2f
2:
	.long	0x3
3:
	.align 8
4:

#pragma once
#include <hip/hip_runtime.h>

#define TOLERANCE 0.05

__device__ 
float oneNorm(const float A[3][3])
{
  float norm = 0.0;
  #pragma unroll 3
  for (int i=0; i<3; i++) 
  {
    float columnAbsSum = fabsf(A[0][i]) + fabsf(A[1][i]) + fabsf(A[2][i]);
    if (columnAbsSum > norm) 
      norm = columnAbsSum;
  }
  return norm;
}

__device__
float infNorm(const float A[3][3])
{
  float norm = 0.0;
  #pragma unroll 3
  for (int i=0; i<3; i++) 
  {
    float rowSum = fabsf(A[i][0]) + fabsf(A[i][1]) + fabsf(A[i][2]);
    if (rowSum > norm) 
      norm = rowSum;
  }
  return norm;
}

// cross product: c = a x b
__device__
void crossProduct(const float* a, const float* b, float* c)
{
  c[0] = a[1] * b[2] - a[2] * b[1];
  c[1] = a[2] * b[0] - a[0] * b[2];
  c[2] = a[0] * b[1] - a[1] * b[0];
}

// Input: M (3x3 mtx)
// Output: Q (3x3 rotation mtx), S (3x3 symmetric mtx)
__device__
void gpuComputePolarDecomposition(float Mk[3][3])
{
  //float Mk[3][3];
  float Ek[3][3];
  float det, M_oneNorm, M_infNorm, E_oneNorm;

  // Mk = M^T
  //for(int i=0; i<3; i++)
  //  for(int j=0; j<3; j++)
  //    Mk[i][j] = M[j][i];

  M_oneNorm = oneNorm(Mk); 
  M_infNorm = infNorm(Mk);

  do 
  {
    float MadjTk[3][3];
 
    // row 2 x row 3
    crossProduct(Mk[1], Mk[2], MadjTk[0]); 
    // row 3 x row 1
    crossProduct(Mk[2], Mk[0], MadjTk[1]);
    // row 1 x row 2
    crossProduct(Mk[0], Mk[1], MadjTk[2]);

    det = Mk[0][0] * MadjTk[0][0] + Mk[0][1] * MadjTk[0][1] + Mk[0][2] * MadjTk[0][2];
    if (det == 0.0) 
    {
      //printf("Warning (polarDecomposition) : zero determinant encountered.\n");
      break;
    }

    float MadjT_one = oneNorm(MadjTk); 
    float MadjT_inf = infNorm(MadjTk);

    float gamma = sqrtf(sqrtf((MadjT_one * MadjT_inf) / (M_oneNorm * M_infNorm)) / fabsf(det));
    float g1 = gamma * 0.5;
    float g2 = 0.5 / (gamma * det);

	#pragma unroll 3
    for(int i=0; i<3; i++)
		#pragma unroll 3
		for(int j=0; j<3; j++)
		{
			Ek[i][j] = Mk[i][j];
			Mk[i][j] = g1 * Mk[i][j] + g2 * MadjTk[i][j];
			Ek[i][j] -= Mk[i][j];
		}

    E_oneNorm = oneNorm(Ek);
    M_oneNorm = oneNorm(Mk);  
    M_infNorm = infNorm(Mk);
  }
  while ( E_oneNorm > M_oneNorm * TOLERANCE );

  // Q = Mk^T 
  //for(int i=0; i<3; i++)
  //  for(int j=0; j<3; j++)
  //    Q[i][j] = Mk[j][i];

  //return (det);
}

	.text
	.file	"GPUPolarDecompose.hip"
	.type	__hip_cuid_,@object             # @__hip_cuid_
	.bss
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               # 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_

	.text
	.p2alignl 7, 3214868480
	.fill 96, 4, 3214868480
	.type	__hip_cuid_,@object             ; @__hip_cuid_
	.section	.bss,"aw",@nobits
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               ; 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_
	.amdgcn_target "amdgcn-amd-amdhsa--gfx1100"
	.amdgpu_metadata
---
amdhsa.kernels:  []
amdhsa.target:   amdgcn-amd-amdhsa--gfx1100
amdhsa.version:
  - 1
  - 2
...

	.end_amdgpu_metadata
