Solution 4.2

C (OpenMP 3.1)

// solves 2-D Laplace equation using a relaxation scheme

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifdef _OPENMP
#include <omp.h>
#endif

int main() {
   
   double *T, *Tnew, *Tmp;
   double tol, var = DBL_MAX, top = 100.0;
   unsigned n, n2, maxIter, i, j, iter = 0;
   int itemsread;
   FILE *fout;

   printf("Enter mesh size, max iterations and tolerance: ");
   itemsread = scanf("%u ,%u ,%lf", &n, &maxIter, &tol);

   if (itemsread!=3) {
      fprintf(stderr, "Input error!\n");
      exit(-1);
   }

#ifdef _OPENMP 
    double startTime = omp_get_wtime();
#else
    time_t startTime = clock();
#endif

   n2 = n+2;
   T = calloc(n2*n2, sizeof(*T));
   Tnew = calloc(n2*n2, sizeof(*T));

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

   // set boundary conditions

   for (i=1; i<=n; i++) {
      T[(n+1)*n2+i] = Tnew[(n+1)*n2+i] = i * top / (n+1);  
      T[i*n2+n+1] = Tnew[i*n2+n+1] = i * top / (n+1);
   }
#pragma omp parallel
 {
   while(var > tol && iter <= maxIter) {
#pragma omp barrier
#pragma omp single 
  {
      ++iter;
      var = 0.0;
  }
#pragma omp for private(j) reduction(max:var)
      for (i=1; i<=n; ++i)
         for (j=1; j<=n; ++j) {         
            Tnew[i*n2+j] = 0.25*( T[(i-1)*n2+j] + T[(i+1)*n2+j] 
                                + T[i*n2+(j-1)] + T[i*n2+(j+1)] );
        
            var = fmax(var, fabs(Tnew[i*n2+j] - T[i*n2+j]));
         }

#pragma omp single nowait
  {
         Tmp=T; T=Tnew; Tnew=Tmp;

         if (iter%100 == 0)
            printf("iter: %8u, variation = %12.4lE\n", iter, var);

  }

   }
 }

#ifdef _OPENMP
    double endTime = omp_get_wtime() - startTime;
#else
    double endTime = (clock() - startTime) / (double) CLOCKS_PER_SEC;
#endif

   printf("Elapsed time (s) = %.2lf\n", endTime);
   printf("Mesh size: %u\n", n);
   printf("Stopped at iteration: %u\n", iter);
   printf("Maximum error: %lE\n", var);

   // saving results to file
   fout = fopen("results", "w");
   if (fout == NULL) {
      perror("results");
      exit(-1);
   }
   for (i=1; i<=n; ++i)
      for (j=1; j<=n; ++j)
         fprintf(fout, "%8u %8u %18.9lE\n", i, j, T[i*n2+j]);
   fclose(fout);
   free(T);
   free(Tnew);
   return 0;
}

C (OpenMP 3.0 or older)

// solves 2-D Laplace equation using a relaxation scheme

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifdef _OPENMP
#include <omp.h>
#endif

int main() {
   
   double *T, *Tnew, *Tmp;
   double tol, var = DBL_MAX, top = 100.0;
   unsigned n, n2, maxIter, i, j, iter = 0;
   int itemsread;
   FILE *fout;

   printf("Enter mesh size, max iterations and tolerance: ");
   itemsread = scanf("%u ,%u ,%lf", &n, &maxIter, &tol);

   if (itemsread!=3) {
      fprintf(stderr, "Input error!\n");
      exit(-1);
   }

#ifdef _OPENMP 
    double startTime = omp_get_wtime();
#else
    time_t startTime = clock();
#endif

   n2 = n+2;
   T = calloc(n2*n2, sizeof(*T));
   Tnew = calloc(n2*n2, sizeof(*T));

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

   // set boundary conditions

   for (i=1; i<=n; i++) {
      T[(n+1)*n2+i] = Tnew[(n+1)*n2+i] = i * top / (n+1);  
      T[i*n2+n+1] = Tnew[i*n2+n+1] = i * top / (n+1);
   }
#pragma omp parallel
 {
   while(var > tol && iter <= maxIter) {
#pragma omp barrier
#pragma omp single 
  {
      ++iter;
      var = 0.0;
  }
   double pvar = 0.0;
#pragma omp for nowait private(j)
      for (i=1; i<=n; ++i)
         for (j=1; j<=n; ++j) {         
            Tnew[i*n2+j] = 0.25*( T[(i-1)*n2+j] + T[(i+1)*n2+j] 
                                + T[i*n2+(j-1)] + T[i*n2+(j+1)] );
        
            pvar = fmax(pvar, fabs(Tnew[i*n2+j] - T[i*n2+j]));
         }
#pragma omp critical
   if (pvar > var) var = pvar;

#pragma omp barrier
#pragma omp single nowait
  {
         Tmp=T; T=Tnew; Tnew=Tmp;

         if (iter%100 == 0)
            printf("iter: %8u, variation = %12.4lE\n", iter, var);

  }

   }
 }

#ifdef _OPENMP
    double endTime = omp_get_wtime() - startTime;
#else
    double endTime = (clock() - startTime) / (double) CLOCKS_PER_SEC;
#endif

   printf("Elapsed time (s) = %.2lf\n", endTime);
   printf("Mesh size: %u\n", n);
   printf("Stopped at iteration: %u\n", iter);
   printf("Maximum error: %lE\n", var);

   // saving results to file
   fout = fopen("results", "w");
   if (fout == NULL) {
      perror("results");
      exit(-1);
   }
   for (i=1; i<=n; ++i)
      for (j=1; j<=n; ++j)
         fprintf(fout, "%8u %8u %18.9lE\n", i, j, T[i*n2+j]);
   fclose(fout);
   free(T);
   free(Tnew);
   return 0;
}

Fortran

program laplace
#ifdef _OPENMP
   use omp_lib
#endif
   implicit none
   integer, parameter                 :: dp=kind(1.d0)
   integer                            :: n, maxIter, i, j, iter = 0
   real (dp), dimension(:,:), pointer :: T, Tnew, Tmp=>null()
   real (dp)                          :: tol, var = 1.d0, top = 100.d0
   integer                            :: ierr
#ifdef _OPENMP
   real(kind(1.d0)) :: startTime, endTime
#else 
   real :: startTime, endTime
#endif

   write(*,*) 'Enter mesh size, max iterations and tollerance:'
   read(*,*,iostat=ierr)  n, maxIter, tol
   if(ierr /= 0)  STOP 'Input error!'

#ifdef _OPENMP
   startTime=omp_get_wtime()
#else 
   call cpu_time(startTime)
#endif

   allocate (T(0:n+1,0:n+1), Tnew(0:n+1,0:n+1),stat=ierr)
   if(ierr/=0) STOP 'T Tnew matrix allocation failed'

   T(0:n,0:n) = 0.d0

   T(n+1,1:n) = (/ (i, i=1,n) /) * (top / (n+1))
    
   T(1:n,n+1) = (/ (i, i=1,n) /) * (top / (n+1))

   Tnew = T
!$omp parallel
   do while (var > tol .and. iter <= maxIter)
!$omp barrier
!$omp single
      iter = iter + 1
      var = 0.d0       
!$omp end single 
!$omp do reduction(max:var)    
      do j = 1, n
         do i = 1, n
            Tnew(i,j) = 0.25d0 * ( T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1) )
            var = max(var, abs( Tnew(i,j) - T(i,j) ))
         end do
      end do
!$omp end do
!$omp single
      Tmp =>T; T =>Tnew; Tnew => Tmp; 

      if( mod(iter,100) == 0 ) write(*,"(a,i8,e12.4)") &
         ' iter, variation:', iter, var
!$omp end single nowait
   end do
!$omp end parallel
#ifdef _OPENMP
   endTime=omp_get_wtime()
#else
   call cpu_time(endTime)
#endif

   write(*,'(/A,F10.4)') ' Elapsed time (s)     =', endTime - startTime
   write(*,*) 'Mesh size            =', n
   write(*,*) 'Stopped at iteration =', iter
   write(*,*) 'The maximum error    =', var
    
   open(10, file='results', action='write', iostat=ierr)
   if(ierr /= 0)  STOP 'Error in opening output file!'
   write(10, "(i8, i8, e18.9)") (( i, j, T(i,j), i=1,n), j=1,n)
   close(10)

   deallocate (T, Tnew)
   nullify(Tmp)

end program laplace