High Performance Scientific Computing   AMath 483/583 Class Notes   Spring Quarter, 2011

#### Previous topic

Rotating particles and Fortran efficiency

#### Next topic

Jacobi iteration using OpenMP with coarse-grain parallel block

# Jacobi iteration using OpenMP with parallel do constructsΒΆ

The code below implements Jacobi iteration for solving the linear system arising from the steady state heat equation with a simple application of parallel do loops using OpenMP.

Compare to:

The code:

```! \$CLASSHG/codes/openmp/jacobi1.f90

program jacobi1
use omp_lib
implicit none
integer :: n, nthreads
real(kind=8), dimension(:), allocatable :: x,u,uold,f
real(kind=8) :: alpha, beta, dx, tol, dumax
real(kind=8), intrinsic :: exp
integer :: i,iter,maxiter

! Specify number of threads to use:
!\$ print "('Using OpenMP with ',i3,' threads')", nthreads

print *, "Input n ... "

! allocate storage for boundary points too:
allocate(x(0:n+1), u(0:n+1), uold(0:n+1), f(0:n+1))

open(unit=20, file="heatsoln.txt", status="unknown")

! grid spacing:
dx = 1.d0 / (n+1.d0)

! boundary conditions:
alpha = 20.d0
beta = 60.d0

!\$omp parallel do
do i=0,n+1
! grid points:
x(i) = i*dx
! source term:
f(i) = 100.*exp(x(i))
! initial guess:
u(i) = alpha + x(i)*(beta-alpha)
enddo

! tolerance and max number of iterations:
tol = 0.1 * dx**2
print *, "Convergence tolerance: tol = ",tol
maxiter = 100000
print *, "Maximum number of iterations: maxiter = ",maxiter

! Jacobi iteratation:

uold = u  ! starting values before updating

do iter=1,maxiter
dumax = 0.d0
!\$omp parallel do reduction(max : dumax)
do i=1,n
u(i) = 0.5d0*(uold(i-1) + uold(i+1) + dx**2*f(i))
dumax = max(dumax, abs(u(i)-uold(i)))
enddo
if (mod(iter,10000)==0) then
print *, iter, dumax
endif
! check for convergence:
if (dumax .lt. tol) exit

!\$omp parallel do
do i=1,n
uold(i) = u(i)   ! for next iteration
enddo
enddo

print *, "Total number of iterations: ",iter

write(20,*) "          x                  u"
do i=0,n+1
write(20,222), x(i), u(i)
enddo
222 format(2e20.10)

print *, "Solution is in heatsoln.txt"

close(20)

end program jacobi1
```