/* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}