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 routines for the QMR algorithm for
C symmetric matrices, using the coupled two-term recurrence variant
C of the Lanczos algorithm without look-ahead.
C
C**********************************************************************
C
SUBROUTINE SSCPX (NDIM,NLEN,NLIM,VECS,TOL,INFO)
C
C Purpose:
C This subroutine uses the QMR algorithm based on the coupled two-
C term variant of the Lanczos process without look-ahead to solve
C linear systems. It runs the algorithm to convergence or until a
C user-specified limit on the number of iterations is reached. It
C is set up to solve symmetric systems starting with identical
C 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 algorithm was first described in the RIACS Technical Report
C 92.15, `An Implementation of the QMR Method Based on Coupled Two-
C Term Recurrences`, June 1992. This implementation does not have
C look-ahead, so it is less robust than the full version.
C
C Parameters:
C For a description of the parameters, see the file "sscpx.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 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 single precision sscpxo(n)
C User-supplied routine, specifies the QMR scaling factors.
C
C Noel M. Nachtigal
C March 30, 1993
C
C**********************************************************************
C
INTRINSIC ABS, FLOAT, AMAX1, SQRT, MAX0
EXTERNAL SLAMCH, SNRM2, SAXPBY, SDOT, SROTG, SSCPXO
REAL SDOT
REAL SLAMCH, SNRM2, SSCPXO
C
INTEGER INFO(4), NDIM, NLEN, NLIM
REAL VECS(NDIM,6)
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, N, RETLBL, TF, TRES, VF
SAVE IERR, N, RETLBL, TF, TRES, VF
REAL DNN, ENN, SCS, SINN, RHSN
SAVE DNN, ENN, SCS, SINN, RHSN
REAL COSN, LNP1N, MAXOMG, OMG, R0, SCPN
SAVE COSN, LNP1N, MAXOMG, OMG, R0, SCPN
REAL SCV, RESN, TMAX, TMIN, TNRM, UCHK, UNRM
SAVE SCV, RESN, TMAX, TMIN, TNRM, UCHK, UNRM
C
C Local variables, transient.
C
INTEGER I, REVCOM
REAL LNN, RHN, RHNM1, RHNP1, RHSNP1, UNM1N, STMP
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 40 returning from AXB, 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.40) THEN
GO TO 40
END IF
END IF
IERR = 1
GO TO 60
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 60
C
C Extract from INFO the output units TF and VF, the true residual
C flag TRES.
C
VF = MAX0(INFO(1),0)
I = VF / 100000
VF = VF - I * 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 Set x_0 = 0 and compute the norm of the initial residual.
C
CALL SAXPBY (NLEN,VECS(1,3),DONE,VECS(1,2),SZERO,VECS(1,3))
CALL SAXPBY (NLEN,VECS(1,1),SZERO,VECS(1,1),SZERO,VECS(1,1))
R0 = SNRM2(NLEN,VECS(1,3),1)
IF ((TOL.GE.DONE).OR.(R0.EQ.SZERO)) GO TO 60
C
C Initialize the variables.
C
N = 1
SCV = R0
ENN = DONE
COSN = DONE
RESN = DONE
SCPN = DONE
SCS = SZERO
SINN = SZERO
LNP1N = SZERO
OMG = SSCPXO(N)
RHSN = OMG * R0
MAXOMG = DONE / OMG
C
C This is one step of the coupled two-term Lanczos algorithm.
C Check whether E_n is nonsingular.
C
20 IF (ENN.EQ.SZERO) THEN
IERR = 8
GO TO 60
END IF
C
C Compute scale factor for the vector w_{n}.
C Check for invariant subspaces, and scale the vectors if needed.
C
IERR = 0
IF (SCPN*SCV.LT.TNRM) IERR = IERR + 16
IF (IERR.NE.0) GO TO 60
DNN = SDOT(NLEN,VECS(1,3),1,VECS(1,3),1) / ( SCV**2 )
IF ((SCV.GE.TMAX).OR.(SCV.LE.TMIN)) THEN
STMP = DONE / SCV
CALL SAXPBY (NLEN,VECS(1,3),STMP,VECS(1,3),SZERO,VECS(1,3))
SCV = DONE
END IF
SCV = DONE / SCV
C
C Build the vector p_n.
C
UNM1N = DNN * LNP1N / ENN
STMP = UNM1N * SCPN / SCV
CALL SAXPBY (NLEN,VECS(1,4),DONE,VECS(1,3),-STMP,VECS(1,4))
SCPN = SCV
C
C Check whether D_n is nonsingular.
C
IF (DNN.EQ.SZERO) THEN
IERR = 8
GO TO 60
END IF
C
C Have the caller carry out AXB, then return here.
C CALL AXB (VECS(1,4),VECS(1,6))
C
INFO(2) = 1
INFO(3) = 4
INFO(4) = 6
RETLBL = 30
RETURN
C
C Compute p_n^T A p_n.
C
30 ENN = SDOT(NLEN,VECS(1,4),1,VECS(1,6),1) * ( SCPN**2 )
C
C Build the vector v_{n+1}.
C
LNN = ENN / DNN
CALL SAXPBY (NLEN,VECS(1,3),DONE,VECS(1,6),-LNN,VECS(1,3))
C
C Compute scale factor for the vector v_{n+1}.
C
SCV = SNRM2(NLEN,VECS(1,3),1)
LNP1N = SCPN * SCV
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
RHN = OMG * LNN
OMG = SSCPXO(N+1)
RHNP1 = OMG * LNP1N
MAXOMG = AMAX1(MAXOMG,DONE/OMG)
C
C Apply the previous rotation.
C
RHNM1 = SINN * RHN
RHN = COSN * 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_i.
C
STMP = RHNM1 * SCS / SCPN
CALL SAXPBY (NLEN,VECS(1,5),DONE,VECS(1,4),-STMP,VECS(1,5))
C
C Compute the new QMR iterate, then scale the search direction.
C
SCS = SCPN / RHN
STMP = SCS * RHSN
CALL SAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),STMP,VECS(1,5))
IF ((ABS(SCS).GE.TMAX).OR.(ABS(SCS).LE.TMIN)) THEN
CALL SAXPBY (NLEN,VECS(1,5),SCS,VECS(1,5),SZERO,VECS(1,5))
SCS = 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 50
C
C Have the caller carry out AXB, then return here.
C CALL AXB (VECS(1,1),VECS(1,6))
C
INFO(2) = 1
INFO(3) = 1
INFO(4) = 6
RETLBL = 40
RETURN
40 CALL SAXPBY (NLEN,VECS(1,6),DONE,VECS(1,2),-DONE,VECS(1,6))
RESN = SNRM2(NLEN,VECS(1,6),1) / R0
UCHK = RESN
C
C Output the convergence history.
C
50 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 60
ELSE IF (IERR.NE.0) THEN
GO TO 60
ELSE IF (UNRM.LT.UCHK/DHUN) THEN
IERR = 4
GO TO 60
ELSE IF (N.GE.NLIM) THEN
IERR = 4
GO TO 60
END IF
C
C Update the running counter.
C
N = N + 1
GO TO 20
C
C Done.
C
60 NLIM = N
RETLBL = 0
INFO(1) = IERR
C
RETURN
END
C
C**********************************************************************
C
REAL FUNCTION SSCPXO (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
SSCPXO = 1.0E0
C
RETURN
END
C
C**********************************************************************