Exercise 2

The following code performs a serial matrix multiplication. Try to parallelize it with OpenMP, acting only on the most important loop. Try also to add the timing functions before and after the loop and print the elapsed time.

C source code

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

   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 source code

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

   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