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 eigenvalue solver, using C the three-term recurrence variant of the look-ahead Lanczos C algorithm. C C********************************************************************** C SUBROUTINE SULAL (NDIM,NLEN,NLIM,MAXN,MAXVW,M,NORM,DWK,IDX,IWK, $ HWK,VECS,INFO) C C Purpose: C This subroutine uses the Lanczos algorithm with look-ahead to C compute eigenvalue estimates. The algorithm was first described C in the RIACS Technical Report 90.45, `An Implementation of the C Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices, Part I`, C by R.W. Freund, M.H. Gutknecht and N.M. Nachtigal, November 1990. C The routine will first run the look-ahead Lanczos algorithm for a C number of steps. It will then prompt the user for the number of C eigenvalues to compute. In addition, the routine will also find C the common eigenvalues in two successive runs. C C Parameters: C For a description of the parameters, see the file "sulal.doc" in C the current directory. C C External routines used: C subroutine saxpby(n,z,a,x,b,y) C Library routine, computes z = a * x + b * y. C single precision sdot(n,dx,incx,dy,incy) C BLAS-1 routine, computes y^H * x. C single precision slamch(ch) C LAPACK routine, computes machine-related constants. C single precision snrm2(n,dx,incx) C BLAS-1 routine, computes the 2-norm of x. C subroutine srandn(n,x,seed) C Library routine, fills x with random numbers. C subroutine sulal1(ndim,nlen,m,maxvw,numchk,n,l,lstar,nl,nlstar,vf, C ierr,adjust,norm1,norma,tmax,tnrm,dwk,iwk,idx,vecs) C Low-level routine, builds the vectors v_{n+1} and w_{n+1}. C subroutine sulalc(maxn,neig,nold,hwk,tf,vf) C Library routine, computes common eigenvalues. C subroutine sulale(maxn,neig,hwk,tf,vf) C Library routine, computes Lanczos eigenvalue approximations. C C Noel M. Nachtigal C October 24, 1990 C C********************************************************************** C INTRINSIC ABS, SQRT, MAX0, MOD EXTERNAL SAXPBY, SDOT, SLAMCH, SNRM2, SULAL1, SULALC, SULALE REAL SDOT, SLAMCH, SNRM2 C INTEGER INFO(4), M, MAXN, MAXVW, NDIM, NLEN, NLIM INTEGER IDX(3,NLIM+2), IWK(M,4) REAL DWK(M,5*M+7), HWK(MAXN,2*MAXN+8), NORM REAL VECS(NDIM,2*MAXVW+4+2) C C Common block variables. C C C Common block SULALX. C REAL NORMA COMMON /SULALX/NORMA C C Miscellaneous parameters. C REAL DONE, DTEN, SZERO PARAMETER (DONE = 1.0E0,DTEN = 1.0E1,SZERO = 0.0E0) C C Local variables, permanent. C INTEGER IERR, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NUMCHK SAVE IERR, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NUMCHK INTEGER RETLBL, TF, VF SAVE RETLBL, TF, VF REAL ADJUST, NORM1, TMAX, TMIN, TNRM SAVE ADJUST, NORM1, TMAX, TMIN, TNRM C C Local variables, transient. C INTEGER I, J, K, NEIG, REVCOM REAL STMP, STMP1, STMP2 C C Local variables. C INTEGER NLAL, NOLD 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 60 returning from AXB, go to label 60 C 2 70 returning from ATXB, go to label 70 C REVCOM = INFO(2) INFO(2) = 0 IF (REVCOM.EQ.0) THEN N = 0 NLAL = 0 MVWBLT = 0 IF (RETLBL.EQ.0) GO TO 10 ELSE IF (REVCOM.EQ.1) THEN IF (RETLBL.EQ.60) THEN GO TO 60 END IF ELSE IF (REVCOM.EQ.2) THEN IF (RETLBL.EQ.70) GO TO 70 END IF IERR = 1 GO TO 110 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 (NLIM.GT.MAXN-2) IERR = 2 IF (M.LT.2*MAXVW+2) IERR = 2 IF (IERR.NE.0) GO TO 110 C C Extract from INFO the output units TF and VF, and the starting C vector flags J and K. C VF = MAX0(INFO(1),0) J = VF / 100000 VF = VF - J * 100000 K = VF / 10000 VF = VF - K * 10000 TF = VF / 100 VF = VF - TF * 100 C C Extract the norms. C NORMA = SZERO NORM1 = ABS(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 = SLAMCH('E') * DTEN TMIN = SQRT(SQRT(SLAMCH('S'))) TMAX = DONE / TMIN C C Set up wrapped indices. The following indices are used: C IDX(1,I) = indices used for IWK and DWK, 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) = 2 + 1 + 0*(MAXVW+2) + MOD(I-1,MAXVW+2) IDX(3,I) = 2 + 1 + 1*(MAXVW+2) + MOD(I-1,MAXVW+2) 20 CONTINUE C C Check whether any starting vectors must be supplied. C IF (K.EQ.0) CALL SRANDN (NLEN,VECS(1,1),1) CALL SAXPBY (NLEN,VECS(1,IDX(2,1)),DONE,VECS(1,1),SZERO,VECS(1, $IDX(2,1))) IF (J.EQ.0) CALL SRANDN (NLEN,VECS(1,2),1) CALL SAXPBY (NLEN,VECS(1,IDX(3,1)),DONE,VECS(1,2),SZERO,VECS(1, $IDX(3,1))) C C Zero out the Lanczos matrix. C DO 40 I = 1, MAXN DO 30 J = 1, MAXN HWK(I,J) = SZERO 30 CONTINUE 40 CONTINUE C C Scale the first pair of Lanczos vectors and check for invariant C subspaces. C STMP1 = SNRM2(NLEN,VECS(1,IDX(2,1)),1) STMP2 = SNRM2(NLEN,VECS(1,IDX(3,1)),1) IF (STMP1.LT.TNRM) IERR = IERR + 16 IF (STMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GO TO 110 DWK(IDX(1,1),5*M+4) = DONE DWK(IDX(1,1),IDX(1,1)) = SDOT(NLEN,VECS(1,IDX(2,1)),1,VECS(1,IDX(3 $,1)),1) / ( STMP1 * STMP2 ) IF ((STMP1.GE.TMAX).OR.(STMP1.LE.TMIN)) THEN STMP = DONE / STMP1 CALL SAXPBY (NLEN,VECS(1,IDX(2,1)),STMP,VECS(1,IDX(2,1)),SZERO, $VECS(1,IDX(2,1))) STMP1 = DONE END IF IF ((STMP2.GE.TMAX).OR.(STMP2.LE.TMIN)) THEN STMP = DONE / STMP2 CALL SAXPBY (NLEN,VECS(1,IDX(3,1)),STMP,VECS(1,IDX(3,1)),SZERO, $VECS(1,IDX(3,1))) STMP2 = DONE END IF DWK(IDX(1,1),5*M+5) = DONE / STMP1 DWK(IDX(1,1),5*M+6) = DONE / STMP2 C C Initialize the counters. C L = 1 N = 1 NLAL = 0 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 This is one step of the coupled two-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 IF (VF.NE.0) WRITE (VF,'(A29)') 'Running look-ahead Lanczos...' 50 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,1)) C INFO(2) = 1 INFO(3) = IDX(2,N) INFO(4) = 1 RETLBL = 60 RETURN C C Have the caller carry out ATXB, then return here. C CALL ATXB (VECS(1,IDX(3,N)),VECS(1,2)) C 60 INFO(2) = 2 INFO(3) = IDX(3,N) INFO(4) = 2 RETLBL = 70 RETURN 70 IF ((VF.NE.0).AND.(N.LE.NMAX)) WRITE (VF,'(A11,I8)') 'Rebuilding:' $, N+1 CALL SULAL1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR,VF,IERR, $ ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,DWK,IWK,IDX,VECS) IF ((IERR.NE.0).AND.(IERR.NE.16).AND.(IERR.NE.32) $ .AND.(IERR.NE.48)) GO TO 110 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 NLAL = N NMAX = MAX0(NMAX,NLAL) C C Extract the column of H. C CALL SAXPBY (N+1,HWK(1,N),SZERO,HWK(1,N),SZERO,HWK(1,N)) DO 80 I = NLSTAR, N+1 HWK(I,N) = DWK(IDX(1,I),2*M+IDX(1,N)) 80 CONTINUE C C Do the next step. C IF ((IERR.NE.0).OR.(N.GE.NLIM)) GO TO 90 N = N + 1 GO TO 50 C C The eigenvalue code starts here. C 90 NOLD = 0 100 NEIG = 0 WRITE (6,'()') WRITE (6,'(A43,I4)') 'Number of Lanczos steps completed : ' $, NLAL WRITE (6,'(A43,$)') 'Number of eigenvalues to compute (0=done): ' READ (5,*) NEIG IF ((NEIG.GT.0).AND.(NEIG.LE.NLAL)) THEN CALL SULALE (MAXN,NEIG,HWK,TF,VF) IF (NEIG.EQ.0) GO TO 100 IF (NOLD.NE.0) CALL SULALC (MAXN,NEIG,NOLD,HWK,TF,VF) NOLD = NEIG CALL SAXPBY (NOLD,HWK(1,2*MAXN+7),DONE,HWK(1,2*MAXN+4),SZERO, $HWK(1,2*MAXN+7)) CALL SAXPBY (NOLD,HWK(1,2*MAXN+8),DONE,HWK(1,2*MAXN+5),SZERO, $HWK(1,2*MAXN+8)) ELSE IF (NEIG.EQ.0) THEN GO TO 110 END IF GO TO 100 C C Done. C 110 RETLBL = 0 NLIM = NLAL INFO(1) = IERR NORM = NORM1 MAXVW = MVWBLT C RETURN END C C********************************************************************** C SUBROUTINE SULALC (MAXN,NEIG,NOLD,HWK,TF,VF) C C Purpose: C This subroutine finds eigenvalues common to the last two runs of C the eigenvalue routine. C C Parameters: C MAXN = the dimensioned size of the array HWK (input). C NEIG = the number of eigs in the current run (input). C NOLD = on entry, the number of eigs in the previous run; on exit C the number of common eigs (input/output). C HWK = single precision work array (input/output). C TF = if not 0, then the routine will output the common eigs to C this unit (input). C VF = if not 0, then the routine will produce verbose messages C on this unit (input). C C External routines used: C subroutine shpsrt(n,wr,wi,wm) C Library routine, sorts wr and wi and wm. C C Noel M. Nachtigal C August 6, 1993 C C********************************************************************** C INTRINSIC ABS EXTERNAL SHPSRT C INTEGER MAXN, NEIG, NOLD, TF, VF REAL HWK(MAXN,2*MAXN+8) C C Local variables. C INTEGER I, J, K REAL DISTM, EIMAG, EMAGN, EREAL, OIMAG, OREAL, OMAGN REAL TOL, TOLM C C Get the number of eigenvalues to find. C 10 WRITE (6,'()') WRITE (6,'(A43,$)') 'Find the common eigenvalues (1=Yes/0=No) ? ' READ (5,*) I IF (I.NE.1) RETURN WRITE (6,'(A43,$)') 'Enter separation tolerance : ' READ (5,*) TOL C C Check for common eigenvalues. C K = 0 DO 30 I = 1, NEIG EREAL = HWK(I,2*MAXN+4) EIMAG = HWK(I,2*MAXN+5) EMAGN = HWK(I,2*MAXN+6) DO 20 J = 1, NOLD OREAL = HWK(J,2*MAXN+7) OIMAG = HWK(J,2*MAXN+8) OMAGN = SQRT(OREAL**2 + OIMAG**2) TOLM = 0.5E0 * TOL * (EMAGN + OMAGN) IF (ABS(EMAGN-OMAGN).GE.TOLM) GO TO 20 DISTM = SQRT((EREAL-OREAL)**2 + (EIMAG-OIMAG)**2) IF (DISTM.GE.TOLM) GO TO 20 K = K + 1 HWK(K,2*MAXN+1) = EREAL HWK(K,2*MAXN+2) = EIMAG HWK(K,2*MAXN+3) = EMAGN GO TO 30 20 CONTINUE 30 CONTINUE C C Sort the final eigenvalues by decreasing real part. C CALL SHPSRT (K,HWK(1,2*MAXN+1),HWK(1,2*MAXN+2),HWK(1,2*MAXN+3)) C C Output the results. C IF (VF.NE.0) WRITE (VF,'(A43,I4)') $ 'Number of common eigenvalues found : ', K IF (TF.NE.0) THEN WRITE (TF,'(A30,E12.7)') 'Common eigenvalues, tolerance: ', TOL WRITE (TF,'(2E25.17)') (HWK(I,2*MAXN+1),HWK(I,2*MAXN+2),I=1,K) WRITE (TF,'()') END IF GO TO 10 C END C C********************************************************************** C SUBROUTINE SULALE (MAXN,NEIG,HWK,TF,VF) C C Purpose: C This subroutine computes the eigenvalues estimates obtained from C the NEIG x NEIG upper left corner of the Lanczos matrix. It uses C the heuristic proposed by Cullum and Willoughby to separate the C genuine eigenvalues from the spurious ones. C C Parameters: C MAXN = the dimensioned size of the array HWK (input). C NEIG = the size of the submatrix whose eigs are computed. On C exit, it is the number of eigs found (input/output). C HWK = single precision work array (input/output). C TF = if not 0, then the routine will output the Lanczos eigs C to this unit (input). C VF = if not 0, then the routine will produce verbose messages C on this unit (input). C C External routines used: C subroutine balanc(nm,n,a,low,igh,scale) C EISPACK routine, balances a matrix. C subroutine hqr(nm,n,low,igh,h,wr,wi,ierr) C EISPACK routine, computes eigenvalues of an upper Hessenberg C matrix. C subroutine saxpby(n,z,a,x,b,y) C Library routine, computes z = a * x + b * y. C subroutine shpsrt(n,wr,wi,wm) C Library routine, sorts wr and wi and wm. C single precision slamch(ch) C LAPACK routine, computes machine-related constants. C C Noel M. Nachtigal C August 6, 1993 C C********************************************************************** C INTRINSIC ABS, SQRT EXTERNAL BALANC, SAXPBY, SHPSRT, SLAMCH, HQR REAL SLAMCH C INTEGER MAXN, NEIG, TF, VF REAL HWK(MAXN,2*MAXN+8) C C Miscellaneous parameters. C REAL DONE, SZERO PARAMETER (DONE = 1.0E0,SZERO = 0.0E0) C C Local variables. C INTEGER I, IERR, J, K, NBEG, NCHK, NEND, NLAL REAL DISTM, EIMAG, EMAGN, EREAL, TOL, TOLM LOGICAL MULTIP C C Compute the separation tolerance. C TOL = SQRT(SLAMCH('E')) C C Extract the NEIG x NEIG matrix. C DO 10 I = 1, NEIG CALL SAXPBY (I+1,HWK(1,MAXN+I),DONE,HWK(1,I),SZERO,HWK(1,MAXN+I $)) 10 CONTINUE C C Balance it and compute its eigenvalues. C IF (VF.NE.0) WRITE (VF,'(A24)') 'Computing eigenvalues...' CALL BALANC (MAXN,NEIG,HWK(1,MAXN+1),I,J,HWK(1,2*MAXN+1)) CALL HQR (MAXN,NEIG,I,J,HWK(1,MAXN+1),HWK(1,2*MAXN+1),HWK(1,2*MAXN $+2),IERR) IF (IERR.NE.0) THEN WRITE (6,'(A23,I2)') 'Error return from HQR: ', IERR RETURN END IF C C Compute the magnitude of the eigenvalues. C DO 20 I = 1, NEIG HWK(I,2*MAXN+3) = SQRT(HWK(I,2*MAXN+1)**2 + HWK(I,2*MAXN+2)**2) 20 CONTINUE C C Sort the eigenvalues in decreasing order of magnitude. C CALL SHPSRT (NEIG,HWK(1,2*MAXN+3),HWK(1,2*MAXN+1),HWK(1,2*MAXN+2)) IF (TF.NE.0) THEN WRITE (TF,'(A33,I4)') 'All Lanczos eigenvalues at step: ', NEIG $ WRITE (TF,'(2E25.17)') (HWK(I,2*MAXN+1),HWK(I,2*MAXN+2),I=1, $NEIG) WRITE (TF,'()') END IF C C Separate the repeated eigenvalues from the simple ones. C I = 1 NEND = 0 NLAL = NEIG NBEG = NEIG + 1 30 IF (I.GT.NLAL) GO TO 70 C C Outer loop, over all eigenvalues. C J = I + 1 EREAL = HWK(I,2*MAXN+1) EIMAG = HWK(I,2*MAXN+2) EMAGN = HWK(I,2*MAXN+3) MULTIP = .FALSE. 40 IF (J.GT.NLAL) GO TO 60 C C Inner loop, find eigenvalues J close to eigenvalue I. "Close" is C defined relative to the average of the two points, so as to be a C symmetric relationship. C TOLM = 0.5E0 * TOL * (EMAGN + HWK(J,2*MAXN+3)) C C If the difference in magnitude between I and J is greater than C the tolerance, skip the remaining eigenvalues, they will be even C further away (the eigenvalues are assumed sorted in some order of C magnitude). C IF (EMAGN-HWK(J,2*MAXN+3).GE.TOLM) GO TO 60 C C Compute the distance between the two eigenvalues. C DISTM = SQRT((EREAL-HWK(J,2*MAXN+1))**2 + (EIMAG-HWK(J,2*MAXN+2)) $**2) C C Check the eigenvalue J. C IF (DISTM.GE.TOLM) THEN J = J + 1 ELSE MULTIP = .TRUE. DO 50 K = J, NLAL-1 HWK(K,2*MAXN+1) = HWK(K+1,2*MAXN+1) HWK(K,2*MAXN+2) = HWK(K+1,2*MAXN+2) HWK(K,2*MAXN+3) = HWK(K+1,2*MAXN+3) 50 CONTINUE NLAL = NLAL - 1 END IF GO TO 40 C C Collect the multiple eigenvalues in ELAL(1:NEND), and the simple C ones in ELAL(NBEG:NEIG). C 60 IF (MULTIP) THEN NEND = NEND + 1 HWK(NEND,2*MAXN+4) = EREAL HWK(NEND,2*MAXN+5) = EIMAG HWK(NEND,2*MAXN+6) = EMAGN ELSE NBEG = NBEG - 1 HWK(NBEG,2*MAXN+4) = EREAL HWK(NBEG,2*MAXN+5) = EIMAG HWK(NBEG,2*MAXN+6) = EMAGN END IF I = I + 1 GO TO 30 70 CONTINUE C C Check whether they are all considered genuine. C IF (NEND.EQ.NLAL) GO TO 120 C C Extract the (NEIG-1) x (NEIG-1) matrix. C NCHK = NEIG - 1 DO 80 I = 1, NCHK CALL SAXPBY (I+1,HWK(1,MAXN+I),DONE,HWK(2,I+1),SZERO,HWK(1,MAXN $+I)) 80 CONTINUE C C Balance it and compute its eigenvalues. C IF (VF.NE.0) WRITE (VF,'(A30)') 'Computing check eigenvalues...' CALL BALANC (MAXN,NCHK,HWK(1,MAXN+1),I,J,HWK(1,2*MAXN+1)) CALL HQR (MAXN,NCHK,I,J,HWK(1,MAXN+1),HWK(1,2*MAXN+1),HWK(1,2*MAXN $+2),IERR) IF (IERR.NE.0) THEN WRITE (6,'(A23,I2)') 'Error return from HQR: ', IERR RETURN END IF C C Compute the magnitude of the check eigenvalues. C DO 90 I = 1, NCHK HWK(I,2*MAXN+3) = SQRT(HWK(I,2*MAXN+1)**2 + HWK(I,2*MAXN+2)**2) 90 CONTINUE C C Sort the eigenvalues in decreasing order of magnitude. C CALL SHPSRT (NCHK,HWK(1,2*MAXN+3),HWK(1,2*MAXN+1),HWK(1,2*MAXN+2)) C C Check the Lanczos eigenvalues that were simple; keep those that C appear only in the Lanczos matrix. C The approach used is brute force; it could be more sophisticated. C DO 110 I = NBEG, NEIG EREAL = HWK(I,2*MAXN+4) EIMAG = HWK(I,2*MAXN+5) EMAGN = HWK(I,2*MAXN+6) DO 100 J = 1, NCHK TOLM = 0.5E0 * TOL * (EMAGN + HWK(J,2*MAXN+3)) IF (ABS(EMAGN-HWK(J,2*MAXN+3)).GE.TOLM) GO TO 100 DISTM = SQRT((EREAL-HWK(J,2*MAXN+1))**2 + (EIMAG-HWK(J, $2*MAXN+2))**2) IF (DISTM.LT.TOLM) GO TO 110 100 CONTINUE C C It appears only in the Lanczos matrix, so it is good. C NEND = NEND + 1 HWK(NEND,2*MAXN+4) = EREAL HWK(NEND,2*MAXN+5) = EIMAG HWK(NEND,2*MAXN+6) = EMAGN 110 CONTINUE C C Sort the final eigenvalues by decreasing real part. C 120 CALL SHPSRT (NEND,HWK(1,2*MAXN+4),HWK(1,2*MAXN+5),HWK(1,2*MAXN+6)) C C Output the results. C IF (VF.NE.0) WRITE (VF,'(A43,I4)') $ 'Number of Lanczos eigenvalues found : ', NEND IF (TF.NE.0) THEN WRITE (TF,'(A30,I4)') 'Filtered eigenvalues at step: ', NEIG WRITE (TF,'(2E25.17)') (HWK(I,2*MAXN+4),HWK(I,2*MAXN+5),I=1, $NEND) WRITE (TF,'()') END IF NEIG = NEND C C Done. C RETURN END C C********************************************************************** C REAL FUNCTION SULALH(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 SULALX. C REAL NORMA COMMON /SULALX/NORMA C C INTEGER I, J C IF ((I.LT.1).OR.(J.LT.1)) THEN SULALH = 0.0E0 ELSE IF (I.EQ.J) THEN SULALH = NORMA ELSE SULALH = 0.0E0 END IF C RETURN END C C********************************************************************** C SUBROUTINE SULAL1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR, $ VF,IERR,ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,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 "sulal.doc" or "suqmr.doc". C C External routines used: 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 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 sqrdc(x,ldx,n,p,qraux,jpvt,work,job) C LINPACK routine, computes the QR factorization of x. C subroutine sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) C LINPACK routine, applies the QR factorization of x. C subroutine ssvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) C LINPACK routine, computes the SVD of x. C subroutine sulal2 (ndim,nlen,m,dwk,iwk,idx,vecs) C Forces closure of an inner block. C single precision sulalh(i,n) C User-supplied routine, computes inner recurrence coefficients. C C Noel M. Nachtigal C July 9, 1993 C C********************************************************************** C INTRINSIC ABS, AMAX1, MAX0 EXTERNAL SAXPBY, SDOT, SLAMCH, SNRM2, SQRDC, SQRSL, SSVDC EXTERNAL SULAL2, SULALH REAL SDOT, SLAMCH, SNRM2, SULALH C INTEGER IERR, L, LSTAR, M, MAXVW, N, NDIM, NL, NLEN, NLSTAR INTEGER IDX(3,3), IWK(M,4), NUMCHK, VF REAL ADJUST, NORM1, NORMA, TMAX, TMIN, TNRM REAL DWK(M,5*M+7), VECS(NDIM,2*MAXVW+4+2) C C Miscellaneous parameters. C REAL DONE, SZERO PARAMETER (DONE = 1.0E0,SZERO = 0.0E0) C C Local variables. C INTEGER I, IJ, J, LBLKSZ, NP1 REAL STMP, STMP1, STMP2, STMP3, STMP4, STMP5 REAL ISTMP1, ISTMP2, ISTMP3, ISTMP4 LOGICAL IBUILT, INNER, RERUN C C Initialize local variables. C NP1 = N + 1 DWK(IDX(1,N),5*M+7) = SZERO RERUN = .FALSE. IBUILT = .FALSE. IERR = 0 LBLKSZ = N - NL + 1 C C Update the norm estimates. C IF (N.LE.NUMCHK) THEN STMP1 = SNRM2(NLEN,VECS(1,1),1) * DWK(IDX(1,N),5*M+5) STMP2 = SNRM2(NLEN,VECS(1,2),1) * DWK(IDX(1,N),5*M+6) NORMA = AMAX1(STMP1,STMP2,NORMA) END IF NORM1 = AMAX1(ADJUST*NORMA,NORM1) C C Clear current column of H. C DO 100 I = 1, M DWK(IDX(1,I),2*M+IDX(1,N)) = SZERO 100 CONTINUE C C Update D^{(n-1)} to D^{(n)}. C DO 120 I = NL, N-1 STMP = SZERO DO 110 J = NL, N-1 STMP = STMP + DWK(IDX(1,I),IDX(1,J)) * DWK(IDX(1,J),2*M+ $IDX(1,N-1)) 110 CONTINUE STMP = ( DWK(IDX(1,I),M+IDX(1,N-1)) - STMP ) / DWK(IDX(1,N) $,2*M+IDX(1,N-1)) DWK(IDX(1,I),IDX(1,N)) = STMP DWK(IDX(1,N),IDX(1,I)) = STMP * DWK(IDX(1,N),5*M+4) / DWK(IDX(1 $,I),5*M+4) 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 STMP = SZERO DO 130 J = NL, I+1 STMP = STMP + DWK(IDX(1,N),IDX(1,J)) * DWK(IDX(1,J),2*M+ $IDX(1,I)) 130 CONTINUE DWK(IDX(1,N),M+IDX(1,I)) = STMP DWK(IDX(1,I),M+IDX(1,N)) = STMP * DWK(IDX(1,I),5*M+4) / $DWK(IDX(1,N),5*M+4) NORMA = AMAX1(NORMA,ABS(DWK(IDX(1,N),M+IDX(1,I))), $ABS(DWK(IDX(1,I),M+IDX(1,N)))) NORM1 = AMAX1(ADJUST*NORMA,NORM1) 140 CONTINUE C C Check whether D_l is nonsingular. C IF (LBLKSZ.EQ.1) THEN INNER = ABS(DWK(IDX(1,NL),IDX(1,NL))).EQ.SZERO ELSE DO 160 I = NL, N DO 150 J = NL, N DWK(I-NL+1,4*M+J-NL+1) = DWK(IDX(1,I),IDX(1,J)) 150 CONTINUE 160 CONTINUE CALL SSVDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+1),DWK(1,5*M $+3),SZERO,0,SZERO,0,DWK(1,5*M+2),0,I) IF (I.NE.0) THEN IERR = -I GOTO 370 END IF STMP = SZERO IF (DWK(1,5*M+1).NE.SZERO) STMP = DWK(LBLKSZ,5*M+1) / DWK(1,5*M $+1) INNER = STMP.LT.(FLOAT(NLEN) * SLAMCH('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,1) = A v_n and VECS(1,2) = A^T w_n. C DWK(IDX(1,NP1),2*M+IDX(1,N)) = SZERO IF (NL.EQ.1) GO TO 180 IF (NL.EQ.NLSTAR+1) THEN DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) = DWK(IDX(1,NLSTAR),M+IDX(1,N)) $ / DWK(IDX(1,NLSTAR),IDX(1,NLSTAR)) STMP = DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),5*M+ $5) / DWK(IDX(1,N),5*M+5) CALL SAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),-STMP,VECS(1,IDX(2, $NLSTAR))) STMP = DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),5*M+ $6) / DWK(IDX(1,N),5*M+6) * DWK(IDX(1,N),5*M+4) / DWK(IDX(1,NLSTAR) $,5*M+4) CALL SAXPBY (NLEN,VECS(1,2),DONE,VECS(1,2),-STMP,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 STMP2 = SZERO DO 170 I = NL-1, NLSTAR, -1 DWK(IDX(1,I),2*M+IDX(1,N)) = DWK(IDX(1,NL-1),M+IDX(1,N)) * $DWK(IDX(1,I),3*M+IDX(1,L-1)) IF (N.NE.NL) GO TO 170 STMP1 = DWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),5*M+5) CALL SAXPBY (NLEN,VECS(1,IDX(2,NL-1)),STMP2,VECS(1,IDX(2,NL- $1)),STMP1,VECS(1,IDX(2,I))) STMP1 = DWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),5*M+6) / $ DWK(IDX(1,I),5*M+4) CALL SAXPBY (NLEN,VECS(1,IDX(3,NL-1)),STMP2,VECS(1,IDX(3,NL- $1)),STMP1,VECS(1,IDX(3,I))) STMP2 = DONE 170 CONTINUE STMP1 = DWK(IDX(1,NL-1),M+IDX(1,N)) / DWK(IDX(1,N),5*M+5) CALL SAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),-STMP1,VECS(1,IDX(2, $NL-1))) STMP1 = DWK(IDX(1,NL-1),M+IDX(1,N)) * DWK(IDX(1,N),5*M+4) / $DWK(IDX(1,N),5*M+6) CALL SAXPBY (NLEN,VECS(1,2),DONE,VECS(1,2),-STMP1,VECS(1,IDX(3, $NL-1))) END IF C C Compute F_{nn}. C 180 DWK(IDX(1,N),M+IDX(1,N)) = SDOT(NLEN,VECS(1,IDX(3,N)),1,VECS(1,1), $1) * DWK(IDX(1,N),5*M+5) * DWK(IDX(1,N),5*M+6) NORMA = AMAX1(NORMA,ABS(DWK(IDX(1,N),M+IDX(1,N)))) NORM1 = AMAX1(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 DWK(IDX(1,NL),2*M+IDX(1,N)) = DWK(IDX(1,NL),M+IDX(1,N)) / $DWK(IDX(1,NL),IDX(1,NL)) ELSE DO 200 J = NL, N DWK(J-NL+1,5*M+1) = DWK(IDX(1,J),M+IDX(1,N)) DO 190 IJ = NL, N DWK(J-NL+1,4*M+IJ-NL+1) = DWK(IDX(1,J),IDX(1,IJ)) 190 CONTINUE 200 CONTINUE CALL SQRDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),0,SZERO,0 $) CALL SQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M $+1), $ SZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),SZERO, $100,J) DO 210 J = NL, N DWK(IDX(1,J),2*M+IDX(1,N)) = DWK(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 STMP = DWK(IDX(1,N),2*M+IDX(1,N)) CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,1),-STMP,VECS(1, $IDX(2,N))) CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,2),-STMP,VECS(1, $IDX(3,N))) DO 220 I = NL, N-1 STMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+5) / $DWK(IDX(1,N),5*M+5) CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1)),- $STMP,VECS(1,IDX(2,I))) STMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+6) / $DWK(IDX(1,N),5*M+6) * DWK(IDX(1,N),5*M+4) / DWK(IDX(1,I),5*M+4) CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)),- $STMP,VECS(1,IDX(3,I))) 220 CONTINUE C C Compute scale factors for the new vectors. C 230 STMP3 = SNRM2(NLEN,VECS(1,IDX(2,NP1)),1) STMP4 = SNRM2(NLEN,VECS(1,IDX(3,NP1)),1) STMP1 = DWK(IDX(1,N),5*M+5) * STMP3 STMP2 = DWK(IDX(1,N),5*M+6) * STMP4 DWK(IDX(1,NP1),2*M+IDX(1,N)) = STMP1 IF (STMP1.LT.TNRM) IERR = IERR + 16 IF (STMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GOTO 360 DWK(IDX(1,NP1),5*M+4) = DWK(IDX(1,N),5*M+4) * STMP1 / STMP2 DWK(IDX(1,NP1),IDX(1,NP1)) = SDOT(NLEN,VECS(1,IDX(3,NP1)),1,VECS(1 $,IDX(2,NP1)),1) / ( STMP3 * STMP4 ) IF ((STMP3.GE.TMAX).OR.(STMP3.LE.TMIN)) THEN STMP = DONE / STMP3 CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),STMP,VECS(1,IDX(2,NP1)), $SZERO,VECS(1,IDX(2,NP1))) STMP3 = DONE END IF IF ((STMP4.GE.TMAX).OR.(STMP4.LE.TMIN)) THEN STMP = DONE / STMP4 CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),STMP,VECS(1,IDX(3,NP1)), $SZERO,VECS(1,IDX(3,NP1))) STMP4 = DONE END IF DWK(IDX(1,NP1),5*M+5) = DONE / STMP3 DWK(IDX(1,NP1),5*M+6) = DONE / STMP4 C C Compute the last column of D_l^{-1}. C IF (LBLKSZ.EQ.1) THEN DWK(IDX(1,NL),3*M+IDX(1,L)) = DONE / DWK(IDX(1,NL),IDX(1,NL)) ELSE DO 240 I = 1, LBLKSZ-1 DWK(I,5*M+1) = SZERO 240 CONTINUE DWK(LBLKSZ,5*M+1) = DONE CALL SQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M $+1), $ SZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),SZERO, $100,J) DO 250 I = NL, N DWK(IDX(1,I),3*M+IDX(1,L)) = DWK(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 STMP1 = SZERO STMP2 = SZERO DO 260 I = IWK(IDX(1,L-1+1),1), N STMP = ABS(DWK(IDX(1,I),2*M+IDX(1,N))) STMP1 = STMP1 + STMP STMP2 = STMP2 + STMP / DWK(IDX(1,I),5*M+4) 260 CONTINUE STMP2 = STMP2 * DWK(IDX(1,N),5*M+4) 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 STMP3 = SZERO STMP4 = SZERO STMP5 = DWK(IDX(1,NP1),IDX(1,NP1)) * DWK(IDX(1,NP1),2*M+IDX(1,N)) $* DWK(IDX(1,N),5*M+4) / DWK(IDX(1,NP1),5*M+4) DO 270 I = NL, N STMP = ABS(STMP5 * DWK(IDX(1,I),3*M+IDX(1,L))) STMP3 = STMP3 + STMP STMP4 = STMP4 + STMP / DWK(IDX(1,I),5*M+4) 270 CONTINUE STMP4 = STMP4 * DWK(IDX(1,NP1),5*M+4) C C Check H_{n_{l-1}:n,n} and H_{n_l:n-1,n+1}. C STMP = AMAX1(STMP1,STMP2,STMP3,STMP4) INNER = STMP.GT.NORM1 IF (.NOT.INNER) GO TO 360 DWK(IDX(1,N),5*M+7) = STMP C C If H_{n_{l-1}:n,n} is bad, build inner vectors. C IF (AMAX1(STMP1,STMP2).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,1) and VECS(1,2). C IBUILT = .TRUE. DO 280 I = NL, N DWK(IDX(1,I),2*M+IDX(1,NP1)) = SULALH(I,N) STMP = DWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),5*M+5) $ / DWK(IDX(1,N),5*M+5) CALL SAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),-STMP,VECS(1,IDX(2,I $))) STMP = DWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),5*M+6) $ / DWK(IDX(1,N),5*M+6) * DWK(IDX(1,N),5*M+4) / DWK(IDX(1,I),5*M+4) CALL SAXPBY (NLEN,VECS(1,2),DONE,VECS(1,2),-STMP,VECS(1,IDX(3,I $))) 280 CONTINUE ISTMP3 = SNRM2(NLEN,VECS(1,1),1) ISTMP4 = SNRM2(NLEN,VECS(1,2),1) ISTMP1 = DWK(IDX(1,N),5*M+5) * ISTMP3 ISTMP2 = DWK(IDX(1,N),5*M+6) * ISTMP4 DWK(IDX(1,NP1),2*M+IDX(1,NP1)) = ISTMP1 IF (ISTMP1.LT.TNRM) IERR = IERR + 16 IF (ISTMP2.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 STMP = SZERO DO 290 J = NLSTAR, NL-1 STMP = STMP + DWK(IDX(1,NL),IDX(1,J)) * DWK(IDX(1,I),2*M+IDX(1, $N)) 290 CONTINUE DO 300 J = NL, N STMP = STMP + DWK(IDX(1,NL),IDX(1,J)) * DWK(IDX(1,I),2*M+IDX(1, $NP1)) 300 CONTINUE STMP1 = SZERO STMP2 = SZERO STMP3 = AMAX1(STMP3,STMP4) STMP = ( DWK(IDX(1,NL),M+IDX(1,N)) - STMP ) / DWK(IDX(1,NP1),2*M+ $IDX(1,NP1)) STMP5 = STMP * DWK(IDX(1,NL),2*M+IDX(1,NL-1)) * DWK(IDX(1,NL-1), $5*M+4) / DWK(IDX(1,NL),5*M+4) DO 310 I = IWK(IDX(1,L-1+1),1), NL-1 STMP = ABS(DWK(IDX(1,I),3*M+IDX(1,L-1)) * STMP5) STMP1 = STMP1 + STMP STMP2 = STMP2 + STMP / DWK(IDX(1,I),5*M+4) 310 CONTINUE STMP2 = STMP2 * DWK(IDX(1,N),5*M+4) 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 = STMP3.GT.AMAX1(STMP1,STMP2) 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 SULAL2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR,ADJUST, $ NORM1,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 SAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,1),SZERO,VECS(1, $IDX(2,NP1))) CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,2),SZERO,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 DWK(IDX(1,I),2*M+IDX(1,N)) = DWK(IDX(1,I),2*M+IDX(1,NP1)) 340 CONTINUE STMP1 = ISTMP1 STMP2 = ISTMP2 STMP3 = ISTMP3 STMP4 = ISTMP4 ELSE DO 350 I = NL, N DWK(IDX(1,I),2*M+IDX(1,N)) = SULALH(I,N) STMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+5) $/ DWK(IDX(1,N),5*M+5) CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1)) $,-STMP,VECS(1,IDX(2,I))) STMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+6) $/ DWK(IDX(1,N),5*M+6) * DWK(IDX(1,N),5*M+4) / DWK(IDX(1,I),5*M+4) CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)) $,-STMP,VECS(1,IDX(3,I))) 350 CONTINUE STMP3 = SNRM2(NLEN,VECS(1,IDX(2,NP1)),1) STMP4 = SNRM2(NLEN,VECS(1,IDX(3,NP1)),1) STMP1 = DWK(IDX(1,N),5*M+5) * STMP3 STMP2 = DWK(IDX(1,N),5*M+6) * STMP4 END IF DWK(IDX(1,NP1),2*M+IDX(1,N)) = STMP1 IF (STMP1.LT.TNRM) IERR = IERR + 16 IF (STMP2.LT.TNRM) IERR = IERR + 32 IF (IERR.NE.0) GOTO 360 DWK(IDX(1,NP1),5*M+4) = DWK(IDX(1,N),5*M+4) * STMP1 / STMP2 DWK(IDX(1,NP1),IDX(1,NP1)) = SDOT(NLEN,VECS(1,IDX(3,NP1)),1,VECS(1 $,IDX(2,NP1)),1) / ( STMP3 * STMP4 ) IF ((STMP3.GE.TMAX).OR.(STMP3.LE.TMIN)) THEN STMP = DONE / STMP3 CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),STMP,VECS(1,IDX(2,NP1)), $SZERO,VECS(1,IDX(2,NP1))) STMP3 = DONE END IF IF ((STMP4.GE.TMAX).OR.(STMP4.LE.TMIN)) THEN STMP = DONE / STMP4 CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),STMP,VECS(1,IDX(3,NP1)), $SZERO,VECS(1,IDX(3,NP1))) STMP4 = DONE END IF DWK(IDX(1,NP1),5*M+5) = DONE / STMP3 DWK(IDX(1,NP1),5*M+6) = DONE / STMP4 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 SULAL2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR, $ ADJUST,NORM1,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 "sulal.doc" or "suqmr.doc". C C External routines used: C subroutine saxpby(n,z,a,x,b,y) C Library routine, computes z = a * x + b * y. C subroutine sqrdc(x,ldx,n,p,qraux,jpvt,work,job) C LINPACK routine, computes the QR factorization of x. C subroutine sqrsl(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 SAXPBY, SQRDC, SQRSL C INTEGER IERR, L, LSTAR, M, N, NDIM, NL, NLEN, NLSTAR INTEGER IDX(3,3), IWK(M,4), VF REAL ADJUST, DWK(M,5*M+7), NORM1, VECS(NDIM,*) C C Miscellaneous parameters. C REAL DONE, SZERO PARAMETER (DONE = 1.0E0,SZERO = 0.0E0) C C Local variables. C INTEGER I, IJ, J, LBLKSZ, NP1 REAL STMP1, STMP2, SCALV, SCALW 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 STMP1 = DWK(IDX(1,J),5*M+7) DO 100 I = NL+1, N STMP2 = DWK(IDX(1,I),5*M+7) IF (STMP2.GT.SZERO) THEN IF ((STMP1.EQ.SZERO).OR.(STMP2.LT.STMP1)) THEN J = I STMP1 = STMP2 END IF END IF 100 CONTINUE IF (STMP1.EQ.SZERO) THEN IF (VF.NE.0) WRITE (VF,'(A47)') $'... no new norm estimates available (aborting).' IERR = 8 RETURN END IF NORM1 = ADJUST * STMP1 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),5*M+7) = SZERO 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 DWK(IDX(1,NL),2*M+IDX(1,NP1)) = DWK(IDX(1,NL),2*M+IDX(1,N)) DWK(IDX(1,NL),2*M+IDX(1,N)) = DWK(IDX(1,NL),M+IDX(1,N)) / $DWK(IDX(1,NL),IDX(1,NL)) ELSE DO 120 J = NL, N DWK(J-NL+1,5*M+1) = DWK(IDX(1,J),M+IDX(1,N)) DO 110 IJ = NL, N DWK(J-NL+1,4*M+IJ-NL+1) = DWK(IDX(1,J),IDX(1,IJ)) 110 CONTINUE 120 CONTINUE CALL SQRDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),0,SZERO,0 $) CALL SQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M $+1), $ SZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),SZERO, $100,J) DO 130 J = NL, N DWK(IDX(1,J),2*M+IDX(1,NP1)) = DWK(IDX(1,J),2*M+IDX(1,N)) DWK(IDX(1,J),2*M+IDX(1,N)) = DWK(J-NL+1,5*M+2) 130 CONTINUE END IF C C Convert inner vectors to regular vectors. C SCALV = DWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),5*M+5) SCALW = DWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),5*M+6) * $DWK(IDX(1,N),5*M+4) / DWK(IDX(1,NP1),5*M+4) DWK(IDX(1,NP1),2*M+IDX(1,N)) = SZERO DO 140 I = NL, N STMP1 = DWK(IDX(1,I),2*M+IDX(1,N)) - DWK(IDX(1,I),2*M+IDX(1,NP1 $)) STMP2 = STMP1 * DWK(IDX(1,I),5*M+5) / SCALV CALL SAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1)),- $STMP2,VECS(1,IDX(2,I))) STMP2 = STMP1 * DWK(IDX(1,I),5*M+6) / SCALW * DWK(IDX(1,N),5*M+ $4) / DWK(IDX(1,I),5*M+4) CALL SAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)),- $STMP2,VECS(1,IDX(3,I))) 140 CONTINUE C RETURN END C C**********************************************************************