Mercurial > hg > CbC > CbC_gcc

comparison libgomp/testsuite/libgomp.fortran/jacobi.f @ 0:a06113de4d67

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children 1830386684a0
comparison
equal deleted inserted replaced
-1:000000000000 0:a06113de4d67
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