program main
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
write(*,*) "Input n,m - grid dimension in x,y direction "
n = 64
m = 64
write(*,*) n, m
write(*,*) "Input alpha - Helmholts constant "
alpha = 0.5
write(*,*) alpha
write(*,*) "Input relax - Successive over-relaxation parameter"
relax = 0.9
write(*,*) relax
write(*,*) "Input tol - error tolerance for iterative solver"
tol = 1.0E-12
write(*,*) tol
write(*,*) "Input mits - Maximum iterations for solver"
mits = 100
write(*,*) mits
call omp_set_num_threads (2)
call driver ()
stop
end
subroutine driver ( )
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
call initialize (n,m,alpha,dx,dy,u,f)
call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits)
call error_check (n,m,alpha,dx,dy,u,f)
return
end
subroutine initialize (n,m,alpha,dx,dy,u,f)
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)
!$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)
implicit none
integer n,m,maxit
double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega
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
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
!$omp parallel
!$omp do
do j=1,m
do i=1,n
uold(i,j) = u(i,j)
enddo
enddo
!$omp do private(resid) reduction(+:error)
do j = 2,m-1
do i = 2,n-1
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
u(i,j) = uold(i,j) - omega * resid
error = error + resid*resid
end do
enddo
!$omp enddo nowait
!$omp end parallel
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
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