OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-old/] [gcc-4.2.2/] [libgomp/] [testsuite/] [libgomp.fortran/] [jacobi.f] - Blame information for rev 867

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 38 julius
* { 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

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.