Subversion Repositories openrisc

* { 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