* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: poisson.cpp,v 1.2 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: poisson.cpp,v $
// Revision 1.2  1995/05/16  16:15:09  dweber
// Final revision for project completion.  Notes:
// 	#ifdef CMF is back to #ifdef CMFORTRAN
//      No equivalent to CMF_descriptors.  Am looking into pointers.
//      Method #6 (and 16) appear to produce incorrect results.
//
// Revision 1.1  1995/04/18  01:02:09  dweber
// Initial revision
//
# Revision 1.1  1995/02/11  20:23:32  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      program main
*
***********************************************************************
*
      integer n
      
c     get the command line arguments
      call parscl()

c     main program loop: read from input.par and execute poisson
      open(unit=27,file='input.par',status='old')
 50   read(27,*,END=100) n
      write(*,*) "Data set: ", n
      call poisson(n,n,n)
      goto 50
100   continue
      
      stop
      end




*
***********************************************************************
*
      subroutine poisson(idim,jdim,kdim)
*
***********************************************************************
*
*
*     This is a sample problem to be used to test the CMStab package.
* It sets up the necessary arrays to solve the poisson equation on the
* unit cube as follows:
*
*       ((d/dx)^2 + (d/dy)^2 + (d/dz)^2) phi = 100.
*
* With the boundary conditions:
*
*       phi(x=0.) = 0.
*       phi(x=1.) = 0.
*       phi(y=0.) = 0.
*       phi(y=1.) = 0.
*       phi(z=0.) = 0.
*       phi(z=1.) = 0.
*
***********************************************************************
*
*     include 'poisson.h'
#include "driver.h"
      integer idim,jdim,kdim
#ifdef CMFORTRAN
      real, array(idim,jdim,kdim) :: phi,rhs,ac,an,as,ae,aw,at,ab
#elif HPF
      real, DIMENSION(idim,jdim,kdim) :: phi,rhs,ac,an,as,ae,aw,at,ab
#endif
      integer i,j,k
      real dx,dy,dz
      character dfile*20
*
***********************************************************************
*
*     Finite differenced equation reads:
*
*       (phi(i+1,j,k) + phi(i-1,j,k)) * dy^2 * dz^2
*     + (phi(i,j+1,k) + phi(i,j-1,k)) * dx^2 * dz^2
*     + (phi(i,j,k+1) + phi(i,j,k-1)) * dx^2 * dy^2
*     - 2 * (dy^2 * dz^2 + dx^2 * dz^2 + dx^2 * dy^2) * phi(i,j) 
*                                          = 100 * (dx^2 * dy^2 * dz^2)
*
*     Let i reference the east-west direction of grid and correspond
* to planes of constant x; let j reference the north-south direction
* and correspond to planes of constant y; let k reference the
* top-bottom direction and correspond to planes of constant z.  Then
* the difference coefficients are:
*
      dx = 1./float(idim-1)
      dy = 1./float(jdim-1)
      dz = 1./float(kdim-1)

      ac=-2.*( dx*dx * dz*dz
     .       + dy*dy * dz*dz
     .       + dx*dx * dy*dy)
      an=dx*dx * dz*dz
      as=dx*dx * dz*dz
      ae=dy*dy * dz*dz
      aw=dy*dy * dz*dz
      at=dx*dx * dy*dy
      ab=dx*dx * dy*dy

      rhs=100.*dx*dx*dy*dy*dz*dz
*
* Now apply Boundary conditions
*

* Do faces

* on x=dx plane
      i = 2
      aw(i,3:jdim-2,3:kdim-2) = 0.

* on x=1.-dx plane
      i = idim-1
      ae(i,3:jdim-2,3:kdim-2) = 0.

* on y=dy plane
      j = 2
      as(3:idim-2,j,3:kdim-2) = 0.

* on y=1.-dy plane
      j = jdim-1
      an(3:idim-2,j,3:kdim-2) = 0.

* on z=dz plane
      k = 2
      ab(3:idim-2,3:jdim-2,k) = 0.

* on z=1.-dz plane
      k = kdim-1
      at(3:idim-2,3:jdim-2,k) = 0.

* Do edges

* along x=dx, y=dy
      i = 2
      j = 2
      aw(i,j,3:kdim-2) = 0.
      as(i,j,3:kdim-2) = 0.

* along x=dx y=1.-dy
      i = 2
      j = jdim-1
      aw(i,j,3:kdim-2) = 0.
      an(i,j,3:kdim-2) = 0.

* along x=1.-dx, y=1.-dy
      i = idim-1
      j = jdim-1
      ae(i,j,3:kdim-2) = 0.
      an(i,j,3:kdim-2) = 0.

* along x=1.-dx, y=dy
      i = idim-1
      j = 2
      ae(i,j,3:kdim-2) = 0.
      as(i,j,3:kdim-2) = 0.


* along y=dy, z=dz
      j = 2
      k = 2
      as(3:idim-2,j,k) = 0.
      ab(3:idim-2,j,k) = 0.

* along y=dy, z=1.-dz
      j = 2
      k = kdim-1
      as(3:idim-2,j,k) = 0.
      at(3:idim-2,j,k) = 0.

* along y=1.-dy, z=1.-dz
      j = jdim-1
      k = kdim-1
      an(3:idim-2,j,k) = 0.
      at(3:idim-2,j,k) = 0.

* along y=1.-dy, z=dz
      j = jdim-1
      k = 2
      an(3:idim-2,j,k) = 0.
      ab(3:idim-2,j,k) = 0.

* along z=dz, x=dx
      k = 2
      i = 2
      ab(i,3:jdim-2,k) = 0.
      aw(i,3:jdim-2,k) = 0.

* along z=dz, x=1.-dx
      k = 2
      i = idim-1
      ab(i,3:jdim-2,k) = 0.
      ae(i,3:jdim-2,k) = 0.

* along z=1.-dz, x=1.-dx
      k = kdim-1
      i = idim-1
      at(i,3:jdim-2,k) = 0.
      ae(i,3:jdim-2,k) = 0.

* along z=1.-dz, x=dx
      k = kdim-1
      i = 2
      at(i,3:jdim-2,k) = 0.
      aw(i,3:jdim-2,k) = 0.

* Do verticies

* at x=dx, y=dy, z=dz
      i = 2
      j = 2
      k = 2
      aw(i,j,k) = 0.
      as(i,j,k) = 0.
      ab(i,j,k) = 0.

* at x=dx, y=dy, z=1.-dz
      i = 2
      j = 2
      k = kdim-1
      aw(i,j,k) = 0.
      as(i,j,k) = 0.
      at(i,j,k) = 0.

* at x=dx, y=1.-dy, z=dz
      i = 2
      j = jdim-1
      k = 2
      aw(i,j,k) = 0.
      an(i,j,k) = 0.
      ab(i,j,k) = 0.

* at x=dx, y=1.-dy, z=1.-dz
      i = 2
      j = jdim-1
      k = kdim-1
      aw(i,j,k) = 0.
      an(i,j,k) = 0.
      at(i,j,k) = 0.

* at x=1.-dx, y=dy, z=dz
      i = idim-1
      j = 2
      k = 2
      ae(i,j,k) = 0.
      as(i,j,k) = 0.
      ab(i,j,k) = 0.

* at x=1.-dz, y=dy, z=1.-dz
      i = idim-1
      j = 2
      k = kdim-1
      ae(i,j,k) = 0.
      as(i,j,k) = 0.
      at(i,j,k) = 0.

* at x=1.-dx, y=1.-dy, z=dz
      i = idim-1
      j = jdim-1
      k = 2
      ae(i,j,k) = 0.
      an(i,j,k) = 0.
      ab(i,j,k) = 0.

* at x=1.-dx, y=1.-dy, z=1.-dz
      i = idim-1
      j = jdim-1
      k = kdim-1
      ae(i,j,k) = 0.
      an(i,j,k) = 0.
      at(i,j,k) = 0.

* Set initial guess as phi=0. everywhere
      phi=0.

*
* Call the a setup routine for the CMStab solver
*
      call driver (ac,an,as,ae,aw,at,ab,phi,rhs,
     .             idim,2,idim-1,jdim,2,jdim-1,jdim,2,jdim-1)

* Output a line of data to a datafile.  Give it a real name.
#ifdef CMFORTRAN
      call datname(dfile, method, iprecon, nodes/4, idim)
#elif HPF
      call datname(dfile, method, iprecon, nodes, idim)
#endif
      open(unit=50, file=dfile, status="UNKNOWN", position="APPEND")
      do i=1,idim
        write(50,*) (i-1)*dx,phi(i,jdim/2,kdim/2)
      enddo
      close(50)

      return
      end
*
***********************************************************************
*



*
***********************************************************************
*
      subroutine parscl()
*
***********************************************************************
*     This subroutine determines if any arguments are to be parsed
* from the commandline.  The arguments currently accepted are:
*
*     -method n     : The method number to test (1-6) (default 1)
*     -preconl
*     -preconr      : Preconditioning (default is none)
*     -iters n      : Number of iterations to test each method
*     -help         : Print out this type of help message
*
*     The variables declarations for these are found in driver.h
*
***********************************************************************
*
      implicit none
#include "driver.h"
      character argument*80, val*80
      integer i
      integer IARGC

*     Default to no preconditioning
      iprecon = 0

*     Default to method 1
      method = 1

*     Default testing iterations to 3
      tstitr = 3

      do 100 i = 1, iargc(), 1
         call getarg(i, argument)

         if (index(argument, "-method") .ne. 0) then
            call getarg(i+1, val)
            read(val, *) method
            if ((method .lt. 1) .or. (method .gt. 6)) then
               write(*,*) "Invalid method number ", method
               call usage()
            endif

         elseif (index(argument, "-preconl") .ne. 0) then
            iprecon = 1

         elseif (index(argument, "-preconr") .ne. 0) then
            iprecon = 2

         elseif (index(argument, "-precons") .ne. 0) then
            iprecon = 3

         elseif (index(argument, "-iters") .ne. 0) then
            call getarg(i+1, val)
            read(val, *) tstitr
            if (tstitr .lt. 1) then
               write(*,*) "Invalid # of interations: ", tstitr
               call usage()
            endif

         elseif (index(argument, "-help") .ne. 0) then
            call usage()

         endif
               
 100  continue

      return
      end


*
***********************************************************************
*
*     Usage:
*     Prints out the correct usage of the poisson program
*
      subroutine usage
      implicit none

      write(*,*) "Usage: poisson [options]"
      write(*,*) "Options:"
      write(*,*) "-help       : Generate this message"
      write(*,*) "-method n   : Choose method #n (1 <= n <= 6)"
      write(*,*) "-iters n    : N test iterations of solver (n > 0)"
      write(*,*) "-preconl    : Left preconditioning"
      write(*,*) "-preconr    : Right preconditioning"
      write(*,*) "-precons    : Split preconditioning"
      write(*,*)

      stop
      end


*
***********************************************************************
*
*     datname
*     Returns a filename for outputting a slice of the data.  It is
* now a unique file name so comparisons can be made between the
* methods and versions, both on the CM-5 and SGI.
*
      subroutine datname(dfile, method, iprecon, nodes, idim)
      implicit none
      character dfile*(*)
      integer method, iprecon, nodes, idim
      integer pcmeth
      character meth*2, nds*3, dims*3

      pcmeth = 10 * iprecon + method
      if (pcmeth .lt. 10) then
         write(meth, 'A1,I1') "0", pcmeth
      else
         write(meth, 'I2') pcmeth
      endif

      if (nodes .lt. 10) then
         write(nds, 'A2,I1') "00", nodes
      elseif (nodes .lt. 100) then
         write(nds, 'A1,I2') "0", nodes
      else
         write(nds, 'I3') nodes
      endif

      if (idim .lt. 10) then
         write(dims, 'A2,I1') "00", idim
      elseif (idim .lt. 100) then
         write(dims, 'A1,I2') "0", idim
      else
         write(dims, 'I3') idim
      endif
      write(dfile, 'A5,A3,A1,A2,A1,A3') "data.",nds,".",meth,".",dims

      return
      end
