C**********************************************************************
C
C     Copyright (C) 1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file "cpyrit.doc" in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
C     This  file  contains  the  routine  for  the  QMR  algorithm  for
C     symmetric matrices,  using the  three-term recurrence  variant of
C     the Lanczos algorithm without look-ahead.
C
C**********************************************************************
C
      SUBROUTINE SSQMX (NDIM,NLEN,NLIM,VECS,TOL,INFO)
C
C     Purpose:
C     This subroutine uses  the QMR algorithm to solve  linear systems.
C     It runs  the algorithm  to convergence or until  a user-specified
C     limit on the number of  iterations is reached.  It is set  up  to
C     solve symmetric systems starting with identical starting vectors.
C
C     The  code is  set up  to solve  the system  A x = b  with initial
C     guess x_0 = 0.  Here  A x = b  denotes the preconditioned system,
C     and  it  is  connected with the  original system as follows.  Let
C     B y = c be the original unpreconditioned system to be solved, and
C     let y_0 be an arbitrary initial guess for its solution.  Then:
C          A x = b, where  A = M_1^{-1} B M_2^{-1},
C                          x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
C     Here M = M_1 M_2 is the preconditioner.
C
C     To recover the final iterate y_n for  the original system B y = c
C     from the final iterate x_n for the preconditioned system A x = b,
C     set
C               y_n = y_0 + M_2^{-1} x_n.
C
C     The implementation does not have look-ahead, so it is less robust
C     than the full version.
C
C     Parameters:
C     For a description of  the parameters, see the file "ssqmx.doc" in
C     the current directory.
C
C     External routines used:
C     single precision slamch(ch)
C        LAPACK routine, computes machine-related constants.
C     single precision snrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine saxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     single precision sdot(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C     single precision zqmxom(n)
C        User-supplied routine, specifies the QMR scaling factors.
C     subroutine srandn(n,x,seed)
C        Library routine, fills x with random numbers.
C     subroutine srotg(a,b,cos,sin)
C        BLAS-1 routine, computes the Givens rotation which rotates the
C        vector [a; b] into [ sqrt(a**2 + b**2); 0 ].
C
C     Noel M. Nachtigal
C     May 25, 1993
C
C**********************************************************************
C
      INTRINSIC ABS, FLOAT, AMAX1, SQRT, MAX0
      INTRINSIC MOD
      EXTERNAL SLAMCH, SNRM2, SAXPBY, SDOT, SRANDN, SROTG, SSQMXO
      REAL SDOT
      REAL SLAMCH, SNRM2, SSQMXO
C
      INTEGER INFO(4), NDIM, NLEN, NLIM
      REAL VECS(NDIM,7)
      REAL TOL
C
C     Miscellaneous parameters.
C
      REAL DHUN, DONE, DTEN, SZERO
      PARAMETER (DHUN = 1.0E2,DONE = 1.0E0,DTEN = 1.0E1,SZERO = 0.0E0)
C
C     Local variables, permanent.
C
      INTEGER IERR, ISN, ISNM1, ISNM2, IVN, IVNM1, IVNP1, IWN, IWNM1
      SAVE    IERR, ISN, ISNM1, ISNM2, IVN, IVNM1, IVNP1, IWN, IWNM1
      INTEGER IWNP1, N, RETLBL, TF, TRES, VF
      SAVE    IWNP1, N, RETLBL, TF, TRES, VF
      REAL DN, DNM1, DNP1, RHSN, SCSN, SCSNM1, SINN, SINNM1
      SAVE           DN, DNM1, DNP1, RHSN, SCSN, SCSNM1, SINN, SINNM1
      REAL COSN, COSNM1, CSIN, CSINP1, MAXOMG, OMGN
      SAVE             COSN, COSNM1, CSIN, CSINP1, MAXOMG, OMGN
      REAL OMGNP1, R0, RESN, RHON, RHONP1, SCVN, SCVNP1
      SAVE             OMGNP1, R0, RESN, RHON, RHONP1, SCVN, SCVNP1
      REAL SCWN, SCWNP1, TMAX, TMIN, TNRM, UCHK, UNRM
      SAVE             SCWN, SCWNP1, TMAX, TMIN, TNRM, UCHK, UNRM
C
C     Local variables, transient.
C
      INTEGER INIT, REVCOM
      REAL RHN, RHNM1, RHNM2, RHNP1, RHSNP1
      REAL STMP, SCVNM1, SCWNM1
C
C     Initialize some of the permanent variables.
C
      DATA RETLBL /0/
C
C     Check the reverse communication flag to see where to branch.
C        REVCOM   RETLBL      Comment
C           0        0    first call, go to label 10
C           1       30    returning from AXB, go to label 30
C           1       60    returning from AXB, go to label 60
C           2       40    returning from ATXB, go to label 40
C
      REVCOM  = INFO(2)
      INFO(2) = 0
      IF (REVCOM.EQ.0) THEN
         N = 0
         IF (RETLBL.EQ.0) GO TO 10
      ELSE IF (REVCOM.EQ.1) THEN
         IF (RETLBL.EQ.30) THEN
            GO TO 30
         ELSE IF (RETLBL.EQ.60) THEN
            GO TO 60
         END IF
      ELSE IF (REVCOM.EQ.2) THEN
         IF (RETLBL.EQ.40) GO TO 40
      END IF
      IERR = 1
      GO TO 80
C
C     Check whether the inputs are valid.
C
 10   IERR = 0
      IF (NDIM.LT.1)        IERR = 2
      IF (NLEN.LT.1)        IERR = 2
      IF (NLIM.LT.1)        IERR = 2
      IF (NLEN.GT.NDIM)     IERR = 2
      IF (IERR.NE.0) GO TO 80
C
C     Extract from INFO the output units TF and VF, the true residual
C     flag TRES, and the left starting vector flag INIT.
C
      VF   = MAX0(INFO(1),0)
      INIT = VF / 100000
      VF   = VF - INIT * 100000
      TRES = VF / 10000
      VF   = VF - TRES * 10000
      TF   = VF / 100
      VF   = VF - TF * 100
C
C     Extract and check the various tolerances.
C
      TNRM = SLAMCH('E') * DTEN
      TMIN = SQRT(SQRT(SLAMCH('S')))
      TMAX = DONE / TMIN
      IF (TOL.LE.SZERO) TOL = SQRT(SLAMCH('E'))
C
C     Start the trace messages and convergence history.
C
      IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') 0, DONE, DONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, DONE, DONE
C
C     Initialize the wrapped indices.
C
      ISNM1 = 5
      ISN   = ISNM1
      IVN   = 3
      IVNP1 = IVN
      IWN   = 3
      IWNP1 = IWN
C
C     Set x_0 = 0 and compute the norm of the initial residual.
C
      CALL SAXPBY (NLEN,VECS(1,IVN),DONE,VECS(1,2),SZERO,VECS(1,IVN))
      CALL SAXPBY (NLEN,VECS(1,1),SZERO,VECS(1,1),SZERO,VECS(1,1))
      R0 = SNRM2(NLEN,VECS(1,IVN),1)
      IF ((TOL.GE.DONE).OR.(R0.EQ.SZERO)) GO TO 80
C
C     Check whether the auxiliary vector must be supplied.
C
CSYM  IF (INIT.EQ.0) CALL SRANDN (NLEN,VECS(1,IWN),1)
C
C     Compute scale factors and check for invariant subspaces.
C
      SCVNP1 = R0
CSYM  SCWNP1 = SNRM2(NLEN,VECS(1,IWN),1)
      SCWNP1 = SCVNP1
      IF (SCVNP1.LT.TNRM) IERR = IERR + 16
      IF (IERR.NE.0) GO TO 80
      DNP1 = SDOT(NLEN,VECS(1,IWN),1,VECS(1,IVN),1) / ( SCVNP1 * SCWNP1
     $)
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         STMP = DONE / SCVNP1
         CALL SAXPBY (NLEN,VECS(1,IVN),STMP,VECS(1,IVN),SZERO,VECS(1,IVN
     $))
         SCVNP1 = DONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         STMP = DONE / SCWNP1
CSYM     CALL SAXPBY (NLEN,VECS(1,IWN),STMP,VECS(1,IWN),SZERO,VECS(1,IWN))
         SCWNP1 = DONE
      END IF
      RHONP1 = SCVNP1
      CSINP1 = SCWNP1
      SCVNP1 = DONE / SCVNP1
      SCWNP1 = DONE / SCWNP1
C
C     Initialize the variables.
C
      N      = 1
      DN     = DONE
      COSN   = DONE
      RESN   = DONE
      COSNM1 = DONE
      OMGN   = SZERO
      SCSN   = SZERO
      SCVN   = SZERO
      SCWN   = SZERO
      SCSNM1 = SZERO
      SINN   = SZERO
      SINNM1 = SZERO
      OMGNP1 = SSQMXO(N)
      RHSN   = OMGNP1 * R0
      MAXOMG = DONE / OMGNP1
C
C     This is one step of the classical Lanczos algorithm.
C
 20   IVNM1 = IVN
      IVN   = IVNP1
      IVNP1 = MOD(N,2) + 3
      IWNM1 = IWN
      IWN   = IWNP1
      IWNP1 = MOD(N,2) + 3
C
C     Check whether D_n is nonsingular.
C
      DNM1 = DN
      DN   = DNP1
      IF (DN.EQ.SZERO) THEN
         IERR = 8
         GO TO 80
      END IF
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,IVN),VECS(1,7))
C
      INFO(2) = 1
      INFO(3) = IVN
      INFO(4) = 7
      RETLBL  = 30
      RETURN
 30   RETLBL = 0
C
C     Compute H_{n-1,n} and build part of the vector v_{n+1}.
C
      SCVNM1 = SCVN
      CSIN   = CSINP1
      SCVN   = SCVNP1
      STMP   = CSIN * DN / DNM1 * SCVNM1 / SCVN
      CALL SAXPBY (NLEN,VECS(1,IVNP1),DONE,VECS(1,7),-STMP,VECS(1,IVNM1)
     $)
C
C     Have the caller carry out ATXB, then return here.
C        CALL ATXB (VECS(1,IWN),VECS(1,7))
C
      INFO(2) = 2
      INFO(3) = IWN
      INFO(4) = 7
      RETLBL  = 40
CSYM  RETURN
 40   RETLBL = 0
C
C     Build part of the vector w_{n+1}.
C
      SCWNM1 = SCWN
      RHON   = RHONP1
      SCWN   = SCWNP1
      STMP   = RHON * DN / DNM1 * SCWNM1 / SCWN
CSYM  CALL SAXPBY (NLEN,VECS(1,IWNP1),DONE,VECS(1,7),-STMP,VECS(1,IWNM1))
C
C     Compute H_{nn} and finish the new vectors.
C
      RHN = SCVN * SCWN * SDOT(NLEN,VECS(1,IWN),1,VECS(1,IVNP1),1) / DN
      CALL SAXPBY (NLEN,VECS(1,IVNP1),DONE,VECS(1,IVNP1),-RHN,VECS(1,IVN
     $))
CSYM  CALL SAXPBY (NLEN,VECS(1,IWNP1),DONE,VECS(1,IWNP1),-RHN,VECS(1,IWN))
C
C     Compute scale factors and check for invariant subspaces.
C
      IERR   = 0
      SCVNP1 = SNRM2(NLEN,VECS(1,IVNP1),1)
CSYM  SCWNP1 = SNRM2(NLEN,VECS(1,IWNP1),1)
      SCWNP1 = SCVNP1
      RHONP1 = SCVN * SCVNP1
      CSINP1 = SCWN * SCWNP1
      RHNP1  = RHONP1
      IF (SCVNP1.LT.TNRM) IERR = IERR + 16
      IF (IERR.NE.0) GO TO 50
      DNP1   = SDOT(NLEN,VECS(1,IWNP1),1,VECS(1,IVNP1),1) / ( SCVNP1 *
     $SCWNP1 )
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         STMP = DONE / SCVNP1
         CALL SAXPBY (NLEN,VECS(1,IVNP1),STMP,VECS(1,IVNP1),SZERO,VECS(1
     $,IVNP1))
         SCVNP1 = DONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         STMP = DONE / SCWNP1
CSYM     CALL SAXPBY (NLEN,VECS(1,IWNP1),STMP,VECS(1,IWNP1),SZERO,VECS(1,IWNP1))
         SCWNP1 = DONE
      END IF
      SCVNP1 = DONE / SCVNP1
      SCWNP1 = DONE / SCWNP1
C
C     The QMR code starts here.
C     Multiply the new column by the previous omegas.
C     Get the next scaling factor omega(i) and update MAXOMG.
C
 50   RHNM1  = CSIN * DN * OMGN / DNM1
      OMGN   = OMGNP1
      RHN    = OMGN * RHN
      OMGNP1 = SSQMXO(N+1)
      RHNP1  = OMGNP1 * RHNP1
      MAXOMG = AMAX1(MAXOMG,DONE/OMGN)
C
C     Apply the previous rotations.
C
      RHNM2  = SINNM1 * RHNM1
      RHNM1  = COSNM1 * RHNM1
      COSNM1 = COSN
      SINNM1 = SINN
      STMP   = RHNM1
      RHNM1  =  COSNM1 * STMP + SINNM1 * RHN
      RHN    = -SINNM1 * STMP + COSNM1 * RHN
C
C     Compute the rotation for the last element (this also applies it).
C
      CALL SROTG (RHN,RHNP1,COSN,SINN)
C
C     Apply the new rotation to the right-hand side vector.
C
      RHSNP1 = -SINN * RHSN
      RHSN   =  COSN * RHSN
C
C     Compute the next search direction s_n.
C
      ISNM2  = ISNM1
      ISNM1  = ISN
      ISN    = MOD(N-1,2) + 5
      STMP   = SCSNM1 * RHNM2 / SCVN
      CALL SAXPBY (NLEN,VECS(1,ISN),DONE,VECS(1,IVN),-STMP,VECS(1,ISNM2)
     $)
      SCSNM1 = SCSN
      STMP   = SCSNM1 * RHNM1 / SCVN
      CALL SAXPBY (NLEN,VECS(1,ISN),DONE,VECS(1,ISN),-STMP,VECS(1,ISNM1)
     $)
      SCSN   = SCVN / RHN
C
C     Compute the new QMR iterate, then scale the search direction.
C
      STMP = SCSN * RHSN
      CALL SAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),STMP,VECS(1,ISN))
      STMP = ABS(SCSN)
      IF ((STMP.GE.TMAX).OR.(STMP.LE.TMIN)) THEN
         CALL SAXPBY (NLEN,VECS(1,ISN),SCSN,VECS(1,ISN),SZERO,VECS(1,ISN
     $))
         SCSN = DONE
      END IF
C
C     Compute the residual norm upper bound.
C     If the scaled upper bound is within one order of magnitude of the
C     target convergence norm, compute the true residual norm.
C
      RHSN = RHSNP1
      UNRM = SQRT(FLOAT(N+1)) * MAXOMG * ABS(RHSNP1) / R0
      UCHK = UNRM
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 70
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,1),VECS(1,7))
C
      INFO(2) = 1
      INFO(3) = 1
      INFO(4) = 7
      RETLBL  = 60
      RETURN
 60   RETLBL = 0
      CALL SAXPBY (NLEN,VECS(1,7),DONE,VECS(1,2),-DONE,VECS(1,7))
      RESN = SNRM2(NLEN,VECS(1,7),1) / R0
      UCHK = RESN
C
C     Output the trace messages and convergence history.
C
 70   IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') N, UNRM, RESN
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') N, UNRM, RESN
C
C     Check for convergence or termination.  Stop if:
C         1. algorithm converged;
C         2. there is an error condition;
C         3. the residual norm upper bound is smaller than the computed
C     residual norm by a factor of at least 100;
C         4. algorithm exceeded the iterations limit.
C
      IF (RESN.LE.TOL) THEN
         IERR = 0
         GO TO 80
      ELSE IF (IERR.NE.0) THEN
         GO TO 80
      ELSE IF (UNRM.LT.UCHK/DHUN) THEN
         IERR = 4
         GO TO 80
      ELSE IF (N.GE.NLIM) THEN
         IERR = 4
         GO TO 80
      END IF
C
C     Update the running counter.
C
      N = N + 1
      GO TO 20
C
C     Done.
C
 80   NLIM    = N
      RETLBL  = 0
      INFO(1) = IERR
C
      RETURN
      END
C
C**********************************************************************
C
      REAL FUNCTION SSQMXO (I)
C
C     Purpose:
C     Returns the scaling parameter OMEGA(I).
C
C     Parameters:
C     I = the index of the parameter OMEGA (input).
C
C     Noel M. Nachtigal
C     March 30, 1993
C
C**********************************************************************
C
      INTEGER I
C
      SSQMXO = 1.0E0
C
      RETURN
      END
C
C**********************************************************************