summaryrefslogtreecommitdiff
path: root/libgomp/testsuite/libgomp.fortran/jacobi.f
diff options
context:
space:
mode:
Diffstat (limited to 'libgomp/testsuite/libgomp.fortran/jacobi.f')
-rw-r--r--libgomp/testsuite/libgomp.fortran/jacobi.f261
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