/*
    -- MAGMA (version 0.3) --
       Univ. of Tennessee, Knoxville
       Univ. of California, Berkeley
       Univ. of Colorado, Denver
       June 2010

    -- Developed by:
       Rajib Nath 
       Stan Tomov

    -- The implementation of these DGEMM routines for the Fermi GPUs are described in
       Nath, R., Tomov, S., Dongarra, J. "An Improved MAGMA GEMM for Fermi GPUs," 
       University of Tennessee Computer Science Technical Report, UT-CS-10-655 
       (also LAPACK working note 227), July 29, 2010. 

       These routines will be included in CUBLAS 3.2.
*/

/*
    blk_M=64 blk_N=64 blk_K=16 nthd_x=64 nthd_y=4
*/

#include <cuda_runtime.h>
#include <cuda.h>

texture<int2,1>  tex_x_double_A;
texture<int2,1>  tex_x_double_B;

static __inline__ __device__ double fetch_x_A(const int& i)
{
  register int2  v = tex1Dfetch(tex_x_double_A, i);
  return __hiloint2double(v.y, v.x);
}
static __inline__ __device__ double fetch_x_B(const int& i)
{
  register int2  v = tex1Dfetch(tex_x_double_B, i);
  return __hiloint2double(v.y, v.x);
}


__global__ void fermiDgemm_v2_kernel_NN(double *C, const double *A, const double *B,  
                                        int m, int n, int k, int lda, int ldb,  
                                        int ldc, double alpha, double beta) 
{
    const  int tx = threadIdx.x;
    const  int ty = threadIdx.y;

    const int iby = blockIdx.y * 64;
    const int ibx = blockIdx.x * 64;
    const int idt = ty * 64 + tx;

    const int res = idt%16;
    const int qot = idt/16;

    __shared__ double  Bb[16][65];
    __shared__ double Abs[64][17];

    double xxA[4];
    double xxB[4];
    
    B+= res+ __mul24(iby + qot * 4, ldb );
    int trackB = res+ __mul24(iby + qot * 4, ldb );

    A += ibx +__mul24( qot, lda) + res ; 
    int trackA =  ibx +__mul24( qot, lda) + res ;

    #pragma unroll
    for(int y=0; y<4; y++)
       Bb[res][qot*4+y] = fetch_x_B( trackB + y * ldb ) ;

    #pragma unroll
    for(int y=0; y<4; y++)
      Abs[res+ y*16][qot] = fetch_x_A(trackA + y*16) ;
    __syncthreads();

    const double *Bend = B + k-16;
   
    double Axs[4];
    double Bxp[4];

    double Cb[16] = {0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0};

     do 
     {
	B += 16;
	A += lda *16  ;
	trackA += 16*lda ; 
	trackB += 16;

        #pragma unroll
        for( int y=0; y<4; y++)
  	  xxB[y] = fetch_x_B( trackB + y*ldb);

        #pragma unroll
        for( int y=0; y<4; y++)
	  xxA[y] = fetch_x_A(trackA + y*16);

        #pragma unroll 
        for( int j1=0;j1<16;j1++){

             #pragma unroll
             for( int y=0; y<4; y++)
                  Bxp[y]= Bb[j1][qot+y*16];

             #pragma unroll
             for( int y=0; y<4; y++)
                  Axs[y] =  Abs[res+y*16][j1] ;

             #pragma unroll 
             for( int x=0; x<4; x++){
                  #pragma unroll 
                  for( int y=0; y<4; y++){
                      Cb[x*4+y]  += Axs[x]*Bxp[y];
		  }
	     }
        }

	__syncthreads();
	#pragma unroll
	for(int y=0; y<4; y++)
	    Abs[res+y*16][qot] = xxA[y]; 

	#pragma unroll
	for(int y=0; y<4; y++)
	    Bb[res][qot*4 + y] = xxB[y];

	__syncthreads();
     }
     while (B < Bend);

     C += res + ibx  + __mul24 (qot +  iby ,ldc);

     #pragma unroll 
     for(int j1=0;j1<16;j1++){

              #pragma unroll
              for( int y=0; y<4; y++)
                  Bxp[y]= Bb[j1][qot + y*16];

              #pragma unroll
              for( int y=0; y<4; y++)
                  Axs[y] =  Abs[res + y*16][j1] ;

              #pragma unroll 
              for( int x=0; x<4; x++){
                  #pragma unroll 
                  for( int y=0;y<4; y++){
                      Cb[x*4 + y]  += Axs[x]*Bxp[y];
		  }
	      }
    }

    #pragma unroll
    for( int y=0;y<4;y++){

       #pragma unroll
       for(int x=0;x<4;x++){
          C[x*16] = alpha*Cb[y+x*4] + beta * C[x*16];
       }
	   
       C += ldc*16;
    }
}

__global__ void fermiDgemm_v2_kernel_TN(double *C, const double *A, const double *B,
                                        int m, int n,  int k,  int lda,  int ldb,  
                                        int ldc, double alpha, double beta) 
{
    const  int tx = threadIdx.x;
    const  int ty = threadIdx.y;

    const int iby = blockIdx.y * 64;
    const int ibx = blockIdx.x * 64;
    const int idt = ty * 64 + tx;

    const int res = idt%16;
    const int qot = idt/16;

    __shared__ double Bb[16][65];
    __shared__ double Abs[64][17];

    double xxA[4];
    double xxB[4];

    B+= res+ __mul24(iby+qot*4, ldb );
    int trackB = res+ __mul24(iby+qot*4, ldb );

    A+= __mul24( ibx + qot, lda )   + res; 
    int trackA = __mul24( ibx + qot, lda ) + res;

    #pragma unroll
    for(int y=0; y<4; y++)
	Bb[res][qot*4+y] =fetch_x_B( trackB + y*ldb );

    #pragma unroll
    for(int y=0; y<4; y++)
	Abs[qot+16*y][res] = fetch_x_A(trackA +  lda*16*y);
    __syncthreads();

    const double *Bend = B + k-16;
   
    double Axs[4];
    double Bxp[4];

    double Cb[16] = {0,0,0,0,    0,0,0,0, 0,0,0,0, 0,0,0,0};

    do 
    {
	B += 16;
	A += 16  ;
	trackA+=16 ; 
	trackB+=16;

        #pragma unroll
        for(int y=0; y<4; y++)
	   xxB[y]=fetch_x_B( trackB + y*ldb);

	#pragma unroll
        for(int y=0; y<4; y++)
	   xxA[y] =fetch_x_A(trackA +  lda*y*16);

	#pragma unroll 
        for(int j1=0;j1<16;j1++){
                  #pragma unroll
                  for(int y=0; y<4; y++)
                      Bxp[y]= Bb[j1][qot+y*16];

                  #pragma unroll
                  for(int y=0; y<4; y++)
                      Axs[y] =  Abs[res+y*16][j1];

                  #pragma unroll 
                  for(int x=0; x<4; x++){
                     #pragma unroll 
                     for(int y=0; y<4; y++){
                         Cb[x*4+y]  += Axs[x]*Bxp[y];
		     }
		  }
	}

	__syncthreads();
	#pragma unroll
	for(int y=0; y<4; y++)
	   Abs[qot+16*y][res] = xxA[y];
 
	#pragma unroll
	for(int y=0; y<4; y++)
	   Bb[res][qot*4+y] =xxB[y];
	__syncthreads();
    } 
    while (B < Bend);

    C += res + ibx  + __mul24 (qot + iby ,ldc);

    #pragma unroll 
    for(int j1=0; j1<16; j1++){
	#pragma unroll
        for(int y=0; y<4; y++)
            Bxp[y]= Bb[j1][qot+y*16];

        #pragma unroll
        for(int y=0; y<4; y++)
            Axs[y] = Abs[res+y*16][j1];

        #pragma unroll 
        for(int x=0; x<4; x++){
            #pragma unroll 
            for(int y=0; y<4; y++){
               Cb[x*4+y] += Axs[x]*Bxp[y];
	    }
	}
    }

    #pragma unroll
    for(int y=0;y<4;y++){
	#pragma unroll
        for(int x=0;x<4;x++)
	   C[x*16] =alpha*Cb[y+x*4] + beta * C[x*16];
	   
	C+=ldc*16;
    }
}


__global__ 
void fermiDgemm_v2_kernel_TT(double *C, const double *A, const double *B,  int m, int n,  
                             int k,  int lda,  int ldb,  int ldc, double alpha, double beta) 
{
    const  int tx = threadIdx.x;
    const  int ty = threadIdx.y;

    const int iby = blockIdx.y * 64;
    const int ibx = blockIdx.x * 64;
    const int idt = ty * 64 + tx;

    const int res = idt% 16;
    const int qot = idt/ 16;

    __shared__ double Bb[16][65];
    __shared__ double Abs[64][17];

    double xxA[4];
    double xxB[4];

    B += iby + tx + __mul24(ty , ldb );
    A += __mul24( ibx + qot , lda )   + res; 

    int trackA = __mul24( ibx + qot, lda ) + res;
    int trackB =  iby+ tx + __mul24(ty , ldb );

    Bb[ty+0*4][tx] = fetch_x_B(trackB+ldb*0);
    Bb[ty+1*4][tx] = fetch_x_B(trackB+ldb*4);
    Bb[ty+2*4][tx] = fetch_x_B(trackB+ldb*8);
    Bb[ty+3*4][tx] = fetch_x_B(trackB+ldb*12);
    Abs[qot +16*0][res] = fetch_x_A(trackA +  lda*16*0) ;
    Abs[qot +16*1][res] = fetch_x_A(trackA +  lda*16*1) ;
    Abs[qot +16*2][res] = fetch_x_A(trackA +  lda*16*2) ;
    Abs[qot +16*3][res] = fetch_x_A(trackA +  lda*16*3) ;
    __syncthreads();

    const double *Bend = B + k*ldb-16*ldb;
   
    double Axs[4];
    double Bxp[4];

    double Cb[16] = {0,0,0,0,    0,0,0,0, 0,0,0,0, 0,0,0,0};

    do 
    {
	B += 16*ldb;
	A += 16  ;
	trackA+=16 ; 
	trackB+=16*ldb;

	xxB[0]=fetch_x_B( trackB + 0*ldb) ;
	xxB[1]=fetch_x_B( trackB + 4*ldb) ;
	xxB[2]=fetch_x_B( trackB + 8*ldb) ;
	xxB[3]=fetch_x_B( trackB + 12*ldb) ;

	xxA[0] =fetch_x_A(trackA +  lda*0*16) ;
	xxA[1]=fetch_x_A(trackA +  lda*1*16 ) ;
	xxA[2]=fetch_x_A(trackA +  lda*2*16 ) ;
	xxA[3]=fetch_x_A(trackA +  lda*3*16 ) ;

        #pragma unroll 
        for( int j1=0;j1<16;j1++){
            #pragma unroll
            for( int y=0;y<4;y++)
                Bxp[y]= Bb[j1][qot + y*16];

            #pragma unroll
            for(int y=0;y<4;y++)
                Axs[y] =  Abs[res + y*16][j1];

            #pragma unroll 
            for( int x=0;x<4;x++){
                #pragma unroll 
                for( int y=0;y<4;y++){
                   Cb[x*4+y]  += Axs[x]*Bxp[y];
		}
	    }
        }
        
	__syncthreads();
	#pragma unroll
	for(int y=0;y<4;y++)
	   Abs[qot+16*y][res] = xxA[y];
 
	#pragma unroll
	for(int y=0;y<4;y++)
	   Bb[ty+y*4][tx] =xxB[y];

	__syncthreads();
    } 
    while (B < Bend);

    C += res + ibx  + __mul24 (qot + iby ,ldc);

    #pragma unroll 
    for( int j1=0;j1<16; j1++){
       #pragma unroll
       for( int y=0;y<4;y++)
           Bxp[y]= Bb[j1][qot+y*16];

       #pragma unroll
       for( int y=0;y<4;y++)
           Axs[y] =  Abs[res+y*16][j1];

       #pragma unroll 
       for( int x=0;x<4;x++){
            #pragma unroll 
            for( int y=0;y<4;y++){
               Cb[x*4+y]  += Axs[x]*Bxp[y];
	    }
	}
    }

    #pragma unroll
    for( int y=0;y<4;y++){
	#pragma unroll
        for(int x=0; x<4; x++)
	   C[x*16] =alpha*Cb[y+x*4] + beta * C[x*16];
		   
	C+=ldc*16;
    }
}

__global__ void fermiDgemm_v2_kernel_NT(double *C, const double *A, const double *B,  
                                        int m, int n,  int k,  int lda,  int ldb,  
                                        int ldc, double alpha, double beta) 
{
    const  int tx = threadIdx.x;
    const  int ty = threadIdx.y;

    const int iby = blockIdx.y * 64;
    const int ibx = blockIdx.x * 64;
    const int idt = ty * 64 + tx;

    const int res = idt%16;
    const int qot = idt/16;

    __shared__ double Bb[16][65];
    __shared__ double Abs[64][17];

    double xxA[4];
    double xxB[4];

    B+= iby+ tx + __mul24(ty , ldb );
    int trackB =  iby+ tx + __mul24(ty , ldb );

    A+= ibx +__mul24( qot, lda) + res ; 
    int trackA =  ibx +__mul24( qot, lda) + res ;

    Bb[ty+0*4][tx] = fetch_x_B(trackB+ldb*0);
    Bb[ty+1*4][tx] = fetch_x_B(trackB+ldb*4);
    Bb[ty+2*4][tx] = fetch_x_B(trackB+ldb*8);
    Bb[ty+3*4][tx] = fetch_x_B(trackB+ldb*12);

    Abs[res+0*16   ][qot] =fetch_x_A(trackA +  0*16) ;
    Abs[res+1*16   ][qot] =fetch_x_A(trackA +  1*16) ;
    Abs[res+2*16   ][qot] =fetch_x_A(trackA +  2*16) ;
    Abs[res+3*16   ][qot] =fetch_x_A(trackA +  3*16) ;
    __syncthreads();

    const double *Bend = B + k*ldb-16*ldb;
   
    double Axs[4];
    double Bxp[4];

    double Cb[16] = {0,0,0,0,    0,0,0,0, 0,0,0,0, 0,0,0,0};

    do 
    {
	B += 16*ldb;
	A += lda *16  ;
	trackA+=16*lda ; 
	trackB+=16*ldb;

	xxB[0]=fetch_x_B( trackB + 0*ldb) ;
	xxB[1]=fetch_x_B( trackB + 4*ldb) ;
	xxB[2]=fetch_x_B( trackB + 8*ldb) ;
	xxB[3]=fetch_x_B( trackB + 12*ldb) ;

	xxA[0]=fetch_x_A(trackA +  0*16) ;
	xxA[1]=fetch_x_A(trackA +  1*16 ) ;
	xxA[2]=fetch_x_A(trackA +  2*16 ) ;
	xxA[3]=fetch_x_A(trackA +  3*16 ) ;

        #pragma unroll 
        for( int j1=0;j1<16;j1++){
            #pragma unroll
            for( int y=0;y<4;y++)
                Bxp[y]= Bb[j1][qot+y*16];

            #pragma unroll
            for(int y=0;y<4;y++)
               Axs[y] = Abs[res+y*16][j1] ;

            #pragma unroll 
            for( int x=0;x<4;x++){
                #pragma unroll 
                for( int y=0;y<4;y++){
                   Cb[x*4+y]  += Axs[x]*Bxp[y];
		}
	    } 
        }
        
	__syncthreads();
	#pragma unroll
	for(  int y=0;y<4;y++)
	    Abs[res+y*16   ][qot] =xxA[y]; 

	#pragma unroll
	for(  int y=0;y<4;y++)
	    Bb[ty+y*4][tx] =xxB[y];
	__syncthreads();
    } 
    while (B < Bend);

    C += res + ibx  + __mul24 (qot + iby ,ldc);

    #pragma unroll 
    for( int j1=0;j1<16;j1++){
        #pragma unroll
        for( int y=0;y<4;y++)
            Bxp[y]= Bb[j1][qot+y*16];

        #pragma unroll
        for( int y=0;y<4;y++)
            Axs[y] =  Abs[res+y*16][j1];

        #pragma unroll 
        for( int x=0; x<4; x++){
            #pragma unroll 
            for( int y=0;y<4;y++){
                 Cb[x*4+y]  += Axs[x]*Bxp[y];
	    }
	}
    }

    #pragma unroll
    for( int y=0;y<4;y++){
	#pragma unroll
        for(int x=0;x<4;x++)
	   C[x*16] =alpha*Cb[y+x*4] + beta * C[x*16];
	   
	C+=ldc*16;
    }
}

#include<stdio.h>
#include "cublas.h"

extern "C" void
magmablas_fermi_dgemm(char TRANSA, char TRANSB, int m , int n , int k , 
                      double alpha, const double *A, int lda, const double *B, 
                      int ldb, double beta, double *C, int ldc ) 
{
   if (m<=0 || n<=0 || k<=0)
     return;

   if( m % (64) !=0 || n% (64)!=0 || k%(64) !=0 )
   {
	printf("Dimension Should Be multiple of %d\n", 64);
	printf("Calling cublasDgemm\n");
	cublasDgemm(TRANSA, TRANSB, m, n, k, alpha, A, lda, B,ldb,
                    beta, C, ldc);
  	return;
   }

   cudaChannelFormatDesc channelDesc = 
       cudaCreateChannelDesc(32, 32, 0, 0, cudaChannelFormatKindSigned);
   cudaError_t  errt = cudaBindTexture(0,tex_x_double_A,A,channelDesc);
   if( errt != cudaSuccess) printf("can not bind to texture \n");
   errt = cudaBindTexture(0,tex_x_double_B,B,channelDesc);
   if( errt != cudaSuccess) printf("can not bind to texture \n");

   dim3 threads( 64, 4 );
   dim3 grid(m/(64)+(m%(64)!=0),n/(64)+(n%(64)!=0));

   if( TRANSB == 'T' || TRANSB == 't') 
	if( TRANSA == 'N' ||  TRANSA == 'n') 
	    fermiDgemm_v2_kernel_NT<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                        ldc, alpha, beta);
	else
            fermiDgemm_v2_kernel_TT<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                        ldc, alpha, beta);
        else
	    if( TRANSA == 'N' || TRANSA == 'n') 
		fermiDgemm_v2_kernel_NN<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                            ldc, alpha, beta);
            else
                fermiDgemm_v2_kernel_TN<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                            ldc, alpha, beta);

   cudaUnbindTexture ( tex_x_double_A ) ;
   cudaUnbindTexture ( tex_x_double_B ) ;
}
