UW AMath High Performance Scientific Computing
 
Coursera Edition

Previous topic

Numerical methods for the Poisson problem

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:

 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