program main
c
c Jacobi iteration in Fortran 90 using OpenMP directives
c
       integer n,maxiters
       common /idat/ n,maxiters

c read values for n and maxiters (not shown),
c the value of n includes the boundaries

       call jacobi()
       stop
       end


       subroutine jacobi()

       integer n,maxiters
       common /idat/ n,maxiters

c allocate storage dynamically for the grids
       double precision grid(n,n), new(n,n)

       integer i,j,iters
       double precision maxdiff,tempdiff

c initialize boundaries of grid and new (not shown)

c initialize interior points to zeroes
!$omp parallel do
!$omp shared(n,grid,new), private(i,j)
      do j = 2,n-1
        do i = 2,n-1
          grid(i,j) = 0.0
          new(i,j) = 0,0
        enddo
      enddo
!$omp end parallel do

c initialize global variables
       iters = 1
       maxdiff = 0.0

c start main computational loop
!$omp parallel
!$omp& shared(n,maxiters,grid,new,iters)
!$omp& private(i,j,tempdiff)
!$omp& reduction(max: maxdiff)

       do while (iters.le.maxiters)

c update points in new
!$omp do
         do j = 2,n-1
           do i = 2,n-1
             new(i,j) = (grid(i-1,j) + grid(i+1,j) +
                 grid(i,j-1) + grid(i,j+1)) * 0.25
           enddo
         enddo
!$omp end do

c update points in grid
!$omp do
         do j = 2,n-1
           do i = 2,n-1
             grid(i,j) = (new(i-1,j) + new(i+1,j) +
                 new(i,j-1) + new(i,j+1)) * 0.25
           enddo
         enddo
!$omp end do

c one process updates the global iteration count
!$omp single
       iters = iters+2
!$omp end single

       enddo
c end of main computational loop

c compute maximum difference into a reduction variable

!$omp do
         do j = 2,n-1
           do i = 2,n-1
             tempdiff = abs(grid(i,j)-new(i1,j))
             maxdiff = max(maxdiff,tempdiff)
               mydiff = tempdiff
           enddo
         enddo
!$omp end parallel do

!$omp end parallel

       return
       end