#include <stdio.h>
#include <math.h>
#include <time.h>
#ifndef N
#define N 500
#endif
float a[N][N],b[N][N],c[N][N];
float ci[N][N],cs[N][N],cd[N][N];

main()
{

/* Petersen and Arbenz book exercise 3.1, variants of matrix
   multiply. Problem. 
   BLAS routines SAXPY and SDOT, using Fortran ordering.

        cc -O3 thisfile.c
        g77 thisfile.o -lblas

                                          wpp 12 Dec. 2003 
*/

   int i,j,k,n=N,one=1;
   float t1,t2,t3,t4,bjk,big,mxer1,mxer2,mxer3,er1,er2,er3,fnm2,ggl();
   void saxpy_(int*,float*,float*,int*,float*,int*);
   float sdot_(int*,float*,int*,float*,int*);
   static float seed=331.0;

   fnm2 = 1.0/((float)N);

   for(j=0;j<N;j++){
      for(i=0;i<N;i++){
          a[j][i] = ggl(&seed);
          b[j][i] = ggl(&seed);
      }
   }

/* outer product procedure */
   t1   = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(j=0;j<N;j++){
      bjk = b[j][0];
      for(i=0;i<N;i++) c[j][i] = a[0][i]*bjk;
      for(k=1;k<N;k++){
         bjk = b[j][k];
#pragma ivdep
         for(i=0;i<N;i++){
            c[j][i] +=  bjk * a[k][i];
         } 
      }
   }
  
   t1   = ((float)clock())/((float) CLOCKS_PER_SEC) - t1;
   printf(" -------------------------------- \n");
   printf(" Compare inner vs. outer methods: \n");
   printf(" -------------------------------- \n");
   printf(" outer product method time t1=%e sec.\n",t1);

/* inner product procedure */
   t2   = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         ci[j][i] = 0.0;
#pragma ivdep
         for(k=0;k<N;k++){
            ci[j][i] +=  a[k][i] * b[j][k];
         }
      }
   }

   t2   = ((float)clock())/((float) CLOCKS_PER_SEC) - t2;
   printf(" inner product method time t2=%e sec.\n",t2);
   er1   = 0.0;
   mxer1 = 0.0;
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         er1 += (ci[j][i] - c[j][i])*(ci[j][i] - c[j][i]);
         big    = fabs(ci[j][i] - c[j][i]);
         mxer1  = (big>mxer1)?big:mxer1;
      }
   }
   er1 = fnm2*sqrt(er1);
   printf(" Error compare inner vs. outer:\n");
   printf(" RMS error for %d x %d matrix mult. = %15.8e\n",N,N,er1);
   printf(" max error for %d x %d matrix mult. = %15.8e\n\n",N,N,mxer1);
   printf(" -------------- \n");
   printf(" BLAS variants:\n");
   printf(" -------------- \n");
/* outer product procedure via SAXPY */
   t3   = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(j=0;j<N;j++){
      bjk = b[j][0];
      for(i=0;i<N;i++) cs[j][i] = a[0][i]*bjk;
      for(k=1;k<N;k++){
         bjk = b[j][k];
         saxpy_( &n, &bjk, &a[k][0], &one, &cs[j][0], &one );
      }
   }
  
   t3   = ((float)clock())/((float) CLOCKS_PER_SEC) - t3;
   printf(" outer product method by SAXPY time t3=%e sec.\n",t3);

/* inner product procedure via SDOT */
   t4   = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         cd[j][i] = sdot_( &n, &a[0][i], &n, &b[j][0],&one );
      }
   }
	
   t4   = ((float)clock())/((float) CLOCKS_PER_SEC) - t4;
   printf(" inner product by SDOT method time  t4=%e sec.\n",t4);
   er2   = 0.0;
   er3   = 0.0;
   mxer2 = 0.0;
   mxer3 = 0.0;
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         er2 += (cs[j][i] - c[j][i])*(cs[j][i] - c[j][i]);
         er3 += (cd[j][i] - c[j][i])*(cd[j][i] - c[j][i]);
         big    = fabs(cs[j][i] - c[j][i]);
         mxer2  = (big>mxer2)?big:mxer2;
         big    = fabs(cd[j][i] - c[j][i]);
         mxer3  = (big>mxer3)?big:mxer3; 
      }
   }
   er2 = fnm2*sqrt(er2);
   er3 = fnm2*sqrt(er3);
   printf(" Error compare SDOT vs. SAXPY variants: \n");
   printf(" RMS err saxpy %d x %d matrix mult. = %15.8e\n",N,N,er2);
   printf(" max err saxpy %d x %d matrix mult. = %15.8e\n",N,N,mxer2);
   printf(" RMS err sdot  %d x %d matrix mult. = %15.8e\n",N,N,er3);
   printf(" max err sdot  %d x %d matrix mult. = %15.8e\n",N,N,mxer3);
}
float ggl(float *ds)
{
/* generate u(0,1) distributed random numbers.
   Seed ds must be saved between calls. ggl is
   essentially the same as the IMSL routine RNUM.

   W. Petersen and M. Troyer, 24 Oct. 2002, ETHZ */

   double t,d2=0.2147483647e10;
   t   = (double) *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = (float) t;
   return((float)((t-1.0e0)/(d2-1.0e0)));
}