Parallel program (user-supplied code)

C=======================================================================
C=======================================================================
C
C       example program poisson:
C               uses Jacobi relaxation to solve the Poisson equation
C
C=======================================================================
C=======================================================================

C=======================================================================
C
C       host process-main program
C
C=======================================================================

        subroutine hostmain
        include 'mesh_uparms.h'
        include 'mesh_parms.h'
        include 'mesh_common.h'

C=======grid size 
C       moved to archetype file mesh_uparms.h
C       integer NX
C       integer NY
C       parameter (NX=10)
C       parameter (NY=10)

        real H
        parameter (H=0.05)
C       how often to check for convergence
        integer NCHECK
        parameter (NCHECK=10)
C       convergence criterion
        real TOL
        parameter (TOL=0.00001)
C       maximum number of steps
        integer MAXSTEPS
        parameter (MAXSTEPS=1000)

        external F
        external G

C=======grid variables (whole array in host)
        real uk(1:NX,1:NY), ukp1(1:NX,1:NY)

        real diff, difflocal, diffmax

C=======initialization
C       interior points
        do i = 2, NX-1
        do j = 2, NY-1
                uk(i,j) = F(i,j,NX,NY,H)
        enddo
        enddo
C       boundary points
        do j = 1, NY
                uk(1,j) = G(1,j,NX,NY,H)
                uk(NX,j) = G(NX,j,NX,NY,H)
        enddo
        do i = 2, NX-1
                uk(i,1) = G(i,1,NX,NY,H)
                uk(i,NY) = G(i,NY,NX,NY,H)
        enddo
C       move to grid
        call mesh_HtoG_host(uk)

C=======main loop (until convergence)
        diffmax = TOL + 1.0

        do k=1,MAXSTEPS

C=======grid computation (compute new values; in grid only)
C=======convergence test (global max)
        if (mod(k,NCHECK) .eq. 0) then
C               (computation in grid only)
                call mesh_merge_real_maxabs(1, difflocal, diffmax)
        endif
C=======grid computation (copy new values to old values; in grid only)

C=======convergence check
        if (diffmax .le. TOL) go to 1000

        enddo

C=======sequential output
1000    continue
        call mesh_GtoH_host(uk)
        print*, 'NX, NY = ', NX, NY
        print*, 'H = ', H
        print*, 'tolerance = ', TOL
        call Fprint
        call Gprint

        if (diffmax .le. TOL) then
                print*, 'convergence occurred in ', k, ' steps'
        else
                print*, 'no convergence in ', MAXSTEPS, ' steps'
                print*, 'no convergence in ', MAXSTEPS, ' steps; ',
     -                  'max. difference ', diffmax
        endif

        do i = 1, NX
                if (NX .gt. 10) print*, ' '
                print 9999, (uk(i,j), j = 1, NY)
        enddo
9999    format(10F8.4)
                
        end

C=======================================================================
C
C       grid process-main program
C
C=======================================================================

        subroutine gridmain
        include 'mesh_uparms.h'
        include 'mesh_parms.h'
        include 'mesh_common.h'

C=======grid size 
C       moved to archetype file mesh_uparms.h
C       integer NX
C       integer NY
C       parameter (NX=10)
C       parameter (NY=10)

        real H
        parameter (H=0.05)
C       how often to check for convergence
        integer NCHECK
        parameter (NCHECK=10)
C       convergence criterion
        real TOL
        parameter (TOL=0.00001)
C       maximum number of steps
        integer MAXSTEPS
        parameter (MAXSTEPS=1000)

        external F
        external G

C=======grid variables (local sections)
        real uk(IXLO:IXHI,IYLO:IYHI), ukp1(1:NXlsize,1:NYlsize)

        real diff, difflocal, diffmax

        logical iempty, jempty

C=======initialization (in host only)
        call mesh_HtoG_grid(uk)

C=======main loop (until convergence)

C       compute loop bounds
        call xintersect(2, NX-1, istart, iend, iempty)
        call yintersect(2, NY-1, jstart, jend, jempty)

        diffmax = TOL + 1.0

        do k=1,MAXSTEPS

C=======grid computation (compute new values)
        call mesh_update_bdry(uk)
        do i = istart, iend
        do j = jstart, jend
                ukp1(i,j) = 0.25*(H*H*F(i,j,NX,NY,H) 
     -                  + uk(i-1,j) + uk(i,j-1) 
     -                  + uk(i+1,j) + uk(i,j+1) )
        enddo
        enddo
C=======convergence test (global max)
        if (mod(k,NCHECK) .eq. 0) then
                difflocal = 0.0
                do i = istart, iend
                do j = jstart, jend
                        diff = abs(ukp1(i,j) - uk(i,j))
                        if (diff .gt. difflocal) difflocal = diff
                enddo
                enddo
                call mesh_merge_real_maxabs(1, difflocal, diffmax)
        endif
C=======grid computation (copy new values to old values)
        do i = istart, iend
        do j = jstart, jend
                uk(i,j) = ukp1(i,j)
        enddo
        enddo

C=======convergence check
        if (diffmax .le. TOL) go to 1000

        enddo

C=======sequential output (I/O in host only)
1000    continue
        call mesh_GtoH_grid(uk)

        end

C=======================================================================

        real function F(i,j,nx,ny,h)
        integer i, j, nx, ny
        real h

        F = 0.0
        end

        subroutine Fprint
        print*, 'F(i,j) = 0.0'
        end

C=======================================================================

        real function G(i,j,nx,ny,h)
        integer i, j, ny, ny
        real h

        G = (i+j)*h
        end

        subroutine Gprint
        print*, 'G(i,j) = (i+j)*H'
        end

C=======================================================================