|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
|
|
High Performance Scientific Computing
AMath 483/583 Class Notes Spring Quarter, 2013 |
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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 | ! $UWHPSC/codes/openmp/jacobi1d_omp1.f90
!
! Jacobi iteration illustrating fine grain parallelism with OpenMP.
!
! Several omp parallel do loops are used. Each time threads will be
! forked and the compiler will decide how to split up the loop.
program jacobi1d_omp1
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
real(kind=8) :: t1,t2
integer :: i,iter,maxiter
! Specify number of threads to use:
nthreads = 2
!$ call omp_set_num_threads(nthreads)
!$ print "('Using OpenMP with ',i3,' threads')", nthreads
print *, "Input n ... "
read *, 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")
call cpu_time(t1)
! 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
call cpu_time(t2)
print '("CPU time = ",f12.8, " seconds")', t2-t1
print *, "Total number of iterations: ",iter
write(20,*) " x u"
do i=0,n+1
write(20,'(2e20.10)'), x(i), u(i)
enddo
print *, "Solution is in heatsoln.txt"
close(20)
end program jacobi1d_omp1
|