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-from-BCG algorithm for C unsymmetric matrices. C C********************************************************************** C SUBROUTINE CUQBG (NDIM,NLEN,NLIM,VECS,TOL,INFO) C C Purpose: C This subroutine uses the QMR-from-BCG algorithm to solve linear C systems. It runs the algorithm to convergence or until a user- C specified iteration limit 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 described in the RIACS Technical Report 91.26, C A Quasi-Minimal Residual Squared Algorithm for Non-Hermitian C Linear Systems, by Freund and Szeto, December 1991. C C Parameters: C For a description of the parameters, see the file "cuqbg.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 C Noel M. Nachtigal C March 9, 1992 C C********************************************************************** C INTRINSIC CABS, FLOAT, SQRT, MAX0 EXTERNAL SLAMCH, SCNRM2, CAXPBY, CDOTU, CRANDN COMPLEX CDOTU REAL SLAMCH, SCNRM2 C INTEGER INFO(4), NLEN, NDIM, 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 ALPHA, RHO SAVE ALPHA, RHO REAL QPRD, R0, RESN, TAU, UNRM, VAR SAVE QPRD, R0, RESN, TAU, UNRM, VAR C C Local variables, transient. C INTEGER INIT, REVCOM COMPLEX BETA, CTMP REAL COS, 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 2 60 returning from ATXB, go to label 60 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 ELSE IF (REVCOM.EQ.2) THEN IF (RETLBL.EQ.60) GO TO 60 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 Check the covergence tolerance. C 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 QPRD = R0 RESN = SONE VAR = SZERO TAU = R0 * R0 RHO = CDOTU(NLEN,VECS(1,7),1,VECS(1,3),1) CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,3),CZERO,VECS(1,4)) CALL CAXPBY (NLEN,VECS(1,8),CONE,VECS(1,7),CZERO,VECS(1,8)) C C This is one step of the QMR-BCG algorithm. C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,4),VECS(1,6)) C 20 INFO(2) = 1 INFO(3) = 4 INFO(4) = 6 RETLBL = 30 RETURN C C Compute \sigma and check for breakdown. C 30 CTMP = CDOTU(NLEN,VECS(1,8),1,VECS(1,6),1) IF (CABS(CTMP).EQ.SZERO) THEN IERR = 8 GO TO 70 END IF C C Compute \alpha and the updated BCG residual. C ALPHA = RHO / CTMP CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,3),-ALPHA,VECS(1,6)) C C Compute the updated BCG residual norm. C STMP = SCNRM2(NLEN,VECS(1,3),1) STMP = STMP * STMP C C Compute the QMR quantities. C BETA = VAR VAR = STMP / TAU COS = SONE / (SONE + VAR ) TAU = STMP * COS QPRD = QPRD * SQRT(SONE - COS) UNRM = QPRD * SQRT(FLOAT(N+1)) C C Compute \hat{p} and the QMR iterate x_n. C CTMP = ALPHA * COS BETA = BETA * COS CALL CAXPBY (NLEN,VECS(1,5),BETA,VECS(1,5),CTMP,VECS(1,4)) CALL CAXPBY (NLEN,VECS(1,1),CONE,VECS(1,1),CONE,VECS(1,5)) C C Compute the updated residual norm. C If the updated residual norm is within one order of magnitude of C the target convergence norm, compute the true residual norm. C IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.STEN).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 CAXPBY (NLEN,VECS(1,6),CONE,VECS(1,2),-CONE,VECS(1,6)) RESN = SCNRM2(NLEN,VECS(1,6),1) / R0 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. the updated residual norm is smaller than the computed C residual norm by a factor of at least 100; C 3. algorithm exceeded the iterations limit. C IF (RESN.LE.TOL) THEN IERR = 0 GO TO 70 ELSE IF (UNRM.LT.RESN/SHUN) THEN IERR = 4 GO TO 70 ELSE IF (N.GE.NLIM) THEN IERR = 4 GO TO 70 END IF C C Check for breakdown due to \rho = 0. C IF (CABS(RHO).EQ.SZERO) THEN IERR = 8 GO TO 70 END IF C C Compute A^T \tilde{q} and \tilde{r}. 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 = 60 RETURN 60 CONTINUE CALL CAXPBY (NLEN,VECS(1,7),CONE,VECS(1,7),-ALPHA,VECS(1,6)) C C Compute \beta and \rho. C CTMP = CDOTU(NLEN,VECS(1,7),1,VECS(1,3),1) BETA = CTMP / RHO RHO = CTMP C C Compute q, \tilde{q}, and A q. C CALL CAXPBY (NLEN,VECS(1,4),BETA,VECS(1,4),CONE,VECS(1,3)) CALL CAXPBY (NLEN,VECS(1,8),BETA,VECS(1,8),CONE,VECS(1,7)) 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**********************************************************************