C********************************************************************** C C Copyright (C) 1991-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 three-term recurrence variant of C the look-ahead Lanczos algorithm. C C********************************************************************** C SUBROUTINE ZUQMR (NDIM,NLEN,NLIM,MAXVW,M,NORM,ZWK,DWK,IDX,IWK, $VECS,TOL,INFO) C C Purpose: C This subroutine uses the QMR algorithm based on the three-term C variant of the look-ahead Lanczos process to solve linear C systems. It runs the QMR algorithm to convergence or until a C user-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 first described in the RIACS Technical Report C 90.46, An Implementation of the Look-Ahead Lanczos Algorithm for C Non-Hermitian Matrices, Part II, by R.W. Freund and N.M. C Nachtigal, November 1990. C C Parameters: C For a description of the parameters, see the file "zuqmr.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,dx,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,dx,incx,dy,incy) C BLAS-1 routine, computes y^H * x. C subroutine zrandn(n,x,seed) C Library routine, fills x with random numbers. C subroutine zrotg(da,db,dcos,dsin) C BLAS-1 routine, computes the Givens rotation which rotates the C vector [a; b] into [ sqrt(a**2 + b**2); 0 ]. C subroutine zuqmr1(ndim,nlen,m,maxvw,numchk,n,l,lstar,nl,nlstar,vf, C ierr,adjust,norm1,norma,tmax,tnrm,zwk,dwk,iwk,idx,vecs) C Low-level routine, builds the vectors v_{n+1} and w_{n+1}. C double precision zuqmro(n) C User-supplied routine, specifies the QMR scaling factors. C C Noel M. Nachtigal C October 24, 1990 C C********************************************************************** C INTRINSIC CDABS, DABS, DBLE, DCMPLX, DCONJG, DMAX1, DSQRT INTRINSIC MAX0, MOD EXTERNAL DLAMCH, DZNRM2, ZAXPBY, ZDOTU, ZUQMRO, ZUQMR1, ZRANDN EXTERNAL ZROTG DOUBLE COMPLEX ZDOTU DOUBLE PRECISION DLAMCH, DZNRM2, ZUQMRO C INTEGER INFO(4), M, MAXVW, NDIM, NLEN, NLIM INTEGER IDX(3,NLIM+2), IWK(M,4) DOUBLE COMPLEX ZWK(M,5*M+7), VECS(NDIM,4*MAXVW+4+4) DOUBLE PRECISION DWK(M,6), NORM, TOL C C Common block variables. C C C Common block ZUQMRX. C DOUBLE PRECISION NORMA COMMON /ZUQMRX/NORMA 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, IEND, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NQMR SAVE IERR, IEND, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NQMR INTEGER NUMCHK, RETLBL, TF, TRES, VF SAVE NUMCHK, RETLBL, TF, TRES, VF DOUBLE PRECISION ADJUST, MAXOMG, NORM1, R0, RESN, TMAX, TMIN SAVE ADJUST, MAXOMG, NORM1, R0, RESN, TMAX, TMIN DOUBLE PRECISION TNRM, UCHK, UNRM SAVE TNRM, UCHK, UNRM C C Local variables, transient. C INTEGER I, IBASE, INIT, ISJ, ISN, J, REVCOM DOUBLE COMPLEX ZTMP, ZTMP1, ZTMP2 DOUBLE PRECISION DTMP, DTMP1, DTMP2 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 40 returning from AXB, go to label 40 C 1 100 returning from AXB, go to label 100 C 2 50 returning from ATXB, go to label 50 C REVCOM = INFO(2) INFO(2) = 0 IF (REVCOM.EQ.0) THEN N = 0 NQMR = 0 MVWBLT = 0 IF (RETLBL.EQ.0) GO TO 10 ELSE IF (REVCOM.EQ.1) THEN IF (RETLBL.EQ.40) THEN GO TO 40 ELSE IF (RETLBL.EQ.100) THEN GO TO 100 END IF ELSE IF (REVCOM.EQ.2) THEN IF (RETLBL.EQ.50) GO TO 50 END IF IERR = 1 GO TO 130 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 (MAXVW.LT.1) IERR = 2 IF (NLEN.GT.NDIM) IERR = 2 IF (M.LT.2*MAXVW+2) IERR = 2 IF (IERR.NE.0) GO TO 130 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 the norms. C NORMA = DZERO NORM1 = DABS(NORM) C C Set the adjustment parameters. C NUMCHK = 25 ADJUST = DTEN C C Extract and check the various tolerances and norm estimates. C TNRM = DLAMCH('E') * DTEN TMIN = DSQRT(DSQRT(DLAMCH('S'))) TMAX = DONE / TMIN IF (TOL.LE.DZERO) TOL = DSQRT(DLAMCH('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 up wrapped indices. The following indices are used: C IDX(1,I) = indices used for the work arrays row-dimensioned M; C IDX(2,I) = indices used for v_i; C IDX(3,I) = indices used for w_i. C DO 20 I = 1, NLIM+2 IDX(1,I) = MOD(I-1,M) + 1 IDX(2,I) = 4 + 1 + 0*(MAXVW+2) + MOD(I-1,MAXVW+2) IDX(3,I) = 4 + 1 + 1*(MAXVW+2) + MOD(I-1,MAXVW+2) 20 CONTINUE C C Set x_0 = 0 and compute the norm of the initial residual. C CALL ZAXPBY (NLEN,VECS(1,IDX(2,1)),ZONE,VECS(1,2),ZZERO,VECS(1,$IDX(2,1))) CALL ZAXPBY (NLEN,VECS(1,1),ZZERO,VECS(1,1),ZZERO,VECS(1,1)) R0 = DZNRM2(NLEN,VECS(1,IDX(2,1)),1) IF ((TOL.GE.DONE).OR.(R0.EQ.DZERO)) GO TO 130 C C Check whether the auxiliary vector must be supplied. C IF (INIT.EQ.0) CALL ZRANDN (NLEN,VECS(1,3),1) CALL ZAXPBY (NLEN,VECS(1,IDX(3,1)),ZONE,VECS(1,3),ZZERO,VECS(1, $IDX(3,1))) C C Scale the first pair of Lanczos vectors and check for invariant C subspaces. C DTMP1 = R0 DTMP2 = DZNRM2(NLEN,VECS(1,IDX(3,1)),1) IF (DTMP1.LT.TNRM) IERR = IERR + 16 IF (DTMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GO TO 130 DWK(IDX(1,1),1) = DONE ZWK(IDX(1,1),IDX(1,1)) = ZDOTU(NLEN,VECS(1,IDX(2,1)),1,VECS(1,$IDX(3,1)),1) / ( DTMP1 * DTMP2 ) IF ((DTMP1.GE.TMAX).OR.(DTMP1.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP1,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(2,1)),ZTMP,VECS(1,IDX(2,1)),ZZERO, $VECS(1,IDX(2,1))) DTMP1 = DONE END IF IF ((DTMP2.GE.TMAX).OR.(DTMP2.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP2,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(3,1)),ZTMP,VECS(1,IDX(3,1)),ZZERO,$VECS(1,IDX(3,1))) DTMP2 = DONE END IF DWK(IDX(1,1),2) = DONE / DTMP1 DWK(IDX(1,1),3) = DONE / DTMP2 C C Initialize the counters. C L = 1 N = 1 NMAX = 0 IWK(IDX(1,0+1),1) = 1 IWK(IDX(1,1+1),1) = 1 MVWBLT = 1 NL = IWK(IDX(1,L+1),1) C C Set up the QMR iteration. C RESN = DONE DWK(IDX(1,1),5) = ZUQMRO(1) ZWK(IDX(1,1),5*M+4) = DWK(IDX(1,1),5) * R0 MAXOMG = DONE / DWK(IDX(1,1),5) C C This is one step of the three-term Lanczos algorithm. C V_n, W_n, D_{n-1}, F_{n-1}, H_{n-1}, and D_{nn} are given. C C Note that the previous step was not necessarily step N-1, as it C could have been a restart. C Except at the first step, the following hold: C NL.LE.N, N-NL.LE.MAXVW-1. C 30 IWK(IDX(1,N),4) = L IWK(IDX(1,N),3) = N C C Build v_{n+1} and w_{n+1}. C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,IDX(2,N)),VECS(1,3)) C INFO(2) = 1 INFO(3) = IDX(2,N) INFO(4) = 3 RETLBL = 40 RETURN C C Have the caller carry out ATXB, then return here. C CALL ATXB (VECS(1,IDX(3,N)),VECS(1,4)) C 40 INFO(2) = 2 INFO(3) = IDX(3,N) INFO(4) = 4 RETLBL = 50 RETURN 50 IF ((VF.NE.0).AND.(N.LE.NMAX)) WRITE (VF,'(A11,I8)') $'Rebuilding:', N+1 CALL ZUQMR1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR,VF,IERR,$ ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,ZWK,DWK,IWK,IDX,VECS $) IF ((IERR.NE.0).AND.(IERR.NE.16).AND.(IERR.NE.32)$ .AND.(IERR.NE.48)) GO TO 130 MVWBLT = MAX0(MVWBLT,NL-IWK(IDX(1,L-1+1),1)) IWK(IDX(1,N),2) = NLSTAR C C Update the counter for steps taken. C NMAX = MAX0(NMAX,N) IF (N.LT.NL-1) GO TO 120 C C The QMR code starts here. C At this point, steps up to NL-1 are guaranteed not to be rebuilt. C No errors are allowed in IERR, other than possibly having found C one or both invariant subspaces, in which case all remaining C iterates are computed. C IEND = NL-1 IF (IERR.NE.0) IEND = N 60 IF (NQMR.GT.IEND-1) GO TO 120 C C Update the QMR iteration counter. C NQMR = NQMR + 1 C C Get the next scaling factor omega(i) and update MAXOMG. C DWK(IDX(1,NQMR+1),5) = ZUQMRO(NQMR+1) MAXOMG = DMAX1(MAXOMG,DONE/DWK(IDX(1,NQMR+1),5)) C C Compute the starting index IBASE for the column of \hat{R}. C IBASE = MAX0(1,IWK(IDX(1,NQMR),2)-1) ZWK(IDX(1,IBASE),5*M+5) = ZZERO C C Multiply the new column by the previous omegas. C DO 70 J = IWK(IDX(1,NQMR),2), NQMR+1 ZWK(IDX(1,J),5*M+5) = DWK(IDX(1,J),5) * ZWK(IDX(1,J),2*M+IDX(1, $NQMR)) 70 CONTINUE C C Apply the previous rotations. C The loop below explicitly implements a call to ZROT: C CALL ZROT (1,ZWK(IDX(1,J-1),5*M+5),1,ZWK(IDX(1,J),5*M+5),1,DWK(IDX(1,J),6),ZWK(IDX(1,J),5*M+6)) C DO 80 J = IBASE+1, NQMR ZTMP1 = ZWK(IDX(1,J),5*M+5) ZTMP2 = ZWK(IDX(1,J-1),5*M+5) ZWK(IDX(1,J-1),5*M+5) = DWK(IDX(1,J),6) * ZTMP2 + ZWK(IDX(1,J),$5*M+6) * ZTMP1 ZWK(IDX(1,J),5*M+5) = DWK(IDX(1,J),6) * ZTMP1 - $DCONJG(ZWK(IDX(1,J),5*M+6)) * ZTMP2 80 CONTINUE C C Compute the rotation for the last element (this also applies it). C CALL ZROTG (ZWK(IDX(1,NQMR),5*M+5),ZWK(IDX(1,NQMR+1),5*M+5),$DWK(IDX(1,NQMR+1),6),ZWK(IDX(1,NQMR+1),5*M+6)) C C Apply the new rotation to the right-hand side vector. C Could be replaced with: C ZWK(IDX(1,NQMR+1),5*M+4) = DZERO C CALL ZROT (1,ZWK(IDX(1,NQMR),5*M+4),1,ZWK(IDX(1,NQMR+1),5*M+4),1,DWK(IDX(1,NQMR+1),6),ZWK(IDX(1,NQMR+1),5*M+6)) C ZWK(IDX(1,NQMR+1),5*M+4) = -DCONJG(ZWK(IDX(1,NQMR+1),5*M+6)) * $ZWK(IDX(1,NQMR),5*M+4) ZWK(IDX(1,NQMR),5*M+4) = DWK(IDX(1,NQMR+1),6) * ZWK(IDX(1,NQMR)$,5*M+4) C C Compute the next search direction s_i. C This is more complicated than it might have to be because storage C for the vectors VECS(1,ISN) is minimized. C ZTMP2 = ZZERO ZTMP = DWK(IDX(1,NQMR),2) ISN = 4 + 1 + 2*(MAXVW+2) + MOD(NQMR-1,2*MAXVW) DO 90 J = IBASE, NQMR-1 ZTMP1 = ZWK(IDX(1,J),5*M+5) * ZWK(IDX(1,J),5*M+7) / ZTMP IF (CDABS(ZTMP1).EQ.DZERO) GO TO 90 ISJ = 4 + 1 + 2*(MAXVW+2) + MOD(J-1,2*MAXVW) CALL ZAXPBY (NLEN,VECS(1,ISN),ZTMP2,VECS(1,ISN),-ZTMP1,VECS(1, $ISJ)) ZTMP2 = ZONE 90 CONTINUE CALL ZAXPBY (NLEN,VECS(1,ISN),ZTMP2,VECS(1,ISN),ZONE,VECS(1,IDX(2,$NQMR))) ZTMP = ZTMP / ZWK(IDX(1,NQMR),5*M+5) C C Compute the new QMR iterate, then scale the search direction. C ZTMP1 = ZTMP * ZWK(IDX(1,NQMR),5*M+4) CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),ZTMP1,VECS(1,ISN)) ZWK(IDX(1,NQMR),5*M+7) = ZTMP DTMP = CDABS(ZTMP) IF ((DTMP.GE.TMAX).OR.(DTMP.LE.TMIN)) THEN ZWK(IDX(1,NQMR),5*M+7) = ZONE CALL ZAXPBY (NLEN,VECS(1,ISN),ZTMP,VECS(1,ISN),ZZERO,VECS(1,ISN $)) 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 UNRM = DSQRT(DBLE(NQMR+1)) * MAXOMG * CDABS(ZWK(IDX(1,NQMR+1),5*M+$4)) / R0 UCHK = UNRM IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 110 C C Have the caller carry out AXB, then return here. C CALL AXB (VECS(1,1),VECS(1,3)) C INFO(2) = 1 INFO(3) = 1 INFO(4) = 3 RETLBL = 100 RETURN 100 CALL ZAXPBY (NLEN,VECS(1,3),ZONE,VECS(1,2),-ZONE,VECS(1,3)) RESN = DZNRM2(NLEN,VECS(1,3),1) / R0 UCHK = RESN C C Output the convergence history. C 110 IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') NQMR, UNRM, RESN IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') NQMR, 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 130 ELSE IF (IERR.NE.0) THEN GO TO 130 ELSE IF (UNRM.LT.UCHK/DHUN) THEN IERR = 4 GO TO 130 ELSE IF (NQMR.GE.NLIM) THEN IERR = 4 GO TO 130 END IF GO TO 60 C C Update the running counter. C 120 IF (IERR.NE.0) GO TO 130 IERR = 4 IF (N.GE.NLIM) GO TO 130 N = N + 1 GO TO 30 C C Done. C 130 RETLBL = 0 NLIM = NQMR INFO(1) = IERR NORM = NORM1 MAXVW = MVWBLT C RETURN END C C********************************************************************** C DOUBLE COMPLEX FUNCTION ZUQMRH(I,J) C C Purpose: C Returns the recurrence coefficients for inner vectors. C C Parameters: C I = row index of the coefficient (input). C J = column index of the coefficient (input). C C Noel M. Nachtigal C July 9, 1993 C C********************************************************************** C C C Common block ZUQMRX. C DOUBLE PRECISION NORMA COMMON /ZUQMRX/NORMA C INTRINSIC DCMPLX C INTEGER I, J C IF ((I.LT.1).OR.(J.LT.1)) THEN ZUQMRH = (0.0D0,0.0D0) ELSE IF (I.EQ.J) THEN ZUQMRH = DCMPLX(NORMA,0.0D0) ELSE ZUQMRH = (0.0D0,0.0D0) END IF C RETURN END C C********************************************************************** C DOUBLE PRECISION FUNCTION ZUQMRO (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 October 7, 1990 C C********************************************************************** C INTEGER I C ZUQMRO = 1.0D0 C RETURN END C C********************************************************************** C SUBROUTINE ZUQMR1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR, $VF,IERR,ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,ZWK,$ DWK,IWK,IDX,VECS) C C Purpose: C This subroutine builds a new pair of vectors v_{n+1} and w_{n+1}. C It is called internally by the Lanczos code, and is not meant to C to be called directly by the user code. C C Parameters: C See descriptions in "zulal.doc" or "zuqmr.doc". 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 zqrdc(x,ldx,n,p,qraux,jpvt,work,job) C LINPACK routine, computes the QR factorization of x. C subroutine zqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) C LINPACK routine, applies the QR factorization of x. C subroutine zsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) C LINPACK routine, computes the SVD of x. C subroutine zuqmr2 (ndim,nlen,m,zwk,dwk,iwk,idx,vecs) C Forces closure of an inner block. C double precision zuqmrh(i,n) C User-supplied routine, computes inner recurrence coefficients. C C Noel M. Nachtigal C July 9, 1993 C c********************************************************************** C INTRINSIC CDABS, DABS, DCMPLX, DMAX1, DREAL, MAX0 EXTERNAL DLAMCH, DZNRM2, ZAXPBY, ZDOTU, ZQRDC, ZQRSL, ZSVDC EXTERNAL ZUQMR2, ZUQMRH DOUBLE COMPLEX ZDOTU, ZUQMRH DOUBLE PRECISION DLAMCH, DZNRM2 C INTEGER IERR, L, LSTAR, M, MAXVW, N, NDIM, NL, NLEN, NLSTAR INTEGER IDX(3,3), IWK(M,4), NUMCHK, VF DOUBLE COMPLEX VECS(NDIM,4*MAXVW+4+4), ZWK(M,5*M+7) DOUBLE PRECISION ADJUST, NORM1, NORMA, TMAX, TMIN, TNRM DOUBLE PRECISION DWK(M,6) C C Miscellaneous parameters. C DOUBLE COMPLEX ZONE, ZZERO PARAMETER (ZONE = (1.0D0,0.0D0),ZZERO = (0.0D0,0.0D0)) DOUBLE PRECISION DONE, DZERO PARAMETER (DONE = 1.0D0,DZERO = 0.0D0) C C Local variables. C INTEGER I, IJ, J, LBLKSZ, NP1 DOUBLE COMPLEX ZTMP, ZTMP1, ZTMP2 DOUBLE PRECISION DTMP, DTMP1, DTMP2, DTMP3, DTMP4 DOUBLE PRECISION IDTMP1, IDTMP2, IDTMP3, IDTMP4 LOGICAL IBUILT, INNER, RERUN C C Initialize local variables. C NP1 = N + 1 DWK(IDX(1,N),4) = DZERO RERUN = .FALSE. IBUILT = .FALSE. IERR = 0 LBLKSZ = N - NL + 1 C C Update the norm estimates. C IF (N.LE.NUMCHK) THEN DTMP1 = DZNRM2(NLEN,VECS(1,3),1) * DWK(IDX(1,N),2) DTMP2 = DZNRM2(NLEN,VECS(1,4),1) * DWK(IDX(1,N),3) NORMA = DMAX1(DTMP1,DTMP2,NORMA) END IF NORM1 = DMAX1(ADJUST*NORMA,NORM1) C C Clear current column of H. C DO 100 I = 1, M ZWK(IDX(1,I),2*M+IDX(1,N)) = ZZERO 100 CONTINUE C C Update D^{(n-1)} to D^{(n)}. C DO 120 I = NL, N-1 ZTMP = ZZERO DO 110 J = NL, N-1 ZTMP = ZTMP + ZWK(IDX(1,I),IDX(1,J)) * ZWK(IDX(1,J),2*M+ $IDX(1,N-1)) 110 CONTINUE ZTMP = ( ZWK(IDX(1,I),M+IDX(1,N-1)) - ZTMP ) / ZWK(IDX(1,N)$,2*M+IDX(1,N-1)) ZWK(IDX(1,I),IDX(1,N)) = ZTMP ZWK(IDX(1,N),IDX(1,I)) = ZTMP * DWK(IDX(1,N),1) / DWK(IDX(1,I), $1) 120 CONTINUE C C Compute l^\star. C LSTAR = MAX0(1,L-1) NLSTAR = IWK(IDX(1,LSTAR+1),1) C C Compute F_{n,1:n-1} and F_{1:n-1,n}. Delay computing F_{nn}. C DO 140 I = MAX0(1,NL-1), N-1 ZTMP = ZZERO DO 130 J = NL, I+1 ZTMP = ZTMP + ZWK(IDX(1,N),IDX(1,J)) * ZWK(IDX(1,J),2*M+$IDX(1,I)) 130 CONTINUE ZWK(IDX(1,N),M+IDX(1,I)) = ZTMP ZWK(IDX(1,I),M+IDX(1,N)) = ZTMP * DWK(IDX(1,I),1) / DWK(IDX(1,N $),1) NORMA = DMAX1(NORMA,CDABS(ZWK(IDX(1,N),M+IDX(1,I))),$CDABS(ZWK(IDX(1,I),M+IDX(1,N)))) NORM1 = DMAX1(ADJUST*NORMA,NORM1) 140 CONTINUE C C Check whether D_l is nonsingular. C IF (LBLKSZ.EQ.1) THEN INNER = CDABS(ZWK(IDX(1,NL),IDX(1,NL))).EQ.DZERO ELSE DO 160 I = NL, N DO 150 J = NL, N ZWK(I-NL+1,4*M+J-NL+1) = ZWK(IDX(1,I),IDX(1,J)) 150 CONTINUE 160 CONTINUE CALL ZSVDC (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+1),ZWK(1,5*M $+3),ZZERO,0,ZZERO,0,ZWK(1,5*M+2),0,I) IF (I.NE.0) THEN IERR = -I GOTO 370 END IF DTMP = DZERO IF (DREAL(ZWK(1,5*M+1)).NE.DZERO) THEN DTMP = DABS( DREAL(ZWK(LBLKSZ,5*M+1)) / DREAL(ZWK(1,5*M+1))) END IF INNER = DTMP.LT.(DBLE(NLEN) * DLAMCH('E')) END IF IF ((VF.NE.0).AND.INNER)$ WRITE (VF,'(A31)') '... moment matrix D is singular' C C Compute H_{n_{l-1}:n_l-1,n}, for n_l > 1. C Also, build the part common to both inner and regular vectors. C Assume that VECS(1,3) = A v_n and VECS(1,4) = A^T w_n. C ZWK(IDX(1,NP1),2*M+IDX(1,N)) = ZZERO IF (NL.EQ.1) GO TO 180 IF (NL.EQ.NLSTAR+1) THEN ZWK(IDX(1,NLSTAR),2*M+IDX(1,N)) = ZWK(IDX(1,NLSTAR),M+IDX(1,N)) $/ ZWK(IDX(1,NLSTAR),IDX(1,NLSTAR)) ZTMP = ZWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),2) /$ DWK(IDX(1,N),2) CALL ZAXPBY (NLEN,VECS(1,3),ZONE,VECS(1,3),-ZTMP,VECS(1,IDX(2, $NLSTAR))) ZTMP = ZWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),3) /$ DWK(IDX(1,N),3) * DWK(IDX(1,N),1) / DWK(IDX(1,NLSTAR),1) CALL ZAXPBY (NLEN,VECS(1,4),ZONE,VECS(1,4),-ZTMP,VECS(1,IDX(3, $NLSTAR))) ELSE IF (NL.GT.NLSTAR+1) THEN C C If a block of size larger than 1 was just closed, combine all its C vectors into one. C ZTMP2 = ZZERO DO 170 I = NL-1, NLSTAR, -1 ZWK(IDX(1,I),2*M+IDX(1,N)) = ZWK(IDX(1,NL-1),M+IDX(1,N)) *$ZWK(IDX(1,I),3*M+IDX(1,L-1)) IF (N.NE.NL) GO TO 170 ZTMP1 = ZWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),2) CALL ZAXPBY (NLEN,VECS(1,IDX(2,NL-1)),ZTMP2,VECS(1,IDX(2,NL- $1)),ZTMP1,VECS(1,IDX(2,I))) ZTMP1 = ZWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),3) /$DWK(IDX(1,I),1) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NL-1)),ZTMP2,VECS(1,IDX(3,NL- $1)),ZTMP1,VECS(1,IDX(3,I))) ZTMP2 = ZONE 170 CONTINUE ZTMP1 = ZWK(IDX(1,NL-1),M+IDX(1,N)) / DWK(IDX(1,N),2) CALL ZAXPBY (NLEN,VECS(1,3),ZONE,VECS(1,3),-ZTMP1,VECS(1,IDX(2,$NL-1))) ZTMP1 = ZWK(IDX(1,NL-1),M+IDX(1,N)) * DWK(IDX(1,N),1) / $DWK(IDX(1,N),3) CALL ZAXPBY (NLEN,VECS(1,4),ZONE,VECS(1,4),-ZTMP1,VECS(1,IDX(3,$NL-1))) END IF C C Compute F_{nn}. C 180 ZWK(IDX(1,N),M+IDX(1,N)) = ZDOTU(NLEN,VECS(1,IDX(3,N)),1,VECS(1,3) $,1) * DWK(IDX(1,N),2) * DWK(IDX(1,N),3) NORMA = DMAX1(NORMA,CDABS(ZWK(IDX(1,N),M+IDX(1,N)))) NORM1 = DMAX1(ADJUST*NORMA,NORM1) IF (INNER) GO TO 320 C C Compute H_{n_l:n,n}. C IWK(IDX(1,L+1+1),1) = NP1 IF (LBLKSZ.EQ.1) THEN ZWK(IDX(1,NL),2*M+IDX(1,N)) = ZWK(IDX(1,NL),M+IDX(1,N)) /$ZWK(IDX(1,NL),IDX(1,NL)) ELSE DO 200 J = NL, N ZWK(J-NL+1,5*M+1) = ZWK(IDX(1,J),M+IDX(1,N)) DO 190 IJ = NL, N ZWK(J-NL+1,4*M+IJ-NL+1) = ZWK(IDX(1,J),IDX(1,IJ)) 190 CONTINUE 200 CONTINUE CALL ZQRDC (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+3),0,ZZERO,0 $) CALL ZQRSL (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+3),ZWK(1,5*M$+1), $ZZERO,ZWK(1,5*M+1),ZWK(1,5*M+2),ZWK(1,5*M+1),ZZERO,$100,J) DO 210 J = NL, N ZWK(IDX(1,J),2*M+IDX(1,N)) = ZWK(J-NL+1,5*M+2) 210 CONTINUE END IF C C Either D_l is nonsingular, or VW is being rerun. For the latter, C just finish building the vectors. For the former, the check for C H_{n_{l-1}:n,n} could be done. However, the look-ahead strategy C requires the smaller of the norm checks for H_{n_{l-1}:n,n} and C H_{n_l:n,n+1}, so both norms are computed first and then checked. C C Build regular vectors. C CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZONE,VECS(1,3),-ZWK(IDX(1,N), $2*M+IDX(1,N)),VECS(1,IDX(2,N))) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,4),-ZWK(IDX(1,N),$2*M+IDX(1,N)),VECS(1,IDX(3,N))) DO 220 I = NL, N-1 ZTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),2) / DWK(IDX(1 $,N),2) CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZONE,VECS(1,IDX(2,NP1)),-$ZTMP,VECS(1,IDX(2,I))) ZTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),3) / DWK(IDX(1 $,N),3) * DWK(IDX(1,N),1) / DWK(IDX(1,I),1) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,IDX(3,NP1)),-$ZTMP,VECS(1,IDX(3,I))) 220 CONTINUE C C Compute scale factors for the new vectors. C 230 DTMP3 = DZNRM2(NLEN,VECS(1,IDX(2,NP1)),1) DTMP4 = DZNRM2(NLEN,VECS(1,IDX(3,NP1)),1) DTMP1 = DWK(IDX(1,N),2) * DTMP3 DTMP2 = DWK(IDX(1,N),3) * DTMP4 ZWK(IDX(1,NP1),2*M+IDX(1,N)) = DCMPLX(DTMP1,DZERO) IF (DTMP1.LT.TNRM) IERR = IERR + 16 IF (DTMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GOTO 360 DWK(IDX(1,NP1),1) = DWK(IDX(1,N),1) * DTMP1 / DTMP2 ZWK(IDX(1,NP1),IDX(1,NP1)) = ZDOTU(NLEN,VECS(1,IDX(3,NP1)),1, $VECS(1,IDX(2,NP1)),1) / ( DTMP3 * DTMP4 ) IF ((DTMP3.GE.TMAX).OR.(DTMP3.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP3,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZTMP,VECS(1,IDX(2,NP1)),$ZZERO,VECS(1,IDX(2,NP1))) DTMP3 = DONE END IF IF ((DTMP4.GE.TMAX).OR.(DTMP4.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP4,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZTMP,VECS(1,IDX(3,NP1)), $ZZERO,VECS(1,IDX(3,NP1))) DTMP4 = DONE END IF DWK(IDX(1,NP1),2) = DONE / DTMP3 DWK(IDX(1,NP1),3) = DONE / DTMP4 C C Compute the last column of D_l^{-1}. C IF (LBLKSZ.EQ.1) THEN ZWK(IDX(1,NL),3*M+IDX(1,L)) = ZONE / ZWK(IDX(1,NL),IDX(1,NL)) ELSE DO 240 I = 1, LBLKSZ-1 ZWK(I,5*M+1) = ZZERO 240 CONTINUE ZWK(LBLKSZ,5*M+1) = ZONE CALL ZQRSL (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+3),ZWK(1,5*M$+1), $ZZERO,ZWK(1,5*M+1),ZWK(1,5*M+2),ZWK(1,5*M+1),ZZERO,$100,J) DO 250 I = NL, N ZWK(IDX(1,I),3*M+IDX(1,L)) = ZWK(I-NL+1,5*M+2) 250 CONTINUE END IF C C If VW is being rerun, then skip to the end. C IF (RERUN) GO TO 360 C C Compute the norm of H_{n_{l-1}:n,n}. C DTMP1 = DZERO DTMP2 = DZERO DO 260 I = IWK(IDX(1,L-1+1),1), N DTMP = CDABS(ZWK(IDX(1,I),2*M+IDX(1,N))) DTMP1 = DTMP1 + DTMP DTMP2 = DTMP2 + DTMP / DWK(IDX(1,I),1) 260 CONTINUE DTMP2 = DTMP2 * DWK(IDX(1,N),1) C C Build the 2nd term for the next step, H_{n_l:n,n+1}, regular. C Compute the norm of H_{n_l:n,n+1}. C DTMP3 = DZERO DTMP4 = DZERO ZTMP = ZWK(IDX(1,NP1),IDX(1,NP1)) * ZWK(IDX(1,NP1),2*M+IDX(1,N)) $* DWK(IDX(1,N),1) / DWK(IDX(1,NP1),1) DO 270 I = NL, N DTMP = CDABS(ZTMP * ZWK(IDX(1,I),3*M+IDX(1,L))) DTMP3 = DTMP3 + DTMP DTMP4 = DTMP4 + DTMP / DWK(IDX(1,I),1) 270 CONTINUE DTMP4 = DTMP4 * DWK(IDX(1,NP1),1) C C Check H_{n_{l-1}:n,n} and H_{n_l:n-1,n+1}. C DTMP = DMAX1(DTMP1,DTMP2,DTMP3,DTMP4) INNER = DTMP.GT.NORM1 IF (.NOT.INNER) GO TO 360 DWK(IDX(1,N),4) = DTMP C C If H_{n_{l-1}:n,n} is bad, build inner vectors. C IF (DMAX1(DTMP1,DTMP2).GT.NORM1) GO TO 320 C C If H_{n_l:n-1,n+1} is bad, check the inner vectors. C This only applies if n_l > 1. C IF (NL.LE.1) GO TO 320 C C Build the inner vectors to compute the 2nd term at the next step. C Get the coefficients H_{n_l:n+1,n} in a temporary location. C Build inner vectors in VECS(1,3) and VECS(1,4). C IBUILT = .TRUE. DO 280 I = NL, N ZWK(IDX(1,I),2*M+IDX(1,NP1)) = ZUQMRH(I,N) ZTMP = ZWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),2) /$DWK(IDX(1,N),2) CALL ZAXPBY (NLEN,VECS(1,3),ZONE,VECS(1,3),-ZTMP,VECS(1,IDX(2,I $))) ZTMP = ZWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),3) /$DWK(IDX(1,N),3) * DWK(IDX(1,N),1) / DWK(IDX(1,I),1) CALL ZAXPBY (NLEN,VECS(1,4),ZONE,VECS(1,4),-ZTMP,VECS(1,IDX(3,I $))) 280 CONTINUE IDTMP3 = DZNRM2(NLEN,VECS(1,3),1) IDTMP4 = DZNRM2(NLEN,VECS(1,4),1) IDTMP1 = DWK(IDX(1,N),2) * IDTMP3 IDTMP2 = DWK(IDX(1,N),3) * IDTMP4 ZWK(IDX(1,NP1),2*M+IDX(1,NP1)) = DCMPLX(IDTMP1,DZERO) IF (IDTMP1.LT.TNRM) IERR = IERR + 16 IF (IDTMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GOTO 330 C C Build the 2nd term for the next step, H_{n_{l-1}:n_l-1,n+1}, C inner. Compute the norm of H_{n_{l-1}:n_l-1,n+1}. C ZTMP = ZZERO DO 290 J = NLSTAR, NL-1 ZTMP = ZTMP + ZWK(IDX(1,NL),IDX(1,J)) * ZWK(IDX(1,I),2*M+IDX(1,$N)) 290 CONTINUE DO 300 J = NL, N ZTMP = ZTMP + ZWK(IDX(1,NL),IDX(1,J)) * ZWK(IDX(1,I),2*M+IDX(1, $NP1)) 300 CONTINUE DTMP1 = DZERO DTMP2 = DZERO DTMP3 = DMAX1(DTMP3,DTMP4) ZTMP = ( ZWK(IDX(1,NL),M+IDX(1,N)) - ZTMP ) / ZWK(IDX(1,NP1),2*M+$IDX(1,NP1)) ZTMP = ZTMP * ZWK(IDX(1,NL),2*M+IDX(1,NL-1)) * DWK(IDX(1,NL-1),1) $/ DWK(IDX(1,NL),1) DO 310 I = IWK(IDX(1,L-1+1),1), NL-1 DTMP = CDABS(ZWK(IDX(1,I),3*M+IDX(1,L-1)) * ZTMP) DTMP1 = DTMP1 + DTMP DTMP2 = DTMP2 + DTMP / DWK(IDX(1,I),1) 310 CONTINUE DTMP2 = DTMP2 * DWK(IDX(1,N),1) C C Compare the inner and regular versions of the 2nd term at the next C step. Build the vector corresponding to the smaller term. C INNER = DTMP3.GT.DMAX1(DTMP1,DTMP2) IF (.NOT.INNER) GO TO 360 C C Build inner vectors. C Check whether the block has to be forced to close. C 320 IF (VF.NE.0) WRITE (VF,'(A7,I5,A9)') 'Vector ',NP1,' is inner' IF (NP1-NL.EQ.MAXVW) THEN CALL ZUQMR2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR,ADJUST,$ NORM1,ZWK,DWK,IWK,IDX,VECS) IF (IERR.NE.0) GO TO 370 LBLKSZ = N - NL + 1 INNER = .FALSE. RERUN = .TRUE. NP1 = N + 1 GO TO 230 END IF C C The temporary vectors contain either just partial inner vectors, C or the completed ones, depending on whether IBUILT is TRUE or C FALSE. In either case, replace the regular vectors. C 330 CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZONE,VECS(1,3),ZZERO,VECS(1, $IDX(2,NP1))) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,4),ZZERO,VECS(1,$IDX(3,NP1))) C C Get the coefficients H_{n_l:n+1,n} and build inner vectors. C IF (IBUILT) THEN DO 340 I = NL, NP1 ZWK(IDX(1,I),2*M+IDX(1,N)) = ZWK(IDX(1,I),2*M+IDX(1,NP1)) 340 CONTINUE DTMP1 = IDTMP1 DTMP2 = IDTMP2 DTMP3 = IDTMP3 DTMP4 = IDTMP4 ELSE DO 350 I = NL, N ZWK(IDX(1,I),2*M+IDX(1,N)) = ZUQMRH(I,N) ZTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),2) / $DWK(IDX(1,N),2) CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZONE,VECS(1,IDX(2,NP1))$,-ZTMP,VECS(1,IDX(2,I))) ZTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),3) / $DWK(IDX(1,N),3) * DWK(IDX(1,N),1) / DWK(IDX(1,I),1) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,IDX(3,NP1))$,-ZTMP,VECS(1,IDX(3,I))) 350 CONTINUE DTMP3 = DZNRM2(NLEN,VECS(1,IDX(2,NP1)),1) DTMP4 = DZNRM2(NLEN,VECS(1,IDX(3,NP1)),1) DTMP1 = DWK(IDX(1,N),2) * DTMP3 DTMP2 = DWK(IDX(1,N),3) * DTMP4 END IF ZWK(IDX(1,NP1),2*M+IDX(1,N)) = DCMPLX(DTMP1,DZERO) IF (DTMP1.LT.TNRM) IERR = IERR + 16 IF (DTMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GOTO 360 DWK(IDX(1,NP1),1) = DWK(IDX(1,N),1) * DTMP1 / DTMP2 ZWK(IDX(1,NP1),IDX(1,NP1)) = ZDOTU(NLEN,VECS(1,IDX(3,NP1)),1, $VECS(1,IDX(2,NP1)),1) / ( DTMP3 * DTMP4 ) IF ((DTMP3.GE.TMAX).OR.(DTMP3.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP3,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZTMP,VECS(1,IDX(2,NP1)),$ZZERO,VECS(1,IDX(2,NP1))) DTMP3 = DONE END IF IF ((DTMP4.GE.TMAX).OR.(DTMP4.LE.TMIN)) THEN ZTMP = DCMPLX(DONE / DTMP4,DZERO) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZTMP,VECS(1,IDX(3,NP1)), $ZZERO,VECS(1,IDX(3,NP1))) DTMP4 = DONE END IF DWK(IDX(1,NP1),2) = DONE / DTMP3 DWK(IDX(1,NP1),3) = DONE / DTMP4 C C If regular vectors were built, update the counters. C 360 IF (.NOT.INNER) THEN L = L + 1 NL = IWK(IDX(1,L+1),1) END IF C 370 RETURN END C C********************************************************************** C SUBROUTINE ZUQMR2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR,$ ADJUST,NORM1,ZWK,DWK,IWK,IDX,VECS) C C Purpose: C This subroutine rebuilds the data for vectors v_{n+1} and w_{n+1} C at the point where a block was forced to close. The routine is C called internally by the look-ahead Lanczos code and is not meant C to be called directly by the user code. C C Parameters: C See descriptions in "zulal.doc" or "zuqmr.doc". C C External routines used: C subroutine zaxpby(n,z,a,x,b,y) C Library routine, computes z = a * x + b * y. C subroutine zqrdc(x,ldx,n,p,qraux,jpvt,work,job) C LINPACK routine, computes the QR factorization of x. C subroutine zqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) C LINPACK routine, applies the QR factorization of x. C C Noel M. Nachtigal C May 31, 1993 C C********************************************************************** C EXTERNAL ZAXPBY, ZQRDC, ZQRSL C INTEGER IERR, L, LSTAR, M, N, NDIM, NL, NLEN, NLSTAR INTEGER IDX(3,3), IWK(M,4), VF DOUBLE COMPLEX VECS(NDIM,*), ZWK(M,5*M+7) DOUBLE PRECISION ADJUST, DWK(M,6), NORM1 C C Miscellaneous parameters. C DOUBLE COMPLEX ZONE, ZZERO PARAMETER (ZONE = (1.0D0,0.0D0),ZZERO = (0.0D0,0.0D0)) DOUBLE PRECISION DZERO PARAMETER (DZERO = 0.0D0) C C Local variables. C INTEGER I, IJ, J, LBLKSZ, NP1 DOUBLE COMPLEX SCALV, SCALW, ZTMP1, ZTMP2 DOUBLE PRECISION DTMP1, DTMP2 C C Find the index of the vector pair with the smallest pass value. C IERR = 0 IF (VF.NE.0) WRITE (VF,'(A20)') 'block did not close:' J = NL DTMP1 = DWK(IDX(1,J),4) DO 100 I = NL+1, N DTMP2 = DWK(IDX(1,I),4) IF (DTMP2.GT.DZERO) THEN IF ((DTMP1.EQ.DZERO).OR.(DTMP2.LT.DTMP1)) THEN J = I DTMP1 = DTMP2 END IF END IF 100 CONTINUE IF (DTMP1.EQ.DZERO) THEN IF (VF.NE.0) WRITE (VF,'(A47)') $'... no new norm estimates available (aborting).' IERR = 8 RETURN END IF NORM1 = ADJUST * DTMP1 IF (VF.NE.0) WRITE (VF,'(A40,I5,E11.4)')$ '... updated norms, restarting from step:', IWK(IDX(1,J),3), $NORM1 IF (IWK(IDX(1,J),3).EQ.N) RETURN N = IWK(IDX(1,J),3) L = IWK(IDX(1,N),4) NL = IWK(IDX(1,L+1),1) C C Initialize local variables. C NP1 = N + 1 DWK(IDX(1,N),4) = DZERO LBLKSZ = N - NL + 1 C C Compute l^\star. C LSTAR = L - 1 NLSTAR = IWK(IDX(1,LSTAR+1),1) C C Compute H_{n_l:n,n}. Save the old coefficients. C IWK(IDX(1,L+1+1),1) = NP1 IF (LBLKSZ.EQ.1) THEN ZWK(IDX(1,NL),2*M+IDX(1,NP1)) = ZWK(IDX(1,NL),2*M+IDX(1,N)) ZWK(IDX(1,NL),2*M+IDX(1,N)) = ZWK(IDX(1,NL),M+IDX(1,N)) /$ZWK(IDX(1,NL),IDX(1,NL)) ELSE DO 120 J = NL, N ZWK(J-NL+1,5*M+1) = ZWK(IDX(1,J),M+IDX(1,N)) DO 110 IJ = NL, N ZWK(J-NL+1,4*M+IJ-NL+1) = ZWK(IDX(1,J),IDX(1,IJ)) 110 CONTINUE 120 CONTINUE CALL ZQRDC (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+3),0,ZZERO,0 $) CALL ZQRSL (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,5*M+3),ZWK(1,5*M$+1), $ZZERO,ZWK(1,5*M+1),ZWK(1,5*M+2),ZWK(1,5*M+1),ZZERO,$100,J) DO 130 J = NL, N ZWK(IDX(1,J),2*M+IDX(1,NP1)) = ZWK(IDX(1,J),2*M+IDX(1,N)) ZWK(IDX(1,J),2*M+IDX(1,N)) = ZWK(J-NL+1,5*M+2) 130 CONTINUE END IF C C Convert inner vectors to regular vectors. C SCALV = ZWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),2) SCALW = ZWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),3) * $DWK(IDX(1,N),1) / DWK(IDX(1,NP1),1) ZWK(IDX(1,NP1),2*M+IDX(1,N)) = ZZERO DO 140 I = NL, N ZTMP1 = ZWK(IDX(1,I),2*M+IDX(1,N)) - ZWK(IDX(1,I),2*M+IDX(1,NP1$)) ZTMP2 = ZTMP1 * DWK(IDX(1,I),2) / SCALV CALL ZAXPBY (NLEN,VECS(1,IDX(2,NP1)),ZONE,VECS(1,IDX(2,NP1)),- $ZTMP2,VECS(1,IDX(2,I))) ZTMP2 = ZTMP1 * DWK(IDX(1,I),3) / SCALW * DWK(IDX(1,N),1) /$DWK(IDX(1,I),1) CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,IDX(3,NP1)),- \$ZTMP2,VECS(1,IDX(3,I))) 140 CONTINUE C RETURN END C C**********************************************************************