* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: cmstab3.cpp,v 1.4 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: cmstab3.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:41  jtowns
# Fix includes to be directives handled by new Makefile with
# preprocessing.
#
# Revision 1.1  1995/02/11  20:22:59  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      subroutine cg_left (matvec,preconl,a,ia,x,b,q,iq,ier)
*
***********************************************************************
*
*     This routine performs the left preconditioned CG (OMIN) iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     z(0) =  v(0) = Q*r(0)
*
*    Compute for n=1,2,3,...
*
*     alpha   = <z(n-1),z(n-1)>/<Q*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)
*     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 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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,z,v,A_v,QA_v
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 A_v(i,j,k,l)  with x(i,j,k,l)
cmf$ align QA_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,z,v,A_v,QA_v
CHPF$ ALIGN r(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:)  WITH x(:,:,:,:)
CHPF$ ALIGN QA_v(:,:,:,:) WITH x(:,:,:,:)
#endif
      integer istop
      real rnorm,znorm,rinfnorm,zinfnorm,bottom
*
***********************************************************************
*
      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,A_v)
      r = b - A_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

 100  iters = iters + 1

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

      x = x + alpha*v
      r = r - alpha*A_v
      z = z - alpha*QA_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)) return

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

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

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine cgnr_left (matvec,preconl,a,ia,x,b,q,iq,ier)
*
***********************************************************************
*
*     This routine performs the left preconditioned CGNR (OMIN)
* iterative algorithm for solving systems of linear equations as
* follows:
*
*     r(0) = b - A*x
*     zz(0) = v(0) = AT*QT*Q*r(0)
*
*    Compute for n=1,2,3,...
*
*     alpha = <zz(n-1),zz(n-1)>/<Q*A*v(n-1),Q*A*v(n-1)>
*     x(n)  = x(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*v(n-1)
*     zz(n) = AT*QT*Q*r(n)
*     beta  = <zz(n),zz(n)>/<zz(n-1),zz(n-1)>
*     v(n)  = zz(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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,z,v,A_v,zz_QA_v
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 A_v(i,j,k,l)     with x(i,j,k,l)
cmf$ align zz_QA_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,z,v,A_v,zz_QA_v
CHPF$ ALIGN r(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN zz_QA_v(:,:,:,:) WITH x(:,:,:,:)
#endif

      integer istop
      real rnorm,rinfnorm,znorm,zinfnorm,bottom
*
***********************************************************************
*
      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,A_v)
      r = b - A_v
      call preconl (0,q,iq,r,z)
      call preconl (1,q,iq,z,v)
      call matvec(1,a,ia,v,zz_QA_v)

      v = zz_QA_v

      rho = sum(zz_QA_v*zz_QA_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)) return

 100  iters = iters + 1

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

      x = x + alpha*v
      r = r - alpha*A_v
      z = z - alpha*zz_QA_v
      call preconl (1,q,iq,z,A_v)
      call matvec (1,a,ia,A_v,zz_QA_v)

      rhop = rho
      rho = sum(zz_QA_v*zz_QA_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)) return

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

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

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicg_left (matvec,preconl,a,ia,x,b,q,iq,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CG  iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     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*v(n-1)>
*     x(n)  = x(n-1) + alpha*v(n-1)
*     r(n)  = r(n-1) - alpha*A*v(n-1)
*     z(n)  = Q*r(n)
*     s(n)  = s(n-1) - alpha*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 the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if <s(n),z(n)> = 0 or
* <w(n),Q*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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r,z,s,v,w,A_v_ATQT_w,QA_v
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 A_v_ATQT_w(i,j,k,l) with x(i,j,k,l)
cmf$ align QA_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,z,s,v,w,A_v_ATQT_w,QA_v
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 A_v_ATQT_w(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN QA_v (:,:,:,:)      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,A_v_ATQT_w)
      r = b - A_v_ATQT_w
      call preconl (0,q,iq,r,z)
      s = z
      v = z
      w = z

      rho = sum(s*z)
      beta = rho

      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

 100  iters = iters + 1

      call matvec (0,a,ia,v,A_v_ATQT_w)
      call preconl (0,q,iq,A_v_ATQT_w,QA_v)
      alpha = sum(w*QA_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

      x = x + alpha*v
      r = r - alpha*A_v_ATQT_w
      z = z - alpha*QA_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)) return

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

      call preconl (1,q,iq,w,QA_v)
      call matvec (1,a,ia,QA_v,A_v_ATQT_w)
      
      s = s - alpha*A_v_ATQT_w

      rhop = rho
      rho = sum(s*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
      v = z + beta*v
      w = s + beta*w

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgs_left (matvec,preconl,a,ia,x,b,q,iq,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned (Bi-)CGS iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     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*v(n-1)>
*     w(n)   = u - alpha*Q*A*v(n-1)
*     u      = u + w(n)
*     r(n)   = r(n-1) - alpha*A*u
*     z(n)   = Q*r(n)
*     x(n)   = x(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 the stopping is criterion satisfied, the process terminates
* and x(n) is the solution of A*x = b; if rho(n) = 0 or
* <z(0),Q*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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
*
***********************************************************************
*
      real alpha,beta,rho,rhop
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .           :: r,z0,u,v,w_QA_v,z_A_u
cmf$ align r(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_QA_v(i,j,k,l) with x(i,j,k,l)
cmf$ align z_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))
     .           :: r,z0,u,v,w_QA_v,z_A_u
CHPF$ ALIGN r(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN z0(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN u(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN w_QA_v(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN z_A_u(:,:,:,:)  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_A_u)
         bottom = sqrt(sum(z_A_u*z_A_u))
      endif

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

 10   call matvec (0,a,ia,x,z_A_u)
      r = b - z_A_u
      call preconl (0,q,iq,r,z0)
      z_A_u = z0
      u = z0
      v = z0

      rho = sum(z0*z_A_u)
      beta = rho
      
      rnorm = 0.
      znorm = sqrt(rho)
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z_A_u,rnorm,znorm,rinfnorm,zinfnorm,
     .             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_u)
      call preconl (0,q,iq,z_A_u,w_QA_v)
      alpha = sum(z0*w_QA_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_QA_v = u - alpha*w_QA_v

      u = u + w_QA_v
      call matvec(0,a,ia,u,z_A_u)

      x = x + alpha*u
      r = r - alpha*z_A_u
      call preconl (0,q,iq,r,z_A_u)

      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r,z_A_u,rnorm,znorm,rinfnorm,zinfnorm,
     .             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(z0*z_A_u)

      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 = z_A_u + beta*w_QA_v
      v = u + beta*w_QA_v + beta*beta*v

      goto 100

      return
      end
*
***********************************************************************
*
      subroutine bicgst_left (matvec,preconl,a,ia,x,b,q,iq,ier)
*
***********************************************************************
*
*     This routine performs the unpreconditioned Bi-CGStab iterative
* algorithm for solving systems of linear equations as follows:
*
*     r(0) = b - A*x
*     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*v(n-1)
*     alpha    = rho(n)/<z(0),Q*A*v(n)>
*     s        = r(n-1) - alpha*A*v(n)
*     t        = Q*s
*     omega(n) = <Q*A*t,t>/<Q*A*t,Q*A*t>
*     x(n)     = x(n-1) + alpha*v(n) + omega(n)*t
*     r(n)     = s - omega(n)*A*t
*     z(n)     = Q*r(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
* <z(0),Q*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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
*
***********************************************************************
*
      real rho,rhop,beta,alpha,omega
#ifdef CMFORTRAN
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: r_s,z_t,z0,v,A_t_A_v,QA_v,QA_t
cmf$ align r_s(i,j,k,l)     with x(i,j,k,l)
cmf$ align z_t(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 A_t_A_v(i,j,k,l) with x(i,j,k,l)
cmf$ align QA_v(i,j,k,l)    with x(i,j,k,l)
cmf$ align QA_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))
     .          :: r_s,z_t,z0,v,A_t_A_v,QA_v,QA_t
CHPF$ ALIGN r_s(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z_t(:,:,:,:)     WITH x(:,:,:,:)
CHPF$ ALIGN z0(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN A_t_A_v(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN QA_v(:,:,:,:)    WITH x(:,:,:,:)
CHPF$ ALIGN QA_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_t)
         bottom = sqrt(sum(z_t*z_t))
      endif

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

 10   call matvec (0,a,ia,x,A_t_A_v)
      r_s = b - A_t_A_v
      call preconl (0,q,iq,r_s,z_t)
      z0 = z_t

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

      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r_s,z_t,rnorm,znorm,rinfnorm,zinfnorm,
     .             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(z0*z_t)

      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 = z_t + beta*v - beta*omega*QA_v

      call matvec(0,a,ia,v,A_t_A_v)
      call preconl (0,q,iq,A_t_A_v,QA_v)
      alpha = sum(z0*QA_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_t_A_v
      z_t = z_t - alpha*QA_v

      call matvec(0,a,ia,z_t,A_t_A_v)
      call preconl (0,q,iq,A_t_A_v,QA_t)
      
      omega = sum(QA_t*z_t)/sum(QA_t*QA_t)

      x = x + alpha*v + omega*z_t
      r_s = r_s - omega*A_t_A_v
      z_t = z_t - omega*QA_t
      
      rnorm = 0.
      znorm = 0.
      rinfnorm = 0.
      zinfnorm = 0.
      call chkcon (r_s,z_t,rnorm,znorm,rinfnorm,zinfnorm,
     .             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_left (matvec,preconl,a,ia,x,b,q,iq,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*v(0)>
*
*  For n = 1,2,3...
*
*     u(n)     = r(n-1) - alpha(n-1)*A*v(n-1)
*     s(n)     = z(n-1) - alpha(n-1)*Q*A*v(n-1)
*     if n is odd:
*          omega  = <Q*A*s(n),s(n)>/<Q*A*s(n),Q*A*s(n)>
*          x(n)   = x(n-1) + alpha(n-1)*v(n-1) + omega*s(n)
*          r(n)   = u(n) - omega*A*s(n)
*          z(n)   = s(n) - omega*Q*A*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*v(n-1))
*     if n is even:
*          denom  = <Q*A*s,Q*A*s>*<s - w,s - w>
*                   - <Q*A*s,s - w>*<Q*A*s,s - w>
*          chi    = - (<Q*A*s,w>*<Q*A*s,s - w>
*                              - <Q*A*s,Q*A*s>*<s-w,w>)/denom
*          eta    = - (<Q*A*s,s - w>*<s - w,w>
*                            - <Q*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)*h(n) + chi*u(n) + eta*A*s(n)
*          z(n)   = (1 - chi)*w(n) + chi*s(n) + eta*Q*A*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*v(n-1))
*     alpha(n) = rho(n)/<z(0),Q*A*v(n)>
*     t(n)     = s(n) - beta*v(n-1)
*     h(n+1)   = u(n) - alpha(n)*A*t(n)
*     w(n+1)   = Q*h(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
* <z(0),Q*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),q(CMF_sizeof_descriptor)
      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
#elif HPF
      real  a(:,:,:,:),  q(:,:,:,:)
      integer ia(:,:,:,:), iq(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, 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
      external preconl
*
***********************************************************************
*
      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,z0,v,vpp,w,t,A_s,A_v,sw_vp_A_t
      real, array(dsize(x,1),dsize(x,2),dsize(x,3),dsize(x,4))
     .          :: z,u,h,QA_s,QA_v,QA_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 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 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)
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 QA_s(i,j,k,l)      with x(i,j,k,l)
cmf$ align QA_v(i,j,k,l)      with x(i,j,k,l)
cmf$ align QA_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,z0,v,vpp,w,t,A_s,A_v,sw_vp_A_t
      real, DIMENSION(size(x,1),size(x,2),size(x,3),size(x,4))
     .     :: z,u,h,QA_s,QA_v,QA_t
CHPF$ ALIGN xpp(:,:,:,:)       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 A_s(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN A_v(:,:,:,:)       WITH x(:,:,:,:)
CHPF$ ALIGN sw_vp_A_t(:,:,:,:) WITH x(:,:,:,:)
CHPF$ ALIGN z(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN u(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN h(:,:,:,:)         WITH x(:,:,:,:)
CHPF$ ALIGN QA_s(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN QA_v(:,:,:,:)      WITH x(:,:,:,:)
CHPF$ ALIGN QA_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,A_v)

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

      call matvec (0,a,ia,v,A_v)
      call preconl (0,q,iq,A_v,QA_v)

      xpp = x
      sw_vp_A_t = v
      vpp = v

      rho = sum(z0*z)
      alpha = rho/sum(z0*QA_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

 100  iters = iters + 1

      u = r - alpha*A_v
      s = z - alpha*QA_v
      call matvec (0,a,ia,s,A_s)
      call preconl (0,q,iq,A_s,QA_s)

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

         x = x + alpha*v + omega*s
         r = u - omega*A_s
         z = s - omega*QA_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)) return

         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_A_t = v
         v = z + beta*sw_vp_A_t - beta*omega*QA_v

      else
         sw_vp_A_t = s - w
         As_As = sum(QA_s*QA_s)
         As_sw = sum(QA_s*sw_vp_A_t)
         sw_sw = sum(sw_vp_A_t*sw_vp_A_t)
         As_w = sum(QA_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)*h + chi*u + eta*A_s
         z = (1. - chi)*w + chi*s + eta*QA_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)) return

         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_A_t = v
         v = z - beta*(1. - chi)*t - beta*chi*sw_vp_A_t - beta*eta*QA_v
         vpp = v

      endif

      t = s - beta*sw_vp_A_t
      sw_vp_A_t = A_s - beta*A_v
      QA_t = QA_s - beta*QA_v

      call matvec (0,a,ia,v,A_v)
      call preconl (0,q,iq,A_v,QA_v)
      alphapp = alpha
      alpha = sum(z0*QA_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_A_t
      w = s - alpha*QA_t

      goto 100

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