/*
   This is an upper-triangularization operation on a 'nearly upper triangular' matrix.
   In order to place rME and rMI in registers for optimal performance, the number of
   non-zero values below the subdiagonal must be known at compile time.  The only
   known way to do this is templating the function. See line 31.

   This is a simplied version that just solves a single matrix of ngr*ngr size.
   The kernel should be launched with # of threads = ngr.
*/
   template<unsigned char NSUBDIAG>
__global__ void kernelDecomposeRegister ( int ngr,
                                          const int * __restrict__ ODL, // array size = ngr
                                          double * __restrict__ MatrixArrayRM, // RM = in row-major order
                                          double *rhs
                                        )

{
   extern __shared__ double sMEM[];

   // require 19 * ngr SMEM storage  (19 kB for ngr = 128)
   // ... one for current row
   // ... one for column[row] of descendants
   // ... one for nAL for each row

   double * __restrict__ sNormal        = sMEM;
   double * __restrict__ sCurrentRow    = sMEM+1;
   double * __restrict__ sCurrentColumn = sCurrentRow + ngr;  // this really only needs to be size 16 (blockDim.x / ngr)
   double * __restrict__ sRhs           = sCurrentColumn + ngr;
   int    * __restrict__ sNAL           = (int *) (sRhs + ngr);

   // require a compile-time known number of descendants to be placed in registers
   double       rME[NSUBDIAG];                // [r]egister [M]atrix [E]lement
   unsigned int rMI[NSUBDIAG];                // [r]egister [M]atrix [I]ndex

   unsigned int tid = threadIdx.x;
   unsigned int tidongr = tid/ngr;
   unsigned int tidmngr = tid%ngr;

   // load first blockDim.x/ngr lines of dense matrix                 // this assumes ngr >= 16   // is this comment correct??   Seems out of date.
   for ( int ireg=0; ireg<NSUBDIAG; ireg++ ){                    // rows per thread is explicit
      rME[ireg] = MatrixArrayRM[tid + ireg*blockDim.x];
      rMI[ireg] = tidongr + ireg;
   }

   // the zeroth column for each rows places column[0] in SMEM        // place the first column into shared memory
   if ( tidmngr == 0 ) {                                              // this is what will be used to compute the row operations
      for ( int ireg=0; ireg<NSUBDIAG; ireg++ ) {
         sCurrentColumn[tidongr + ireg] = rME[ireg];
      }
   }

   // compute nAL for each row from ODL and place in SMEM  // nAL = number of elements below a given diagona.
   if ( tid < ngr-1 ) {   // this can be made an else clause from the previous if to avoid divergence.
      sNAL[tid] = ODL[tid+1] - ODL[tid];
   }

   // place the RHS in SMEM
   if ( tid < ngr ) {
      sRhs[tid] = rhs[tid];
   }

   __syncthreads();

   // normalize first line of dense matrix - just first ngr threads
   if ( tid == 0 ) {
      sRhs[0] /= sCurrentColumn[0];
      MatrixArrayRM[0] = 1.0;
   }
   else if ( tid < ngr ) {
      rME[0] /= sCurrentColumn[0];
      // place result in SMEM
      sCurrentRow[tid] = rME[0];
      // write to GMEM
      MatrixArrayRM[tid] = rME[0];
   } 

   __syncthreads();

   for ( int irow=0; irow<ngr-1 ; irow++ ) {
      for ( int ireg=0; ireg<NSUBDIAG; ireg++ ) {
         if ( rMI[ireg] == irow ) {                                       // if this row is the diagonal, 
            rMI[ireg] += NSUBDIAG;
            if ( rMI[ireg] < ngr ) {                                       //   if we haven't progressed past the end of the matrix                
               rME[ireg] = MatrixArrayRM[(tidmngr) + rMI[ireg]*ngr];	 //     then read in the next row.
            }
         }

         if ( rMI[ireg] - irow <= sNAL[irow] ) {                          // check if this row is within the number of descendants of the current diagonal
            // if so, then subtract irow 
            if ( tidmngr >= irow ) {
               rME[ireg] -= sCurrentColumn[rMI[ireg]] * sCurrentRow[tidmngr];
            }
            // update RHS
            if ( tidmngr == 0 ) {
               sRhs[rMI[ireg]] -= sCurrentColumn[rMI[ireg]] * sRhs[irow];
            }
         }

         if ( rMI[ireg] == irow+1 && tidmngr == irow+1 ) {               // if this will be the next eliminated (i.e. next iteration this will be the diagonal row)
            *sNormal = 1.0/rME[ireg];                                     //   then cache diagonal value so the row can be normalized
         }

      }
      __syncthreads();

      for ( int ireg=0; ireg < NSUBDIAG; ireg++ )
      {
         if ( rMI[ireg] == irow+1 ) {                                    // if this thread has row irow+1 cached (i.e. the next row to eliminate)
            // write out
            if ( tidmngr > irow ) {                                       // if the column is in the upper triangular
               // normalize
               rME[ireg] *= *sNormal;                                      // normalize
               sCurrentRow[tidmngr] = rME[ireg];                           // write to shared
               MatrixArrayRM[(tidmngr)+(irow+1)*ngr] = rME[ireg];          // also write to GMEM (as it is the fully upper-triangularized result)
            }
            if ( tidmngr == 0 ) {                                         // use the zeroth thread in the row to normalize the RHS
               sRhs[rMI[ireg]] *= *sNormal;
            }
         }

         if ( tidmngr == irow+1 ) {
            if (rMI[ireg] < ngr) {
               sCurrentColumn[rMI[ireg]] = rME[ireg];
            }
         }

      }

      // syncthreads ------------------------------------------------
      __syncthreads();    

   }

   // output updated rhs
   if ( tid < ngr ) {
      rhs[tid] = sRhs[tid];
   }
}

	.file	"tmpxft_003aea7c_00000000-6_UTMatrixOperation.cudafe1.cpp"
	.text
#APP

#NO_APP
	.type	_ZL26__cudaUnregisterBinaryUtilv, @function
_ZL26__cudaUnregisterBinaryUtilv:
.LFB2012:
	.cfi_startproc
	endbr64
	movq	_ZL20__cudaFatCubinHandle(%rip), %rdi
	jmp	__cudaUnregisterFatBinary@PLT
	.cfi_endproc
.LFE2012:
	.size	_ZL26__cudaUnregisterBinaryUtilv, .-_ZL26__cudaUnregisterBinaryUtilv
	.section	.text.startup,"ax",@progbits
	.type	_ZL24__sti____cudaRegisterAllv, @function
_ZL24__sti____cudaRegisterAllv:
.LFB2035:
	.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
.LFE2035:
	.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>
/*
   This is an upper-triangularization operation on a 'nearly upper triangular' matrix.
   In order to place rME and rMI in registers for optimal performance, the number of
   non-zero values below the subdiagonal must be known at compile time.  The only
   known way to do this is templating the function. See line 31.

   This is a simplied version that just solves a single matrix of ngr*ngr size.
   The kernel should be launched with # of threads = ngr.
*/
   template<unsigned char NSUBDIAG>
__global__ void kernelDecomposeRegister ( int ngr,
                                          const int * __restrict__ ODL, // array size = ngr
                                          double * __restrict__ MatrixArrayRM, // RM = in row-major order
                                          double *rhs
                                        )

{
   extern __shared__ double sMEM[];

   // require 19 * ngr SMEM storage  (19 kB for ngr = 128)
   // ... one for current row
   // ... one for column[row] of descendants
   // ... one for nAL for each row

   double * __restrict__ sNormal        = sMEM;
   double * __restrict__ sCurrentRow    = sMEM+1;
   double * __restrict__ sCurrentColumn = sCurrentRow + ngr;  // this really only needs to be size 16 (blockDim.x / ngr)
   double * __restrict__ sRhs           = sCurrentColumn + ngr;
   int    * __restrict__ sNAL           = (int *) (sRhs + ngr);

   // require a compile-time known number of descendants to be placed in registers
   double       rME[NSUBDIAG];                // [r]egister [M]atrix [E]lement
   unsigned int rMI[NSUBDIAG];                // [r]egister [M]atrix [I]ndex

   unsigned int tid = threadIdx.x;
   unsigned int tidongr = tid/ngr;
   unsigned int tidmngr = tid%ngr;

   // load first blockDim.x/ngr lines of dense matrix                 // this assumes ngr >= 16   // is this comment correct??   Seems out of date.
   for ( int ireg=0; ireg<NSUBDIAG; ireg++ ){                    // rows per thread is explicit
      rME[ireg] = MatrixArrayRM[tid + ireg*blockDim.x];
      rMI[ireg] = tidongr + ireg;
   }

   // the zeroth column for each rows places column[0] in SMEM        // place the first column into shared memory
   if ( tidmngr == 0 ) {                                              // this is what will be used to compute the row operations
      for ( int ireg=0; ireg<NSUBDIAG; ireg++ ) {
         sCurrentColumn[tidongr + ireg] = rME[ireg];
      }
   }

   // compute nAL for each row from ODL and place in SMEM  // nAL = number of elements below a given diagona.
   if ( tid < ngr-1 ) {   // this can be made an else clause from the previous if to avoid divergence.
      sNAL[tid] = ODL[tid+1] - ODL[tid];
   }

   // place the RHS in SMEM
   if ( tid < ngr ) {
      sRhs[tid] = rhs[tid];
   }

   __syncthreads();

   // normalize first line of dense matrix - just first ngr threads
   if ( tid == 0 ) {
      sRhs[0] /= sCurrentColumn[0];
      MatrixArrayRM[0] = 1.0;
   }
   else if ( tid < ngr ) {
      rME[0] /= sCurrentColumn[0];
      // place result in SMEM
      sCurrentRow[tid] = rME[0];
      // write to GMEM
      MatrixArrayRM[tid] = rME[0];
   } 

   __syncthreads();

   for ( int irow=0; irow<ngr-1 ; irow++ ) {
      for ( int ireg=0; ireg<NSUBDIAG; ireg++ ) {
         if ( rMI[ireg] == irow ) {                                       // if this row is the diagonal, 
            rMI[ireg] += NSUBDIAG;
            if ( rMI[ireg] < ngr ) {                                       //   if we haven't progressed past the end of the matrix                
               rME[ireg] = MatrixArrayRM[(tidmngr) + rMI[ireg]*ngr];	 //     then read in the next row.
            }
         }

         if ( rMI[ireg] - irow <= sNAL[irow] ) {                          // check if this row is within the number of descendants of the current diagonal
            // if so, then subtract irow 
            if ( tidmngr >= irow ) {
               rME[ireg] -= sCurrentColumn[rMI[ireg]] * sCurrentRow[tidmngr];
            }
            // update RHS
            if ( tidmngr == 0 ) {
               sRhs[rMI[ireg]] -= sCurrentColumn[rMI[ireg]] * sRhs[irow];
            }
         }

         if ( rMI[ireg] == irow+1 && tidmngr == irow+1 ) {               // if this will be the next eliminated (i.e. next iteration this will be the diagonal row)
            *sNormal = 1.0/rME[ireg];                                     //   then cache diagonal value so the row can be normalized
         }

      }
      __syncthreads();

      for ( int ireg=0; ireg < NSUBDIAG; ireg++ )
      {
         if ( rMI[ireg] == irow+1 ) {                                    // if this thread has row irow+1 cached (i.e. the next row to eliminate)
            // write out
            if ( tidmngr > irow ) {                                       // if the column is in the upper triangular
               // normalize
               rME[ireg] *= *sNormal;                                      // normalize
               sCurrentRow[tidmngr] = rME[ireg];                           // write to shared
               MatrixArrayRM[(tidmngr)+(irow+1)*ngr] = rME[ireg];          // also write to GMEM (as it is the fully upper-triangularized result)
            }
            if ( tidmngr == 0 ) {                                         // use the zeroth thread in the row to normalize the RHS
               sRhs[rMI[ireg]] *= *sNormal;
            }
         }

         if ( tidmngr == irow+1 ) {
            if (rMI[ireg] < ngr) {
               sCurrentColumn[rMI[ireg]] = rME[ireg];
            }
         }

      }

      // syncthreads ------------------------------------------------
      __syncthreads();    

   }

   // output updated rhs
   if ( tid < ngr ) {
      rhs[tid] = sRhs[tid];
   }
}

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