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 CSQMX (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 "csqmx.doc" in
C     the current directory.
C
C     External routines used:
C     real slamch(ch)
C        LAPACK routine, computes machine-related constants.
C     real scnrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine caxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     real cdotu(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C     real zqmxom(n)
C        User-supplied routine, specifies the QMR scaling factors.
C     subroutine crandn(n,x,seed)
C        Library routine, fills x with random numbers.
C     subroutine crotg(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 CABS, ABS, FLOAT, CMPLX, CONJG, AMAX1, SQRT, MAX0
      INTRINSIC MOD
      EXTERNAL SLAMCH, SCNRM2, CAXPBY, CDOTU, CRANDN, CROTG, CSQMXO
      COMPLEX CDOTU
      REAL SLAMCH, SCNRM2, CSQMXO
C
      INTEGER INFO(4), NDIM, NLEN, NLIM
      COMPLEX VECS(NDIM,7)
      REAL TOL
C
C     Miscellaneous parameters.
C
      COMPLEX CONE, CZERO
      PARAMETER (CONE = (1.0E0,0.0E0),CZERO = (0.0E0,0.0E0))
      REAL SHUN, SONE, STEN, SZERO
      PARAMETER (SHUN = 1.0E2,SONE = 1.0E0,STEN = 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
      COMPLEX 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
      COMPLEX RHN, RHNM1, RHNM2, RHNP1, RHSNP1, CTMP
      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') * STEN
      TMIN = SQRT(SQRT(SLAMCH('S')))
      TMAX = SONE / 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, SONE, SONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, SONE, SONE
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 CAXPBY (NLEN,VECS(1,IVN),CONE,VECS(1,2),CZERO,VECS(1,IVN))
      CALL CAXPBY (NLEN,VECS(1,1),CZERO,VECS(1,1),CZERO,VECS(1,1))
      R0 = SCNRM2(NLEN,VECS(1,IVN),1)
      IF ((TOL.GE.SONE).OR.(R0.EQ.SZERO)) GO TO 80
C
C     Check whether the auxiliary vector must be supplied.
C
CSYM  IF (INIT.EQ.0) CALL CRANDN (NLEN,VECS(1,IWN),1)
C
C     Compute scale factors and check for invariant subspaces.
C
      SCVNP1 = R0
CSYM  SCWNP1 = SCNRM2(NLEN,VECS(1,IWN),1)
      SCWNP1 = SCVNP1
      IF (SCVNP1.LT.TNRM) IERR = IERR + 16
      IF (IERR.NE.0) GO TO 80
      DNP1 = CDOTU(NLEN,VECS(1,IWN),1,VECS(1,IVN),1) / ( SCVNP1 * SCWNP1
     $ )
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / SCVNP1,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IVN),CTMP,VECS(1,IVN),CZERO,VECS(1,IVN
     $))
         SCVNP1 = SONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / SCWNP1,SZERO)
CSYM     CALL CAXPBY (NLEN,VECS(1,IWN),CTMP,VECS(1,IWN),CZERO,VECS(1,IWN))
         SCWNP1 = SONE
      END IF
      RHONP1 = SCVNP1
      CSINP1 = SCWNP1
      SCVNP1 = SONE / SCVNP1
      SCWNP1 = SONE / SCWNP1
C
C     Initialize the variables.
C
      N      = 1
      DN     = CONE
      COSN   = SONE
      RESN   = SONE
      COSNM1 = SONE
      OMGN   = SZERO
      SCSN   = SZERO
      SCVN   = SZERO
      SCWN   = SZERO
      SCSNM1 = SZERO
      SINN   = CZERO
      SINNM1 = CZERO
      OMGNP1 = CSQMXO(N)
      RHSN   = OMGNP1 * R0
      MAXOMG = SONE / 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 (CABS(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
      CTMP   = CSIN * DN / DNM1 * SCVNM1 / SCVN
      CALL CAXPBY (NLEN,VECS(1,IVNP1),CONE,VECS(1,7),-CTMP,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
      CTMP   = RHON * DN / DNM1 * SCWNM1 / SCWN
CSYM  CALL CAXPBY (NLEN,VECS(1,IWNP1),CONE,VECS(1,7),-CTMP,VECS(1,IWNM1))
C
C     Compute H_{nn} and finish the new vectors.
C
      RHN = SCVN * SCWN * CDOTU(NLEN,VECS(1,IWN),1,VECS(1,IVNP1),1) / DN
      CALL CAXPBY (NLEN,VECS(1,IVNP1),CONE,VECS(1,IVNP1),-RHN,VECS(1,IVN
     $))
CSYM  CALL CAXPBY (NLEN,VECS(1,IWNP1),CONE,VECS(1,IWNP1),-RHN,VECS(1,IWN))
C
C     Compute scale factors and check for invariant subspaces.
C
      IERR   = 0
      SCVNP1 = SCNRM2(NLEN,VECS(1,IVNP1),1)
CSYM  SCWNP1 = SCNRM2(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   = CDOTU(NLEN,VECS(1,IWNP1),1,VECS(1,IVNP1),1) / ( SCVNP1 * 
     $SCWNP1 )
      IF ((SCVNP1.GE.TMAX).OR.(SCVNP1.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / SCVNP1,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IVNP1),CTMP,VECS(1,IVNP1),CZERO,VECS(1
     $,IVNP1))
         SCVNP1 = SONE
      END IF
      IF ((SCWNP1.GE.TMAX).OR.(SCWNP1.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / SCWNP1,SZERO)
CSYM     CALL CAXPBY (NLEN,VECS(1,IWNP1),CTMP,VECS(1,IWNP1),CZERO,VECS(1,IWNP1))
         SCWNP1 = SONE
      END IF
      SCVNP1 = SONE / SCVNP1
      SCWNP1 = SONE / 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 = CSQMXO(N+1)
      RHNP1  = OMGNP1 * RHNP1
      MAXOMG = AMAX1(MAXOMG,SONE/OMGN)
C
C     Apply the previous rotations.
C
      RHNM2  = SINNM1 * RHNM1
      RHNM1  = COSNM1 * RHNM1
      COSNM1 = COSN
      SINNM1 = SINN
      CTMP   = RHNM1
      RHNM1  =  COSNM1 * CTMP + SINNM1 * RHN
      RHN    = -CONJG(SINNM1) * CTMP + COSNM1 * RHN
C
C     Compute the rotation for the last element (this also applies it).
C
      CALL CROTG (RHN,RHNP1,COSN,SINN)
C
C     Apply the new rotation to the right-hand side vector.
C
      RHSNP1 = -CONJG(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
      CTMP   = SCSNM1 * RHNM2 / SCVN
      CALL CAXPBY (NLEN,VECS(1,ISN),CONE,VECS(1,IVN),-CTMP,VECS(1,ISNM2)
     $)
      SCSNM1 = SCSN
      CTMP   = SCSNM1 * RHNM1 / SCVN
      CALL CAXPBY (NLEN,VECS(1,ISN),CONE,VECS(1,ISN),-CTMP,VECS(1,ISNM1)
     $)
      SCSN   = SCVN / RHN
C
C     Compute the new QMR iterate, then scale the search direction.
C
      CTMP = SCSN * RHSN
      CALL CAXPBY (NLEN,VECS(1,1),CONE,VECS(1,1),CTMP,VECS(1,ISN))
      STMP = CABS(SCSN)
      IF ((STMP.GE.TMAX).OR.(STMP.LE.TMIN)) THEN
         CALL CAXPBY (NLEN,VECS(1,ISN),SCSN,VECS(1,ISN),CZERO,VECS(1,ISN
     $))
         SCSN = CONE
      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 * CABS(RHSNP1) / R0
      UCHK = UNRM
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.STEN).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 CAXPBY (NLEN,VECS(1,7),CONE,VECS(1,2),-CONE,VECS(1,7))
      RESN = SCNRM2(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/SHUN) 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 CSQMXO (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
      CSQMXO = 1.0E0
C
      RETURN
      END
C
C**********************************************************************