UW AMath 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

This Page

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:
   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")

    ! 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