* -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*
* $Id: cmstab1.cpp,v 1.5 1995/05/16 16:15:09 dweber Exp dweber $
* /*
* $Log: cmstab1.cpp,v $
# Revision 1.5  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.4  1995/04/22  02:49:50  dweber
# Changed the timing routines to use etime.
# etime declared as real*4 so the type wouldn't be promoted.
#
# 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:24  jtowns
# Fix includes to be directives handled by new Makefile with
# preprocessing.
#
# Revision 1.1  1995/02/11  20:22:34  jtowns
# Initial revision
#
* */
*
***********************************************************************
*
      subroutine cmstabdrv (matvec,preconl,preconr,a,ia,x,b,
     .                q,iq,p,ip,rparm,iparm,ier)
*
***********************************************************************
*
*     cmstab is the interface to CMStab, consisting of six conjugate-
* gradient-based iterative algorithms written to exploit the SIMD
* programming paradigm available on the CM-5....
*
***********************************************************************
*
* Version History:
*
*     V1.00:  First public release
*
***********************************************************************
*
* Calling Parameters:
*
*
*     matvec:     subroutine matvec (job,a,ia,x,y)
*            cmf$ layout a(..., :news, ...)
*            cmf$ layout ia(..., :news, ...)
*            cmf$ layout x(:serial,:news,:news,:news)
*            cmf$ layout y(:serial,:news,:news,:news)
*          This is a user written subroutine that must have the above
*          specified parameter list and layout directives.  The
*          subroutine must perform the matrix vector products specified
*          by the integer parameter job:
*               job = 0  =>  y = A*x
*               job = 1  =>  y = AT*x
*          where A is the coefficient matrix of the system, and AT is
*          the transpose of A.  A is represtented by the arrays a and
*          ia (see below).  All the routines in CMStab require calls
*          with job = 0; the routines for the CGNR and Bi-CG
*          algorithms also require calls to matvec with job = 1.  The
*          parameters x and y are real arrays with the same layout as
*          the array x passed to cmstabdrv.  Only the parameter y may
*          be changed in the subroutine.  The matvec routine may 
*          actually have any name the user chooses, but must be
*          declared external in the user's calling routine.
*
*     preconl:    subroutine preconl (job,q,iq,x,y)
*            cmf$ layout q(..., :news, ...)
*            cmf$ layout iq(..., :news, ...)
*            cmf$ layout x(:serial,:news,:news,:news)
*            cmf$ layout y(:serial,:news,:news,:news)
*          If left- or split-preconditioning is specified by the user,
*          then the user must provide a subroutine having the above
*          specified parameter list and layout directives.  If neither
*          left- nor split-preconditioning is specified, this is a 
*          dummy argument.  The subroutine must perform the matrix 
*          vector products specified by the integer parameter job:
*               job = 0  =>  y = Q*x
*               job = 1  =>  y = QT*x
*          where Q is the left  preconditioning matrix and QT is the
*          transpose of Q.  Q is represented by the arrays q and iq.
*          If left preconditioning is specified, all the routines make
*          calls to preconl with job = 0.  The routines for the CGNR
*          and Bi-CG algorithms also make calls to preconl with
*          job = 1.  The parameters x and y are real arrays with the
*          same layout as the array x passed to cmstabdrv.  Only the
*          parameter y may be changed in the soubroutine.  The preconl
*          routine may actually have any name the user chooses, but
*          must be declared external in the user's calling routine
*
*     preconr:    subroutine preconr (job,p,ip,x,y)
*            cmf$ layout p(..., :news, ...)
*            cmf$ layout ip(..., :news, ...)
*            cmf$ layout x(:serial,:news,:news,:news)
*            cmf$ layout y(:serial,:news,:news,:news)
*          If split-preconditioning is specified by the user, then the
*          user must provide a subroutine having the above specified
*          parameter list and layout directives.  If split-
*          preconditioning is not specified, this is a dummy argument
*          The subroutine must perform the matrix vector products
*          specified by the integer parameter job:
*               job = 0  =>  y = P*x
*               job = 1  =>  y = PT*x
*               job = 2  =>  y = inv(P)*x
*          where P is the right preconditioning matrix, PT is the
*          transpose of P, and inv(P) is the matrix inverse of P.  P
*          is represented by the arrays p and ip. If split-
*          preconditioning is specified, all the routines make calls
*          to preconr with job = 0 and job = 2.  The routines for the
*          CGNR and Bi-CG algorithms also make calls to preconr with
*          job = 1.  The parameters x and y are real arrays with the
*          same layout as the array x passed to cmstabdrv.  Only the
*          parameter y may be changed in the soubroutine.  The preconr
*          routine may actually have any name the user chooses, but
*          must be declared external in the user's calling routine
*
*     a:          real a(..., :news, ...)
*          This is a real CM array of rank 4 containing the coefficient
*          array; the length of each axis is arbirary, but must be at
*          least 1.  a is only accessed by calls to the subroutine
*          matvec.  The storage format for the array is completely
*          specified by the user, limited only by the given shape and
*          may make use of the integer array ia for additional
*          information.
*
*     ia:         integer ia(..., :news, ...)
*          This is an integer CM array of rank 4 containing the users
*          integer work array usually containing additional information
*          about the coefficient array, a, needed by matvec.  The
*          length of each axis is arbirary, but must be at least 1.
*          ia is only accessed by calls to the subroutine matvec.
*
*     x:          real x(:serial,:news,:news,:news)
*          This is a real CM array of rank 4 containing the initial
*          guess, x0, to the solution of the linear system A x = b; the
*          length of each axis is arbirary, but must be at least 1. The
*          storage format for this array is completely specified by the
*          user, limited only by the given shape.  On exit, x contains
*          the final approximate solution to A x = b.
*
*     b:          real b(:serial,:news,:news,:news)
*          This is a real CM array of rank 4 containing the right-hand
*          side vector of the linear system A x = b and must be stored
*          in the same format as x.
*
*     q:          real q(..., :news, ...)
*          This is a real CM array of rank 4 containing the users
*          left-preconditioning array Q; the length of each axis is
*          arbirary, but must be at least 1.  The storage format for
*          this array is completely specified by the user, limited only
*          by the given shape.  q is only accessed by calls to the
*          subroutine preconl.
*
*     iq:         integer iq(..., :news, ...)
*          This is an integer CM array of rank 4 containing the users
*          integer work array usually containing additional information
*          about the left-preconditioning array, q, needed by preconl.
*          The length of each axis is arbitrary, but must be at least
*          1.  iq is only accessed by calls to the subroutine preconl.
*
*     p:          real p(..., :news, ...)
*          This is a real CM array of rank 4 containing the users
*          right-preconditioning array p; the length of each axis is
*          arbirary, but must be at least 1.  The storage format for
*          this array is completely specified by the user, limited only
*          by the given shape.p is only accessed by calls to the
*          subroutine preconr.
*
*     ip:         integer ip(..., :news, ...)
*          This is an integer CM array of rank 4 containing the users
*          integer work array usually containing additional information
*          about the right-preconditioning array, p, needed by preconr.
*          The length of each axis is arbitrary, but must be at least 1.
*          ip is only accessed by calls to the subroutine preconr.
*
*     rparm:      real rparm(20)
*           cmf$ layout rparm(:serial)
*          This is a front end (FE) array of real input parameters as
*          follows:
*
*          rparm(1) = zeta
*               This is the value used in the stopping test (see
*               iparm(1) = ntest below).
*
*          rparm(2) = stptst
*               On output, this contains the last computed value of
*               the stopping test.
*
*          rparm(3) = elapsed
*               On output, this is the elapsed (wall clock) time
*               during the iterative solution process in seconds, if
*               the value of timer (see iparm(11) = timer below) is
*               positive.  This time includes time spent in the user
*               subroutines matvec, preconl and preconr.  If timer is
*               negative, the value of elapsed is unchanged.
*
*          rparm(4) = busy
*               On output, this is the CM busy time during the
*               iterative solve, i.e., the time during the iterative
*               solution process that the CM was actually performing
*               work, if the value of timer (see iparm(11) = timer
*               below) is positive.  This time includes time spent in
*               the user subroutines matvec, preconl and preconr.  If
*               timer is negative, the value of busy is unchanged.
*
*          rparm(5) through rparm(20) are reserved for expansion.
*
*     iparm:      integer iparm(20)
*           cmf$ layout iparm(:serial)
*          This is a front end (FE) array of integer input parameters as
*          follows:
*
*          iparm(1) = ntest
*               The stopping test number, in the range 1 to 6,
*               specifying which stopping test is to be used to
*               terminate iteration as follows:
*                    ntest = 1      ||r||          <= zeta
*                    ntest = 2      ||r||/||b||    <= zeta
*                    ntest = 3      ||z||          <= zeta
*                    ntest = 4      ||z||/||Q*b||  <= zeta
*                    ntest = 5      ainf(r)        <= zeta
*                    ntest = 6      ainf(z)        <= zeta
*               where r = residual = b - A*x
*                     b = right-hand side of linear system A*x = b
*                     z = preconditioned residual = Q*(b - A*x)
*                     ||r|| = the Euclidean norm (2-norm) of r
*                     ainf(r) = the infinity norm (max-norm) of r
*
*          iparm(2) = itmax
*               This is the maximum number of iterations allowed.
*
*          iparm(3) = iters
*               On output, this is the actual number of iterations
*               performed by the solver.
*
*          iparm(4) = level
*               This is the level and type of information output by the
*               solver (if iounit > 0).
*                    level = 0     Output error messages only
*                    level = 1     Output input parameters
*                    level = 2     Output stptst at each iteration
*                    level = 3     Output both level = 1 and level = 2
*                                       information
*                     
*          iparm(5) = iounit
*               This is the FORTRAN unit number used for output.
*
*          iparm(6) = iprecon
*               This indicates the type of preconditioning to be
*               performed.
*                    iprecon = 0     No preconditioning
*                    iprecon = 1     Left preconditioning
*                    iprecon = 3     Split preconditioning
*
*          iparm(7) = method
*               This indicates the iterative method to be used.
*                    method = 1     CG      Conjugate Gradients on A
*                    method = 2     CGNR    Conjugate Gradients on AT*A
*                    method = 3     BICG    Bi-Conjugate Gradients on A
*                    method = 4     BICGS   (Bi-)Conjugate Gradients
*                                                Squared on A
*                    method = 5     BICGST  Bi-CGStab on A
*                    method = 6     BICGST2 Bi-CGStab2 on A
*               Where A is the coefficient matrix in the linear system
*               A*x = b.  And AT is the transpose of A.
*
*          iparm(8) = irstrt
*               This indicates whether or not to restart the iterative
*               solver in the case of breakdown for the Bi-CG, Bi-CGS,
*               Bi-CGStab and Bi-CGStab2 algorithms.
*                    irstrt = 0     Do not restart iterative method
*                    irstrt = 1     Allow up to maxrstrt restarts
*                                        of iterative method
*
*          iparm(9) = maxrstrt
*                If irstrt = 1, this indicates the maximum number of
*                restarts allowed for the iterative solver.
*
*          iparm(10) = restarts
*                On output, this is the actual number of restarts of
*                the iterative solver.
*
*          iparm(11) = timer
*                If timer is a non-negative integer, CMStab will return
*                timing information in rparm(3) = elapsed and
*                rparm(4) = busy.  CMStab will use the value of
*                timer as the timer to use.
*
*          iparm(12) through iparm(20) are reserved for expansion.
*
*     ier:        integer
*          This is an error flag.  On output a non-zero value indicates
*          a fatal or warning message was encountered as follows.
*               ier = 0000     No error detected
*
*             Warning Messages:
*               ier = -20      irstrt = 1 but maxrstrt = 0
*               ier = -21      irstrt = 0 but maxrstrt > 0
*
*             Fatal Error Messsages:
*               ier = 1010     Invalid value for ntest
*               ier = 1011     Requested test involving left 
*                                   preconditioning for a method that
*                                   does not involve left
*                                   preconditioning
*               ier = 1020     Non-positive value for itmax
*               ier = 1040     Invalid value for level
*               ier = 1050     Non-positive value for iounit
*               ier = 1060     Invalid value for iprecon
*               ier = 1070     Invalid value for method
*               ier = 1080     Invalid value for irstrt
*               ier = 1090     Non-positive value for maxrstrt
*               ier = 1410     Non-positive value for zeta
*
*               ier = 4001     Failure to converge in itmax iterations
*               ier = 4010     Breakdown in computing direction vector
*                              magnitude in iterative solver; divisor
*                              is zero
*               ier = 4011     Breakdown in computing direction vector
*                              magnitude in iterative solver; divisor
*                              is zero and maxrstrt reached
*               ier = 4020     Iterative solver has stagnated
*
***********************************************************************
*
*     Note on memory usage:  CM Fortran supports the dynamic memory
* allocation extension of Fortran 90.  The subroutines in this package
* will allocate their own workspace arrays as necessary.  The table
* below indicates the amount of memory that will be required for the
* CMStab code to run.  This does not include additional storage
* required by the users routines (matvec, preconl and preconr).
*
* Unpreconditioned systems:
*
*  CG:         real = size(a) + 5*size(x)
*           integer = size(ia)
*
*  CGNR:       real = size(a) + 5*size(x) 
*           integer = size(ia)
*
*  BiCG:       real = size(a) + 7*size(x)
*           integer = size(ia)
*
*  BiCGS:      real = size(a) + 8*size(x)
*           integer = size(ia)
*
*  BiCGStab:   real = size(a) + 7*size(x)
*           integer = size(ia)
*
*  BiCGStab2:  real = size(a) + 13*size(x)
*           integer = size(ia)
*
* Left preconditioned systems:
*
*  CG:         real = size(a) + size(q) + 7*size(x)
*           integer = size(ia) + size(iq)
*
*  CGNR:       real = size(a) + size(q) + 7*size(x)
*           integer = size(ia) + size(iq)
*
*  BiCG:       real = size(a) + size(q) + 9*size(x)
*           integer = size(ia) + size(iq)
*
*  BiCGS:      real = size(a) + size(q) + 8*size(x)
*           integer = size(ia) + size(iq)
*
*  BiCGStab:   real = size(a) + size(q) + 9*size(x)
*           integer = size(ia) + size(iq)
*
*  BiCGStab2:  real = size(a) + size(q) + 19*size(x)
*           integer = size(ia) + size(iq)
*
* Split preconditioned systems:
*
*  CG:         real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
*  CGNR:       real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
*  BiCG:       real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
*  BiCGS:      real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
*  BiCGStab:   real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
*  BiCGStab2:  real = size(a) + size(q) + size(p) +  *size(x)
*           integer = size(ia) + size(iq) + size(ip)
*
***********************************************************************
*
      implicit none

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

      integer a(CMF_sizeof_descriptor),q(CMF_sizeof_descriptor),
     .     p(CMF_sizeof_descriptor)
cmf$ layout a(:serial)
cmf$ layout q(:serial)
cmf$ layout p(:serial)

      integer ia(CMF_sizeof_descriptor),iq(CMF_sizeof_descriptor),
     .     ip(CMF_sizeof_descriptor)
cmf$ layout ia(:serial)
cmf$ layout iq(:serial)
cmf$ layout ip(:serial)
#elif HPF
      real a(:,:,:,:), q(:,:,:,:), p(:,:,:,:)
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)

      integer ia(:,:,:,:), iq(:,:,:,:), ip(:,:,:,:)
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, BLOCK, BLOCK, BLOCK)
#endif

      real x(:,:,:,:), b(:,:,:,:), rparm(20)
      integer iparm(20),ier

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



*
***********************************************************************
*
* Input parameters
*
      real zeta,stptst
      real*4 elapsed, busy, start, tarray(2)

      integer ntest,itmax,iters,level,iounit,iprecon,method,
     .     irstrt,maxrstrt,restarts,timer
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /mthdints/ iprecon,method
      common /rstints/ irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /timing/ timer, elapsed, busy


*
***********************************************************************
*
      external matvec,preconl,preconr

#ifdef CMFORTRAN
      double precision cm_timer_read_cm_busy, cm_timer_read_elapsed
      external cm_timer_read_cm_busy, cm_timer_read_elapsed
#elif HPF
      real*4 etime
      external etime
#endif

*
***********************************************************************
*
      ier = 0

      call code(0,rparm,iparm)

      if ((level.eq.1).or.(level.eq.3)) call output(2)

      call chkparm(ier)
      if (ier.gt.0) return

      if (timer.ge.0) then
#ifdef CMFORTRAN
        call cm_timer_clear(timer)
        call cm_timer_start(timer)
#elif HPF
        start = etime(tarray)
#endif
      endif

*
*
* Select preconditioning and method
*
*

*
*  No preconditioning
*
      if (iprecon .eq. 0) then

*  Conjugate Gradient Method
         if (method .eq. 1) then
            call cg (matvec,a,ia,x,b,ier)

*  Conjugate Gradient Method Applied to the Normal Equations
         elseif (method .eq. 2) then
            call cgnr (matvec,a,ia,x,b,ier)

*  Bi-Conjugate Gradient Method
         elseif (method .eq. 3) then
            call bicg (matvec,a,ia,x,b,ier)

*  Bi-Conjugate Gradient Squared Method
         elseif (method .eq. 4) then
            call bicgs (matvec,a,ia,x,b,ier)

*  BiCGStab Method
         elseif (method .eq. 5) then
            call bicgst (matvec,a,ia,x,b,ier)

*  BiCGStab2 Method
         elseif (method .eq. 6) then
            call bicgst2 (matvec,a,ia,x,b,ier)

*  Invalid Method Selection
         else
            ier=1070

         endif

*
*  Left Preconditioning
*
      elseif (iprecon .eq. 1) then

*  Conjugate Gradient Method
         if (method .eq. 1) then
            call cg_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Conjugate Gradient Method Applied to the Normal Equations
         elseif (method .eq. 2) then
            call cgnr_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Bi-Conjugate Gradient Method
         elseif (method .eq. 3) then
            call bicg_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Bi-Conjugate Gradient Squared Method
         elseif (method .eq. 4) then
            call bicgs_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Bi-CGStab Method
         elseif (method .eq. 5) then
            call bicgst_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Bi-CGStab2 Method
         elseif (method .eq. 6) then
            call bicgst2_left (matvec,preconl,a,ia,x,b,q,iq,ier)

*  Invalid Method Selection
         else
            ier=1070

         endif

*
*  Split Preconditioning
*
      elseif (iprecon .eq. 3) then

*  Conjugate Gradient Method
         if (method .eq. 1) then
            call cg_split (matvec,preconl,preconr,a,ia,x,b,
     .                     q,iq,p,ip,ier)

*  Conjugate Gradient Method Applied to the Normal Equations
         elseif (method .eq. 2) then
            call cgnr_split (matvec,preconl,preconr,a,ia,x,b,
     .                       q,iq,p,ip,ier)

*  Bi-Conjugate Gradient Method
         elseif (method .eq. 3) then
            call bicg_split (matvec,preconl,preconr,a,ia,x,b,
     .                       q,iq,p,ip,ier)

*  Bi-Conjugate Gradient Squared Method
         elseif (method .eq. 4) then
            call bicgs_split (matvec,preconl,preconr,a,ia,x,b,
     .                        q,iq,p,ip,ier)

*  Bi-CGStab Method
         elseif (method .eq. 5) then
            call bicgst_split (matvec,preconl,preconr,a,ia,x,b,
     .                         q,iq,p,ip,ier)

*  Bi-CGStab2 Method
         elseif (method .eq. 6) then
            call bicgst2_split (matvec,preconl,preconr,a,ia,x,b,
     .                          q,iq,p,ip,ier)

         else
            ier=1070

         endif

*
*  Invalid Preconditioning Selected
*
      else 
         ier=1060

      endif


      if (timer.ge.0) then
#ifdef CMFORTRAN
        call cm_timer_stop(timer)
        elapsed = cm_timer_read_elapsed(timer)
        busy = cm_timer_read_cm_busy(timer)
#elif HPF
        elapsed = etime(tarray) - start
        busy = elapsed
#endif
      endif

      if (ier.ne.0) call prnerr(ier)

      call code(1,rparm,iparm)

      return
      end
*
***********************************************************************
*
      subroutine dfault (rparm,iparm)
*
***********************************************************************
*
*     This routine sets up the default values for the elements of the
* rparm and iparm arrays.
*
***********************************************************************
*
* CALLING PARAMETERS:
*
*
*     iparm:      integer (20)
*          This is an array of integer input parameters as specified
*          in the comments of the routine cmstab.
*
*     rparm:      real(20)
*          This is an array of real input parameters as specified
*          in the comments of the routine cmstab.
*
***********************************************************************
*
      implicit none

      real rparm(20)
      integer iparm(20)
#ifdef CMFORTRAN
cmf$ layout rparm(:serial)
cmf$ layout iparm(:serial)
#elif HPF
CHPF$ DISTRIBUTE rparm(*)
CHPF$ DISTRIBUTE iparm(*)
#endif

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

* Choose ||r||/||b|| <= zeta stopping test
      iparm(1) = 2
* Set itmax = 500
      iparm(2) = 500
* Get all the output available	
      iparm(4) = 3
* Output to the terminal
      iparm(5) = 6
* No preconditioning	
      iparm(6) = 0
* Choose simple CG algorithm
      iparm(7) = 1
* Don't allow restarts (doesn't matter with CG algorithm)
      iparm(8) = 0
* Set maxrstrt = 0
      iparm(9) = 0
* Use timer 0 for timing information
      iparm(11) = 0
* iparm(12) through iparm(20) unused, just set to zero
      iparm(12) = 0
      iparm(13) = 0
      iparm(14) = 0
      iparm(15) = 0
      iparm(16) = 0
      iparm(17) = 0
      iparm(18) = 0
      iparm(19) = 0
      iparm(20) = 0

* Set default stopping tolerance to 1.e-6
      rparm(1) = 1.e-6
* Set timing information to 0.
      rparm(3) = 0.
      rparm(4) = 0.
* rparm(5) through rparm(20) unused, just set to zero
      rparm(5) = 0.
      rparm(6) = 0.
      rparm(7) = 0.
      rparm(8) = 0.
      rparm(9) = 0.
      rparm(10) = 0.
      rparm(11) = 0.
      rparm(12) = 0.
      rparm(13) = 0.
      rparm(14) = 0.
      rparm(15) = 0.
      rparm(16) = 0.
      rparm(17) = 0.
      rparm(18) = 0.
      rparm(19) = 0.
      rparm(20) = 0.

      return
      end
*
***********************************************************************
*
      subroutine code(job,rparm,iparm)
*
***********************************************************************
*
*     This routine encodes/decodes the input/output parameters with
* respect to variables in common blocks.
*
***********************************************************************
*
      implicit none

      real rparm(20)
      integer job, iparm(20)
#ifdef CMFORTRAN
cmf$ layout rparm(:serial)
cmf$ layout iparm(:serial)
#elif HPF
CHPF$ DISTRIBUTE rparm(*)
CHPF$ DISTRIBUTE iparm(*)
#endif

      real zeta,stptst
      real*4 elapsed, busy
      integer ntest,itmax,iters,level,iounit,iprecon,method,
     .     irstrt,maxrstrt,restarts,timer
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /mthdints/ iprecon,method
      common /rstints/ irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
      common /timing/ timer,elapsed,busy
*
***********************************************************************
*
* Decode input parameters
      if (job .eq. 0) then
         ntest      = iparm(1)
         itmax      = iparm(2)
         iters      = iparm(3)
         level      = iparm(4)
         iounit     = iparm(5)
         iprecon    = iparm(6)
         method     = iparm(7)
         irstrt     = iparm(8)
         maxrstrt   = iparm(9)
         restarts   = iparm(10)
         timer      = iparm(11)

         zeta       = rparm(1)
         stptst     = rparm(2)
         elapsed    = rparm(3)
         busy       = rparm(4)

* Encode output parameters
      elseif (job .eq. 1) then
         iparm(1)  = ntest
         iparm(2)  = itmax
         iparm(3)  = iters
         iparm(4)  = level
         iparm(5)  = iounit
         iparm(6)  = iprecon
         iparm(7)  = method
         iparm(8)  = irstrt
         iparm(9)  = maxrstrt
         iparm(10) = restarts
         iparm(11) = timer

         rparm(1)  = zeta
         rparm(2)  = stptst
         rparm(3)  = elapsed
         rparm(4)  = busy

      endif

      return
      end
*
***********************************************************************
*
      subroutine chkparm(ier)
*
***********************************************************************
*
*     This routine checks for valid values of the input parameters.
*
***********************************************************************
*
* CALLING PARAMETERS:
*
*
*     ier:       Integer
*          This is an error flag.  On output a non-zero value indicates
*          a fatal or warning message was encountered.  The possible
*          values of ier that this routine can return are:
*               ier = 0000     No error detected
*
*             Warning Messages:
*               ier = -20      irstrt = 1 but maxrstrt = 0
*               ier = -21      irstrt = 0 but maxrstrt > 0
*
*             Fatal Error Messsages:
*               ier = 1010     Invalid value for ntest
*               ier = 1011     Requested test involving left 
*                                   preconditioning for a method that
*                                   does not involve left
*                                   preconditioning
*               ier = 1020     Non-positive value for itmax
*               ier = 1040     Invalid value for level
*               ier = 1050     Invalid value for iounit
*               ier = 1060     Invalid value for iprecon
*               ier = 1070     Invalid value for method
*               ier = 1080     Invalid value for irstrt
*               ier = 1090     Non-positive value for maxrstrt
*               ier = 1410     Non-positive value for zeta
*
***********************************************************************
*
      implicit none

      integer ier

      real zeta,stptst
      integer ntest,itmax,iters,level,iounit,iprecon,method,
     .     irstrt,maxrstrt,restarts
      common /ints/ ntest,itmax,iters
      common /ioints/ level,iounit
      common /mthdints/ iprecon,method
      common /rstints/ irstrt,maxrstrt,restarts
      common /stop_test/ zeta,stptst
*
***********************************************************************
*

*
* Fatal errors
*
      if ((ntest.lt.1).or.(ntest.gt.6)) then
        ier=1010
        call prnerr(ier)
      endif
      if (((iprecon.eq.0).or.(iprecon.eq.2)).and.(ntest.eq.4)) then
        ier=1011
        call prnerr(ier)
      endif
      if (itmax.lt.1) then
        ier=1020
        call prnerr(ier)
      endif
      if ((level.lt.0).or.(level.gt.3)) then
        ier=1040
        call prnerr(ier)
      endif
      if ((iounit.lt.0).or.(iounit.gt.99)) then
        ier=1050
        call prnerr(ier)
      endif
      if ((iprecon.lt.0).or.(iprecon.gt.3)) then
        ier=1060
        call prnerr(ier)
      endif
      if ((method.lt.1).or.(method.gt.6)) then
        ier=1070
        call prnerr(ier)
      endif
      if ((irstrt.lt.0).or.(irstrt.gt.1)) then
        ier=1080
        call prnerr(ier)
      endif
      if (maxrstrt.lt.0) then
        ier=1090
        call prnerr(ier)
      endif

      if (zeta.le.0.)  then
        ier=1410
        call prnerr(ier)
      endif

      if (ier.gt.0) return

*
* Warnings
*
      if ( (irstrt.eq.1).and.(maxrstrt.eq.0) ) then
        ier=-20
        call prnerr(ier)
      endif
      if ( (irstrt.eq.0).and.(maxrstrt.ne.0) ) then
        ier=-21
        call prnerr(ier)
      endif

      return
      end
*
***********************************************************************
*
      subroutine prnerr (ier)
*
***********************************************************************
*
*     This routine prints error messages for errors encountered
*     during the execution of CMStab as specfied by ier.
*
***********************************************************************
*
* CALLING PARAMETERS:
*
*
*     ier:       Integer
*          This is the error flag.  On output a non-zero value
*          indicates a fatal or warning message was encountered as
*          specified in the comments of the routine cmstab.
*
***********************************************************************
*
      implicit none

      integer ier

      integer level,iounit,iprecon,method
      common /ioints/ level,iounit
      common /mthdints/ iprecon,method
*
***********************************************************************
*
      write (iounit,*)
      write (iounit,*) '****************************************'
      write (iounit,*) '*'
      write (iounit,*)

*
* Warning (non-fatal) Error Messages
*
      if (ier .eq. -20) then
         write (iounit,*) '  irstrt = 1 but maxrstrt = 0'

      elseif (ier .eq. -21) then
         write (iounit,*) '  irstrt = 0 but maxrstrt > 0'

*
* Fatal Error Messages
*
      elseif (ier .eq. 1010) then
         write (iounit,*) '  Invalid value for ntest'

      elseif (ier .eq. 1011) then
         write (iounit,*) '  Requested test involving left'
         write (iounit,*) '     preconditioning for a method that'
         write (iounit,*) '     does not involve left preconditioning'

      elseif (ier .eq. 1020) then
         write (iounit,*) '  Non-positive value for itmax'

      elseif (ier .eq. 1040) then
         write (iounit,*) '  Invalid value for level'

      elseif (ier .eq. 1050) then
         write (iounit,*) '  Invalid value for iounit'

      elseif (ier .eq. 1060) then
         write (iounit,*) '  Invalid value for iprecon'

      elseif (ier .eq. 1070) then
         write (iounit,*) '  Invalid value for method'

      elseif (ier .eq. 1080) then
         write (iounit,*) '  Invalid value for irstrt'

      elseif (ier .eq. 1090) then
         write (iounit,*) '  Non-positive value for maxrstrt'

      elseif (ier .eq. 1410) then
         write (iounit,*) '  Non-positive value for zeta'

      elseif (ier .eq. 4001) then
         write (iounit,*) '  Failure to converge in itmax iterations'

      elseif (ier .eq. 4010) then
         write (iounit,*) '  Breakdown in computing direction vector'
         write (iounit,*) '  magnitude in iterative solver; divisor'
         write (iounit,*) '  is zero'

      elseif (ier .eq. 4011) then
         write (iounit,*) '  Breakdown in computing direction vector'
         write (iounit,*) '  magnitude in iterative solver; divisor'
         write (iounit,*) '  is zero and maxrstrt reached'

      elseif (ier .eq. 4020) then
         write (iounit,*) '  Iterative solver has stagnated'

      endif

      write (iounit,*)
      write (iounit,*) '*'
      write (iounit,*) '****************************************'

      return
      end
*
***********************************************************************
*
      subroutine output(job)
*
***********************************************************************
*
*     This subroutine is the output utility for user information
* (i.e., level = 1, 2, or 3).
*
***********************************************************************
*
* CALLING PARAMETERS:
*
*
*     job :      Integer
*          This specifies the type of output requested.
*               job = -1    Output fact that solver was restarted
*               job = 0     Output the header for printing stptst for
*                           each iteration
*               job = 1     Output stoptst for the current iteration
*               job = 2     Output the level = 1 (input paramaters)
*                           information
*
***********************************************************************
*
      implicit none

      integer job

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

      character*10 itmethod(6)
      data itmethod(1) /'CG'/,
     .     itmethod(2) /'CGNR'/,
     .     itmethod(3) /'Bi-CG'/,
     .     itmethod(4) /'Bi-CGS'/,
     .     itmethod(5) /'Bi-CGStab'/,
     .     itmethod(6) /'Bi-CGStab2'/

      character*5 precontype(0:3)
      data precontype(0) /'None'/,
     .     precontype(1) /'Left'/,
     .     precontype(2) /'Right'/,
     .     precontype(3) /'Split'/

      character*13 stoptest(6)
      data stoptest(1) /'||r||'/,
     .     stoptest(2) /'||r||/||b||'/,
     .     stoptest(3) /'||z||'/,
     .     stoptest(4) /'||z||/||Q*b||'/,
     .     stoptest(5) /'ainf(r)'/,
     .     stoptest(6) /'ainf(z)'/

      integer iparm(20)
      character*10 iparms(10)
      data iparms(1) /'ntest'/,
     .     iparms(2) /'itmax'/,
     .     iparms(3) /'iters'/,
     .     iparms(4) /'level'/,
     .     iparms(5) /'iounit'/,
     .     iparms(6) /'iprecon'/,
     .     iparms(7) /'method'/,
     .     iparms(8) /'irstrt'/,
     .     iparms(9) /'maxrstrt'/,
     .     iparms(10) /'restarts'/

      real rparm(20)
      character*10 rparms(2)
      data rparms(1) /'zeta'/,
     .     rparms(2) /'stptst'/

      integer i
*
***********************************************************************
*

*
* Output notification that solver was restarted
*
      if (job .eq. -1) then
         write(iounit,*) '***'
         write(iounit,*) 'Iterative solver restarted; restarts = ',
     .        restarts
         write(iounit,*) '***'

*
* Output header for print stptst for each iteration
*
      elseif (job .eq. 0) then
 10      format(X,'  Method          : ',A)
 20      format(X,'  Preconditioning : ',A)
 30      format(X,'  Stopping Test   : ',A)
 40      format(X,'  zeta            : ',E12.5)
 50      format(X,'     iteration          stptst')
 55      format(X,5X,I9,7X,E12.5)
         write(iounit,*)
         write(iounit,10) itmethod(method)
         write(iounit,20) precontype(iprecon)
         write(iounit,30) stoptest(ntest)
         write(iounit,40) zeta
         write(iounit,*)
         write(iounit,50)

*
* Output stptst for current iteration
*
      elseif (job .eq. 1) then
         write(iounit,55)iters,stptst

*
* Output the input parameters
*
      elseif (job .eq. 2) then
         call code(1,rparm,iparm)
 110     format(X,'     iparm(',I2,') =  ',I12,3X,'(',A10,')')
 120     format(X,'     rparm(',I2,') =  ',E12.5,3X,'(',A10,')')
         write(iounit,*) ' Input parameters '
         do 500 i=1,10
 500        write(iounit,110) i,iparm(i),iparms(i)
            write(iounit,*)
            do 510 i=1,2
 510           write(iounit,120) i,rparm(i),rparms(i)

            endif

            return
            end
*
***********************************************************************
*
      subroutine chkcon (r,z,rnorm,znorm,rinfnorm,zinfnorm,
     .                   bottom,istop,ier)
*
***********************************************************************
*
*     This routine checks the convergence criterion.  If the stopping
* test is satisfied, it signals the calling program that the iterative
* solver has converged.
*
***********************************************************************
*
* CALLING PARAMETERS:
*
*
*     r:
*          This is the array of unpreconditioned resuduals.  If
*          necessary, chkcon will use it in computing stptst.
*
*     z:
*          This is the array of preconditioned resuduals.  If
*          necessary, chkcon will use it in computing stptst.
*
*     rnorm:      real
*          This is the Euclidean norm of the unpreconditioned
*          residuals.  If it is needed as part of the stopping test and
*          has not already been computed, it will be computed by
*          chkcon.
*
*     znorm:      real
*          This is the Euclidean norm of the preconditioned residuals.
*          If it is needed as part of the stopping test and has not
*          already been computed, it will be computed by chkcon.
*
*     rinfnorm:   real
*          This is the infinity norm of the unpreconditioned residuals.
*          If it is needed as part of the stopping test and has not
*          already been computed, it will be computed by chkcon.
*
*     zinfnorm:   real
*          This is the infinity norm of the preconditioned residuals.
*          If it is needed as part of the stopping test and has not
*          already been computed, it will be computed by chkcon.
*
*     bottom:     real
*          This is the denominator of the stopping test and is computed
*          by the calling routine.
*
*     istop:      integer
*          On output, this tells the calling routine whether or not the
*          stopping test is satisfied.
*               istop = 0     Stopping test not satisfied
*               istop = 1     Stopping test satisfied
*
*     ier:        integer
*          This is the error flag.  On output a non-zero value
*          indicates a fatal or warning message was encountered as
*          specified in the comments of the routine cmstab.
*
***********************************************************************
*
      implicit none

      real r(:,:,:,:),z(:,:,:,:)
#ifdef CMFORTRAN
cmf$ layout r(:serial,:news,:news,:news)
cmf$ layout z(:serial,:news,:news,:news)
#elif HPF
CHPF$ DISTRIBUTE r(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE z(*, BLOCK, BLOCK, BLOCK)
#endif

      real rnorm,znorm,rinfnorm,zinfnorm,bottom
      integer istop,ier

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

      real top,mean,stddev
      real stptst_1,stptst_2,stptst_3,stptst_4,stptst_5
      save stptst_1,stptst_2,stptst_3,stptst_4,stptst_5
*
***********************************************************************
*

* ||r|| <= zeta  OR  ||r||/||b|| <= zeta
      if ((ntest .eq. 1) .or. (ntest .eq. 2)) then
         if (rnorm.eq.0.) rnorm = sqrt(sum(r*r))
         top = rnorm

* ||z|| <= zeta  OR  ||z||/||Q*b|| <= zeta
      elseif ((ntest .eq. 3) .or. (ntest .eq. 4)) then
         if (znorm.eq.0.) znorm = sqrt(sum(z*z))
         top = znorm

* ainf(r) <= zeta
      elseif (ntest .eq. 5) then
         if (rinfnorm.eq.0.) rinfnorm = maxval(abs(r))
         top = rinfnorm

* ainf(z) <= zeta
      elseif (ntest .eq. 6) then
         if (zinfnorm.eq.0.) zinfnorm = maxval(abs(z))
         top = zinfnorm

      else
         ier = 1010
         return

      endif

      stptst = top/bottom

      istop = 0
      if (stptst.le.zeta) then
        istop = 1
        return
      endif

* Test for stagnation
      if (iters.eq.0) then
        stptst_1 = 0.
        stptst_2 = 0.
        stptst_3 = 0.
        stptst_4 = 0.
        stptst_5 = 0.
      endif

      mean = (stptst_1 + stptst_2 + stptst_3 + stptst_4 + stptst_5)/5.
      if (abs(mean - stptst)/stptst.le.zeta) then
        stddev = sqrt(( (stptst_1 - mean)*(stptst_1 - mean)
     .                 +(stptst_2 - mean)*(stptst_2 - mean)
     .                 +(stptst_3 - mean)*(stptst_3 - mean)
     .                 +(stptst_4 - mean)*(stptst_4 - mean)
     .                 +(stptst_5 - mean)*(stptst_5 - mean) )*0.25)
        if (stddev/mean.le.zeta) ier = 4020
      endif

      stptst_5 = stptst_4
      stptst_4 = stptst_3
      stptst_3 = stptst_2
      stptst_2 = stptst_1
      stptst_1 = stptst

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