#include <stdio.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define N 500
float a[N][N],b[N][N],c[N][N];
float cc[N][N];
main()
{
   int i,j,k;
   float seed,t1,t2,bjk,error,fnm2,ggl(float*);

   fnm2 = 1.0/((float)N);
   fnm2 = fnm2*fnm2;

   seed = 331.0;          /* seed for ggl */

   for(j=0;j<N;j++){
      for(i=0;i<N;i++){
          a[j][i] = ggl(&seed);
          b[j][i] = ggl(&seed);
      }
   }
   t1   = ((float)clock())/((float) CLOCKS_PER_SEC);
#pragma omp parallel shared(a,b,c) private(i,j,k,bjk)
#pragma omp parallel for
   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];
         for(i=0;i<N;i++){
            c[j][i] += a[k][i]*bjk;
         } 
      }
   }
   t1   = ((float)clock())/((float) CLOCKS_PER_SEC) - t1;
   printf(" t1=%e\n",t1);
   t2   = ((float)clock())/((float) CLOCKS_PER_SEC);
#pragma omp shared(a,b,c) private(i,j,k)
#pragma omp parallel for
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         cc[j][i] = 0.0;
         for(k=0;k<N;k++){
            cc[j][i] += a[k][i]*b[j][k]; 
         }
      }
   }
   t2   = ((float)clock())/((float) CLOCKS_PER_SEC) - t2;
   printf(" t2=%e\n",t2);
   error = 0.0;
   for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
         error += (cc[j][i] - c[j][i])*(cc[j][i] - c[j][i]);
      }
   }
   error = sqrt(error*fnm2);
   printf(" RMS error for %d x %d matrix mult. = %15.8e\n",
   N,N,error);
}
#include <math.h>
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:
   a modification of a fortran version from
   I. Vattulainen, Tampere Univ. of Technology,
   Finland, 1992 */
                                                                                
   double t,d2=0.2147483647e10;
   t   = (float) *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = (float) t;
   return((float) ((t-1.0e0)/(d2-1.0e0)));
}