* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: usrsubs.cpp,v 1.7 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: usrsubs.cpp,v $
# Revision 1.7  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.5  1995/04/18  01:02:09  dweber
# Isolated CMF & HPF specific features with #ifdef statements.
# Tried to update all HPF directives to be HPF specific - not CMF holdovers.
# Deleted the CMF select(var) - case(val) statements, and am using
#         if-then-elseif statements.
# Created architecture types to match the architecture.  This will need
#         some work as HPF has no concept similar to the CMF_architecture
#         functions.
# Since HPF doesn't define a type like CMF_sizeof_descriptor, I had to
#         hardcode the array layouts for HPF.  This may change.
#
# Revision 1.4  1995/03/03  05:57:17  jtowns
# Enable use of matirx transpose for matvec and preconditioning.  This
# is quite inefficient currently since we must generate the transpose
# evertim each routine is called since we cannot use a SAVE statement.
# This needs to be fixed!!
#
# Revision 1.3  1995/02/27  23:18:13  jtowns
# Missed changing the #include in usrmv1.  Fixed to include correct
# header.
#
# Revision 1.2  1995/02/22  22:16:59  jtowns
# Fix includes to be directoives handled by new Makefile with
# preprocessing.
#
# Revision 1.1  1995/02/11  20:23:38  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      subroutine usrmv (job,a,ia,x,y)
*
***********************************************************************
*
#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif
      
      integer job, ia(3,1,1,1), n
#ifdef CMFORTRAN
cmf$ layout ia(:serial,:news,:news,:news)
#endif
#ifdef HPF
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      integer a(CMF_sizeof_descriptor),x(CMF_sizeof_descriptor),
     .     y(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout x(:serial)
cmf$ layout y(:serial)
#elif HPF
      real a(:,:,:,:), x(:,:,:,:), y(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE y(*, BLOCK, BLOCK, BLOCK)
#endif

*
***********************************************************************
*
      n= ia(1,1,1,1)

      call usrmv1(job,a,n,x,y)

      return
      end
*
***********************************************************************
*
      subroutine usrmv1 (job,a,n,x,y)
*
***********************************************************************
*     include 'poisson.h'

#ifdef CMFORTRAN
*  Use pshift CMSSL routine for commnications
#include "cm/cmssl-cmf.h"
#endif

      integer job, n

#ifdef CMFORTRAN
      real, array(7,n,n,n) :: a
      real, array(1,n,n,n) :: x,y
cmf$ layout a(:serial,:news,:news,:news)
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout y(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(7,n,n,n) :: a
      real, DIMENSION(1,n,n,n) :: x,y
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE y(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      real, array(1,n,n,n) :: xn,xs,xe,xw,xt,xb
cmf$ layout xn(:serial,:news,:news,:news)
cmf$ layout xs(:serial,:news,:news,:news)
cmf$ layout xe(:serial,:news,:news,:news)
cmf$ layout xw(:serial,:news,:news,:news)
cmf$ layout xt(:serial,:news,:news,:news)
cmf$ layout xb(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(1,n,n,n) :: xn,xs,xe,xw,xt,xb
CHPF$ DISTRIBUTE xn(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xs(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xe(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xw(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xt(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xb(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      real, array(7,n,n,n) :: aT
cmf$ layout aT(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(7,n,n,n) :: aT
CHPF$ DISTRIBUTE aT(*, BLOCK, BLOCK, BLOCK)
#endif

      integer setup,ier
      save setup

      integer ntest,itmax,iters
      common /ints/ ntest,itmax,iters
*
***********************************************************************
*

* On first call set up
#ifdef CMFORTRAN
      if (iters.eq.0) then
         setup = pshift_setup (6,x,ier,
     .        CMSSL_CSHIFT,2, 1,
     .        CMSSL_CSHIFT,2,-1,
     .        CMSSL_CSHIFT,3, 1,
     .        CMSSL_CSHIFT,3,-1,
     .        CMSSL_CSHIFT,4, 1,
     .        CMSSL_CSHIFT,4,-1)
      endif
#endif

      if (job .eq. 0) then
#ifdef CMFORTRAN
         call pshift(6,setup,ier,
     .        CMSSL_CSHIFT,xe,x,2, 1,
     .        CMSSL_CSHIFT,xw,x,2,-1,
     .        CMSSL_CSHIFT,xn,x,3, 1,
     .        CMSSL_CSHIFT,xs,x,3,-1,
     .        CMSSL_CSHIFT,xt,x,4, 1,
     .        CMSSL_CSHIFT,xb,x,4,-1)
#elif HPF
         xe(:,:,:,:) = cshift(x(:,:,:,:), 1, 2)
         xw(:,:,:,:) = cshift(x(:,:,:,:), -1, 2)
         xn(:,:,:,:) = cshift(x(:,:,:,:), 1, 3)
         xs(:,:,:,:) = cshift(x(:,:,:,:), -1, 3)
         xt(:,:,:,:) = cshift(x(:,:,:,:), 1, 4)
         xb(:,:,:,:) = cshift(x(:,:,:,:), -1, 4)
#endif
         y(1,:,:,:) = x(1,:,:,:) * a(1,:,:,:)
     .        +xn(1,:,:,:) * a(2,:,:,:)
     .        +xs(1,:,:,:) * a(3,:,:,:)
     .        +xe(1,:,:,:) * a(4,:,:,:)
     .        +xw(1,:,:,:) * a(5,:,:,:)
     .        +xt(1,:,:,:) * a(6,:,:,:)
     .        +xb(1,:,:,:) * a(7,:,:,:)

      elseif (job .eq. 1) then
*         if (iters.le.1) then

         aT(1,:,:,:)= a(1,:,:,:)
#ifdef CMFORTRAN
*        Syntax for CMF: eoshift(array, dim, shift [, boundary])
         aT(2,:,:,:) = eoshift(a(2,:,:,:), 2, -1)
         aT(3,:,:,:) = eoshift(a(3,:,:,:), 2,  1)
         aT(4,:,:,:) = eoshift(a(4,:,:,:), 1, -1)
         aT(5,:,:,:) = eoshift(a(5,:,:,:), 1,  1)
         aT(6,:,:,:) = eoshift(a(6,:,:,:), 3, -1)
         aT(7,:,:,:) = eoshift(a(7,:,:,:), 3,  1)
#elif defined( HPF ) && defined( NO_EOSHIFT )
         aT(2,:,:,:) = cshift(a(2,:,:,:), -1, 2)
         aT(3,:,:,:) = cshift(a(3,:,:,:),  1, 2)
         aT(4,:,:,:) = cshift(a(4,:,:,:), -1, 1)
         aT(5,:,:,:) = cshift(a(5,:,:,:),  1, 1)
         aT(6,:,:,:) = cshift(a(6,:,:,:), -1, 3)
         aT(7,:,:,:) = cshift(a(7,:,:,:),  1, 3)
#elif defined( HPF )
*        Syntax for HPF: eoshift(array, shift [, boundary][, dim])
         aT(2,:,:,:) = eoshift(a(2,:,:,:), -1, 2)
         aT(3,:,:,:) = eoshift(a(3,:,:,:),  1, 2)
         aT(4,:,:,:) = eoshift(a(4,:,:,:), -1, 1)
         aT(5,:,:,:) = eoshift(a(5,:,:,:),  1, 1)
         aT(6,:,:,:) = eoshift(a(6,:,:,:), -1, 3)
         aT(7,:,:,:) = eoshift(a(7,:,:,:),  1, 3)
#endif

*         endif
#ifdef CMFORTRAN
         call pshift(6,setup,ier,
     .        CMSSL_CSHIFT,xe,x,2, 1,
     .        CMSSL_CSHIFT,xw,x,2,-1,
     .        CMSSL_CSHIFT,xn,x,3, 1,
     .        CMSSL_CSHIFT,xs,x,3,-1,
     .        CMSSL_CSHIFT,xt,x,4, 1,
     .        CMSSL_CSHIFT,xb,x,4,-1)
#elif HPF
         xe(:,:,:,:) = cshift(x(:,:,:,:), 1, 2)
         xw(:,:,:,:) = cshift(x(:,:,:,:), -1, 2)
         xn(:,:,:,:) = cshift(x(:,:,:,:), 1, 3)
         xs(:,:,:,:) = cshift(x(:,:,:,:), -1, 3)
         xt(:,:,:,:) = cshift(x(:,:,:,:), 1, 4)
         xb(:,:,:,:) = cshift(x(:,:,:,:), -1, 4)
#endif


         y(1,:,:,:) = x(1,:,:,:) * aT(1,:,:,:)
     .        +xn(1,:,:,:) * aT(3,:,:,:)
     .        +xs(1,:,:,:) * aT(2,:,:,:)
     .        +xe(1,:,:,:) * aT(5,:,:,:)
     .        +xw(1,:,:,:) * aT(4,:,:,:)
     .        +xt(1,:,:,:) * aT(7,:,:,:)
     .        +xb(1,:,:,:) * aT(6,:,:,:)


      else
         print*, 'usrmv called with job = ',job
         print*, 'which is not supported by usrmv'

      endif

      return
      end
*
***********************************************************************
*
      subroutine neumann_left (job,q,iq,x,y)
*
***********************************************************************
*
#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer job, n
#ifdef CMFORTRAN
      integer iq(3,1,1,1)
cmf$ layout iq(:serial,:news,:news,:news)
#elif HPF
      integer iq(:,:,:,:)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      integer q(CMF_sizeof_descriptor),x(CMF_sizeof_descriptor),
     .     y(CMF_sizeof_descriptor)
cmf$ layout q(:serial)
cmf$ layout x(:serial)
cmf$ layout y(:serial)
#elif HPF
      real q(:,:,:,:), x(:,:,:,:), y(:,:,:,:)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE y(*, BLOCK, BLOCK, BLOCK)
#endif

*
***********************************************************************
*

      n= iq(1,1,1,1)

      call neumann_left1(job,q,n,x,y)

      return
      end
*
***********************************************************************
*
      subroutine neumann_left1 (job,q,n,x,y)
*
***********************************************************************
*
*      Perform Neumann Polynomial Preconditioning by the following
* derivation (assumes that q and iq are identical to a and ia):
*
*     A = D - C
*
*     D = diagonal of A
*     C = off diagonal elements of -A    
*
*     A^-1 = (D - C)^-1 
*          = D^-1*(I - C*D^-1)^-1  
*          = D^-1*(I + C*D^-1 + C*D^-1*C*D^-1 + C*D^-1*C*D^-1*C*D^-1..)
*
*     Q = A^-1
*
*     If we also assume that the problem has been diagonally scaled,
* then D = I, and so,
*
*     Q = I + C + C*C + C*C*C + C*C*C*C + ...
*
*     Now A^-1 can be approximated by a truncation of this series.
*
*     Lets use a third order approximation to the series.
*
***********************************************************************
*
      integer m,n
      parameter (m=3)
*
***********************************************************************
*  Use grid sparse CMSSL routines
*
*    a <=> west
*    b <=> south
*    c <=> center
*    d <=> east
*    e <=> north
*    t <=> top
*    b <=> bottom
*
***********************************************************************
*
*     include 'poisson.h'
#ifdef CMFORTRAN
#include "cm/cmssl-cmf.h"
#endif

#ifdef CMFORTRAN
      real, array(7,n,n,n) :: q
      real, array(1,n,n,n) :: x,y
cmf$ layout q(:serial,:news,:news,:news)
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout y(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(:,:,:,:) :: q
      real, DIMENSION(:,:,:,:) :: x,y
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE y(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      real, array(1,n,n,n) :: xn,xs,xe,xw,xt,xb,zz
cmf$ layout xn(:serial,:news,:news,:news)
cmf$ layout xs(:serial,:news,:news,:news)
cmf$ layout xe(:serial,:news,:news,:news)
cmf$ layout xw(:serial,:news,:news,:news)
cmf$ layout xt(:serial,:news,:news,:news)
cmf$ layout xb(:serial,:news,:news,:news)
cmf$ layout zz(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .             :: xn,xs,xe,xw,xt,xb,zz
CHPF$ DISTRIBUTE xn(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xs(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xe(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xw(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xt(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE xb(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE zz(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      real, array(7,n,n,n) :: qT
cmf$ layout qT(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(7,n,n,n) :: qT
CHPF$ DISTRIBUTE qT(*, BLOCK, BLOCK, BLOCK)
#endif

      integer ntest,itmax,iters
      common /ints/ ntest,itmax,iters

      integer setup, ier
      save setup

      integer axes(3)
#ifdef CMFORTRAN
cmf$ layout axes(:serial)
#elif HPF
CHPF$ DISTRIBUTE axes(*)
#endif
*
***********************************************************************
*

* On first call set up
#ifdef CMFORTRAN
      if (iters.eq.0) then
        setup = pshift_setup (6,x,ier,
     .                     CMSSL_CSHIFT,2, 1,
     .                     CMSSL_CSHIFT,2,-1,
     .                     CMSSL_CSHIFT,3, 1,
     .                     CMSSL_CSHIFT,3,-1,
     .                     CMSSL_CSHIFT,4, 1,
     .                     CMSSL_CSHIFT,4,-1)
      endif
#endif

      if (job .eq. 0) then
         y = x
         zz = x

         do 100 i=1,m
#ifdef CMFORTRAN
            call pshift(6,setup,ier,
     .           CMSSL_CSHIFT,xe,zz,2, 1,
     .           CMSSL_CSHIFT,xw,zz,2,-1,
     .           CMSSL_CSHIFT,xn,zz,3, 1,
     .           CMSSL_CSHIFT,xs,zz,3,-1,
     .           CMSSL_CSHIFT,xt,zz,4, 1,
     .           CMSSL_CSHIFT,xb,zz,4,-1)
#elif HPF
            xe(:,:,:,:) = cshift(zz(:,:,:,:), 1, 2)
            xw(:,:,:,:) = cshift(zz(:,:,:,:), -1, 2)
            xn(:,:,:,:) = cshift(zz(:,:,:,:), 1, 3)
            xs(:,:,:,:) = cshift(zz(:,:,:,:), -1, 3)
            xt(:,:,:,:) = cshift(zz(:,:,:,:), 1, 4)
            xb(:,:,:,:) = cshift(zz(:,:,:,:), -1, 4)
#endif
            zz(1,:,:,:) =
     .           -xn(1,:,:,:) * q(2,:,:,:)
     .           -xs(1,:,:,:) * q(3,:,:,:)
     .           -xe(1,:,:,:) * q(4,:,:,:)
     .           -xw(1,:,:,:) * q(5,:,:,:)
     .           -xt(1,:,:,:) * q(6,:,:,:)
     .           -xb(1,:,:,:) * q(7,:,:,:)

            y = y + zz
 100     continue

      elseif (job .eq. 1) then
*          if (iters.le.1) then

         qT(1,:,:,:)= q(1,:,:,:)
#ifdef CMFORTRAN
*        Syntax for CMF: eoshift(array, dim, shift [, boundary])
         qT(2,:,:,:) = eoshift(q(2,:,:,:), 2, -1, 0.0)
         qT(3,:,:,:) = eoshift(q(3,:,:,:), 2,  1, 0.0)
         qT(4,:,:,:) = eoshift(q(4,:,:,:), 1, -1, 0.0)
         qT(5,:,:,:) = eoshift(q(5,:,:,:), 1,  1, 0.0)
         qT(6,:,:,:) = eoshift(q(6,:,:,:), 3, -1, 0.0)
         qT(7,:,:,:) = eoshift(q(7,:,:,:), 3,  1, 0.0)
#elif defined( HPF ) && defined( NO_EOSHIFT )
         qT(2,:,:,:) = cshift(q(2,:,:,:), -1, 2)
         qT(3,:,:,:) = cshift(q(3,:,:,:),  1, 2)
         qT(4,:,:,:) = cshift(q(4,:,:,:), -1, 1)
         qT(5,:,:,:) = cshift(q(5,:,:,:),  1, 1)
         qT(6,:,:,:) = cshift(q(6,:,:,:), -1, 3)
         qT(7,:,:,:) = cshift(q(7,:,:,:),  1, 3)
#elif HPF
*        Syntax for HPF: eoshift(array, shift [, boundary][, dim])
         qT(2,:,:,:) = eoshift(q(2,:,:,:), -1, 0.0, 2)
         qT(3,:,:,:) = eoshift(q(3,:,:,:),  1, 0.0, 2)
         qT(4,:,:,:) = eoshift(q(4,:,:,:), -1, 0.0, 1)
         qT(5,:,:,:) = eoshift(q(5,:,:,:),  1, 0.0, 1)
         qT(6,:,:,:) = eoshift(q(6,:,:,:), -1, 0.0, 3)
         qT(7,:,:,:) = eoshift(q(7,:,:,:),  1, 0.0, 3)
#endif         
*          endif

         do 200 i=1,m

#ifdef CMFORTRAN
            call pshift(6,setup,ier,
     .           CMSSL_CSHIFT,xe,x,2, 1,
     .           CMSSL_CSHIFT,xw,x,2,-1,
     .           CMSSL_CSHIFT,xn,x,3, 1,
     .           CMSSL_CSHIFT,xs,x,3,-1,
     .           CMSSL_CSHIFT,xt,x,4, 1,
     .           CMSSL_CSHIFT,xb,x,4,-1)
#elif HPF
            xe(:,:,:,:) = cshift(x(:,:,:,:), 1, 2)
            xw(:,:,:,:) = cshift(x(:,:,:,:), -1, 2)
            xn(:,:,:,:) = cshift(x(:,:,:,:), 1, 3)
            xs(:,:,:,:) = cshift(x(:,:,:,:), -1, 3)
            xt(:,:,:,:) = cshift(x(:,:,:,:), 1, 4)
            xb(:,:,:,:) = cshift(x(:,:,:,:), -1, 4)
#endif

            zz(1,:,:,:) = 
     .           -xn(1,:,:,:) * qT(3,:,:,:)
     .           -xs(1,:,:,:) * qT(2,:,:,:)
     .           -xe(1,:,:,:) * qT(5,:,:,:)
     .           -xw(1,:,:,:) * qT(4,:,:,:)
     .           -xt(1,:,:,:) * qT(7,:,:,:)
     .           -xb(1,:,:,:) * qT(6,:,:,:)

            y = y + zz
 200     continue

      else
         print*, 'preconl called with job = ',job
         print*, 'which is not supported by preconl1'

      endif

      return
      end
*
***********************************************************************
*
      subroutine preconr
*
***********************************************************************
*
* This is a dummy routine
*
***********************************************************************
*
      return
      end
*
***********************************************************************
*
