#include <stdio.h>
#include <stdlib.h>

double ddot_(int *L, const double *A, int *incx, const double *B, int *incy);

int main(  int argc, char *argv[] )
{ 

   double *A, *B, *Bt, *C, *tC;
   int L, M, N, L0, L1;

   int i, j, k, t, st, incr[2] = { 1, 1};;
   double v, rmin, rmax;

      fprintf(stdout,"Input dimensions of matrices  L, M, N:\n");
      fscanf(stdin,"%d, %d, %d",&L, &M, &N);
      fprintf(stdout," Matrices dimensions are: %d, %d, %d\n",L,M,N);

         A = (double*) malloc(sizeof((double)1.0)*L*M);
         B = (double*) malloc(sizeof((double)1.0)*N*M);
         C = (double*) malloc(sizeof((double)1.0)*L*N);
         if ( A == NULL || B == NULL || C == NULL ) {
            fprintf(stderr,"Error allocating matrices \n");
      	    return(-1);
         }


   /*  Initializations */

 	   for ( i = 0; i < L*M; i++ ) A[i] = 0.0;
   	   for ( i = 0; i < N*M; i++ ) B[i] = 0.0;
   	   for ( i = 0; i < L*N; i++ ) C[i] = 0.0;

   	   /* Row-wise store */
  	   for ( i = 0; i < L; i++ ) {
   	   	   for ( k = 0; k < M; k++ ) {
   	   	   	   A[ i*M + k] = (double)(i-k) / (double)(i+k+1);
   	   	   }
   	   }
   	   
   	   /* Row-wise store */
   	  for ( j = 0; j < N; j++ ) {
   	  	  for ( k = 0; k < M; k++ ) {
   	  	  	  B[ k*N + j ] = 0.01 * (double)(k-j) / (double)(k+j+1);
   	  	  }
   	  }

   /*  Compute multiplication using BLAS Fortran DDOT function 
	   First of all transpose B matrix
	*/
	
      Bt = (double*) malloc(sizeof((double)1.0)*N*M);
	// Compute Bt = B transpose:
		.	.	.
		
	// Call ddot_ for computing matrix multiply:
		.	.	.
	 
      /* Compare results */ 
      tC = (double*) malloc(sizeof((double)1.0)*L*N);
      if ( tC == NULL ) {
         fprintf(stderr," Error allocating matrix tC \n");
      }

      for ( i = 0; i < L; i++ ) {
         for ( k = 0; k < M; k++ ) {
            A[ i*M + k] = (double)(i-k) / (double)(i+k+1);
         }
      }
   
	// Compute matrix multiply using for loops (result in tC):
		.	.	.
		
      rmin = rmax = abs(tC[0] - C[0]);
      for ( j = 0; j < N; j++ ) {
         for ( i = 0; i < L; i++ ) {
            v = abs(tC[ i*N + j ] - C[ i*N + j ]);
            if ( v < rmin ) {
               rmin = v;
            }
            if ( v > rmax ) {
               rmax = v;
            }
         }
      }

      fprintf(stderr," Check results: rmin, rmax = %f, %f\n",rmin,rmax);
     	   	   
  
   return(0);
}