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

    -- Developed by:
       Stan Tomov
       Rajib Nath

    -- The implementation of these SGEMM routines for the Fermi GPUs is based 
       on "auto-tuning" the algorithm 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=96 blk_N=96 blk_K=16 nthd_x=64 nthd_y=4
*/

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

texture<float,1>  tex_x_float_A;
texture<float,1>  tex_x_float_B;

static __inline__ __device__ float fetch_x_A(const int& i, const float * x)
{
  return tex1Dfetch(tex_x_float_A, i);
}
static __inline__ __device__ float fetch_x_B(const int& i, const float * x)
{
  return tex1Dfetch(tex_x_float_B, i);
}


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

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

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

    __shared__ float Bb[16][97];
    __shared__ float Abs[96][17];

    float xxA[6];
    float xxB[6];
    
    B+= res + __mul24(iby + qot * 6, ldb );
    A+= ibx + __mul24( qot, lda) + res ; 

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

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

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

    __syncthreads();

    const float *Bend = B + k-16;
   
    float Axs[6];
    float Bxp[6];

    float Cb[36] = {0,0,0,0,0,0, 0,0,0,0,0,0,  0,0,0,0,0,0, 0,0,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 j1=0;j1<8;j1++){

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

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

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

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

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

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

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

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

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

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

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

	__syncthreads();

    } 
    while (B < Bend);

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

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

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

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

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

//========================================================================

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

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

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

    __shared__ float Bb[16][97];
    __shared__ float Abs[96][17];

    float xxA[6];
    float xxB[6];

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

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

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

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

    __syncthreads();

    const float *Bend = B + k-16;
   
    float Axs[6];
    float Bxp[6];

    float Cb[36] = {0,0,0,0,0,0, 0,0,0,0,0,0,  0,0,0,0,0,0, 0,0,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<6; y++)
	   xxB[y] = fetch_x_B( trackB + y*ldb, B);

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

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

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

            #pragma unroll 
            for( int x=0;x<6; x++){
                #pragma unroll 
                for(int y=0; y<6; y++)
                    Cb[x*6 + y] += Axs[x]*Bxp[y];
	    }
        }
           
	__syncthreads();
	#pragma unroll
	for(int y=0; y<6; y++)
	   Abs[qot + 16*y][res] = xxA[y]; 
	
	#pragma unroll
	for( int y=0; y<6; y++)
	   Bb[res][qot*6 + 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<6; y++)
            Bxp[y] = Bb[j1][qot + y*16];

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

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

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

//========================================================================

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

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

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

    __shared__ float Bb[16][97];
    __shared__ float Abs[96][17];

    float xxA[6];
    float xxB[6];

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

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

    #pragma unroll
    for(int y=0; y<6; y++)
      Bb[qot][res+16*y] = fetch_x_B(trackB+16*y, B);

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

    __syncthreads();

    const float *Bend = B + k*ldb - 16*ldb;
   
    float Axs[6];
    float Bxp[6];
    
    float Cb[36] = {0,0,0,0,0,0, 0,0,0,0,0,0,  0,0,0,0,0,0, 0,0,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;

	#pragma unroll
        for( int y=0; y<6; y++)
           xxB[y] = fetch_x_B(trackB + 16*y, B);
		
        #pragma unroll
        for( int y=0; y<6; y++)
	   xxA[y] = fetch_x_A(trackA + lda*y*16, A);

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

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

           #pragma unroll 
           for( int x=0; x<6; x++){
               #pragma unroll 
               for( int y=0;y<6;y++){
                  Cb[x*6+y] += Axs[x]*Bxp[y];
	       }
	    }
        }
            
	__syncthreads();
	#pragma unroll
	for( int y=0; y<6; y++)
	   Abs[qot + 16*y][res] = xxA[y];
 
	#pragma unroll
	for( int y=0; y<6; y++)
           Bb[qot][res+y*16] = 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<6; y++)
            Bxp[y]= Bb[j1][qot + y*16];

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

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

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


//========================================================================

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

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

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

    __shared__ float Bb[16][97];
    __shared__ float Abs[96][17];

    float xxA[6];
    float xxB[6];

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

    #pragma unroll
    for(int y=0; y<6; y++)
        Bb[qot][res+16*y] = fetch_x_B(trackB+16*y, B);

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

    __syncthreads();

    const float *Bend = B + k*ldb - 16*ldb;
   
    float Axs[6];
    float Bxp[6];

    float Cb[36] = {0,0,0,0,0,0, 0,0,0,0,0,0,  0,0,0,0,0,0, 0,0,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;

        #pragma unroll
        for( int y=0; y<6; y++)
          xxB[y] = fetch_x_B( trackB + 16*y, B);

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

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

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

	#pragma unroll
	for( int y=0; y<6; y++)
            Bb[qot][res+y*16] = 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<6; y++)
                      Bxp[y] = Bb[j1][qot + y*16];

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

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

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

//=================================================================================

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

extern "C" void
magmablas_fermi_sgemm(char TRANSA, char TRANSB, int m , int n , int k , float alpha, 
                      const float *A, int lda, const float *B, int ldb, float beta, 
                      float *C, int ldc ) 
{

        if (m<=0 || n<=0 || k<=0)
           return;

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

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

        errt = cudaBindTexture(0,tex_x_float_B,B,channelDesc);
	if( errt != cudaSuccess) printf("can not bind to texture \n");

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

        if( TRANSB == 'T' || TRANSB == 't') 
	  if( TRANSA == 'N' ||  TRANSA == 'n') 
           fermiSgemm_v2_kernel_NT<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                       ldc, alpha, beta);
	  else
           fermiSgemm_v2_kernel_TT<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                       ldc, alpha, beta);
        else
	  if( TRANSA == 'N' || TRANSA == 'n') 
           fermiSgemm_v2_kernel_NN<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                       ldc, alpha, beta);
          else
           fermiSgemm_v2_kernel_TN<<< grid, threads>>>(C, A, B, m, n, k, lda, ldb, 
                                                       ldc, alpha, beta);

 	cudaUnbindTexture ( tex_x_float_A ) ;
	cudaUnbindTexture ( tex_x_float_B ) ;
}



