* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: cmstab5.cpp,v 1.4 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: cmstab5.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:21  jtowns
# Fix includes to be directives handled by new Makefile with
# preprocessing.
#
# Revision 1.1  1995/02/11  20:23:07  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      subroutine cg_split (matvec,preconl,preconr,a,ia,x,b,
     .                     q,iq,p,ip,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned CG (OMIN) iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     y(0) = inv(P)*x
*     z(0) = v(0) = Q*r(0)
*
*    Compute for n=0,1,2,...
*
*     alpha   = <z(n-1),z(n-1)>/<Q*A*P*v(n-1),v(n-1)>
*     y(n)  = y(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*P*v(n-1)
*     z(n)  = Q*r(n)
*     beta    = <z(n),z(n)>/<z(n-1),z(n-1)>
*     v(n)  = z(n) + beta*v(n-1)
*
*     If r(n) <= zeta, the process terminates and x(n) = P*y(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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: y,r,z,v,QAP_v,AP_v
cmf$ align y(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 z(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 QAP_v(i,j,k,l) with x(i,j,k,l)
cmf$ align AP_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))
     .          :: y,r,z,v,QAP_v,AP_v
CHPF$ ALIGN y(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)  WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      stptst = 0.

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then
         call preconl (0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

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

      call matvec (0,a,ia,x,AP_v)
      r = b - AP_v
      call preconl (0,q,iq,r,z)
      v = z

      rho = sum(z*z)

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

      call preconr(2,p,ip,x,y)
      
 100  iters = iters + 1

      call preconr (0,p,ip,v,QAP_v)
      call matvec (0,a,ia,QAP_v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)
      alpha = rho/sum(QAP_v*v)

      y = y + alpha*v
      r = r - alpha*AP_v
      z = z - alpha*QAP_v

      rhop = rho
      rho = sum(z*z)

      rnorm = 0.
      znorm = sqrt(rho)
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) then
         call preconr (0,p,ip,y,x)
         return
      endif

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

      if (iters.ge.itmax) then
         ier = 4001
         call preconr (0,p,ip,y,x)
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine cgnr_split (matvec,preconl,preconr,a,ia,x,b,
     .                       q,iq,p,ip,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned CGNR (OMIN) iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     y(0) = inv(P)*x(0)
*     z(0) = Q*r(0)
*     zz(0) = v(0) = PT*AT*QT*z(0)
*
*    Compute for n=1,2,3,...
*
*     alpha = <zz(n-1),zz(n-1)>/<Q*A*P*v(n-1),Q*A*P*v(n-1)>
*     y(n)  = y(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*P*v(n-1)
*     z(n)  = Q*r(n)
*     zz(n) = PT*AT*QT*z(n)
*     beta  = <zz(n),zz(n)>/<zz(n-1),zz(n-1)>
*     v(n)  = zz(n) + beta*v(n-1)
*
*     If z(n) <= zeta, the process terminates and x(n) = P*y(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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      real alpha,beta,rho,rhop

#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: y,r,z,v,zz,AP_v,QAP_v
cmf$ align y(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 z(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 zz(i,j,k,l)    with x(i,j,k,l)
cmf$ align AP_v(i,j,k,l)  with x(i,j,k,l)
cmf$ align QAP_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))
     .          :: y,r,z,v,zz,AP_v,QAP_v
CHPF$ ALIGN y(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN zz(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      stptst = 0.

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then
         call preconl (0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

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

      call matvec (0,a,ia,x,AP_v)
      r = b - AP_v
      call preconl (0,q,iq,r,z)

      call preconl (1,q,iq,z,AP_v)
      call matvec (1,a,ia,AP_v,QAP_v)
      call preconr (1,p,ip,QAP_v,zz)
      v = zz

      rho = sum(zz*zz)

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

      call preconr(2,p,ip,x,y)

 100  iters = iters + 1

      call preconr (0,p,ip,v,QAP_v)
      call matvec(0,a,ia,QAP_v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)
      alpha = rho/sum(QAP_v*QAP_v)

      y = y + alpha*v
      r = r - alpha*AP_v
      z = z - alpha*QAP_v

      call preconl (1,q,iq,z,QAP_v)
      call matvec (1,a,ia,QAP_v,AP_v)
      call preconr (1,p,ip,AP_v,zz)
      
      rhop = rho
      rho = sum(zz*zz)

      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) then
         call preconr (0,p,ip,y,x)
         return
      endif

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

      if (iters.ge.itmax) then
         ier = 4001
         call preconr (0,p,ip,y,x)
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicg_split (matvec,preconl,preconr,a,ia,x,b,
     .                       q,iq,p,ip,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CG iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     y(0) = inv(P)*x(0)
*     z(0) =  s(0) = v(0) = w(0) = Q*r(0)
*     beta = <s(0),z(0)>
*     
*    Compute for n=1,2,3,...
*
*     alpha = <s(n-1),z(n-1)>/<w(n-1),Q*A*P*v(n-1)>
*     y(n)  = y(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*P*v(n-1)
*     z(n)  = Q*r(n)
*     s(n)  = s(n-1) - alpha*PT*AT*QT*w(n-1)
*     beta  = <s(n),z(n)>/<s(n-1),z(n-1)>
*     v(n)  = z(n) + beta*v(n-1)
*     w(n)  = s(n) + beta*w(n-1)
*
*     If r(n) <= zeta, the process terminates and x(n) = P*y(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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: y,r,z,s,v,w,AP_v,QAP_v,PTATQT_w
cmf$ align y(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 z(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 AP_v(i,j,k,l)     with x(i,j,k,l)
cmf$ align QAP_v(i,j,k,l)    with x(i,j,k,l)
cmf$ align PTATQT_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))
     .          :: y,r,z,s,v,w,AP_v,QAP_v,PTATQT_w
CHPF$ ALIGN y(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN s(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN PTATQT_w(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      stptst = 0.
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then
         call preconl (0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

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

 10   call matvec(0,a,ia,x,AP_v)
      r = b - AP_v
      call preconl (0,q,iq,r,z)
      s = z
      v = z
      w = z

      rho = sum(s*z)

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

      call preconr(2,p,ip,x,y)
      
 100  iters = iters + 1
      
      call preconr (0,p,ip,v,QAP_v)
      call matvec(0,a,ia,QAP_v,AP_v)
      call matvec(0,a,ia,v,AP_v)

      call preconl (0,q,iq,AP_v,QAP_v)
      alpha = sum(w*QAP_v)

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

      alpha = rho/alpha

      y = y + alpha*v
      r = r - alpha*AP_V
      z = z - alpha*QAP_v

      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) then
         call preconr (0,p,ip,y,x)
         return
      endif

      call preconl (1,q,iq,w,QAP_v)
      call matvec(1,a,ia,QAP_v,AP_v)
      call preconr (1,p,ip,AP_v,PTATQT_w)
      s = s - alpha*PTATQT_w

      rhop = rho
      rho = sum(s*r)

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

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

      if (iters.ge.itmax) then
         ier = 4001
         call preconr (0,p,ip,y,x)
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgs_split (matvec,preconl,preconr,a,ia,x,b,
     .                        q,iq,p,ip,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned (Bi-)CGS iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     y(0) = inv(P)*x(0)
*     z(0) =  u = v(0) = Q*r(0)
*     rho(0) = <z(0),z(0)>
*     
*    Compute for n=1,2,3,...
*
*     alpha  = rho(n-1)/<z(0),Q*A*P*v(n-1)>
*     w(n)   = u - alpha*Q*A*P*v(n-1)
*     u      = u + w(n)
*     r(n)   = r(n-1) - alpha*A*P*u
*     z(n)   = Q*r(n)
*     y(n)   = y(n-1) + alpha*u
*     rho(n) = <z(0),z(n)>
*     beta   = rho(n)/rho(n-1)
*     u      = z(n) + beta*w(n)
*     v(n)   = u + beta*w(n) + beta^2*v(n-1)
*
*     If r(n) <= zeta, the process terminates and x(n) = P*y(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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .           :: y,r,z0,z,u,v,w,AP_v,QAP_v,AP_u
cmf$ align y(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 z(i,j,k,l)     with x(i,j,k,l)
cmf$ align z0(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 AP_v(i,j,k,l)  with x(i,j,k,l)
cmf$ align QAP_v(i,j,k,l) with x(i,j,k,l)
cmf$ align AP_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))
     .           :: y,r,z0,z,u,v,w,AP_v,QAP_v,AP_u
CHPF$ ALIGN y(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z0(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN u(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN AP_u(:,:,:,:)  WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      stptst = 0.
      restarts = 0

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then 
         call preconl(0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

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

 10   call matvec(0,a,ia,x,AP_v)
      r = b - AP_v
      call preconl(0,q,ip,r,z)
      z0 = z
      u = z
      v = z

      rho = sum(z0*z)
      
      rnorm = 0.
      znorm = sqrt(rho)
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) return
      
      call preconr (2,p,ip,x,y)

 100  iters = iters + 1

      call preconr (0,p,ip,v,QAP_v)
      call matvec(0,a,ia,QAP_v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)
      alpha = sum(z0*QAP_v)

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

      alpha = rho/alpha
      w = u - alpha*QAP_v

      u = u + w
      call preconr(0,p,ip,u,QAP_v)
      call matvec(0,a,ia,QAP_v,AP_u)

      y = y + alpha*u
      r = r - alpha*AP_u
      call preconl(0,q,iq,r,z)

      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) then
         call preconr (0,p,ip,y,x)
         return
      endif

      rhop = rho
      rho = sum(z0*z)

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

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

      if (iters.ge.itmax) then
         ier = 4001
         call preconr (0,p,ip,y,x)
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgst_split (matvec,preconl,preconr,a,ia,x,b,
     .                         q,iq,p,ip,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CGStab iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     y(0) = inv(P)*x(0)
*     z(0) = Q*r(0)
*     rho(0) = alpha = omega(0) = 1.
*     v(0) = 0.
*
*    Compute for n=1,2,3,...
*
*     rho(n)   = <z(0),z(n-1)>
*     beta     = ( rho(n)/rho(n-1) )*( alpha/omega(n-1) )
*     v(n)     = z(n-1) + beta*v(n-1) - beta*omega(n-1)*Q*A*P*v(n-1)
*     alpha    = rho(n)/<z(0),Q*A*P*v(n)>
*     t        = r(n-1) - alpha*A*P*v(n)
*     s        = Q*t
*     omega(n) = <Q*A*P*s,s>/<Q*A*P*s,Q*A*P*s>
*     y(n)     = y(n-1) + alpha*v(n) + omega(n)*s
*     r(n)     = t - omega(n)*A*P*s
*     z(n)     = Q*r(n)
*
*     If r(n) <= zeta, the process terminates and x(n) = P*y(n) is the
* solution of A*x = b; if r(n) <> 0 but rho(n) = 0 or
* <r(0),A*v(n)> = 0 or <A*s,A*s> = 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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      real rho,rhop,beta,alpha,omega

#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: y,r,z,z0,v,s,t,AP_s,QAP_s,AP_v,QAP_v
cmf$ align y(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 z(i,j,k,l)     with x(i,j,k,l)
cmf$ align z0(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 s(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 AP_s(i,j,k,l)  with x(i,j,k,l)
cmf$ align QAP_s(i,j,k,l) with x(i,j,k,l)
cmf$ align AP_v(i,j,k,l)  with x(i,j,k,l)
cmf$ align QAP_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))
     .          :: y,r,z,z0,v,s,t,AP_s,QAP_s,AP_v,QAP_v
CHPF$ ALIGN y(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z0(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN s(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN t(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN AP_s(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN QAP_s(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:) WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      restarts = 0
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      stptst = 0.

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then
         call preconl (0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

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

 10   call matvec(0,a,ia,x,AP_v)
      r = b - AP_v
      call preconl (0,q,iq,r,z)
      z0 = z
      s = z

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

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

      call preconr (2,p,ip,x,y)

 100  iters = iters + 1

      rhop = rho
      rho = sum(z0*z)

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

      beta = (rho/rhop)*(alpha/omega)      
      v = z + beta*v - beta*omega*QAP_v

      call preconr (0,p,ip,v,QAP_v)
      call matvec(0,a,ia,QAP_v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)
      alpha = sum(z0*QAP_v)

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

      alpha = rho/alpha

      t = r - alpha*AP_v
      s = z - alpha*QAP_v
      call preconr (0,p,ip,s,QAP_s)
      call matvec(0,a,ia,QAP_s,AP_s)
      call preconl (0,q,iq,AP_s,QAP_s)

      omega = sum(QAP_s*s)/sum(QAP_s*QAP_s)

      y = y + alpha*v + omega*s
      r = t - omega*AP_s
      z = s - omega*QAP_s
      
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .             bottom,istop,ier)
      if (level.ge.2) call output(1)
      if ((istop.eq.1).or.(ier.gt.0)) then
         call preconr (0,p,ip,y,x)
         return
      endif

      if (iters.ge.itmax) then
         ier = 4001
         call preconr (0,p,ip,y,x)
         return
      endif

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgst2_split (matvec,preconl,preconr,a,ia,x,b,
     .                          q,iq,p,ip,ier)
*
**********************************************************************
*
*     This routine performs the left preconditioned Bi-CGStab2
* iterative algorithm for solving systems of linear equations as
* follows:
*
*     r(0) = b - A*x
*     z(0) = v(0) = Q*r(0)
*     rho(0) = <z(0),z(0)>
*     alpha(0) = rho(0)/<z(0),Q*A*P*v(0)>
*
*  For n = 1,2,3...
*
*     u(n)     = r(n-1) - alpha(n-1)*A*P*v(n-1)
*     s(n)     = z(n-1) - alpha(n-1)*Q*A*P*v(n-1)
*     if n is odd:
*          omega  = <Q*A*P*s(n),s(n)>/<Q*A*P*s(n),Q*A*P*s(n)>
*          y(n)   = y(n-1) + alpha(n-1)*v(n-1) + omega*s(n)
*          r(n)   = u(n) - omega*A*P*s(n)
*          z(n)   = s(n) - omega*Q*A*P*s(n-1)
*          rho(n) = <z(0),z(n)>
*          beta     = - alpha(n-1)*rho(n)/(omega*rho(n-1))
*          v(n)   = z(n) - beta*(v(n-1) - omega*Q*A*P*v(n-1))
*     if n is even:
*          denom  = <Q*A*P*s,Q*A*P*s>*<s - w,s - w>
*                   - <Q*A*P*s,s - w>*<Q*A*P*s,s - w>
*          chi    = - (<Q*A*P*s,w>*<Q*A*P*s,s - w>
*                              - <Q*A*P*s,Q*A*P*s>*<s-w,w>)/denom
*          eta    = - (<Q*A*P*s,s - w>*<s - w,w>
*                            - <Q*A*P*s,w>*<s - w,s - w>)/denom
*          y(n)   = (1 - chi)*(y(n-2) - alpha(n-2)*v(n-2) 
*                                + alpha(n-1)*t(n-1))
*                     + chi*(y(n-1) + alpha(n-1)*v(n-1)) - eta*s(n)
*          r(n)   = (1 - chi)*h(n) + chi*u(n) + eta*A*P*s(n)
*          z(n)   = (1 - chi)*w(n) + chi*s(n) + eta*Q*A*P*s(s)
*          rho(n) = <z(0),z(n)>
*          beta   = ( rho(n)/(rho(n-1) )*( alpha(n-1)/eta )
*          v(n)   = z(n) - beta*((1 - chi)*t(n-1) 
*                                    + chi*v(n-1) + eta*Q*A*P*v(n-1))
*     alpha(n) = rho(n)/<z(0),Q*A*P*v(n)>
*     t(n)     = s(n) - beta*v(n-1)
*     h(n+1)   = u(n) - alpha(n)*A*P*t(n)
*     w(n+1)   = Q*h(n+1)
*
*     If the stopping is criterion satisfied, the process terminates
* and x(n)=P*y(n) is the solution of A*x = b; if rho(n) = 0. or
* <z(0),Q*A*P*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), q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:),  p(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, 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,preconl,preconr
*
***********************************************************************
*
      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))
     .          :: y,ypp,r,s,z0,v,vpp,w,t,AP_s,AP_v,sw_vp_AP_t
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: z,u,h,QAP_s,QAP_v,QAP_t
cmf$ align y(i,j,k,l)          with x(i,j,k,l)
cmf$ align ypp(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 z0(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 AP_s(i,j,k,l)       with x(i,j,k,l)
cmf$ align AP_v(i,j,k,l)       with x(i,j,k,l)
cmf$ align sw_vp_AP_t(i,j,k,l) with x(i,j,k,l)
cmf$ align z(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 h(i,j,k,l)          with x(i,j,k,l)
cmf$ align QAP_s(i,j,k,l)      with x(i,j,k,l)
cmf$ align QAP_v(i,j,k,l)      with x(i,j,k,l)
cmf$ align QAP_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))
     .          :: y,ypp,r,s,z0,v,vpp,w,t,AP_s,AP_v,sw_vp_AP_t
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .          :: z,u,h,QAP_s,QAP_v,QAP_t
CHPF$ ALIGN y(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN ypp(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN r(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN s(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN z0(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN vpp(:,:,:,:)        WITH x(:,:,:,:)
CHPF$ ALIGN w(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN t(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN AP_s(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN AP_v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN sw_vp_AP_t(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN u(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN h(:,:,:,:)          WITH x(:,:,:,:)
CHPF$ ALIGN QAP_s(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN QAP_v(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN QAP_t(:,:,:,:)      WITH x(:,:,:,:)
#endif

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

      bottom = 1.
      if (ntest.eq.2) bottom = sqrt(sum(b*b))
      if (ntest.eq.4) then
         call preconl (0,q,iq,b,z)
         bottom = sqrt(sum(z*z))
      endif

      iters = 0
      if (level.ge.2) call output(0)
 10   call matvec(0,a,ia,x,AP_v)

      r = b - AP_v
      call preconl (0,q,iq,r,z0)
      z = z0
      v = z0

      call preconr(0,p,ip,v,QAP_v)
      call matvec (0,a,ia,v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)

      sw_vp_AP_t = v
      vpp = v

      rho = sum(z0*z)
      alpha = rho/sum(z0*QAP_v)

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

      call preconr(2,p,ip,x,y)
      ypp = y

 100  iters = iters + 1

      u = r - alpha*AP_v
      s = z - alpha*QAP_v
      call preconr(0,p,ip,s,QAP_s)
      call matvec (0,a,ia,s,AP_s)
      call preconl (0,q,iq,AP_s,QAP_s)

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

         y = y + alpha*v + omega*s
         r = u - omega*AP_s
         z = s - omega*QAP_s

         rnorm = 0.
         znorm = 0.
         rinfnorm = 0.
         zinfnorm = 0.
         call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .                bottom,istop,ier)
         if (level.ge.2) call output(1)
         if ((istop.eq.1).or.(ier.gt.0)) then
            call preconr (0,p,ip,y,x)
            return
         endif

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

         rhop = rho
         rho = sum(z0*z)

         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_AP_t = v
         v = z + beta*sw_vp_AP_t - beta*omega*QAP_v


      else


         sw_vp_AP_t = s - w
         As_As = sum(QAP_s*QAP_s)
         As_sw = sum(QAP_s*sw_vp_AP_t)
         sw_sw = sum(sw_vp_AP_t*sw_vp_AP_t)
         As_w = sum(QAP_s*w)
         sw_w = sum(sw_vp_AP_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

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

         r = (1. - chi)*h + chi*u + eta*AP_s
         z = (1. - chi)*w + chi*s + eta*QAP_s

         rnorm = 0.
         znorm = 0.
         rinfnorm = 0.
         zinfnorm = 0.
         call chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .                bottom,istop,ier)
         if (level.ge.2) call output(1)
         if ((istop.eq.1).or.(ier.gt.0)) then
            call preconr (0,p,ip,y,x)
            return
         endif

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

         rhop = rho
         rho = sum(z0*z)
         
         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_AP_t = v
         v = z - beta*(1. - chi)*t - beta*chi*sw_vp_AP_t
     .         - beta*eta*QAP_v
         vpp = v

      endif

      t = s - beta*sw_vp_AP_t
      sw_vp_AP_t = AP_s - beta*AP_v
      QAP_t = QAP_s - beta*QAP_v

      call preconr(0,p,ip,v,QAP_v)
      call matvec (0,a,ia,v,AP_v)
      call preconl (0,q,iq,AP_v,QAP_v)
      alphapp = alpha
      alpha = sum(z0*QAP_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
      h = u - alpha*sw_vp_AP_t
      w = s - alpha*QAP_t

      goto 100

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