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 ZSLAL (NDIM,NLEN,NLIM,MAXN,MAXVW,M,NORM,ZWK,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 "zslal.doc" in
C     the current directory.
C
C     External routines used:
C     subroutine daxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
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 zslal1(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     subroutine zslalc(maxn,neig,nold,hwk,tf,vf)
C        Library routine, computes common eigenvalues.
C     subroutine zslale(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 DABS, DIMAG, DREAL, DSQRT, MAX0, MOD
      EXTERNAL DAXPBY, DLAMCH, DZNRM2, ZAXPBY, ZDOTU, ZRANDN, ZSLAL1
      EXTERNAL ZSLALC, ZSLALE
      DOUBLE COMPLEX ZDOTU
      DOUBLE PRECISION DLAMCH, DZNRM2
C
      INTEGER INFO(4), M, MAXN, MAXVW, NDIM, NLEN, NLIM
      INTEGER IDX(3,NLIM+2), IWK(M,4)
      DOUBLE COMPLEX VECS(NDIM,1*MAXVW+2+1), ZWK(M,5*M+3)
      DOUBLE PRECISION DWK(M,4), HWK(MAXN,4*MAXN+8), NORM
C
C     Common block variables.
C
C
C     Common block ZSLALX.
C
      DOUBLE PRECISION NORMA
      COMMON /ZSLALX/NORMA
C
C     Miscellaneous parameters.
C
      DOUBLE COMPLEX ZONE, ZZERO
      PARAMETER (ZONE = (1.0D0,0.0D0),ZZERO = (0.0D0,0.0D0))
      DOUBLE PRECISION DONE, DTEN, DZERO
      PARAMETER (DONE = 1.0D0,DTEN = 1.0D1,DZERO = 0.0D0)
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
      DOUBLE PRECISION ADJUST, NORM1, TMAX, TMIN, TNRM
      SAVE             ADJUST, NORM1, TMAX, TMIN, TNRM
C
C     Local variables, transient.
C
      INTEGER I, J, K, NEIG, REVCOM
      DOUBLE COMPLEX ZTMP
      DOUBLE PRECISION DTMP1, DTMP2
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 = 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
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) = 1 + 1 + 0*(MAXVW+2) + MOD(I-1,MAXVW+2)
         IDX(3,I) = IDX(2,I)
 20   CONTINUE
C
C     Check whether any starting vectors must be supplied.
C
      IF (K.EQ.0) CALL ZRANDN (NLEN,VECS(1,1),1)
      CALL ZAXPBY (NLEN,VECS(1,IDX(2,1)),ZONE,VECS(1,1),ZZERO,VECS(1,
     $IDX(2,1)))
CSYM  IF (J.EQ.0) CALL ZRANDN (NLEN,VECS(1,1),1)
CSYM  CALL ZAXPBY (NLEN,VECS(1,IDX(3,1)),ZONE,VECS(1,1),ZZERO,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) = DZERO
            HWK(I,MAXN+J) = DZERO
 30      CONTINUE
 40   CONTINUE
C
C     Scale the first pair of Lanczos vectors and check for invariant
C     subspaces.
C
      DTMP1 = DZNRM2(NLEN,VECS(1,IDX(2,1)),1)
CSYM  DTMP2 = DZNRM2(NLEN,VECS(1,IDX(3,1)),1)
      DTMP2 = DTMP1
      IF (DTMP1.LT.TNRM) IERR = IERR + 16
      IF (DTMP2.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 110
      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)
CSYM     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
      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,1))
C
 60   INFO(2) = 2
      INFO(3) = IDX(3,N)
      INFO(4) = 2
      RETLBL  = 70
CSYM  RETURN
 70   IF ((VF.NE.0).AND.(N.LE.NMAX)) WRITE (VF,'(A11,I8)') 'Rebuilding:'
     $, N+1
      CALL ZSLAL1 (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 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 DAXPBY (N+1,HWK(1,N),DZERO,HWK(1,N),DZERO,HWK(1,N))
      CALL DAXPBY (N+1,HWK(1,MAXN+N),DZERO,HWK(1,MAXN+N),DZERO,HWK(1,
     $MAXN+N))
      DO 80 I = NLSTAR, N+1
         ZTMP       = ZWK(IDX(1,I),2*M+IDX(1,N))
         HWK(I,N) = DREAL(ZTMP)
         HWK(I,MAXN+N) = DIMAG(ZTMP)
 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 ZSLALE (MAXN,NEIG,HWK,TF,VF)
         IF (NEIG.EQ.0) GO TO 100
         IF (NOLD.NE.0) CALL ZSLALC (MAXN,NEIG,NOLD,HWK,TF,VF)
         NOLD = NEIG
         CALL DAXPBY (NOLD,HWK(1,4*MAXN+7),DONE,HWK(1,4*MAXN+4),DZERO,
     $HWK(1,4*MAXN+7))
         CALL DAXPBY (NOLD,HWK(1,4*MAXN+8),DONE,HWK(1,4*MAXN+5),DZERO,
     $HWK(1,4*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 ZSLALC (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   = double 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 dhpsrt(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 DABS
      EXTERNAL DHPSRT
C
      INTEGER MAXN, NEIG, NOLD, TF, VF
      DOUBLE PRECISION HWK(MAXN,4*MAXN+8)
C
C     Local variables.
C
      INTEGER I, J, K
      DOUBLE PRECISION DISTM, EIMAG, EMAGN, EREAL, 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,4*MAXN+4)
         EIMAG = HWK(I,4*MAXN+5)
         EMAGN = HWK(I,4*MAXN+6)
         DO 20 J = 1, NOLD
            DISTM = DSQRT(HWK(J,4*MAXN+7)**2 + HWK(J,4*MAXN+8)**2)
            TOLM  = 0.5D0 * TOL * (EMAGN + DISTM)
            IF (DABS(EMAGN-DISTM).GE.TOLM) GO TO 20
            DISTM = DSQRT((EREAL-HWK(J,4*MAXN+7))**2 + (EIMAG-HWK(J,
     $4*MAXN+8))**2)
            IF (DISTM.GE.TOLM) GO TO 20
            K        = K + 1
            HWK(K,4*MAXN+1) = EREAL
            HWK(K,4*MAXN+2) = EIMAG
            HWK(K,4*MAXN+3) = EMAGN
            GO TO 30
 20      CONTINUE
 30   CONTINUE
C
C     Sort the final eigenvalues by decreasing real part.
C
      CALL DHPSRT (K,HWK(1,4*MAXN+1),HWK(1,4*MAXN+2),HWK(1,4*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,4*MAXN+1),HWK(I,4*MAXN+2),I=1,K)
         WRITE (TF,'()')
      END IF
      GO TO 10
C
      END
C
C**********************************************************************
C
      SUBROUTINE ZSLALE (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   = double 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 cbal(nm,n,ar,ai,low,igh,scale)
C        EISPACK routine, balances a matrix.
C     subroutine comqr(nm,n,low,igh,hr,hi,wr,wi,ierr)
C        EISPACK routine, computes eigenvalues of an upper Hessenberg
C        matrix.
C     subroutine daxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     subroutine dhpsrt(n,wr,wi,wm)
C        Library routine, sorts wr and wi and wm.
C     double precision dlamch(ch)
C        LAPACK routine, computes machine-related constants.
C
C     Noel M. Nachtigal
C     August 6, 1993
C
C**********************************************************************
C
      INTRINSIC DABS, DSQRT
      EXTERNAL CBAL, COMQR, DAXPBY, DHPSRT, DLAMCH
      DOUBLE PRECISION DLAMCH
C
      INTEGER MAXN, NEIG, TF, VF
      DOUBLE PRECISION HWK(MAXN,4*MAXN+8)
C
C     Miscellaneous parameters.
C
      DOUBLE PRECISION DONE, DZERO
      PARAMETER (DONE = 1.0D0,DZERO = 0.0D0)
C
C     Local variables.
C
      INTEGER I, IERR, J, K, NBEG, NCHK, NEND, NLAL
      DOUBLE PRECISION DISTM, EIMAG, EMAGN, EREAL, TOL, TOLM
      LOGICAL MULTIP
C
C     Compute the separation tolerance.
C
      TOL = DSQRT(DLAMCH('E'))
C
C     Extract the NEIG x NEIG matrix.
C
      DO 10 I = 1, NEIG
         CALL DAXPBY (I+1,HWK(1,2*MAXN+I),DONE,HWK(1,I),DZERO,HWK(1,
     $2*MAXN+I))
         CALL DAXPBY (I+1,HWK(1,3*MAXN+I),DONE,HWK(1,MAXN+I),DZERO,HWK(1
     $,3*MAXN+I))
 10   CONTINUE
C
C     Balance it and compute its eigenvalues.
C
      IF (VF.NE.0) WRITE (VF,'(A24)') 'Computing eigenvalues...'
      CALL CBAL (MAXN,NEIG,HWK(1,2*MAXN+1),HWK(1,3*MAXN+1),I,J,HWK(1,
     $4*MAXN+1))
      CALL COMQR (MAXN,NEIG,I,J,HWK(1,2*MAXN+1),HWK(1,3*MAXN+1),HWK(1,
     $4*MAXN+1),HWK(1,4*MAXN+2),IERR)
      IF (IERR.NE.0) THEN
         WRITE (6,'(A25,I2)') 'Error return from COMQR: ', IERR
         RETURN
      END IF
C
C     Compute the magnitude of the eigenvalues.
C
      DO 20 I = 1, NEIG
         HWK(I,4*MAXN+3) = DSQRT(HWK(I,4*MAXN+1)**2 + HWK(I,4*MAXN+2)**2
     $)
 20   CONTINUE
C
C     Sort the eigenvalues in decreasing order of magnitude.
C
      CALL DHPSRT (NEIG,HWK(1,4*MAXN+3),HWK(1,4*MAXN+1),HWK(1,4*MAXN+2))
      IF (TF.NE.0) THEN
         WRITE (TF,'(A33,I4)') 'All Lanczos eigenvalues at step: ', NEIG
         WRITE (TF,'(2E25.17)') (HWK(I,4*MAXN+1),HWK(I,4*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,4*MAXN+1)
      EIMAG  = HWK(I,4*MAXN+2)
      EMAGN  = HWK(I,4*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.5D0 * TOL * (EMAGN + HWK(J,4*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,4*MAXN+3).GE.TOLM) GO TO 60
C
C     Compute the distance between the two eigenvalues.
C
      DISTM = DSQRT((EREAL-HWK(J,4*MAXN+1))**2 + (EIMAG-HWK(J,4*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,4*MAXN+1) = HWK(K+1,4*MAXN+1)
            HWK(K,4*MAXN+2) = HWK(K+1,4*MAXN+2)
            HWK(K,4*MAXN+3) = HWK(K+1,4*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,4*MAXN+4) = EREAL
         HWK(NEND,4*MAXN+5) = EIMAG
         HWK(NEND,4*MAXN+6) = EMAGN
      ELSE
         NBEG = NBEG - 1
         HWK(NBEG,4*MAXN+4) = EREAL
         HWK(NBEG,4*MAXN+5) = EIMAG
         HWK(NBEG,4*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 DAXPBY (I+1,HWK(1,2*MAXN+I),DONE,HWK(2,I+1),DZERO,HWK(1,
     $2*MAXN+I))
         CALL DAXPBY (I+1,HWK(1,3*MAXN+I),DONE,HWK(2,MAXN+I+1),DZERO,
     $HWK(1,3*MAXN+I))
 80   CONTINUE
C
C     Balance it and compute its eigenvalues.
C
      IF (VF.NE.0) WRITE (VF,'(A30)') 'Computing check eigenvalues...'
      CALL CBAL (MAXN,NCHK,HWK(1,2*MAXN+1),HWK(1,3*MAXN+1),I,J,HWK(1,
     $4*MAXN+1))
      CALL COMQR (MAXN,NCHK,I,J,HWK(1,2*MAXN+1),HWK(1,3*MAXN+1),HWK(1,
     $4*MAXN+1),HWK(1,4*MAXN+2),IERR)
      IF (IERR.NE.0) THEN
         WRITE (6,'(A25,I2)') 'Error return from COMQR: ', IERR
         RETURN
      END IF
C
C     Compute the magnitude of the check eigenvalues.
C
      DO 90 I = 1, NCHK
         HWK(I,4*MAXN+3) = DSQRT(HWK(I,4*MAXN+1)**2 + HWK(I,4*MAXN+2)**2
     $)
 90   CONTINUE
C
C     Sort the eigenvalues in decreasing order of magnitude.
C
      CALL DHPSRT (NCHK,HWK(1,4*MAXN+3),HWK(1,4*MAXN+1),HWK(1,4*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,4*MAXN+4)
         EIMAG = HWK(I,4*MAXN+5)
         EMAGN = HWK(I,4*MAXN+6)
         DO 100 J = 1, NCHK
            TOLM = 0.5D0 * TOL * (EMAGN + HWK(J,4*MAXN+3))
            IF (DABS(EMAGN-HWK(J,4*MAXN+3)).GE.TOLM) GO TO 100
            DISTM = DSQRT((EREAL-HWK(J,4*MAXN+1))**2 + (EIMAG-HWK(J,
     $4*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,4*MAXN+4) = EREAL
         HWK(NEND,4*MAXN+5) = EIMAG
         HWK(NEND,4*MAXN+6) = EMAGN
 110  CONTINUE
C
C     Sort the final eigenvalues by decreasing real part.
C
 120  CALL DHPSRT (NEND,HWK(1,4*MAXN+4),HWK(1,4*MAXN+5),HWK(1,4*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,4*MAXN+4),HWK(I,4*MAXN+5),I=1,
     $NEND)
         WRITE (TF,'()')
      END IF
      NEIG = NEND
C
C     Done.
C
      RETURN
      END
C
C**********************************************************************
C
      DOUBLE COMPLEX FUNCTION ZSLALH(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 ZSLALX.
C
      DOUBLE PRECISION NORMA
      COMMON /ZSLALX/NORMA
C
      INTRINSIC DCMPLX
C
      INTEGER I, J
C
      IF ((I.LT.1).OR.(J.LT.1)) THEN
         ZSLALH = (0.0D0,0.0D0)
      ELSE IF (I.EQ.J) THEN
         ZSLALH = DCMPLX(NORMA,0.0D0)
      ELSE
         ZSLALH = (0.0D0,0.0D0)
      END IF
C
      RETURN
      END
C
C**********************************************************************
C
      SUBROUTINE ZSLAL1 (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 "zsqmr.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 zslal2 (ndim,nlen,m,zwk,dwk,iwk,idx,vecs)
C        Forces closure of an inner block.
C     double precision zslalh(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 ZSLAL2, ZSLALH
      DOUBLE COMPLEX ZDOTU, ZSLALH
      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,1*MAXVW+2+1), ZWK(M,5*M+3)
      DOUBLE PRECISION ADJUST, NORM1, NORMA, TMAX, TMIN, TNRM
      DOUBLE PRECISION DWK(M,4)
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,1),1) * DWK(IDX(1,N),2)
CSYM     DTMP2 = DZNRM2(NLEN,VECS(1,1),1) * DWK(IDX(1,N),3)
         DTMP2 = DTMP1
         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,1) = A v_n and VECS(1,1) = 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,1),ZONE,VECS(1,1),-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)
CSYM     CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),-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)
CSYM        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,1),ZONE,VECS(1,1),-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)
CSYM     CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),-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,1)
     $,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,1),-ZWK(IDX(1,N),
     $2*M+IDX(1,N)),VECS(1,IDX(2,N)))
CSYM  CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,1),-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)
CSYM     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)
CSYM  DTMP4      = DZNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
      DTMP4      = DTMP3
      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)
CSYM     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,1) and VECS(1,1).
C
      IBUILT = .TRUE.
      DO 280 I = NL, N
         ZWK(IDX(1,I),2*M+IDX(1,NP1)) = ZSLALH(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,1),ZONE,VECS(1,1),-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)
CSYM     CALL ZAXPBY (NLEN,VECS(1,1),ZONE,VECS(1,1),-ZTMP,VECS(1,IDX(3,I)))
 280  CONTINUE
      IDTMP3       = DZNRM2(NLEN,VECS(1,1),1)
CSYM  IDTMP4       = DZNRM2(NLEN,VECS(1,1),1)
      IDTMP4       = IDTMP3
      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 ZSLAL2 (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,1),ZZERO,VECS(1,
     $IDX(2,NP1)))
CSYM  CALL ZAXPBY (NLEN,VECS(1,IDX(3,NP1)),ZONE,VECS(1,1),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)) = ZSLALH(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)
CSYM        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)
CSYM     DTMP4 = DZNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
         DTMP4 = DTMP3
         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)
CSYM     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 ZSLAL2 (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 "zsqmr.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+3)
      DOUBLE PRECISION ADJUST, DWK(M,4), 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)
CSYM     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**********************************************************************