Solution 2

C

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

int main(int argc, char **argv)
{
   int n;
   int i, j, k;

   if(argc != 2) {
      fprintf(stderr,"Usage: %s  matrix size\n", argv[0]);
      exit(EXIT_FAILURE);
    }

   n = atoi(argv[1]);
   if ( n > 0) {
     printf("Matrix size is %d\n",n);
   }
   else {
     fprintf(stderr,"Error, matrix size is %d \n", n);
     exit(EXIT_FAILURE);
   }

   double (*a)[n] =  malloc(sizeof(double[n][n]));
   double (*b)[n] =  malloc(sizeof(double[n][n]));
   double (*c)[n] =  malloc(sizeof(double[n][n]));

   if ( a == NULL || b == NULL || c == NULL) {
       fprintf(stderr, "Not enough memory!\n");
       exit(EXIT_FAILURE);
    }

   for (i=0; i<n; i++) 
      for (j=0; j<n; j++) {
         a[i][j] = ((double)rand())/((double)RAND_MAX);
         b[i][j] = ((double)rand())/((double)RAND_MAX);
         c[i][j] = 0.0;
      }

#pragma omp parallel for private(k,j)
   for (i=0; i<n; i++) 
      for (k=0; k<n; k++)  
         for (j=0; j<n; j++) 
            c[i][j] += a[i][k]*b[k][j];

   //check a random element
   i = rand()%n;
   j = rand()%n;
   double d = 0.0;
   for (k=0; k<n; k++)
      d += a[i][k]*b[k][j];

   printf("Check on a random element: %18.9lE\n", fabs(d-c[i][j]));

   return 0;

}

FORTRAN

Program matrix_matrix_prod
   implicit none
   integer :: n
   real(kind(1.d0)), dimension(:,:), allocatable   :: a, b, c
   real(kind(1.d0)) :: d
   integer            :: i, j, k, ierr
   character(len=128) :: command
   character(len=80) :: arg

   call get_command_argument(0,command)
   if (command_argument_count() /= 1) then
     write(0,*) 'Usage:', trim(command), '   matrix size'
     stop
   else
     call get_command_argument(1,arg)
     read(arg,*) n
   endif
   if (n > 0 ) then
     write(*,*) 'Matrix size is ', n
   else
     write(0,*) "Error, matrix size is ", n
   endif

   allocate(a(n,n),b(n,n),c(n,n),stat=ierr)

   if(ierr/=0) STOP 'a,b,c matrix allocation failed'

   call random_number(a)
   call random_number(b)
   c = 0.d0

!$omp parallel do
   do j=1, n
      do k=1, n
         do i=1, n
            c(i,j) = c(i,j) + a(i,k)*b(k,j)
         end do
      end do
   end do

   call random_number(d)
   i = int( d*n+1)
   call random_number(d)
   j = int( d*n+1)
   d = 0.d0
   do k=1, n
         d = d + a(i,k)*b(k,j)
   end do

   write(*,*) "Check on a random element:" , abs(d-c(i,j))


end program matrix_matrix_prod