/* icc -tpp2 -parallel -par-threshold0 -par-report3 -openmp -lm */

#include <stdio.h>
#include <time.h>
#include <omp.h>
#include <math.h>

#define N 900

float a[N][N],b[N][N],c[N][N];

int main()
{
   int i,j,k;
   float seed,t,ggl(float*);

   int num_threads=omp_get_max_threads();
   printf(" \n number of processors: %d\n",num_threads);
   
   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);
      }
   
   t   = (float)clock()/(float) CLOCKS_PER_SEC;
   
   for(j=0;j<N;j++)
      for(i=0;i<N;i++)
      {  c[j][i] = 0.0;
         for(k=0;k<N;k++)
            c[j][i] += b[k][i]*a[j][k]; 
      }
   
   t   = (float)clock()/(float) CLOCKS_PER_SEC - t;
   
   printf(" c[555][555]=%15.8e\n",c[555][555]);
   printf("\n CPU time=%e seconds\n",t);
   
   return 0;
}

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)));
}