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 parallized 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 initialzed.
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* Initilize 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
262