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

first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children 1830386684a0
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libgomp/testsuite/libgomp.fortran/jacobi.f	Fri Jul 17 14:47:48 2009 +0900
@@ -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