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 TFQMR algorithm. C C********************************************************************** C SUBROUTINE ZUTFX (NDIM,NLEN,NLIM,VECS,TOL,INFO) C C Purpose: C This subroutine uses the TFQMR 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. 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 91.18, A Transpose-Free Quasi-Minimal Residual Algorithm for C Non-Hermitian Linear Systems, by Roland Freund, September 1991, C which subsequently appeared in SIAM J. Sci. Comput., 14 (1993), C pp. 470--482. C C Parameters: C For a description of the parameters, see the file "zutfx.doc" in C the current directory. C C External routines used: C double precision dlamch(ch) C LAPACK routine, computes machine-related constants. C double precision dznrm2(n,x,incx) C BLAS-1 routine, computes the 2-norm of x. C subroutine zaxpby(n,z,a,x,b,y) C Library routine, computes z = a * x + b * y. C double precision zdotu(n,x,incx,y,incy) C BLAS-1 routine, computes y^H * x. C subroutine zrandn(n,x,seed) C Library routine, fills x with random numbers. C C Noel M. Nachtigal C April 13, 1993 C C********************************************************************** C INTRINSIC CDABS, DBLE, DSQRT, MAX0 EXTERNAL DLAMCH, DZNRM2, ZAXPBY, ZDOTU, ZRANDN DOUBLE COMPLEX ZDOTU DOUBLE PRECISION DLAMCH, DZNRM2 C INTEGER INFO(4), NDIM, NLEN, NLIM DOUBLE COMPLEX VECS(NDIM,9) DOUBLE PRECISION TOL C C Miscellaneous parameters. C DOUBLE COMPLEX ZONE, ZZERO PARAMETER (ZONE = (1.0D0,0.0D0),ZZERO = (0.0D0,0.0D0)) DOUBLE PRECISION DHUN, DONE, DTEN, DZERO PARAMETER (DHUN = 1.0D2,DONE = 1.0D0,DTEN = 1.0D1,DZERO = 0.0D0) C C Local variables, permanent. C INTEGER IERR, N, RETLBL, TF, TRES, VF SAVE IERR, N, RETLBL, TF, TRES, VF DOUBLE COMPLEX ALPHA, BETA, ETA, RHO SAVE ALPHA, BETA, ETA, RHO DOUBLE PRECISION COS, VAR, R0, RESN, TAU, UCHK, UNRM SAVE COS, VAR, R0, RESN, TAU, UCHK, UNRM C C Local variables, transient. C INTEGER INIT, REVCOM DOUBLE COMPLEX ZTMP DOUBLE PRECISION DTMP 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 1 60 returning from AXB, go to label 60 C 1 70 returning from AXB, go to label 70 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 ELSE IF (RETLBL.EQ.60) THEN GO TO 60 ELSE IF (RETLBL.EQ.70) THEN GO TO 70 END IF END IF IERR = 1 GO TO 90 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 90 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 convergence tolerance. C IF (TOL.LE.DZERO) TOL = DSQRT(DLAMCH('E')) C C Start the trace messages and convergence history. C IF (VF.NE.0) WRITE (VF,'(2I8,2E11.4)') 0, 0, DONE, DONE IF (TF.NE.0) WRITE (TF,'(2I8,2E11.4)') 0, 0, DONE, DONE C C Set x_0 = 0 and compute the norm of the initial residual. C CALL ZAXPBY (NLEN,VECS(1,5),ZONE,VECS(1,2),ZZERO,VECS(1,5)) CALL ZAXPBY (NLEN,VECS(1,1),ZZERO,VECS(1,1),ZZERO,VECS(1,1)) R0 = DZNRM2(NLEN,VECS(1,5),1) IF ((TOL.GE.DONE).OR.(R0.EQ.DZERO)) GO TO 90 C C Check whether the auxiliary vector must be supplied. C IF (INIT.EQ.0) CALL ZRANDN (NLEN,VECS(1,3),1) C C Initialize the variables. C N = 1 RESN = DONE RHO = ZONE VAR = DZERO ETA = ZZERO TAU = R0 * R0 IERR = 8 CALL ZAXPBY (NLEN,VECS(1,8),ZZERO,VECS(1,8),ZZERO,VECS(1,8)) CALL ZAXPBY (NLEN,VECS(1,4),ZZERO,VECS(1,4),ZZERO,VECS(1,4)) CALL ZAXPBY (NLEN,VECS(1,6),ZZERO,VECS(1,6),ZZERO,VECS(1,6)) C C This is one step of the TFQMR algorithm. C Compute \beta_{n-1} and \rho_{n-1}. C 20 ZTMP = ZDOTU(NLEN,VECS(1,3),1,VECS(1,5),1) BETA = ZTMP / RHO RHO = ZTMP C C Compute y_{2n-1}, v_{n-1}, and A y_{2n-1}. C CALL ZAXPBY (NLEN,VECS(1,4),BETA,VECS(1,4),ZONE,VECS(1,8)) CALL ZAXPBY (NLEN,VECS(1,6),ZONE,VECS(1,5),BETA,VECS(1,6)) C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,6),VECS(1,9)) C INFO(2) = 1 INFO(3) = 6 INFO(4) = 9 RETLBL = 30 RETURN 30 CALL ZAXPBY (NLEN,VECS(1,4),BETA,VECS(1,4),ZONE,VECS(1,9)) C C Compute \sigma{n-1} and check for breakdowns. C ZTMP = ZDOTU(NLEN,VECS(1,3),1,VECS(1,4),1) IF ((CDABS(ZTMP).EQ.DZERO).OR.(CDABS(RHO).EQ.DZERO)) THEN IERR = 8 GO TO 90 END IF C C Compute \alpha_{n-1}, d_{2n-1} and w_{2n}. C ALPHA = RHO / ZTMP ZTMP = VAR * ETA / ALPHA CALL ZAXPBY (NLEN,VECS(1,7),ZONE,VECS(1,6),ZTMP,VECS(1,7)) CALL ZAXPBY (NLEN,VECS(1,5),ZONE,VECS(1,5),-ALPHA,VECS(1,9)) C C Compute \varepsilon_{2n-1}^2, \eta_{2n-1}^2, c_{2n-1}^2, and C \tau_{2n-1}^2. C DTMP = DZNRM2(NLEN,VECS(1,5),1) DTMP = DTMP * DTMP VAR = DTMP / TAU COS = DONE / ( DONE + VAR ) TAU = DTMP * COS ETA = ALPHA * COS C C Compute x_{2n-1} and the upper bound for its residual norm. C CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),ETA,VECS(1,7)) 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 UNRM = DSQRT(DBLE(2*N) * TAU) / R0 UCHK = UNRM IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN)) GO TO 50 C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,1),VECS(1,9)) C INFO(2) = 1 INFO(3) = 1 INFO(4) = 9 RETLBL = 40 RETURN 40 CALL ZAXPBY (NLEN,VECS(1,9),ZONE,VECS(1,2),-ZONE,VECS(1,9)) RESN = DZNRM2(NLEN,VECS(1,9),1) / R0 UCHK = RESN C C Output the trace messages and convergence history. C 50 IF (VF.NE.0) WRITE (VF,'(2I8,2E11.4)') N, 2*N-1, UNRM, RESN IF (TF.NE.0) WRITE (TF,'(2I8,2E11.4)') N, 2*N-1, UNRM, RESN C C Check for convergence or termination. Stop if: C 1. algorithm converged; C 2. the residual norm upper bound is smaller than the computed C residual norm by a factor of at least 100. C IF (RESN.LE.TOL) THEN IERR = 0 GO TO 90 ELSE IF (UNRM.LT.UCHK/DHUN) THEN IERR = 4 GO TO 90 END IF C C Compute y_{2n}, A y_{2n}, d_{2n}, and w_{2n+1}. C CALL ZAXPBY (NLEN,VECS(1,6),ZONE,VECS(1,6),-ALPHA,VECS(1,4)) ZTMP = VAR * COS CALL ZAXPBY (NLEN,VECS(1,7),ZONE,VECS(1,6),ZTMP,VECS(1,7)) C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,6),VECS(1,8)) C INFO(2) = 1 INFO(3) = 6 INFO(4) = 8 RETLBL = 60 RETURN 60 CALL ZAXPBY (NLEN,VECS(1,5),ZONE,VECS(1,5),-ALPHA,VECS(1,8)) C C Compute \varepsilon_{2n}^2, \eta_{2n}^2, c_{2n}^2, and C \tau_{2n}^2. C DTMP = DZNRM2(NLEN,VECS(1,5),1) DTMP = DTMP * DTMP VAR = DTMP / TAU COS = DONE / ( DONE + VAR ) TAU = DTMP * COS ETA = ALPHA * COS C C Compute x_{2n}. C CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),ETA,VECS(1,7)) 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 UNRM = DSQRT(DBLE(2*N+1) * TAU) / R0 UCHK = UNRM IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 80 C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,1),VECS(1,9)) C INFO(2) = 1 INFO(3) = 1 INFO(4) = 9 RETLBL = 70 RETURN 70 CALL ZAXPBY (NLEN,VECS(1,9),ZONE,VECS(1,2),-ZONE,VECS(1,9)) RESN = DZNRM2(NLEN,VECS(1,9),1) / R0 UCHK = UNRM C C Output the trace messages and convergence history. C 80 IF (VF.NE.0) WRITE (VF,'(2I8,2E11.4)') N, 2*N, UNRM, RESN IF (TF.NE.0) WRITE (TF,'(2I8,2E11.4)') N, 2*N, UNRM, RESN C C Check for convergence or termination. Stop if: C 1. algorithm converged; C 2. the residual norm upper bound 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 90 ELSE IF (UNRM.LT.UCHK/DHUN) THEN IERR = 4 GO TO 90 ELSE IF (N.GE.NLIM) THEN IERR = 4 GO TO 90 END IF C C Update the running counter. C N = N + 1 GO TO 20 C C Done. C 90 NLIM = N RETLBL = 0 INFO(1) = IERR C RETURN END C C**********************************************************************