diff options
Diffstat (limited to 'libgomp/testsuite/libgomp.fortran/jacobi.f')
| -rw-r--r-- | libgomp/testsuite/libgomp.fortran/jacobi.f | 261 |
1 files changed, 261 insertions, 0 deletions
diff --git a/libgomp/testsuite/libgomp.fortran/jacobi.f b/libgomp/testsuite/libgomp.fortran/jacobi.f new file mode 100644 index 000000000..b27e20f27 --- /dev/null +++ b/libgomp/testsuite/libgomp.fortran/jacobi.f @@ -0,0 +1,261 @@ +* { dg-do run } + + program main +************************************************************ +* program to solve a finite difference +* discretization of Helmholtz equation : +* (d2/dx2)u + (d2/dy2)u - alpha u = f +* using Jacobi iterative method. +* +* Modified: Sanjiv Shah, Kuck and Associates, Inc. (KAI), 1998 +* Author: Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998 +* +* Directives are used in this code to achieve paralleism. +* All do loops are parallized with default 'static' scheduling. +* +* Input : n - grid dimension in x direction +* m - grid dimension in y direction +* alpha - Helmholtz constant (always greater than 0.0) +* tol - error tolerance for iterative solver +* relax - Successice over relaxation parameter +* mits - Maximum iterations for iterative solver +* +* On output +* : u(n,m) - Dependent variable (solutions) +* : f(n,m) - Right hand side function +************************************************************* + implicit none + + integer n,m,mits,mtemp + include "omp_lib.h" + double precision tol,relax,alpha + + common /idat/ n,m,mits,mtemp + common /fdat/tol,alpha,relax +* +* Read info +* + write(*,*) "Input n,m - grid dimension in x,y direction " + n = 64 + m = 64 +* read(5,*) n,m + write(*,*) n, m + write(*,*) "Input alpha - Helmholts constant " + alpha = 0.5 +* read(5,*) alpha + write(*,*) alpha + write(*,*) "Input relax - Successive over-relaxation parameter" + relax = 0.9 +* read(5,*) relax + write(*,*) relax + write(*,*) "Input tol - error tolerance for iterative solver" + tol = 1.0E-12 +* read(5,*) tol + write(*,*) tol + write(*,*) "Input mits - Maximum iterations for solver" + mits = 100 +* read(5,*) mits + write(*,*) mits + + call omp_set_num_threads (2) + +* +* Calls a driver routine +* + call driver () + + stop + end + + subroutine driver ( ) +************************************************************* +* Subroutine driver () +* This is where the arrays are allocated and initialzed. +* +* Working varaibles/arrays +* dx - grid spacing in x direction +* dy - grid spacing in y direction +************************************************************* + implicit none + + integer n,m,mits,mtemp + double precision tol,relax,alpha + + common /idat/ n,m,mits,mtemp + common /fdat/tol,alpha,relax + + double precision u(n,m),f(n,m),dx,dy + +* Initialize data + + call initialize (n,m,alpha,dx,dy,u,f) + +* Solve Helmholtz equation + + call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits) + +* Check error between exact solution + + call error_check (n,m,alpha,dx,dy,u,f) + + return + end + + subroutine initialize (n,m,alpha,dx,dy,u,f) +****************************************************** +* Initializes data +* Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2) +* +****************************************************** + implicit none + + integer n,m + double precision u(n,m),f(n,m),dx,dy,alpha + + integer i,j, xx,yy + double precision PI + parameter (PI=3.1415926) + + dx = 2.0 / (n-1) + dy = 2.0 / (m-1) + +* Initilize initial condition and RHS + +!$omp parallel do private(xx,yy) + do j = 1,m + do i = 1,n + xx = -1.0 + dx * dble(i-1) ! -1 < x < 1 + yy = -1.0 + dy * dble(j-1) ! -1 < y < 1 + u(i,j) = 0.0 + f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy) + & - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy) + enddo + enddo +!$omp end parallel do + + return + end + + subroutine jacobi (n,m,dx,dy,alpha,omega,u,f,tol,maxit) +****************************************************************** +* Subroutine HelmholtzJ +* Solves poisson equation on rectangular grid assuming : +* (1) Uniform discretization in each direction, and +* (2) Dirichlect boundary conditions +* +* Jacobi method is used in this routine +* +* Input : n,m Number of grid points in the X/Y directions +* dx,dy Grid spacing in the X/Y directions +* alpha Helmholtz eqn. coefficient +* omega Relaxation factor +* f(n,m) Right hand side function +* u(n,m) Dependent variable/Solution +* tol Tolerance for iterative solver +* maxit Maximum number of iterations +* +* Output : u(n,m) - Solution +***************************************************************** + implicit none + integer n,m,maxit + double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega +* +* Local variables +* + integer i,j,k,k_local + double precision error,resid,rsum,ax,ay,b + double precision error_local, uold(n,m) + + real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2 + real te1,te2 + real second + external second +* +* Initialize coefficients + ax = 1.0/(dx*dx) ! X-direction coef + ay = 1.0/(dy*dy) ! Y-direction coef + b = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff + + error = 10.0 * tol + k = 1 + + do while (k.le.maxit .and. error.gt. tol) + + error = 0.0 + +* Copy new solution into old +!$omp parallel + +!$omp do + do j=1,m + do i=1,n + uold(i,j) = u(i,j) + enddo + enddo + +* Compute stencil, residual, & update + +!$omp do private(resid) reduction(+:error) + do j = 2,m-1 + do i = 2,n-1 +* Evaluate residual + resid = (ax*(uold(i-1,j) + uold(i+1,j)) + & + ay*(uold(i,j-1) + uold(i,j+1)) + & + b * uold(i,j) - f(i,j))/b +* Update solution + u(i,j) = uold(i,j) - omega * resid +* Accumulate residual error + error = error + resid*resid + end do + enddo +!$omp enddo nowait + +!$omp end parallel + +* Error check + + k = k + 1 + + error = sqrt(error)/dble(n*m) +* + enddo ! End iteration loop +* + print *, 'Total Number of Iterations ', k + print *, 'Residual ', error + + return + end + + subroutine error_check (n,m,alpha,dx,dy,u,f) + implicit none +************************************************************ +* Checks error between numerical and exact solution +* +************************************************************ + + integer n,m + double precision u(n,m),f(n,m),dx,dy,alpha + + integer i,j + double precision xx,yy,temp,error + + dx = 2.0 / (n-1) + dy = 2.0 / (m-1) + error = 0.0 + +!$omp parallel do private(xx,yy,temp) reduction(+:error) + do j = 1,m + do i = 1,n + xx = -1.0d0 + dx * dble(i-1) + yy = -1.0d0 + dy * dble(j-1) + temp = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy) + error = error + temp*temp + enddo + enddo + + error = sqrt(error)/dble(n*m) + + print *, 'Solution Error : ',error + + return + end |