* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: cmstab2.cpp,v 1.4 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: cmstab2.cpp,v $
# Revision 1.4  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.3  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.2  1995/02/22  22:17:09  jtowns
# Fix includes to be directives handled by new Makefile with
# preprocessing.
#
# Revision 1.1  1995/02/11  20:22:51  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      subroutine cg (matvec,a,ia,x,b,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned CG (OMIN) iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = v(0) = b - A*x
*
* Compute for n=1,2,3,...
*
*     alpha = <r(n-1),r(n-1)>/<A*v(n-1),v(n-1)>
*     x(n)  = x(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*v(n-1)
*     beta  = <r(n),r(n)>/<r(n-1),r(n-1)>
*     v(n)  = r(n) + beta*v(n-1)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b.
*
***********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif
      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit

      external matvec
*
***********************************************************************
*
*
      real alpha,beta,rho,rhop

#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,v,A_v
cmf$ align r(i,j,k,l)   with x(i,j,k,l)
cmf$ align v(i,j,k,l)   with x(i,j,k,l)
cmf$ align A_v(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: r,v,A_v
CHPF$ ALIGN r(:,:,:,:)   WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)   WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:) WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

      call matvec (0,a,ia,x,A_v)
      r = b - A_v
      v = r

      rho = sum(r*r)

      rnorm = sqrt(rho)
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .     bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

 100  iters = iters + 1

      call matvec (0,a,ia,v,A_v)
      alpha = rho/sum(A_v*v)

      x = x + alpha*v
      r = r - alpha*A_v

      rhop = rho
      rho = sum(r*r)

      rnorm = sqrt(rho)
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

      if (iters.ge.itmax) then
         ier = 4001
         return
      endif

      beta = rho/rhop
      v = r + beta*v

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine cgnr (matvec,a,ia,x,b,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned CGNR (OMIN) iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     z(0) = v(0) = AT*r(0)
*
* Compute for n=1,2,3,...
*
*     alpha = <z(n-1),z(n-1)>/<A*v(n-1),A*v(n-1)>
*     x(n)  = x(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*v(n-1)
*     z(n)  = AT*r(n)
*     beta  = <z(n),z(n)>/<z(n-1),z(n-1)>
*     v(n)  = z(n) + beta*v(n-1)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b.
*
***********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit

      external matvec
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,v,z_A_v
cmf$ align r(i,j,k,l)   with x(i,j,k,l)
cmf$ align v(i,j,k,l)   with x(i,j,k,l)
cmf$ align z_A_v(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: r,v,z_A_v
CHPF$ ALIGN r(:,:,:,:)   WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)   WITH x(:,:,:,:)
CHPF$ ALIGN z_A_v(:,:,:,:) WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

      call matvec (0,a,ia,x,z_A_v)
      r = b - z_A_v
      call matvec (1,a,ia,r,z_A_v)
      v = z_A_v

      rho = sum(z_A_v*z_A_v)

      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return
      
 100  iters = iters + 1

      call matvec(0,a,ia,v,z_A_v)
      alpha = rho/sum(z_A_v*z_A_v)

      x = x + alpha*v
      r = r - alpha*z_A_v
      call matvec(1,a,ia,r,z_A_v)

      rhop = rho
      rho = sum(z_A_v*z_A_v)

      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

      if (iters.ge.itmax) then
         ier = 4001
         return
      endif

      beta = rho/rhop
      v = z_A_v + beta*v

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicg (matvec,a,ia,x,b,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CG iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = s(0) = v(0) = w(0) = b - A*x
*     beta = rho(0) = <s(0),r(0)>
*     
*    Compute for n=1,2,3,...
*
*     alpha  = rho(n-1)/<w(n-1),A*v(n-1)>
*     x(n)   = x(n-1) + alpha*v(n-1)
*     r(n)   = r(n-1) - alpha*A*v(n-1)
*     s(n)   = s(n-1) - alpha*AT*w(n-1)
*     rho(n) = <s(n),r(n)>
*     beta   = rho(n)/rho(n-1)
*     v(n)   = r(n) + beta*v(n-1)
*     w(n)   = s(n) + beta*w(n-1)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if <s(n),r(n)> = 0 or
* <w(n),A*v(n)> = 0., the algorithm breaks down and will exit if either
* irstrt = 0, or irstrt = 1 and restarts = maxrstrt.
*
***********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit,irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /rstints/ irstrt,maxrstrt,restarts

      external matvec
*
***********************************************************************
*
      real alpha,beta,rho,rhop

#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,s,v,w,A_v_AT_w
cmf$ align r(i,j,k,l)        with x(i,j,k,l)
cmf$ align s(i,j,k,l)        with x(i,j,k,l)
cmf$ align v(i,j,k,l)        with x(i,j,k,l)
cmf$ align w(i,j,k,l)        with x(i,j,k,l)
cmf$ align A_v_AT_w(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: r,s,v,w,A_v_AT_w
CHPF$ ALIGN r(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN s(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN A_v_AT_w(:,:,:,:) WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

 10   call matvec(0,a,ia,x,A_v_AT_w)
      r = b - A_v_AT_w
      s = r
      v = r
      w = r

      rho = sum(s*r)

      rnorm = sqrt(rho)
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .     bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return
      
 100  iters = iters + 1
      
      call matvec(0,a,ia,v,A_v_AT_w)
      alpha = sum(w*A_v_AT_w)

      if (alpha.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      alpha = rho/alpha

      x = x + alpha*v
      r = r - alpha*A_v_AT_w

      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

      if (iters.ge.itmax) then
         ier = 4001
         return
      endif

      call matvec(1,a,ia,w,A_v_AT_w)
      s = s - alpha*A_v_AT_w

      rhop = rho
      rho = sum(s*r)

      if (rho.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      beta = rho/rhop
      v = r + beta*v
      w = s + beta*w

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgs (matvec,a,ia,x,b,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned (Bi-)CGS iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = u = v(0) = b - A*x
*     rho(0) = <r(0),r(0)>
*     
*    Compute for n=1,2,3,...
*
*     alpha  = rho(n-1)/<r(0),A*v(n-1)>
*     w(n)   = u - alpha*A*v(n-1)
*     u      = u + w(n)
*     r(n)   = r(n-1) - alpha*A*u
*     x(n)   = x(n-1) + alpha*u
*     rho(n) = <r(0),r(n)>
*     beta   = rho(n)/rho(n-1)
*     u      = r(n) + beta*w(n)
*     v(n)   = u + beta*w(n) + beta^2*v(n-1)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if rho(n) = 0 or
* <r(0),A*v(n)> = 0., the algorithm breaks down and will exit if either
* irstrt = 0, or irstrt = 1 and restarts = maxrstrt.
*
***********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit,irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /rstints/ irstrt,maxrstrt,restarts

      external matvec
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .           :: r0,r,u,v,w,A_v_A_u
cmf$ align r0(i,j,k,l)      with x(i,j,k,l)
cmf$ align r(i,j,k,l)       with x(i,j,k,l)
cmf$ align u(i,j,k,l)       with x(i,j,k,l)
cmf$ align v(i,j,k,l)       with x(i,j,k,l)
cmf$ align w(i,j,k,l)       with x(i,j,k,l)
cmf$ align A_v_A_u(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .           :: r0,r,u,v,w,A_v_A_u
CHPF$ ALIGN r0(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN u(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN A_v_A_u(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

 10   call matvec(0,a,ia,x,A_v_A_u)
      r = b - A_v_A_u
      r0 = r
      u = r
      v = r

      rho = sum(r0*r)
      
      rnorm = sqrt(rho)
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return
      
 100  iters = iters + 1

      call matvec(0,a,ia,v,A_v_A_u)
      alpha = sum(r0*A_v_A_u)

      if (alpha.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      alpha = rho/alpha

      w = u - alpha*A_v_A_u

      u = u + w
      call matvec(0,a,ia,u,A_v_A_u)

      x = x + alpha*u
      r = r - alpha*A_v_A_u

      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

      if (iters.ge.itmax) then
         ier = 4001
         return
      endif

      rhop = rho
      rho = sum(r0*r)

      if (rho.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      beta = rho/rhop
      u = r + beta*w
      v = u + beta*w + beta*beta*v

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgst (matvec,a,ia,x,b,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CGStab iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     rho(0) = alpha = omega(0) = 1.
*     v(0) = 0.
*
*    Compute for n=1,2,3,...
*
*     rho(n)   = <r(0),r(n-1)>
*     beta     = ( rho(n)/rho(n-1) )*( alpha/omega(n-1) )
*     v(n)     = r(n-1) + beta*v(n-1) - beta*omega(n-1)*A*v(n-1)
*     alpha    = rho(n)/<r(0),A*v(n)>
*     s        = r(n-1) - alpha*A*v(n)
*     omega(n) = <A*s,s>/<A*s,A*s>
*     x(n)     = x(n-1) + alpha*v(n) + omega(n)*s
*     r(n)     = s - omega(n)*A*s
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if rho(n) = 0 or
* <r(0),A*v(n)> = 0., the algorithm breaks down and will exit if either
* irstrt = 0, or irstrt = 1 and restarts = maxrstrt.
*
***********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit,irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /rstints/ irstrt,maxrstrt,restarts

      external matvec
*
***********************************************************************
*
      real rho,rhop,beta,alpha,omega
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r_s,r0,v,A_s,A_v
cmf$ align r_s(i,j,k,l) with x(i,j,k,l)
cmf$ align r0(i,j,k,l)  with x(i,j,k,l)
cmf$ align v(i,j,k,l)   with x(i,j,k,l)
cmf$ align A_s(i,j,k,l) with x(i,j,k,l)
cmf$ align A_v(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: r_s,r0,v,A_s,A_v
CHPF$ ALIGN r_s(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN r0(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)   WITH x(:,:,:,:)
CHPF$ ALIGN A_s(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

 10   call matvec(0,a,ia,x,A_v)
      r0 = b - A_v
      r_s = r0

      rho = 1.
      alpha = 1.
      omega = 1.
      A_v = 0.
      v = 0.

      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r_s,r_s,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

 100  iters = iters + 1

      rhop = rho
      rho = sum(r0*r_s)

      if (rho.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      beta = (rho/rhop)*(alpha/omega)      

      v = r_s + beta*v - beta*omega*A_v

      call matvec(0,a,ia,v,A_v)
      alpha = sum(r0*A_v)

      if (alpha.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      alpha = rho/alpha

      r_s = r_s - alpha*A_v

      call matvec(0,a,ia,r_s,A_s)
      omega = sum(A_s*r_s)/sum(A_s*A_s)

      x = x + alpha*v + omega*r_s
      r_s = r_s - omega*A_s
      
      rnorm = 0.
      rinfnorm = 0.
      call chkcon (r_s,r_s,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

      if (iters.ge.itmax) then
         ier = 4001
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgst2 (matvec,a,ia,x,b,ier)
*
**********************************************************************
*
*     This routine performs the unpreconditioned Bi-CGStab2 iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = v(0) = b - A*x
*     rho(0) = <r(0),r(0)>
*     alpha(0) = rho(0)/<r(0),A*v(0)>
*
*  For n = 1,2,3...
*
*     s(n)     = r(n-1) - alpha(n-1)*A*v(n-1)
*     if n is odd:
*          omega  = <A*s(n),s(n)>/<A*s(n),A*s(n)>
*          x(n)   = x(n-1) + alpha(n-1)*v(n-1) + omega*s(n)
*          r(n)   = s(n) - omega*A*s(n)
*          rho(n) = <r(0),r(n)>
*          beta     = - ( rho(n)/rho(n-1) ) * ( alpha(n-1)/omega )
*          v(n)   = r(n) - beta*(v(n-1) - omega*A*v(n-1))
*     if n is even:
*          denom  = <A*s,A*s>*<s - w,s - w>
*                   - <A*s,s - w>*<A*s,s - w>
*          chi    = (<A*s,w>*<A*s,s - w>
*                              - <A*s,A*s>*<s-w,w>)/denom
*          eta    = (<A*s,s - w>*<s - w,w>
*                            - <A*s,w>*<s - w,s - w>)/denom
*          x(n)   = (1 - chi)*(x(n-2) - alpha(n-2)*v(n-2) 
*                                + alpha(n-1)*t(n-1))
*                     + chi*(x(n-1) + alpha(n-1)*v(n-1)) - eta*s(n)
*          r(n)   = (1 - chi)*w(n) + chi*s(n) + eta*A*s(n)
*          rho(n) = <r(0),r(n)>
*          beta   = ( rho(n)/(rho(n-1) )*( alpha(n-1)/eta )
*          v(n)   = r(n) - beta*((1 - chi)*t(n-1) 
*                                    + chi*v(n-1) + eta*A*v(n-1))
*     alpha(n) = rho(n)/<r(0),A*v(n)>
*     t(n)     = s(n) - beta*v(n-1)
*     w(n+1)   = s(n) - alpha(n)*A*t(n)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if rho(n) = 0 or
* <r(0),A*v(n)> = 0., the algorithm breaks down and will exit if either
* irstrt = 0, or irstrt = 1 and restarts = maxrstrt.
*
**********************************************************************
*
      implicit none

#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer ier
      real x(:,:,:,:),b(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

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

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit,irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /rstints/ irstrt,maxrstrt,restarts

      external matvec
*
***********************************************************************
*
      real rho,rhop,beta,alpha,alphapp,omega,denom,chi,eta,
     .     As_As,As_w,As_sw,sw_w,sw_sw
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: xpp,r,s,r0,v,vpp,w,t,A_s,A_v,sw_vp_A_t
cmf$ align xpp(i,j,k,l)       with x(i,j,k,l)
cmf$ align r(i,j,k,l)         with x(i,j,k,l)
cmf$ align s(i,j,k,l)         with x(i,j,k,l)
cmf$ align r0(i,j,k,l)        with x(i,j,k,l)
cmf$ align v(i,j,k,l)         with x(i,j,k,l)
cmf$ align vpp(i,j,k,l)       with x(i,j,k,l)
cmf$ align w(i,j,k,l)         with x(i,j,k,l)
cmf$ align t(i,j,k,l)         with x(i,j,k,l)
cmf$ align A_s(i,j,k,l)       with x(i,j,k,l)
cmf$ align A_v(i,j,k,l)       with x(i,j,k,l)
cmf$ align sw_vp_A_t(i,j,k,l) with x(i,j,k,l)
#elif HPF
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: xpp,r,s,r0,v,vpp,w,t,A_s,A_v,sw_vp_A_t
CHPF$ ALIGN xpp(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN s(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN r0(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN vpp(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN t(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN A_s(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN sw_vp_A_t(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,rinfnorm,bottom
*
***********************************************************************
*
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))

      iters = 0
      if (level.ge.2) call output(0)

 10   call matvec(0,a,ia,x,A_v)
      r = b - A_v
      r0 = r
      v = r

      call matvec (0,a,ia,v,A_v)

      xpp = x
      sw_vp_A_t = v
      vpp = v

      rho = sum(r0*r)
      alpha = rho/sum(r0*A_v)

      rnorm = rho
      rinfnorm = 0.
      call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return

 100  iters = iters + 1

      s = r - alpha*A_v
      call matvec (0,a,ia,s,A_s)

      if (mod(iters,2).eq.1) then
         omega = sum(A_s*s)/sum(A_s*A_s)

         x = x + alpha*v + omega*s
         r = s - omega*A_s

         rnorm = 0.
         rinfnorm = 0.
         call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .                bottom,istop,ier)
         if (level.ge.2) call output(1)
         if ((istop.eq.1).or.(ier.gt.0)) return

         if (iters.ge.itmax) then
            ier = 4001
            return
         endif

         rhop = rho
         rho = sum(r0*r)

         if (rho.eq.0.) then
            if (irstrt.eq.0) then
               ier = 4010
               return
            else
               if (restarts.le.maxrstrt) then
                  restarts = restarts + 1
                  call output(-1)
                  goto 10
               else
                  ier = 4011
                  return
               endif
            endif
         endif

         beta = (rho/rhop)*(alpha/omega)

         sw_vp_A_t = v
         v = r + beta*sw_vp_A_t - beta*omega*A_v


      else


         sw_vp_A_t = s - w
         As_As = sum(A_s*A_s)
         As_sw = sum(A_s*sw_vp_A_t)
         sw_sw = sum(sw_vp_A_t*sw_vp_A_t)
         As_w = sum(A_s*w)
         sw_w = sum(sw_vp_A_t*w)

         denom = As_As*sw_sw - As_sw*As_sw
         chi = (As_w*As_sw - As_As*sw_w)/denom
         eta = (As_sw*sw_w - As_w*sw_sw)/denom

         x = (1. - chi)*xpp - (1. - chi)*alphapp*vpp + (1. - chi)*alpha*t
     .        + chi*x + chi*alpha*v - eta*s
         xpp = x

         r = (1. - chi)*w + chi*s + eta*A_s

         rnorm = 0.
         rinfnorm = 0.
         call chkcon (r,r,rnorm,rnorm,rinfnorm,rinfnorm,
     .                bottom,istop,ier)
         if (level.ge.2) call output(1)
         if ((istop.eq.1).or.(ier.gt.0)) return

         if (iters.ge.itmax) then
            ier = 4001
            return
         endif

         rhop = rho
         rho = sum(r0*r)
         
         if (rho.eq.0.) then
            if (irstrt.eq.0) then
               ier = 4010
               return
            else
               if (restarts.le.maxrstrt) then
                  restarts = restarts + 1
                  call output(-1)
                  goto 10
               else
                  ier = 4011
                  return
               endif
            endif
         endif

         beta = (rho/rhop)*(alpha/eta)
         sw_vp_A_t = v
         v = r - beta*(1. - chi)*t - beta*chi*sw_vp_A_t - beta*eta*A_v
         vpp = v

      endif

      t = s - beta*sw_vp_A_t
      sw_vp_A_t = A_s - beta*A_v

      call matvec (0,a,ia,v,A_v)
      alphapp = alpha
      alpha = sum(r0*A_v)

      if (alpha.eq.0.) then
         if (irstrt.eq.0) then
            ier = 4010
            return
         else
            if (restarts.le.maxrstrt) then
               restarts = restarts + 1
               call output(-1)
               goto 10
            else
               ier = 4011
               return
            endif
         endif
      endif

      alpha = rho/alpha

      w = s - alpha*sw_vp_A_t

      goto 100

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