//$Id: mycudamath.cu,v 1.2 2010/05/15 16:24:57 afs Exp $

__device__ void choldcU(float* a, int* pn, float* y)  
{
	int n =  *pn;
	int i,j,k;
	float sum;
	unsigned int ij, ik, jk, ii, ji, nk;
	//
	for (i=0;i<n;i++) 
		for (j=0;j<n;j++) {
			ij=i+n*j;
			y[ij] = a[ij]; 
		}
	//
	for (i=0;i<n;i++) {
		for (j=i;j<n;j++) {
			ij=i+n*j;
			for (sum=y[ij],k=i-1;k>=0;k--) {
				nk=n*k;
		 		ik=i+nk; jk=j+nk;
				sum -= y[ik]*y[jk];
			}
			ii=i+n*i;
			if (i == j) {
				if (sum <= 0.0f)  //	A, with rounding errors, is not positive-deﬁnite.
					return;
					// throw("\n*** Cholesky failed\n")
				y[ii]=sqrt(sum);
			}
			else {
				ji=j+n*i;
				y[ji]=sum/y[ii];
			}
		}
	}
	for (i=0;i<n;i++) // output as Upper
		for (j=0;j<i;j++) {
			ji=j+n*i;
			ij=i+n*j;
			y[ji] = y[ij];
			y[ij] = 0.0f;
		}

/*
	int l;
	for (i=0;i<n;i++){ // by columns
		for (j=0;j<n;j++) {
			l=i+n*j;
			printf("%e ", y[l]);
		}
		printf("\n");
	}
	printf("----------\n");
*/
}

//----------------------------------------------------------------------
// Input as L^T
__device__ void elsolveU(float* el, int* b, float* y, int* pn)
{
// Solve L y = b, where L is the lower triangular matrix in the stored Cholesky decomposi-
// tion. b[0..n-1] is input as the right-hand side vector. The solution vector is returned in
// y[0..n-1].
	int n = *pn;
	int i,k;
	int ki, ii;
	float sum;
	// if (b.size() != n || y.size() != n) throw("bad lengths");
	for (i=n-1;i>=0;i--) { // Solve L^T x = y.
		for (sum=b[i],k=i+1;k<n;k++){
			ki=k*n+i;
			sum -= el[ki]*y[k];
		}
		ii=i*n+i;
		y[i]=sum/el[ii];
	}
}


//----------------------------------------------------------------------

__device__ void mdgbacksolve(float* a, int* pn, float* y)
{
	int n =  *pn;
	// __shared__ float el[25];
	float el[25];
  choldcU(a,pn,el);
	//
	// __shared__ int b[5]; 
	int b[5]; 

	for(int i=0; i < n; i++) b[i]= 0;
	for(int i=0; i < n; i++) {
		b[i]=1;
		elsolveU(el, b , y+i*n, pn); 
		b[i]=0;
	}

}

//----------------------------------------------------------------------

__device__ void mtcrossp(float *a, float *b, float *c, int *pn)
{
	int n = *pn;
	float sum;
	int ik, kj, ij, nk;
	//
  for (int i=0; i<n; i++)
		for (int j=0; j<n; j++) {
		  sum = 0.0f;
		  for (int k=0; k<n; k++) {
				nk=n*k;
				ik=i+nk;
				kj=nk+j;
				sum += a[ik] * b[kj];
			}
			ij=i+n*j;
			c[ij] = sum;
		}
}

//----------------------------------------------------------------------
// product of matrix (mxn) x vector (n) 

__device__ void mvprod(float *a, float *b, float *c, int *pm, int *pn)
{
	int m = *pm;
	int n = *pn;
	float sum;
	int ik;
	//
  for (int i=0; i<m; i++) {
		sum = 0.0f;
		for (int k=0; k<n; k++) {
			ik=i+m*k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}

//----------------------------------------------------------------------
/*
// product of matrix (mxn) x vector (m) 

__device__ void mvprodm(float *a, float *b, float *c, int *pm)
{
	int m = *pm;
	float sum;
	int ik;
	//
  for (int i=0; i<m; i++) {
		sum = 0.0f;
		for (int k=0; k<m; k++) {
			ik=i+m*k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}

*/

//----------------------------------------------------------------------
//
// product of vectors 1xm X mx1
// 
__device__ void vprod(float *a, float *b, float *c, int *pm)
{
	int m = *pm;
	float sum;
	//
	sum = 0.0f;
  for (int i=0; i<m; i++) {
		sum += a[i] * b[i];
	}
  *c = sum;
}

//----------------------------------------------------------------------


// t(matrix)*vector - with R-order (by column)
__device__ void mvtcrossp(float *a, float *b, float *c, int *pm, int* pn)
{
	int m = *pm;
	int n = *pn;
	float sum;
	int ik, mi;
	//
  for (int i=0; i<n; i++) {
		sum = 0.0f;
		mi=m*i;
		for (int k=0; k<m; k++) {
			ik=mi+k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}


//----------------------------------------------------------------------

	.file	"tmpxft_00246587_00000000-6_mycudamath.cudafe1.cpp"
	.text
#APP

#NO_APP
	.type	_ZL26__cudaUnregisterBinaryUtilv, @function
_ZL26__cudaUnregisterBinaryUtilv:
.LFB2017:
	.cfi_startproc
	endbr64
	movq	_ZL20__cudaFatCubinHandle(%rip), %rdi
	jmp	__cudaUnregisterFatBinary@PLT
	.cfi_endproc
.LFE2017:
	.size	_ZL26__cudaUnregisterBinaryUtilv, .-_ZL26__cudaUnregisterBinaryUtilv
	.globl	_Z7choldcUPfPiS_
	.type	_Z7choldcUPfPiS_, @function
_Z7choldcUPfPiS_:
.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	_Z7choldcUPfPiS_, .-_Z7choldcUPfPiS_
	.globl	_Z8elsolveUPfPiS_S0_
	.type	_Z8elsolveUPfPiS_S0_, @function
_Z8elsolveUPfPiS_S0_:
.LFB2009:
	.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
.LFE2009:
	.size	_Z8elsolveUPfPiS_S0_, .-_Z8elsolveUPfPiS_S0_
	.globl	_Z12mdgbacksolvePfPiS_
	.type	_Z12mdgbacksolvePfPiS_, @function
_Z12mdgbacksolvePfPiS_:
.LFB2156:
	.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
.LFE2156:
	.size	_Z12mdgbacksolvePfPiS_, .-_Z12mdgbacksolvePfPiS_
	.globl	_Z8mtcrosspPfS_S_Pi
	.type	_Z8mtcrosspPfS_S_Pi, @function
_Z8mtcrosspPfS_S_Pi:
.LFB2158:
	.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
.LFE2158:
	.size	_Z8mtcrosspPfS_S_Pi, .-_Z8mtcrosspPfS_S_Pi
	.globl	_Z6mvprodPfS_S_PiS0_
	.type	_Z6mvprodPfS_S_PiS0_, @function
_Z6mvprodPfS_S_PiS0_:
.LFB2012:
	.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
.LFE2012:
	.size	_Z6mvprodPfS_S_PiS0_, .-_Z6mvprodPfS_S_PiS0_
	.globl	_Z5vprodPfS_S_Pi
	.type	_Z5vprodPfS_S_Pi, @function
_Z5vprodPfS_S_Pi:
.LFB2160:
	.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
.LFE2160:
	.size	_Z5vprodPfS_S_Pi, .-_Z5vprodPfS_S_Pi
	.globl	_Z9mvtcrosspPfS_S_PiS0_
	.type	_Z9mvtcrosspPfS_S_PiS0_, @function
_Z9mvtcrosspPfS_S_PiS0_:
.LFB2162:
	.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
.LFE2162:
	.size	_Z9mvtcrosspPfS_S_PiS0_, .-_Z9mvtcrosspPfS_S_PiS0_
	.section	.text.startup,"ax",@progbits
	.type	_ZL24__sti____cudaRegisterAllv, @function
_ZL24__sti____cudaRegisterAllv:
.LFB2040:
	.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
.LFE2040:
	.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:


#include <hip/hip_runtime.h>
//$Id: mycudamath.cu,v 1.2 2010/05/15 16:24:57 afs Exp $

__device__ void choldcU(float* a, int* pn, float* y)  
{
	int n =  *pn;
	int i,j,k;
	float sum;
	unsigned int ij, ik, jk, ii, ji, nk;
	//
	for (i=0;i<n;i++) 
		for (j=0;j<n;j++) {
			ij=i+n*j;
			y[ij] = a[ij]; 
		}
	//
	for (i=0;i<n;i++) {
		for (j=i;j<n;j++) {
			ij=i+n*j;
			for (sum=y[ij],k=i-1;k>=0;k--) {
				nk=n*k;
		 		ik=i+nk; jk=j+nk;
				sum -= y[ik]*y[jk];
			}
			ii=i+n*i;
			if (i == j) {
				if (sum <= 0.0f)  //	A, with rounding errors, is not positive-deﬁnite.
					return;
					// throw("\n*** Cholesky failed\n")
				y[ii]=sqrt(sum);
			}
			else {
				ji=j+n*i;
				y[ji]=sum/y[ii];
			}
		}
	}
	for (i=0;i<n;i++) // output as Upper
		for (j=0;j<i;j++) {
			ji=j+n*i;
			ij=i+n*j;
			y[ji] = y[ij];
			y[ij] = 0.0f;
		}

/*
	int l;
	for (i=0;i<n;i++){ // by columns
		for (j=0;j<n;j++) {
			l=i+n*j;
			printf("%e ", y[l]);
		}
		printf("\n");
	}
	printf("----------\n");
*/
}

//----------------------------------------------------------------------
// Input as L^T
__device__ void elsolveU(float* el, int* b, float* y, int* pn)
{
// Solve L y = b, where L is the lower triangular matrix in the stored Cholesky decomposi-
// tion. b[0..n-1] is input as the right-hand side vector. The solution vector is returned in
// y[0..n-1].
	int n = *pn;
	int i,k;
	int ki, ii;
	float sum;
	// if (b.size() != n || y.size() != n) throw("bad lengths");
	for (i=n-1;i>=0;i--) { // Solve L^T x = y.
		for (sum=b[i],k=i+1;k<n;k++){
			ki=k*n+i;
			sum -= el[ki]*y[k];
		}
		ii=i*n+i;
		y[i]=sum/el[ii];
	}
}


//----------------------------------------------------------------------

__device__ void mdgbacksolve(float* a, int* pn, float* y)
{
	int n =  *pn;
	// __shared__ float el[25];
	float el[25];
  choldcU(a,pn,el);
	//
	// __shared__ int b[5]; 
	int b[5]; 

	for(int i=0; i < n; i++) b[i]= 0;
	for(int i=0; i < n; i++) {
		b[i]=1;
		elsolveU(el, b , y+i*n, pn); 
		b[i]=0;
	}

}

//----------------------------------------------------------------------

__device__ void mtcrossp(float *a, float *b, float *c, int *pn)
{
	int n = *pn;
	float sum;
	int ik, kj, ij, nk;
	//
  for (int i=0; i<n; i++)
		for (int j=0; j<n; j++) {
		  sum = 0.0f;
		  for (int k=0; k<n; k++) {
				nk=n*k;
				ik=i+nk;
				kj=nk+j;
				sum += a[ik] * b[kj];
			}
			ij=i+n*j;
			c[ij] = sum;
		}
}

//----------------------------------------------------------------------
// product of matrix (mxn) x vector (n) 

__device__ void mvprod(float *a, float *b, float *c, int *pm, int *pn)
{
	int m = *pm;
	int n = *pn;
	float sum;
	int ik;
	//
  for (int i=0; i<m; i++) {
		sum = 0.0f;
		for (int k=0; k<n; k++) {
			ik=i+m*k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}

//----------------------------------------------------------------------
/*
// product of matrix (mxn) x vector (m) 

__device__ void mvprodm(float *a, float *b, float *c, int *pm)
{
	int m = *pm;
	float sum;
	int ik;
	//
  for (int i=0; i<m; i++) {
		sum = 0.0f;
		for (int k=0; k<m; k++) {
			ik=i+m*k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}

*/

//----------------------------------------------------------------------
//
// product of vectors 1xm X mx1
// 
__device__ void vprod(float *a, float *b, float *c, int *pm)
{
	int m = *pm;
	float sum;
	//
	sum = 0.0f;
  for (int i=0; i<m; i++) {
		sum += a[i] * b[i];
	}
  *c = sum;
}

//----------------------------------------------------------------------


// t(matrix)*vector - with R-order (by column)
__device__ void mvtcrossp(float *a, float *b, float *c, int *pm, int* pn)
{
	int m = *pm;
	int n = *pn;
	float sum;
	int ik, mi;
	//
  for (int i=0; i<n; i++) {
		sum = 0.0f;
		mi=m*i;
		for (int k=0; k<m; k++) {
			ik=mi+k;
			sum += a[ik] * b[k];
		}
		c[i] = sum;
	}
}


//----------------------------------------------------------------------

	.text
	.file	"mycudamath.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
