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 unsymmetric matrices, using the coupled two-term recurrence C variant of the Lanczos algorithm without look-ahead. C C********************************************************************** C SUBROUTINE CUCPX (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. 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 "cucpx.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 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 real cucpxo(n) C User-supplied routine, specifies the QMR scaling factors. C C Noel M. Nachtigal C March 30, 1993 C C********************************************************************** C INTRINSIC CABS, ABS, FLOAT, CMPLX, CONJG, AMAX1, SQRT, MAX0 EXTERNAL SLAMCH, SCNRM2, CAXPBY, CDOTU, CRANDN, CROTG, CUCPXO COMPLEX CDOTU REAL SLAMCH, SCNRM2, CUCPXO C INTEGER INFO(4), NDIM, NLEN, NLIM COMPLEX VECS(NDIM,8) 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, N, RETLBL, TF, TRES, VF SAVE IERR, N, RETLBL, TF, TRES, VF COMPLEX DNN, ENN, SCS, SINN, RHSN SAVE DNN, ENN, SCS, SINN, RHSN REAL COSN, GAMN, LNP1N, MAXOMG, OMG, R0, SCPN, SCQN SAVE COSN, GAMN, LNP1N, MAXOMG, OMG, R0, SCPN, SCQN REAL SCV, RESN, TMAX, TMIN, TNRM, UCHK, UNRM SAVE SCV, RESN, TMAX, TMIN, TNRM, UCHK, UNRM C C Local variables, transient. C INTEGER INIT, REVCOM COMPLEX LNN, RHN, RHNM1, RHNP1, RHSNP1, UNM1N, CTMP REAL GAMNM1, SCW 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 50 returning from AXB, go to label 50 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.50) THEN GO TO 50 END IF ELSE IF (REVCOM.EQ.2) THEN IF (RETLBL.EQ.40) GO TO 40 END IF IERR = 1 GO TO 70 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 70 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 Set x_0 = 0 and compute the norm of the initial residual. C CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,2),CZERO,VECS(1,3)) CALL CAXPBY (NLEN,VECS(1,1),CZERO,VECS(1,1),CZERO,VECS(1,1)) R0 = SCNRM2(NLEN,VECS(1,3),1) IF ((TOL.GE.SONE).OR.(R0.EQ.SZERO)) GO TO 70 C C Check whether the auxiliary vector must be supplied. C IF (INIT.EQ.0) CALL CRANDN (NLEN,VECS(1,7),1) C C Initialize the variables. C N = 1 SCV = R0 ENN = CONE COSN = SONE GAMN = SONE RESN = SONE SCPN = SONE SCQN = SONE SCS = CZERO SINN = CZERO LNP1N = SZERO OMG = CUCPXO(N) RHSN = OMG * R0 MAXOMG = SONE / OMG C C This is one step of the coupled two-term Lanczos algorithm. C Check whether E_n is nonsingular. C 20 IF (CABS(ENN).EQ.SZERO) THEN IERR = 8 GO TO 70 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 SCW = SCNRM2(NLEN,VECS(1,7),1) IF (SCPN*SCV.LT.TNRM) IERR = IERR + 16 IF (SCQN*SCW.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GO TO 70 GAMNM1 = GAMN GAMN = GAMN * SCPN / SCQN * SCV / SCW DNN = CDOTU(NLEN,VECS(1,3),1,VECS(1,7),1) / ( SCV * SCW ) IF ((SCV.GE.TMAX).OR.(SCV.LE.TMIN)) THEN CTMP = CMPLX(SONE / SCV,SZERO) CALL CAXPBY (NLEN,VECS(1,3),CTMP,VECS(1,3),CZERO,VECS(1,3)) SCV = SONE END IF IF ((SCW.GE.TMAX).OR.(SCW.LE.TMIN)) THEN CTMP = CMPLX(SONE / SCW,SZERO) CALL CAXPBY (NLEN,VECS(1,7),CTMP,VECS(1,7),CZERO,VECS(1,7)) SCW = SONE END IF SCV = SONE / SCV SCW = SONE / SCW C C Build the vectors p_n and q_n. C UNM1N = DNN * LNP1N * GAMNM1 / ( GAMN * ENN ) CTMP = UNM1N * SCPN / SCV CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,3),-CTMP,VECS(1,4)) CTMP = UNM1N * SCQN / SCW * GAMN / GAMNM1 CALL CAXPBY (NLEN,VECS(1,8),CONE,VECS(1,7),-CTMP,VECS(1,8)) SCPN = SCV SCQN = SCW C C Check whether D_n is nonsingular. C IF (CABS(DNN).EQ.SZERO) THEN IERR = 8 GO TO 70 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 q_n^T A p_n. C 30 ENN = SCPN * SCQN * CDOTU(NLEN,VECS(1,8),1,VECS(1,6),1) C C Build the vector v_{n+1}. C LNN = ENN / DNN CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,6),-LNN,VECS(1,3)) C C Have the caller carry out ATXB, then return here. C CALL ATXB (VECS(1,8),VECS(1,6)) C INFO(2) = 2 INFO(3) = 8 INFO(4) = 6 RETLBL = 40 RETURN C C Build the vector w_{n+1}. C 40 LNN = ENN / DNN CALL CAXPBY (NLEN,VECS(1,7),CONE,VECS(1,6),-LNN,VECS(1,7)) C C Compute scale factor for the vector v_{n+1}. C SCV = SCNRM2(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 = CUCPXO(N+1) RHNP1 = OMG * LNP1N MAXOMG = AMAX1(MAXOMG,SONE/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 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_i. C CTMP = RHNM1 * SCS / SCPN CALL CAXPBY (NLEN,VECS(1,5),CONE,VECS(1,4),-CTMP,VECS(1,5)) C C Compute the new QMR iterate, then scale the search direction. C SCS = SCPN / RHN CTMP = SCS * RHSN CALL CAXPBY (NLEN,VECS(1,1),CONE,VECS(1,1),CTMP,VECS(1,5)) IF ((CABS(SCS).GE.TMAX).OR.(CABS(SCS).LE.TMIN)) THEN CALL CAXPBY (NLEN,VECS(1,5),SCS,VECS(1,5),CZERO,VECS(1,5)) SCS = 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 60 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 = 50 RETURN 50 CALL CAXPBY (NLEN,VECS(1,6),CONE,VECS(1,2),-CONE,VECS(1,6)) RESN = SCNRM2(NLEN,VECS(1,6),1) / R0 UCHK = RESN C C Output the convergence history. C 60 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 70 ELSE IF (IERR.NE.0) THEN GO TO 70 ELSE IF (UNRM.LT.UCHK/SHUN) THEN IERR = 4 GO TO 70 ELSE IF (N.GE.NLIM) THEN IERR = 4 GO TO 70 END IF C C Update the running counter. C N = N + 1 GO TO 20 C C Done. C 70 NLIM = N RETLBL = 0 INFO(1) = IERR C RETURN END C C********************************************************************** C REAL FUNCTION CUCPXO (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 CUCPXO = 1.0E0 C RETURN END C C**********************************************************************