libgomp/testsuite/libgomp.fortran/jacobi.f

   1 * { dg-do run }
   2
   3       program main
   4 ************************************************************
   5 * program to solve a finite difference
   6 * discretization of Helmholtz equation :
   7 * (d2/dx2)u + (d2/dy2)u - alpha u = f
   8 * using Jacobi iterative method.
   9 *
  10 * Modified: Sanjiv Shah,       Kuck and Associates, Inc. (KAI), 1998
  11 * Author:   Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998
  12 *
  13 * Directives are used in this code to achieve paralleism.
  14 * All do loops are parallelized with default 'static' scheduling.
  15 *
  16 * Input :  n - grid dimension in x direction
  17 *          m - grid dimension in y direction
  18 *          alpha - Helmholtz constant (always greater than 0.0)
  19 *          tol   - error tolerance for iterative solver
  20 *          relax - Successice over relaxation parameter
  21 *          mits  - Maximum iterations for iterative solver
  22 *
  23 * On output
  24 *       : u(n,m) - Dependent variable (solutions)
  25 *       : f(n,m) - Right hand side function
  26 *************************************************************
  27       implicit none
  28
  29       integer n,m,mits,mtemp
  30       include "omp_lib.h"
  31       double precision tol,relax,alpha
  32
  33       common /idat/ n,m,mits,mtemp
  34       common /fdat/tol,alpha,relax
  35 *
  36 * Read info
  37 *
  38       write(*,*) "Input n,m - grid dimension in x,y direction "
  39       n = 64
  40       m = 64
  41 *     read(5,*) n,m
  42       write(*,*) n, m
  43       write(*,*) "Input alpha - Helmholts constant "
  44       alpha = 0.5
  45 *     read(5,*) alpha
  46       write(*,*) alpha
  47       write(*,*) "Input relax - Successive over-relaxation parameter"
  48       relax = 0.9
  49 *     read(5,*) relax
  50       write(*,*) relax
  51       write(*,*) "Input tol - error tolerance for iterative solver"
  52       tol = 1.0E-12
  53 *     read(5,*) tol
  54       write(*,*) tol
  55       write(*,*) "Input mits - Maximum iterations for solver"
  56       mits = 100
  57 *     read(5,*) mits
  58       write(*,*) mits
  59
  60       call omp_set_num_threads (2)
  61
  62 *
  63 * Calls a driver routine
  64 *
  65       call driver ()
  66
  67       stop
  68       end
  69
  70       subroutine driver ( )
  71 *************************************************************
  72 * Subroutine driver ()
  73 * This is where the arrays are allocated and initialized.
  74 *
  75 * Working varaibles/arrays
  76 *     dx  - grid spacing in x direction
  77 *     dy  - grid spacing in y direction
  78 *************************************************************
  79       implicit none
  80
  81       integer n,m,mits,mtemp
  82       double precision tol,relax,alpha
  83
  84       common /idat/ n,m,mits,mtemp
  85       common /fdat/tol,alpha,relax
  86
  87       double precision u(n,m),f(n,m),dx,dy
  88
  89 * Initialize data
  90
  91       call initialize (n,m,alpha,dx,dy,u,f)
  92
  93 * Solve Helmholtz equation
  94
  95       call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits)
  96
  97 * Check error between exact solution
  98
  99       call  error_check (n,m,alpha,dx,dy,u,f)
 100
 101       return
 102       end
 103
 104       subroutine initialize (n,m,alpha,dx,dy,u,f)
 105 ******************************************************
 106 * Initializes data
 107 * Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2)
 108 *
 109 ******************************************************
 110       implicit none
 111
 112       integer n,m
 113       double precision u(n,m),f(n,m),dx,dy,alpha
 114
 115       integer i,j, xx,yy
 116       double precision PI
 117       parameter (PI=3.1415926)
 118
 119       dx = 2.0 / (n-1)
 120       dy = 2.0 / (m-1)
 121
 122 * Initialize initial condition and RHS
 123
 124 !$omp parallel do private(xx,yy)
 125       do j = 1,m
 126          do i = 1,n
 127             xx = -1.0 + dx * dble(i-1)        ! -1 < x < 1
 128             yy = -1.0 + dy * dble(j-1)        ! -1 < y < 1
 129             u(i,j) = 0.0
 130             f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy)
 131      &           - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy)
 132          enddo
 133       enddo
 134 !$omp end parallel do
 135
 136       return
 137       end
 138
 139       subroutine jacobi (n,m,dx,dy,alpha,omega,u,f,tol,maxit)
 140 ******************************************************************
 141 * Subroutine HelmholtzJ
 142 * Solves poisson equation on rectangular grid assuming :
 143 * (1) Uniform discretization in each direction, and
 144 * (2) Dirichlect boundary conditions
 145 *
 146 * Jacobi method is used in this routine
 147 *
 148 * Input : n,m   Number of grid points in the X/Y directions
 149 *         dx,dy Grid spacing in the X/Y directions
 150 *         alpha Helmholtz eqn. coefficient
 151 *         omega Relaxation factor
 152 *         f(n,m) Right hand side function
 153 *         u(n,m) Dependent variable/Solution
 154 *         tol    Tolerance for iterative solver
 155 *         maxit  Maximum number of iterations
 156 *
 157 * Output : u(n,m) - Solution
 158 *****************************************************************
 159       implicit none
 160       integer n,m,maxit
 161       double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega
 162 *
 163 * Local variables
 164 *
 165       integer i,j,k,k_local
 166       double precision error,resid,rsum,ax,ay,b
 167       double precision error_local, uold(n,m)
 168
 169       real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2
 170       real te1,te2
 171       real second
 172       external second
 173 *
 174 * Initialize coefficients
 175       ax = 1.0/(dx*dx) ! X-direction coef
 176       ay = 1.0/(dy*dy) ! Y-direction coef
 177       b  = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff
 178
 179       error = 10.0 * tol
 180       k = 1
 181
 182       do while (k.le.maxit .and. error.gt. tol)
 183
 184          error = 0.0
 185
 186 * Copy new solution into old
 187 !$omp parallel
 188
 189 !$omp do
 190          do j=1,m
 191             do i=1,n
 192                uold(i,j) = u(i,j)
 193             enddo
 194          enddo
 195
 196 * Compute stencil, residual, & update
 197
 198 !$omp do private(resid) reduction(+:error)
 199          do j = 2,m-1
 200             do i = 2,n-1
 201 *     Evaluate residual
 202                resid = (ax*(uold(i-1,j) + uold(i+1,j))
 203      &                + ay*(uold(i,j-1) + uold(i,j+1))
 204      &                 + b * uold(i,j) - f(i,j))/b
 205 * Update solution
 206                u(i,j) = uold(i,j) - omega * resid
 207 * Accumulate residual error
 208                error = error + resid*resid
 209             end do
 210          enddo
 211 !$omp enddo nowait
 212
 213 !$omp end parallel
 214
 215 * Error check
 216
 217          k = k + 1
 218
 219          error = sqrt(error)/dble(n*m)
 220 *
 221       enddo                     ! End iteration loop
 222 *
 223       print *, 'Total Number of Iterations ', k
 224       print *, 'Residual                   ', error
 225
 226       return
 227       end
 228
 229       subroutine error_check (n,m,alpha,dx,dy,u,f)
 230       implicit none
 231 ************************************************************
 232 * Checks error between numerical and exact solution
 233 *
 234 ************************************************************
 235
 236       integer n,m
 237       double precision u(n,m),f(n,m),dx,dy,alpha
 238
 239       integer i,j
 240       double precision xx,yy,temp,error
 241
 242       dx = 2.0 / (n-1)
 243       dy = 2.0 / (m-1)
 244       error = 0.0
 245
 246 !$omp parallel do private(xx,yy,temp) reduction(+:error)
 247       do j = 1,m
 248          do i = 1,n
 249             xx = -1.0d0 + dx * dble(i-1)
 250             yy = -1.0d0 + dy * dble(j-1)
 251             temp  = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy)
 252             error = error + temp*temp
 253          enddo
 254       enddo
 255
 256       error = sqrt(error)/dble(n*m)
 257
 258       print *, 'Solution Error : ',error
 259
 260       return
 261       end