C**********************************************************************
C
C     Copyright (C) 1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file "cpyrit.doc" in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
C     This  file  contains  the  routines  for  the  QMR  algorithm for
C     unsymmetric  matrices,  using  the  coupled  two-term  recurrence
C     variant of the look-ahead Lanczos algorithm.
C
C**********************************************************************
C
      SUBROUTINE CUCPL (NDIM,NLEN,NLIM,MAXPQ,MAXVW,M,MVEC,NORMS,ZWK,
     $                  DWK,IDX,IWK,VECS,TOL,INFO)
C
C     Purpose:
C     This subroutine uses the QMR algorithm based  on the coupled two-
C     term  variant of  the look-ahead Lanczos  process to solve linear
C     systems.  It runs the algorithm  to convergence or  until a user-
C     specified limit on the number of iterations is reached.
C
C     The  code is  set up  to solve  the system  A x = b  with initial
C     guess x_0 = 0.  Here  A x = b  denotes the preconditioned system,
C     and  it  is  connected with the  original system as follows.  Let
C     B y = c be the original unpreconditioned system to be solved, and
C     let y_0 be an arbitrary initial guess for its solution.  Then:
C          A x = b, where  A = M_1^{-1} B M_2^{-1},
C                          x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
C     Here M = M_1 M_2 is the preconditioner.
C
C     To recover the final iterate y_n for  the original system B y = c
C     from the final iterate x_n for the preconditioned system A x = b,
C     set
C               y_n = y_0 + M_2^{-1} x_n.
C
C     The algorithm  was first described  in the RIACS Technical Report
C     92.15, `An Implementation of the QMR Method Based on Coupled Two-
C     Term Recurrences`,  June 1992.  A good reference  for the details
C     of  this implementation  is the  RIACS  Technical  Report  92.19,
C     `Implementation  details of the coupled  QMR algorithm`,  by R.W.
C     Freund and N.M. Nachtigal, October 1992.
C
C     Parameters:
C     For a description of  the parameters, see the file "cucpl.doc" in
C     the current directory.
C
C     External routines used:
C     real slamch(ch)
C        LAPACK routine, computes machine-related constants.
C     real scnrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine caxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     real cdotu(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C     subroutine cqrdc(x,ldx,n,p,qraux,jpvt,work,job)
C        LINPACK routine, computes the QR factorization of x.
C     subroutine cqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
C        LINPACK routine, applies the QR factorization of x.
C     subroutine crandn(n,x,seed)
C        Library routine, fill x with random numbers.
C     subroutine crotg(a,b,cos,sin)
C        BLAS-1 routine, computes the Givens rotation which rotates the
C        vector [a; b] into [ sqrt(a**2 + b**2); 0 ].
C     subroutine csvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
C        LINPACK routine, computes the SVD of x.
C     subroutine cucpl1(ndim,nlen,m,n,k,kstar,l,lstar,mk,mkstar,nl,
C        nlstar,vf,ierr,adjust,norm1,norm2,zwk,dwk,iwk,idx,vecs)
C        Low-level routine, rebuilds vectors PQ in CPL.
C     subroutine cucpl2(ndim,nlen,m,n,k,kstar,l,lstar,mk,mkstar,nl,
C        nlstar,vf,ierr,adjust,norm1,norm2,zwk,dwk,iwk,idx,vecs)
C        Low-level routine, rebuilds vectors VW in CPL.
C     real function cucpll(i,n)
C        User-supplied routine, computes inner recurrence coefficients.
C     real cucplom(n)
C        User-supplied routine, specifies the QMR scaling factors.
C     real function cucplu(i,n)
C        User-supplied routine, computes inner recurrence coefficients.
C
C     Noel M. Nachtigal
C     March 1, 1992
C
C**********************************************************************
C
      INTRINSIC CABS, ABS, FLOAT, CMPLX, CONJG, AMAX1, REAL, SQRT
      INTRINSIC MAX0, MIN0, MOD
      EXTERNAL SLAMCH, SCNRM2, CAXPBY, CDOTU, CQRDC, CQRSL, CRANDN
      EXTERNAL CROTG, CSVDC 
      EXTERNAL CUCPL1, CUCPL2, CUCPLL, CUCPLO, CUCPLU
      COMPLEX CDOTU, CUCPLL, CUCPLU
      REAL SLAMCH, SCNRM2, CUCPLO
C
      INTEGER INFO(4), M, MAXPQ, MAXVW, MVEC, NDIM, NLEN, NLIM
      INTEGER IDX(6,NLIM+2), IWK(M,13)
      COMPLEX VECS(NDIM,5*MVEC+4-1), ZWK(M,8*M+7)
      REAL DWK(M,11), NORMS(2), TOL
C
C     Common block variables.
C
C
C     Common block CUCPLX.
C
      REAL NORMA
      COMMON /CUCPLX/NORMA
C
C     Miscellaneous parameters.
C
      COMPLEX CONE, CZERO
      PARAMETER (CONE = (1.0E0,0.0E0),CZERO = (0.0E0,0.0E0))
      REAL SHUN, SONE, STEN, SZERO
      PARAMETER (SHUN = 1.0E2,SONE = 1.0E0,STEN = 1.0E1,SZERO = 0.0E0)
C
C     Local variables, permanent.
C
      LOGICAL IBUILT, INNER, RERUN
      SAVE    IBUILT, INNER, RERUN
      INTEGER IEND, IERR, K, KBLKSZ, KSTAR, L, LBLKSZ, LSTAR, MK
      SAVE    IEND, IERR, K, KBLKSZ, KSTAR, L, LBLKSZ, LSTAR, MK
      INTEGER MKMAX, MKSTAR, MPQBLT, MVWBLT, N, NL, NLMAX, NLSTAR, NMAX
      SAVE    MKMAX, MKSTAR, MPQBLT, MVWBLT, N, NL, NLMAX, NLSTAR, NMAX
      INTEGER NP1, NQMR, NUMCHK, PQBEG, RETLBL, TF, TRES, VF, VWBEG
      SAVE    NP1, NQMR, NUMCHK, PQBEG, RETLBL, TF, TRES, VF, VWBEG
      REAL ADJUST, MAXOMG, NORM1, NORM2, TMAX, TMIN, TNRM
      SAVE             ADJUST, MAXOMG, NORM1, NORM2, TMAX, TMIN, TNRM
      REAL R0, RESN, UCHK, UNRM
      SAVE             R0, RESN, UCHK, UNRM
C
C     Local variables, transient.
C
      INTEGER I, IBASE, IBLKSZ, IJ, INIT, J, KEND, LEND, MI, MIP1
      INTEGER NI, NIP1, REVCOM
      COMPLEX CTMP, CTMP1, CTMP2
      REAL STMP, STMP1, STMP2, STMP3, STMP4
      REAL ISTMP1, ISTMP2, ISTMP3, ISTMP4
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      200    returning from AXB, go to label 200
C           1      350    returning from AXB, go to label 350
C           1      720    returning from AXB, go to label 720
C           2      370    returning from ATXB, go to label 370
C
      REVCOM  = INFO(2)
      INFO(2) = 0
      IF (REVCOM.EQ.0) THEN
         N      = 0
         NQMR   = 0
         MPQBLT = 0
         MVWBLT = 0
         IF (RETLBL.EQ.0) GO TO 10
      ELSE IF (REVCOM.EQ.1) THEN
         IF (RETLBL.EQ.200) THEN
            GO TO 200
         ELSE IF (RETLBL.EQ.350) THEN
            GO TO 350
         ELSE IF (RETLBL.EQ.720) THEN
            GO TO 720
         END IF
      ELSE IF (REVCOM.EQ.2) THEN
         IF (RETLBL.EQ.370) GO TO 370
      END IF
      IERR = 1
      GO TO 750
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 (MAXPQ.LT.1)         IERR = 2
      IF (MAXVW.LT.1)         IERR = 2
      IF (NLEN.GT.NDIM)       IERR = 2
      IF (MVEC.LT.MAXPQ+1)    IERR = 2
      IF (MVEC.LT.MAXVW+1)    IERR = 2
      IF (M.LT.MAXVW+MAXPQ+2) IERR = 2
      IF (IERR.NE.0) GO TO 750
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 = SZERO
      NORM1 = ABS(NORMS(1))
      NORM2 = ABS(NORMS(2))
C
C     Set the adjustment parameters.
C
      NUMCHK = 25
      ADJUST = STEN
C
C     Extract and check the various tolerances and norm estimates.
C
      TNRM = SLAMCH('E') * STEN
      TMIN = SQRT(SQRT(SLAMCH('S')))
      TMAX = SONE / TMIN
      IF (TOL.LE.SZERO) TOL = SQRT(SLAMCH('E'))
C
C     Start the trace messages and convergence history.
C
      IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') 0, SONE, SONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, SONE, SONE
C
C     Set 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        IDX(4,I) = indices used for p_i;
C        IDX(5,I) = indices used for q_i;
C        IDX(6,I) = indices used for s_i.
C
      DO 20 I = 1, NLIM+2
         IDX(1,I) = MOD(I-1,M) + 1
         IDX(2,I) = 4 + 1 + 0*MVEC + MOD(I-1,MVEC)
         IDX(3,I) = 4 + 1 + 1*MVEC + MOD(I-1,MVEC)
         IDX(4,I) = 4 + 1 + 2*MVEC + MOD(I-1,MVEC)
         IDX(5,I) = 4 + 1 + 3*MVEC + MOD(I-1,MVEC)
         IDX(6,I) = 4 + 1 + 4*MVEC + MOD(I-1,MVEC-1)
 20   CONTINUE
C
C     Set x_0 = 0 and compute the norm of the initial residual.
C
      CALL CAXPBY (NLEN,VECS(1,IDX(2,1)),CONE,VECS(1,2),CZERO,VECS(1,
     $IDX(2,1)))
      CALL CAXPBY (NLEN,VECS(1,1),CZERO,VECS(1,1),CZERO,VECS(1,1))
      R0 = SCNRM2(NLEN,VECS(1,IDX(2,1)),1)
      IF ((TOL.GE.SONE).OR.(R0.EQ.SZERO)) GO TO 750
C
C     Check whether the auxiliary vector must be supplied.
C
      IF (INIT.EQ.0) CALL CRANDN (NLEN,VECS(1,3),1)
      CALL CAXPBY (NLEN,VECS(1,IDX(3,1)),CONE,VECS(1,3),CZERO,VECS(1,
     $IDX(3,1)))
C
C     Scale the first pair of Lanczos vectors and check for invariant
C     subspaces.
C
      STMP1 = R0
      STMP2 = SCNRM2(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 750
      DWK(IDX(1,1),3)   = SONE
      ZWK(IDX(1,1),IDX(1,1)) = CDOTU(NLEN,VECS(1,IDX(2,1)),1,VECS(1,
     $IDX(3,1)),1) / ( STMP1 * STMP2 )
      IF ((STMP1.GE.TMAX).OR.(STMP1.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP1,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(2,1)),CTMP,VECS(1,IDX(2,1)),CZERO,
     $VECS(1,IDX(2,1)))
         STMP1 = SONE
      END IF
      IF ((STMP2.GE.TMAX).OR.(STMP2.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP2,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,1)),CTMP,VECS(1,IDX(3,1)),CZERO,
     $VECS(1,IDX(3,1)))
         STMP2 = SONE
      END IF
      DWK(IDX(1,1),4) = SONE / STMP1
      DWK(IDX(1,1),5) = SONE / STMP2
C
C     Initialize the counters.
C
      K      = 0
      L      = 1
      N      = 1
      NMAX   = 0
      KSTAR  = 1
      LSTAR  = 1
      MKMAX  = 0
      NLMAX  = 1
      PQBEG  = 1
      VWBEG  = 1
      IWK(IDX(1,0+1),2) = 1
      IWK(IDX(1,1+1),2) = 1
      IWK(IDX(1,0+1),1) = 1
      IWK(IDX(1,1+1),1) = 1
      MPQBLT = 1
      MVWBLT = 1
      MK     = IWK(IDX(1,K+1),1)
      NL     = IWK(IDX(1,L+1),2)
      IWK(IDX(1,N),3) = NL
      NLSTAR = IWK(IDX(1,LSTAR+1),2)
      IWK(IDX(1,N),5) = NLSTAR
C
C     Set up the QMR iteration.
C
      RESN   = SONE
      DWK(IDX(1,1),1) = CUCPLO(1)
      ZWK(IDX(1,1),8*M+4) = DWK(IDX(1,1),1) * R0
      MAXOMG = SONE / DWK(IDX(1,1),1)
C
C     This is one step of the coupled two-term Lanczos algorithm.
C     V_n,  W_n,  P_{n-1}, Q_{n-1}, D_{n-1}, E_{n-1}, F_{n-1}, L_{n-1},
C     U_{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, MK.LE.N-1, N-NL.LE.MAXVW-1, N-MK.LE.MAXPQ.
C
C     The references in the comments  to step numbers correspond to the
C     description in Section 4 of the report `Implementation details of
C     the coupled QMR algorithm`, by Freund and Nachtigal, RIACS report
C     92.19, October 1992.
C
 30   IWK(IDX(1,N),9)   = L
      IWK(IDX(1,N),7)   = N
      IWK(IDX(1,N),8)  = K
      IWK(IDX(1,N),12)  = K
      IWK(IDX(1,N),10) = KSTAR
      IWK(IDX(1,N),13) = KSTAR
      IWK(IDX(1,N),11)  = LSTAR
C
C     Build p_n and q_n.
C     There are three possible cases:
C       N > NMAX   : the PQ sequence is being built new.
C       N > MKMAX  : the PQ sequence is being rebuilt, build as if new.
C       N <= MKMAX : the VW sequence is being rebuilt, just rerun PQ.
C     The first two are identical in terms of code.
C     The code is inlined because of the matrix multiplications.
C
C     Initialize Lanczos step variables.
C
      NP1        = N + 1
      DWK(IDX(1,N),11) = SZERO
      KBLKSZ     = N - MK
      IBUILT     = .FALSE.
      IERR       = 0
C
C     Clear current column of U.
C
      DO 40 I = 1, M
         ZWK(IDX(1,I),6*M+IDX(1,N)) = CZERO
 40   CONTINUE
      ZWK(IDX(1,N),6*M+IDX(1,N)) = CONE
C
C     The first step is different, deal with it here.
C
      IF (N.EQ.1) THEN
         NP1     = 2
         MKSTAR  = 1
         DWK(IDX(1,1),9) = SONE
         DWK(IDX(1,1),10) = SONE
         DWK(IDX(1,1),7)  = DWK(IDX(1,1),4)
         DWK(IDX(1,1),8)  = DWK(IDX(1,1),5)
         INNER   = .FALSE.
         CALL CAXPBY (NLEN,VECS(1,IDX(4,1)),CONE,VECS(1,IDX(2,1)),CZERO,
     $VECS(1,IDX(4,1)))
         CALL CAXPBY (NLEN,VECS(1,IDX(5,1)),CONE,VECS(1,IDX(3,1)),CZERO,
     $VECS(1,IDX(5,1)))
         GO TO 340
      END IF
C
C     Step 1:
C     Update D^{(n-1)} to D^{(n)}.
C
      DO 60 I = NL, N-1
         CTMP = CZERO
         DO 50 J = MAX0(NLSTAR,IWK(IDX(1,I),3)), N-1
            CTMP = CTMP + ZWK(IDX(1,I),IDX(1,J)) * ZWK(IDX(1,J),2*M+
     $IDX(1,N-1))
 50      CONTINUE
         CTMP     = ( ZWK(IDX(1,I),M+IDX(1,N-1)) - CTMP ) / ZWK(IDX(1,N)
     $,2*M+IDX(1,N-1))
         ZWK(IDX(1,I),IDX(1,N)) = CTMP
         ZWK(IDX(1,N),IDX(1,I)) = CTMP * DWK(IDX(1,N),3) / DWK(IDX(1,I),
     $3)
 60   CONTINUE
C
C     Step 2:
C     Compute k^\star.
C
      I = KSTAR
      DO 70 J = I+1, K
         IF (IWK(IDX(1,J+1),1).LE.NL-1) KSTAR = J
 70   CONTINUE
      MKSTAR = IWK(IDX(1,KSTAR+1),1)
C
C     Check memory allocation  (MKSTAR.GE.PQBEG).  Update  the range of
C     valid indices:  the new vectors  p_n and q_n  will overwrite  the
C     vectors at location N-MVEC.
C
      IF (MKSTAR.LT.PQBEG) THEN
         IERR = 64
         GO TO 750
      END IF
      PQBEG = MAX0(PQBEG,N-MVEC+1)
C
C     Step 3:
C     Compute F_{n,1:n-1} = F_{n,m_{k^\star}:n-1}.
C
      DO 90 I = MKSTAR, N-1
         CTMP = CZERO
         DO 80 J = MAX0(NL,IWK(IDX(1,I),5)), I+1
            CTMP = CTMP + ZWK(IDX(1,N),IDX(1,J)) * ZWK(IDX(1,J),2*M+
     $IDX(1,I))
 80      CONTINUE
         ZWK(IDX(1,N),M+IDX(1,I)) = CTMP
         STMP     = CABS(CTMP) / DWK(IDX(1,I),9)
         NORMA    = AMAX1(NORMA,STMP)
         NORM1    = AMAX1(ADJUST*NORMA,NORM1)
         NORM2    = AMAX1(ADJUST*NORMA,NORM2)
 90   CONTINUE
C
C     Step 4:
C     Check whether E_k is nonsingular.
C
      IF (KBLKSZ.EQ.1) THEN
         INNER = CABS(ZWK(IDX(1,MK),5*M+IDX(1,MK))).EQ.SZERO
      ELSE
         DO 110 I = MK, N-1
            DO 100 J = MK, N-1
               ZWK(I-MK+1,4*M+J-MK+1) = ZWK(IDX(1,I),5*M+IDX(1,J))
 100        CONTINUE
 110     CONTINUE
         CALL CSVDC (ZWK(1,4*M+1),M,KBLKSZ,KBLKSZ,ZWK(1,8*M+1),ZWK(1,8*M
     $+3),CZERO,0,CZERO,0,ZWK(1,8*M+2),0,I)
         IF (I.NE.0) THEN
            IERR = -I
            GO TO 750
         END IF
         STMP = SZERO
         IF (REAL(ZWK(1,8*M+1)).NE.SZERO) THEN
            STMP = ABS( REAL(ZWK(KBLKSZ,8*M+1)) / REAL(ZWK(1,8*M+1)))
         END IF
         INNER = STMP.LT.(FLOAT(NLEN) * SLAMCH('E'))
      END IF
      KEND = K
      IF (INNER) KEND = K - 1
C
C     Step 5:
C     Compute U_{m_i:m_{i+1}-1,n}, for i = k^\star, ..., kend, kend = k
C     or kend = k-1.
C
      IF (KEND.EQ.K) IWK(IDX(1,K+1+1),1) = N
      DO 150 I = KSTAR, KEND
         MI     = IWK(IDX(1,I+1),1)
         MIP1   = IWK(IDX(1,I+1+1),1)
         IBLKSZ = MIP1 - MI
         IF (IBLKSZ.EQ.1) THEN
            ZWK(IDX(1,MI),6*M+IDX(1,N)) = ZWK(IDX(1,N),M+IDX(1,MI)) / 
     $ZWK(IDX(1,MI),5*M+IDX(1,MI)) * DWK(IDX(1,MI),3) / DWK(IDX(1,N),3)
         ELSE
            DO 130 J = MI, MIP1-1
               ZWK(J-MI+1,8*M+1) = ZWK(IDX(1,N),M+IDX(1,J)) * DWK(IDX(1,
     $J),3) / DWK(IDX(1,N),3)
               DO 120 IJ = MI, MIP1-1
                  ZWK(J-MI+1,4*M+IJ-MI+1) = ZWK(IDX(1,J),5*M+IDX(1,IJ))
 120           CONTINUE
 130        CONTINUE
            CALL CQRDC (ZWK(1,4*M+1),M,IBLKSZ,IBLKSZ,ZWK(1,8*M+3),0,
     $CZERO,0)
            CALL CQRSL (ZWK(1,4*M+1),M,IBLKSZ,IBLKSZ,ZWK(1,8*M+3),ZWK(1,
     $8*M+1),
     $                  CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),
     $CZERO,100,J)
            DO 140 J = MI, MIP1-1
               ZWK(IDX(1,J),6*M+IDX(1,N)) = ZWK(J-MI+1,8*M+2)
 140        CONTINUE
         END IF
 150  CONTINUE
C
C     Step 6:
C     Build the part common to both inner and regular vectors.
C
      DWK(IDX(1,N),7) = DWK(IDX(1,N),4)
      DWK(IDX(1,N),8) = DWK(IDX(1,N),5)
      CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,IDX(2,N)),CZERO,VECS(1,3))
      CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,IDX(3,N)),CZERO,VECS(1,4))
      DO 160 I = MKSTAR, MK-1
         CTMP = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),7) / DWK(IDX(1
     $,N),7)
         CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,3),-CTMP,VECS(1,IDX(4,I
     $)))
         CTMP = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),8) / DWK(IDX(1
     $,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,4),-CTMP,VECS(1,IDX(5,I
     $)))
 160  CONTINUE
C
C     Check whether PQ is being rerun.
C
      RERUN = (N.LE.NMAX).AND.(N.LE.MKMAX)
      IF (RERUN) THEN
C
C     Check whether E_k has become singular (this should not happen).
C
         IF ((IWK(IDX(1,K+1+1),1).EQ.N).AND.INNER) THEN
            IERR = 128
            GO TO 750
         END IF
C
C     Determine whether this was an inner vector or not.
C
         INNER = IWK(IDX(1,K+1+1),1).NE.N
         IF (VF.NE.0) THEN
            IF (INNER) THEN
               WRITE (VF,'(A17,I5,A14)') 'Rerunning vector ',N,
     $' (PQ) as inner'
            ELSE
               WRITE (VF,'(A17,I5,A16)') 'Rerunning vector ',N,
     $' (PQ) as regular'
            END IF
         END IF
C
C     If building an inner vector, set the coefficients U_{m_k:n-1,n}.
C
         IF (INNER) THEN
            DO 170 I = MK, N-1
               ZWK(IDX(1,I),6*M+IDX(1,N)) = CUCPLU(I,N)
 170        CONTINUE
         END IF
      ELSE
C
C     PQ is not being rerun, determine whether it is being rebuilt.
C
         IF ((VF.NE.0).AND.(N.LE.NMAX)) THEN
            WRITE (VF,'(A15,I8)') 'Rebuilding P&Q:', N
         END IF
C
C     If E_k is singular, build an inner vector.
C
         IF (INNER) THEN
            IF (VF.NE.0) WRITE (VF,'(A31)') 
     $'... moment matrix E is singular'
            GO TO 300
         END IF
      END IF
C
C     Either E_k is nonsingular, or PQ is being rerun.  For the latter,
C     just finish building  the vectors.  For the former,  the check in
C     Step 7 (G_{m_{k-1}:n-1,n-1})  could be done.  However,  the look-
C     ahead strategy requires  the smaller of the checks  in Step 7 and
C     Step 9, in case the norm  estimates need updating.  Hence, Step 7
C     and Step 9  are switched,  and their  checks  are later performed
C     together.
C
C     Leave the common part of the vectors in the temporary vectors.
C
      CALL CAXPBY (NLEN,VECS(1,IDX(4,N)),CONE,VECS(1,3),CZERO,VECS(1,
     $IDX(4,N)))
      CALL CAXPBY (NLEN,VECS(1,IDX(5,N)),CONE,VECS(1,4),CZERO,VECS(1,
     $IDX(5,N)))
C
C     Step 8:
C     Build regular vectors and compute A p_n and q_n^T A p_n.
C
      DO 180 I = MK, N-1
         CTMP = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),7) / DWK(IDX(1
     $,N),7)
         CALL CAXPBY (NLEN,VECS(1,IDX(4,N)),CONE,VECS(1,IDX(4,N)),-CTMP,
     $VECS(1,IDX(4,I)))
         CTMP = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),8) / DWK(IDX(1
     $,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,IDX(5,N)),CONE,VECS(1,IDX(5,N)),-CTMP,
     $VECS(1,IDX(5,I)))
 180  CONTINUE
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,IDX(4,N)),VECS(1,IDX(2,NP1)))
C
 190  INFO(2) = 1
      INFO(3) = IDX(4,N)
      INFO(4) = IDX(2,NP1)
      RETLBL  = 200
      RETURN
 200  ZWK(IDX(1,N),5*M+IDX(1,N)) = CDOTU(NLEN,VECS(1,IDX(5,N)),1,VECS(1,
     $IDX(2,NP1)),1) * DWK(IDX(1,N),7) * DWK(IDX(1,N),8)
C
C     Step 9:
C     Compute the last column of E_k^{-1}.
C
      IF (KBLKSZ.EQ.1) THEN
         ZWK(IDX(1,MK),7*M+IDX(1,K)) = CONE / ZWK(IDX(1,MK),5*M+IDX(1,MK
     $))
      ELSE
         DO 210 I = 1, KBLKSZ-1
            ZWK(I,8*M+1) = CZERO
 210     CONTINUE
         ZWK(KBLKSZ,8*M+1) = CONE
         CALL CQRSL (ZWK(1,4*M+1),M,KBLKSZ,KBLKSZ,ZWK(1,8*M+3),ZWK(1,8*M
     $+1),
     $               CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),CZERO,
     $100,J)
         DO 220 I = MK, N-1
            ZWK(IDX(1,I),7*M+IDX(1,K)) = ZWK(I-MK+1,8*M+2)
 220     CONTINUE
      END IF
C
C     Check for zero norms.
C
      DWK(IDX(1,N),9) = SCNRM2(NLEN,VECS(1,IDX(4,N)),1) * DWK(IDX(1,N),7
     $)
      DWK(IDX(1,N),10) = SCNRM2(NLEN,VECS(1,IDX(5,N)),1) * DWK(IDX(1,N),
     $8)
      IF (DWK(IDX(1,N),9).EQ.SZERO) GO TO 360
      IF (DWK(IDX(1,N),10).EQ.SZERO) GO TO 360
C
C     If PQ is being rerun, then skip to the end.
C
      IF (RERUN) GO TO 360
C
C     Step 7:
C     Build the full column G_{m_{k-1}:n-1,n-1} and compute its norm.
C
      STMP1 = SZERO
      STMP2 = SZERO
      DO 240 I = IWK(IDX(1,K-1+1),1), N-1
         CTMP = SZERO
         DO 230 J = MAX0(I,NLSTAR), N
            CTMP = CTMP + ZWK(IDX(1,I),6*M+IDX(1,J)) * ZWK(IDX(1,J),2*M+
     $IDX(1,N-1))
 230     CONTINUE
         STMP  = CABS(CTMP)
         STMP1 = STMP1 + STMP * DWK(IDX(1,I),9)
         STMP2 = STMP2 + STMP * DWK(IDX(1,I),10) / DWK(IDX(1,I),3)
 240  CONTINUE
      STMP1 = STMP1 / DWK(IDX(1,N-1),9)
      STMP2 = STMP2 * DWK(IDX(1,N-1),3) / DWK(IDX(1,N-1),10)
C
C     Step 9:
C     Build the 2nd term for the next step, G_{m_k:n-1,n}, regular.
C     Compute the norm of G_{m_k:n-1,n}.
C
      STMP3 = SZERO
      STMP4 = SZERO
      CTMP  = ZWK(IDX(1,N),5*M+IDX(1,N)) * ZWK(IDX(1,N),2*M+IDX(1,N-1)) 
     $* DWK(IDX(1,N-1),3) / DWK(IDX(1,N),3)
      DO 250 I = MK, N-1
         STMP  = CABS(CTMP * ZWK(IDX(1,I),7*M+IDX(1,K)))
         STMP3 = STMP3 + STMP * DWK(IDX(1,I),9)
         STMP4 = STMP4 + STMP * DWK(IDX(1,I),10) / DWK(IDX(1,I),3)
 250  CONTINUE
      STMP3 = STMP3 / DWK(IDX(1,N),9)
      STMP4 = STMP4 * DWK(IDX(1,N),3) / DWK(IDX(1,N),10)
C
C     Steps 7 and 9:
C     Check G_{m_{k-1}:n-1,n-1} and G_{m_k:n-1,n}.
C
      STMP  = AMAX1(STMP1,STMP2,STMP3,STMP4)
      INNER = STMP.GT.NORM1
      IF (.NOT.INNER) GO TO 360
      DWK(IDX(1,N),11) = STMP
C
C     If G_{m_{k-1}:n-1,n-1} is bad, build inner vectors.
C
      IF (AMAX1(STMP1,STMP2).GT.NORM1) GO TO 300
C
C     If G_{m_k:n-1,n} is bad, check the inner vectors.
C     This only applies if m_k > 1.
C
      IF (MK.LE.1) GO TO 300
C
C     Build the inner vectors to compute the 2nd term at the next step.
C     Get the coefficients U_{m_k:n-1,n} in a temporary location.
C     Build the inner vectors, their norms are needed.
C
      IBUILT = .TRUE.
      DO 260 I = MK, N-1
         ZWK(IDX(1,I),6*M+IDX(1,NP1)) = CUCPLU(I,N)
         CTMP       = ZWK(IDX(1,I),6*M+IDX(1,NP1)) * DWK(IDX(1,I),7) / 
     $DWK(IDX(1,N),7)
         CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,3),-CTMP,VECS(1,IDX(4,I
     $)))
         CTMP       = ZWK(IDX(1,I),6*M+IDX(1,NP1)) * DWK(IDX(1,I),8) / 
     $DWK(IDX(1,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,4),-CTMP,VECS(1,IDX(5,I
     $)))
 260  CONTINUE
      DWK(IDX(1,NP1),9) = SCNRM2(NLEN,VECS(1,3),1) * DWK(IDX(1,N),7)
      DWK(IDX(1,NP1),10) = SCNRM2(NLEN,VECS(1,4),1) * DWK(IDX(1,N),8)
      IF (DWK(IDX(1,NP1),9).EQ.SZERO) GO TO 310
      IF (DWK(IDX(1,NP1),10).EQ.SZERO) GO TO 310
C
C     Build the 2nd term for the next step, G_{m_{k-1}:m_k-1,n}, inner.
C     Compute the norm of G_{m_{k-1}:m_k-1,n}.
C
      CTMP = CZERO
      DO 270 J = MKSTAR, MK-1
         CTMP = CTMP + ZWK(IDX(1,J),6*M+IDX(1,N)) * ZWK(IDX(1,J),5*M+
     $IDX(1,MK)) / DWK(IDX(1,J),3)
 270  CONTINUE
      DO 280 J = MK, N-1
         CTMP = CTMP + ZWK(IDX(1,J),6*M+IDX(1,NP1)) * ZWK(IDX(1,J),5*M+
     $IDX(1,MK)) / DWK(IDX(1,J),3)
 280  CONTINUE
      STMP1 = SZERO
      STMP2 = SZERO
      STMP3 = AMAX1(STMP3,STMP4)
      CTMP  = ZWK(IDX(1,N),M+IDX(1,MK)) - CTMP * DWK(IDX(1,N),3)
      CTMP  = CTMP * ZWK(IDX(1,MK),2*M+IDX(1,MK-1)) * DWK(IDX(1,MK-1),3)
     $ / DWK(IDX(1,MK),3)
      DO 290 I = IWK(IDX(1,K-1+1),1), MK-1
         STMP  = CABS(CTMP * ZWK(IDX(1,I),7*M+IDX(1,K-1)))
         STMP1 = STMP1 + STMP * DWK(IDX(1,I),9)
         STMP2 = STMP2 + STMP * DWK(IDX(1,I),10) / DWK(IDX(1,I),3)
 290  CONTINUE
      STMP1 = STMP1 / DWK(IDX(1,NP1),9)
      STMP2 = STMP2 * DWK(IDX(1,N),3) / DWK(IDX(1,NP1),10)
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 P&Q block has to be forced to close.
C
 300  IF (VF.NE.0) WRITE (VF,'(A7,I5,A14)') 'Vector ',N,' (PQ) is inner'
      IF (N-MK.EQ.MAXPQ) THEN
         CALL CUCPL1 (NDIM,NLEN,M,N,K,KSTAR,L,LSTAR,MK,MKSTAR,NL,NLSTAR,
     $                VF,IERR,ADJUST,NORM1,NORM2,ZWK,DWK,IWK,IDX,VECS)
         IF (IERR.NE.0) GO TO 750
         INNER  = .FALSE.
         KBLKSZ = N - MK
         RERUN  = .TRUE.
         NP1    = N + 1
         GO TO 190
      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
 310  CALL CAXPBY (NLEN,VECS(1,IDX(4,N)),CONE,VECS(1,3),CZERO,VECS(1,
     $IDX(4,N)))
      CALL CAXPBY (NLEN,VECS(1,IDX(5,N)),CONE,VECS(1,4),CZERO,VECS(1,
     $IDX(5,N)))
C
C     Step 11:
C     Get the coefficients U_{m_k:n-1,n} and build inner vectors.
C
      IF (IBUILT) THEN
         DO 320 I = MK, N-1
            ZWK(IDX(1,I),6*M+IDX(1,N)) = ZWK(IDX(1,I),6*M+IDX(1,NP1))
 320     CONTINUE
         DWK(IDX(1,N),9) = DWK(IDX(1,NP1),9)
         DWK(IDX(1,N),10) = DWK(IDX(1,NP1),10)
      ELSE
         DO 330 I = MK, N-1
            ZWK(IDX(1,I),6*M+IDX(1,N)) = CUCPLU(I,N)
            CTMP     = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),7) / 
     $DWK(IDX(1,N),7)
            CALL CAXPBY (NLEN,VECS(1,IDX(4,N)),CONE,VECS(1,IDX(4,N)),-
     $CTMP,VECS(1,IDX(4,I)))
            CTMP     = ZWK(IDX(1,I),6*M+IDX(1,N)) * DWK(IDX(1,I),8) / 
     $DWK(IDX(1,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
            CALL CAXPBY (NLEN,VECS(1,IDX(5,N)),CONE,VECS(1,IDX(5,N)),-
     $CTMP,VECS(1,IDX(5,I)))
 330     CONTINUE
         DWK(IDX(1,N),9) = SCNRM2(NLEN,VECS(1,IDX(4,N)),1) * DWK(IDX(1,N
     $),7)
         DWK(IDX(1,N),10) = SCNRM2(NLEN,VECS(1,IDX(5,N)),1) * DWK(IDX(1,
     $N),8)
      END IF
C
C     Step 11:
C     Compute A p_n and q_n^T A p_n.
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,IDX(4,N)),VECS(1,IDX(2,NP1)))
C
 340  INFO(2) = 1
      INFO(3) = IDX(4,N)
      INFO(4) = IDX(2,NP1)
      RETLBL  = 350
      RETURN
 350  ZWK(IDX(1,N),5*M+IDX(1,N)) = CDOTU(NLEN,VECS(1,IDX(5,N)),1,VECS(1,
     $IDX(2,NP1)),1) * DWK(IDX(1,N),7) * DWK(IDX(1,N),8)
C
C     Step 10:
C     If regular vectors were built, update the counters.
C
 360  IF (.NOT.INNER) THEN
         K  = K + 1
         MK = IWK(IDX(1,K+1),1)
      END IF
C
C     Update counters.
C
      MPQBLT = MAX0(MPQBLT,MK-IWK(IDX(1,K-1+1),1))
      MKMAX  = MAX0(MK,MKMAX)
      IWK(IDX(1,N),6) = MKSTAR
      IWK(IDX(1,N),4) = MK
C
C     Check for invariant subspaces.
C
      IF (DWK(IDX(1,N),9).LT.TNRM) IERR = IERR + 16
      IF (DWK(IDX(1,N),10).LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 750
C
C     Step 13:
C     Compute A^T q_n.
C     Have the caller carry out ATXB, then return here.
C        CALL ATXB (VECS(1,IDX(5,N)),VECS(1,IDX(3,NP1)))
C
      INFO(2) = 2
      INFO(3) = IDX(5,N)
      INFO(4) = IDX(3,NP1)
      RETLBL  = 370
      RETURN
C
C     Update the norm estimates.
C
 370  IF (N.LE.NUMCHK) THEN
         STMP1 = SCNRM2(NLEN,VECS(1,IDX(2,NP1)),1) * DWK(IDX(1,N),7) / 
     $DWK(IDX(1,N),9)
         STMP2 = SCNRM2(NLEN,VECS(1,IDX(3,NP1)),1) * DWK(IDX(1,N),8) / 
     $DWK(IDX(1,N),10)
         NORMA = AMAX1(STMP1,STMP2,NORMA)
      END IF
      STMP  = CABS(ZWK(IDX(1,N),5*M+IDX(1,N))) / ( DWK(IDX(1,N),9) * 
     $DWK(IDX(1,N),10) )
      NORMA = AMAX1(STMP,NORMA)
      NORM1 = AMAX1(ADJUST*NORMA,NORM1)
      NORM2 = AMAX1(ADJUST*NORMA,NORM2)
C
C     Build v_{n+1} and w_{n+1}.
C     There are three cases:
C       N > NMAX   : the VW sequence is being built new.
C       N >= NLMAX : the VW sequence is being rebuilt, build as if new.
C       N < NLMAX  : the PQ sequence is being rebuilt, just rerun VW.
C     The first two are identical in terms of code.
C
      IWK(IDX(1,N),12)   = K
      IWK(IDX(1,N),13)  = KSTAR
      DWK(IDX(1,N),6) = SZERO
      IBUILT     = .FALSE.
      LBLKSZ     = N - NL + 1
C
C     Clear current column of L.
C
      DO 380 I = 1, M
         ZWK(IDX(1,I),2*M+IDX(1,N)) = CZERO
 380  CONTINUE
C
C     Step 14:
C     Update E^{(n-1)} to E^{(n)}.
C
      DO 400 I = MK, N-1
         CTMP = CZERO
         DO 390 J = MAX0(MKSTAR,IWK(IDX(1,I),4)), N-1
            CTMP = CTMP + ZWK(IDX(1,J),6*M+IDX(1,N)) * ZWK(IDX(1,J),5*M+
     $IDX(1,I)) / DWK(IDX(1,J),3)
 390     CONTINUE
         CTMP     = ZWK(IDX(1,N),M+IDX(1,I)) - CTMP * DWK(IDX(1,N),3)
         ZWK(IDX(1,N),5*M+IDX(1,I)) = CTMP
         ZWK(IDX(1,I),5*M+IDX(1,N)) = CTMP * DWK(IDX(1,I),3) / DWK(IDX(1
     $,N),3)
         STMP1    = CABS(ZWK(IDX(1,N),5*M+IDX(1,I))) / (DWK(IDX(1,N),10)
     $ * DWK(IDX(1,I),9))
         STMP2    = CABS(ZWK(IDX(1,I),5*M+IDX(1,N))) / (DWK(IDX(1,I),10)
     $ * DWK(IDX(1,N),9))
         NORMA    = AMAX1(NORMA,STMP1,STMP2)
         NORM1    = AMAX1(ADJUST*NORMA,NORM1)
         NORM2    = AMAX1(ADJUST*NORMA,NORM2)
 400  CONTINUE
C
C     Step 15:
C     Compute l^\star.
C
      I = LSTAR
      DO 410 J = I+1, L
         IF (IWK(IDX(1,J+1),2).LE.MK) LSTAR = J
 410  CONTINUE
      NLSTAR = IWK(IDX(1,LSTAR+1),2)
C
C     Check memory allocation  (NLSTAR.GE.VWBEG).  Update  the range of
C     valid indices: the new vectors v_{n+1} and w_{n+1} will overwrite
C     the vectors at location N+1-MVEC.
C
      IF (NLSTAR.LT.VWBEG) THEN
         IERR = 64
         GO TO 750
      END IF
      VWBEG = MAX0(VWBEG,NP1-MVEC+1)
C
C     Step 16:
C     Compute F_{1:n,n} = F_{n_{l^\star}:n,n}.
C
      DO 430 I = NLSTAR, N
         CTMP = CZERO
         DO 420 J = MAX0(MK,IWK(IDX(1,I),6)), I
            CTMP = CTMP + ZWK(IDX(1,J),6*M+IDX(1,I)) * ZWK(IDX(1,J),5*M+
     $IDX(1,N)) / DWK(IDX(1,J),3)
 420     CONTINUE
         ZWK(IDX(1,I),M+IDX(1,N)) = CTMP * DWK(IDX(1,I),3)
         STMP     = CABS(CTMP) / DWK(IDX(1,N),9)
         NORMA    = AMAX1(NORMA,STMP)
         NORM1    = AMAX1(ADJUST*NORMA,NORM1)
         NORM2    = AMAX1(ADJUST*NORMA,NORM2)
 430  CONTINUE
C
C     Step 17:
C     Check whether D_l is nonsingular.
C
      IF (LBLKSZ.EQ.1) THEN
         INNER = CABS(ZWK(IDX(1,NL),IDX(1,NL))).EQ.SZERO
      ELSE
         DO 450 I = NL, N
            DO 440 J = NL, N
               ZWK(I-NL+1,4*M+J-NL+1) = ZWK(IDX(1,I),IDX(1,J))
 440        CONTINUE
 450     CONTINUE
         CALL CSVDC (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,8*M+1),ZWK(1,8*M
     $+3),CZERO,0,CZERO,0,ZWK(1,8*M+2),0,I)
         IF (I.NE.0) THEN
            IERR = -I
            GO TO 750
         END IF
         STMP = SZERO
         IF (REAL(ZWK(1,8*M+1)).NE.SZERO) THEN
            STMP = ABS( REAL(ZWK(LBLKSZ,8*M+1)) / REAL(ZWK(1,8*M+1)))
         END IF
         INNER = STMP.LT.(FLOAT(NLEN) * SLAMCH('E'))
      END IF
      LEND = L
      IF (INNER) LEND = L - 1
C
C     Step 18:
C     Compute L_{n_i:n_{i+1}-1,n}, for i = l^\star, ..., lend, lend = l
C     or lend = l-1.
C
      IF (LEND.EQ.L) IWK(IDX(1,L+1+1),2) = NP1
      DO 490 I = LSTAR, LEND
         NI     = IWK(IDX(1,I+1),2)
         NIP1   = IWK(IDX(1,I+1+1),2)
         IBLKSZ = NIP1 - NI
         IF (IBLKSZ.EQ.1) THEN
            ZWK(IDX(1,NI),2*M+IDX(1,N)) = ZWK(IDX(1,NI),M+IDX(1,N)) / 
     $ZWK(IDX(1,NI),IDX(1,NI))
         ELSE
            DO 470 J = NI, NIP1-1
               ZWK(J-NI+1,8*M+1) = ZWK(IDX(1,J),M+IDX(1,N))
               DO 460 IJ = NI, NIP1-1
                  ZWK(J-NI+1,4*M+IJ-NI+1) = ZWK(IDX(1,J),IDX(1,IJ))
 460           CONTINUE
 470        CONTINUE
            CALL CQRDC (ZWK(1,4*M+1),M,IBLKSZ,IBLKSZ,ZWK(1,8*M+3),0,
     $CZERO,0)
            CALL CQRSL (ZWK(1,4*M+1),M,IBLKSZ,IBLKSZ,ZWK(1,8*M+3),ZWK(1,
     $8*M+1),
     $                  CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),
     $CZERO,100,J)
            DO 480 J = NI, NIP1-1
               ZWK(IDX(1,J),2*M+IDX(1,N)) = ZWK(J-NI+1,8*M+2)
 480        CONTINUE
         END IF
 490  CONTINUE
C
C     Step 19:
C     Build the part common to both inner and regular vectors.
C     Assume that V(NP1) = A p_n and W(NP1) = A^T q_n.
C
      ZWK(IDX(1,NP1),2*M+IDX(1,N)) = CZERO
      DO 500 I = NLSTAR, NL-1
         CTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),4) / DWK(IDX(1
     $,N),7)
         CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CONE,VECS(1,IDX(2,NP1)),-
     $CTMP,VECS(1,IDX(2,I)))
         CTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5) / DWK(IDX(1
     $,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CONE,VECS(1,IDX(3,NP1)),-
     $CTMP,VECS(1,IDX(3,I)))
 500  CONTINUE
C
C     Check whether VW is being rerun.
C
      RERUN = (N.LE.NMAX).AND.(N.LT.NLMAX)
      IF (RERUN) THEN
C
C     Check whether D_l has become singular (this should not happen).
C
         IF ((IWK(IDX(1,L+1+1),2).EQ.NP1).AND.INNER) THEN
            IERR = 128
            GO TO 750
         END IF
C
C     Determine whether this was an inner vector or not.
C
         INNER = IWK(IDX(1,L+1+1),2).NE.NP1
         IF (VF.NE.0) THEN
            IF (INNER) THEN
               WRITE (VF,'(A17,I5,A14)') 'Rerunning vector ',NP1,
     $' (VW) as inner'
            ELSE
               WRITE (VF,'(A17,I5,A16)') 'Rerunning vector ',NP1,
     $' (VW) as regular'
            END IF
         END IF
C
C     If building an inner vector, set the coefficients L_{n_l:n,n}.
C
         IF (INNER) THEN
            DO 510 I = NL, N
               ZWK(IDX(1,I),2*M+IDX(1,N))  = CUCPLL(I,N)
 510        CONTINUE
         END IF
      ELSE
C
C     VW is not being rerun, determine whether it is being rebuilt.
C
         IF ((VF.NE.0).AND.(N.LE.NMAX)) THEN
            WRITE (VF,'(A15,I8)') 'Rebuilding V&W:', NP1
         END IF
C
C     If E_k is singular, build an inner vector.
C
         IF (INNER) THEN
            IF (VF.NE.0) WRITE (VF,'(A31)') 
     $'... moment matrix D is singular'
            GO TO 630
         END IF
      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 in
C     Step 20  (H_{n_{l-1}:n,n}) could be done. However, the look-ahead
C     strategy  requires  the smaller of the checks in Step 20 and Step
C     22, in case the norm estimates need updating.  Hence, Step 20 and
C     Step 22  are switched,  and  their  checks  are  later  performed
C     together.
C
C     Save the common part of the vectors.
C
      CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,IDX(2,NP1)),CZERO,VECS(1,3
     $))
      CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,IDX(3,NP1)),CZERO,VECS(1,4
     $))
C
C     Step 21.
C     Build regular vectors.
C
      DO 520 I = NL, N
         CTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),4) / DWK(IDX(1
     $,N),7)
         CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CONE,VECS(1,IDX(2,NP1)),-
     $CTMP,VECS(1,IDX(2,I)))
         CTMP = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5) / DWK(IDX(1
     $,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CONE,VECS(1,IDX(3,NP1)),-
     $CTMP,VECS(1,IDX(3,I)))
 520  CONTINUE
C
C     Compute scale factors for the new vectors.
C
 530  STMP3      = SCNRM2(NLEN,VECS(1,IDX(2,NP1)),1)
      STMP4      = SCNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
      STMP1      = DWK(IDX(1,N),7) * STMP3
      STMP2      = DWK(IDX(1,N),8) * STMP4
      ZWK(IDX(1,NP1),2*M+IDX(1,N)) = CMPLX(STMP1,SZERO)
      IF (STMP1.LT.TNRM) IERR = IERR + 16
      IF (STMP2.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 670
      DWK(IDX(1,NP1),3)     = DWK(IDX(1,N),3) * STMP1 / STMP2
      ZWK(IDX(1,NP1),IDX(1,NP1)) = CDOTU(NLEN,VECS(1,IDX(3,NP1)),1,
     $VECS(1,IDX(2,NP1)),1) / ( STMP3 * STMP4 )
      IF ((STMP3.GE.TMAX).OR.(STMP3.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP3,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CTMP,VECS(1,IDX(2,NP1)),
     $CZERO,VECS(1,IDX(2,NP1)))
         STMP3 = SONE
      END IF
      IF ((STMP4.GE.TMAX).OR.(STMP4.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP4,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CTMP,VECS(1,IDX(3,NP1)),
     $CZERO,VECS(1,IDX(3,NP1)))
         STMP4 = SONE
      END IF
      DWK(IDX(1,NP1),4) = SONE / STMP3
      DWK(IDX(1,NP1),5) = SONE / STMP4
C
C     Step 22:
C     Compute the last column of D_l^{-1}.
C
      IF (LBLKSZ.EQ.1) THEN
         ZWK(IDX(1,NL),3*M+IDX(1,L)) = CONE / ZWK(IDX(1,NL),IDX(1,NL))
      ELSE
         DO 540 I = 1, LBLKSZ-1
            ZWK(I,8*M+1) = CZERO
 540     CONTINUE
         ZWK(LBLKSZ,8*M+1) = CONE
         CALL CQRSL (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,8*M+3),ZWK(1,8*M
     $+1),
     $               CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),CZERO,
     $100,J)
         DO 550 I = NL, N
            ZWK(IDX(1,I),3*M+IDX(1,L)) = ZWK(I-NL+1,8*M+2)
 550     CONTINUE
      END IF
C
C     If VW is being rerun, then skip to the end.
C
      IF (RERUN) GO TO 670
C
C     Step 20:
C     Build the full column H_{n_{l-1}:n,n} and compute its norm.
C
      STMP1 = SZERO
      STMP2 = SZERO
      DO 570 I = IWK(IDX(1,L-1+1),2), N
         CTMP = CZERO
         DO 560 J = MAX0(1,I-1), N
            CTMP = CTMP + ZWK(IDX(1,I),2*M+IDX(1,J)) * ZWK(IDX(1,J),6*M+
     $IDX(1,N))
 560     CONTINUE
         STMP  = CABS(CTMP)
         STMP1 = STMP1 + STMP
         STMP2 = STMP2 + STMP / DWK(IDX(1,I),3)
 570  CONTINUE
      STMP2 = STMP2 * DWK(IDX(1,N),3)
C
C     Step 22:
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
      CTMP  = ZWK(IDX(1,NP1),IDX(1,NP1)) * ZWK(IDX(1,NP1),2*M+IDX(1,N)) 
     $* DWK(IDX(1,N),3) / DWK(IDX(1,NP1),3)
      DO 580 I = NL, N
         STMP  = CABS(CTMP * ZWK(IDX(1,I),3*M+IDX(1,L)))
         STMP3 = STMP3 + STMP
         STMP4 = STMP4 + STMP / DWK(IDX(1,I),3)
 580  CONTINUE
      STMP4 = STMP4 * DWK(IDX(1,NP1),3)
C
C     Steps 20 and 22:
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.NORM2
      IF (.NOT.INNER) GO TO 670
      DWK(IDX(1,N),6) = STMP
C
C     If H_{n_{l-1}:n,n} is bad, build inner vectors.
C
      IF (AMAX1(STMP1,STMP2).GT.NORM2) GO TO 630
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 630
C
C     Build the inner vectors to compute the 2nd term at the next step.
C     Get the coefficients L_{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 590 I = NL, N
         ZWK(IDX(1,I),2*M+IDX(1,NP1)) = CUCPLL(I,N)
         CTMP       = ZWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),4) / 
     $DWK(IDX(1,N),7)
         CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,3),-CTMP,VECS(1,IDX(2,I
     $)))
         CTMP       = ZWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),5) / 
     $DWK(IDX(1,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,4),CONE,VECS(1,4),-CTMP,VECS(1,IDX(3,I
     $)))
 590  CONTINUE
      ISTMP3       = SCNRM2(NLEN,VECS(1,3),1)
      ISTMP4       = SCNRM2(NLEN,VECS(1,4),1)
      ISTMP1       = DWK(IDX(1,N),7) * ISTMP3
      ISTMP2       = DWK(IDX(1,N),8) * ISTMP4
      ZWK(IDX(1,NP1),2*M+IDX(1,NP1)) = CMPLX(ISTMP1,SZERO)
      IF (ISTMP1.LT.TNRM) IERR = IERR + 16
      IF (ISTMP2.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GO TO 640
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
      CTMP = CZERO
      DO 600 J = NLSTAR, NL-1
         CTMP = CTMP + ZWK(IDX(1,NL),IDX(1,J)) * ZWK(IDX(1,I),2*M+IDX(1,
     $N))
 600  CONTINUE
      DO 610 J = NL, N
         CTMP = CTMP + ZWK(IDX(1,NL),IDX(1,J)) * ZWK(IDX(1,I),2*M+IDX(1,
     $NP1))
 610  CONTINUE
      STMP1 = SZERO
      STMP2 = SZERO
      STMP3 = AMAX1(STMP3,STMP4)
      CTMP  = ( ZWK(IDX(1,NL),M+IDX(1,N)) - CTMP ) / ZWK(IDX(1,NP1),2*M+
     $IDX(1,NP1))
      CTMP  = CTMP * ZWK(IDX(1,NL),2*M+IDX(1,NL-1)) * DWK(IDX(1,NL-1),3)
     $ / DWK(IDX(1,NL),3)
      DO 620 I = IWK(IDX(1,L-1+1),2), NL-1
         STMP  = CABS(ZWK(IDX(1,I),3*M+IDX(1,L-1)) * CTMP)
         STMP1 = STMP1 + STMP
         STMP2 = STMP2 + STMP / DWK(IDX(1,I),3)
 620  CONTINUE
      STMP2 = STMP2 * DWK(IDX(1,N),3)
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 670
C
C     Build inner vectors.
C     Check whether the V&W block has to be forced to close.
C
 630  IF (VF.NE.0) WRITE (VF,'(A7,I5,A14)') 'Vector ',NP1,
     $' (VW) is inner'
      IF (NP1-NL.EQ.MAXVW) THEN
         CALL CUCPL2 (NDIM,NLEN,M,N,K,KSTAR,L,LSTAR,MK,MKSTAR,NL,NLSTAR,
     $                VF,IERR,ADJUST,NORM1,NORM2,ZWK,DWK,IWK,IDX,VECS)
         IF (IERR.NE.0) GO TO 750
         LBLKSZ = N - NL + 1
         INNER  = .FALSE.
         RERUN  = .TRUE.
         NP1    = N + 1
         GO TO 530
      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
 640  CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CONE,VECS(1,3),CZERO,VECS(1,
     $IDX(2,NP1)))
      CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CONE,VECS(1,4),CZERO,VECS(1,
     $IDX(3,NP1)))
C
C     Step 24:
C     Get the coefficients L_{n_l:n+1,n} and build inner vectors.
C
      IF (IBUILT) THEN
         DO 650 I = NL, NP1
            ZWK(IDX(1,I),2*M+IDX(1,N)) = ZWK(IDX(1,I),2*M+IDX(1,NP1))
 650     CONTINUE
         STMP1 = ISTMP1
         STMP2 = ISTMP2
         STMP3 = ISTMP3
         STMP4 = ISTMP4
      ELSE
         DO 660 I = NL, N
            ZWK(IDX(1,I),2*M+IDX(1,N))  = CUCPLL(I,N)
            CTMP      = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),4) / 
     $DWK(IDX(1,N),7)
            CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CONE,VECS(1,IDX(2,NP1))
     $,-CTMP,VECS(1,IDX(2,I)))
            CTMP      = ZWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5) / 
     $DWK(IDX(1,N),8) * DWK(IDX(1,N),3) / DWK(IDX(1,I),3)
            CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CONE,VECS(1,IDX(3,NP1))
     $,-CTMP,VECS(1,IDX(3,I)))
 660     CONTINUE
         STMP3      = SCNRM2(NLEN,VECS(1,IDX(2,NP1)),1)
         STMP4      = SCNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
         STMP1      = DWK(IDX(1,N),7) * STMP3
         STMP2      = DWK(IDX(1,N),8) * STMP4
         ZWK(IDX(1,NP1),2*M+IDX(1,N)) = CMPLX(STMP1,SZERO)
         IF (STMP1.LT.TNRM) IERR = IERR + 16
         IF (STMP2.LT.TNRM) IERR = IERR + 32
      END IF
      IF (IERR.NE.0) GO TO 670
      DWK(IDX(1,NP1),3)     = DWK(IDX(1,N),3) * STMP1 / STMP2
      ZWK(IDX(1,NP1),IDX(1,NP1)) = CDOTU(NLEN,VECS(1,IDX(3,NP1)),1,
     $VECS(1,IDX(2,NP1)),1) / ( STMP3 * STMP4 )
      IF ((STMP3.GE.TMAX).OR.(STMP3.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP3,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CTMP,VECS(1,IDX(2,NP1)),
     $CZERO,VECS(1,IDX(2,NP1)))
         STMP3 = SONE
      END IF
      IF ((STMP4.GE.TMAX).OR.(STMP4.LE.TMIN)) THEN
         CTMP = CMPLX(SONE / STMP4,SZERO)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CTMP,VECS(1,IDX(3,NP1)),
     $CZERO,VECS(1,IDX(3,NP1)))
         STMP4 = SONE
      END IF
      DWK(IDX(1,NP1),4) = SONE / STMP3
      DWK(IDX(1,NP1),5) = SONE / STMP4
C
C     Step 23:
C     If regular vectors were built, update the counters.
C
 670  IF (.NOT.INNER) THEN
         L  = L + 1
         NL = IWK(IDX(1,L+1),2)
      END IF
C
C     Update counters.
C
      MVWBLT = MAX0(MVWBLT,NL-IWK(IDX(1,L-1+1),2))
      NLMAX  = MAX0(NL,NLMAX)
      IWK(IDX(1,N),5) = NLSTAR
      IWK(IDX(1,N),3) = NL
C
C     Update the counter for steps taken.
C
      NMAX = MAX0(NMAX,N)
      IF ((N.LT.MKMAX).OR.(N.LT.NLMAX-1)) GO TO 740
C
C     The QMR code starts here.
C     At this point, (N.GE.MKMAX).AND.(N.GE.NLMAX-1), so that  steps up
C     to MIN(MKMAX,NLMAX-1) are guaranteed not to be rebuilt.  Also, no
C     errors are allowed in IERR, other than  possibly having found one
C     or both invariant subspaces, in which case all remaining iterates
C     are computed.
C
      IEND = MIN0(MKMAX,NLMAX-1)
      IF (IERR.NE.0) IEND = N
 680  IF (NQMR.GT.IEND-1) GO TO 740
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),1) = CUCPLO(NQMR+1)
      MAXOMG      = AMAX1(MAXOMG,SONE/DWK(IDX(1,NQMR+1),1))
C
C     Compute the starting index IBASE for the column of \hat{R}.
C
      IBASE       = MAX0(1,IWK(IDX(1,NQMR),5)-1)
      ZWK(IDX(1,IBASE),8*M+5) = CZERO
C
C     Multiply the new column by the previous omegas.
C
      DO 690 J = IWK(IDX(1,NQMR),5), NQMR+1
         ZWK(IDX(1,J),8*M+5) = DWK(IDX(1,J),1) * ZWK(IDX(1,J),2*M+IDX(1,
     $NQMR))
 690  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),8*M+5),1,ZWK(IDX(1,J),8*M+5),1,DWK(IDX(1,J),2),ZWK(IDX(1,J),8*M+6))
C
      DO 700 J = IBASE+1, NQMR
         CTMP1     = ZWK(IDX(1,J),8*M+5)
         CTMP2     = ZWK(IDX(1,J-1),8*M+5)
         ZWK(IDX(1,J-1),8*M+5) = DWK(IDX(1,J),2) * CTMP2 + ZWK(IDX(1,J),
     $8*M+6) * CTMP1
         ZWK(IDX(1,J),8*M+5)   = DWK(IDX(1,J),2) * CTMP1 - 
     $CONJG(ZWK(IDX(1,J),8*M+6)) * CTMP2
 700  CONTINUE
C
C     Compute the rotation for the last element (this also applies it).
C
      CALL CROTG (ZWK(IDX(1,NQMR),8*M+5),ZWK(IDX(1,NQMR+1),8*M+5),
     $DWK(IDX(1,NQMR+1),2),ZWK(IDX(1,NQMR+1),8*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),8*M+4) = SZERO
C        CALL ZROT (1,ZWK(IDX(1,NQMR),8*M+4),1,ZWK(IDX(1,NQMR+1),8*M+4),1,DWK(IDX(1,NQMR+1),2),ZWK(IDX(1,NQMR+1),8*M+6))
C
      ZWK(IDX(1,NQMR+1),8*M+4) = -CONJG(ZWK(IDX(1,NQMR+1),8*M+6)) * 
     $ZWK(IDX(1,NQMR),8*M+4)
      ZWK(IDX(1,NQMR),8*M+4)   =  DWK(IDX(1,NQMR+1),2) * ZWK(IDX(1,NQMR)
     $,8*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,IDX(6,NQMR)) is minimized.
C
      CTMP2 = CZERO
      CTMP  = DWK(IDX(1,NQMR),7)
      DO 710 J = IBASE, NQMR-1
         CTMP1 = ZWK(IDX(1,J),8*M+5) * ZWK(IDX(1,J),8*M+7) / CTMP
         IF (CABS(CTMP1).EQ.SZERO) GO TO 710
         CALL CAXPBY (NLEN,VECS(1,IDX(6,NQMR)),CTMP2,VECS(1,IDX(6,NQMR))
     $,-CTMP1,VECS(1,IDX(6,J)))
         CTMP2 = CONE
 710  CONTINUE
      CALL CAXPBY (NLEN,VECS(1,IDX(6,NQMR)),CTMP2,VECS(1,IDX(6,NQMR)),
     $CONE,VECS(1,IDX(4,NQMR)))
      CTMP = CTMP / ZWK(IDX(1,NQMR),8*M+5)
C
C     Compute the new QMR iterate, then scale the search direction.
C
      CTMP1 = CTMP * ZWK(IDX(1,NQMR),8*M+4)
      CALL CAXPBY (NLEN,VECS(1,1),CONE,VECS(1,1),CTMP1,VECS(1,IDX(6,NQMR
     $)))
      ZWK(IDX(1,NQMR),8*M+7) = CTMP
      STMP      = CABS(CTMP)
      IF ((STMP.GE.TMAX).OR.(STMP.LE.TMIN)) THEN
         ZWK(IDX(1,NQMR),8*M+7) = CONE
         CALL CAXPBY (NLEN,VECS(1,IDX(6,NQMR)),CTMP,VECS(1,IDX(6,NQMR)),
     $CZERO,VECS(1,IDX(6,NQMR)))
      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 = SQRT(FLOAT(NQMR+1)) * MAXOMG * CABS(ZWK(IDX(1,NQMR+1),8*M+4
     $)) / R0
      UCHK = UNRM
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.STEN).AND.(N.LT.NLIM)) GO TO 730
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  = 720
      RETURN
 720  CALL CAXPBY (NLEN,VECS(1,3),CONE,VECS(1,2),-CONE,VECS(1,3))
      RESN = SCNRM2(NLEN,VECS(1,3),1) / R0
      UCHK = RESN
C
C     Output the convergence history.
C
 730  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 750
      ELSE IF (IERR.NE.0) THEN
         GO TO 750
      ELSE IF (UNRM.LT.UCHK/SHUN) THEN
         IERR = 4
         GO TO 750
      ELSE IF (NQMR.GE.NLIM) THEN
         IERR = 4
         GO TO 750
      END IF
      GO TO 680
C
C     Update the running counter.
C
 740  IF (IERR.NE.0) GO TO 750
      IERR = 4
      IF (N.GE.NLIM) GO TO 750
      N = N + 1
      GO TO 30
C
C     Done.
C
 750  RETLBL   = 0
      NLIM     = NQMR
      INFO(1)  = IERR
      MAXPQ    = MPQBLT
      MAXVW    = MVWBLT
      NORMS(1) = NORM1
      NORMS(2) = NORM2
C
      RETURN
      END
C
C**********************************************************************
C
      SUBROUTINE CUCPL1 (NDIM,NLEN,M,N,K,KSTAR,L,LSTAR,MK,MKSTAR,NL,
     $                   NLSTAR,VF,IERR,ADJUST,NORM1,NORM2,ZWK,DWK,IWK,
     $                   IDX,VECS)
C
C     Purpose:
C     This subroutine rebuilds the data  for vectors p_n and q_n at the
C     point where a block was  forced to close.  The routine  is called
C     internally by the coupled Lanczos code,  and is  not meant  to be
C     called directly by the user code.
C
C     Parameters:
C        See descriptions in the main routine CUCPL.
C
C     External routines used:
C     real scnrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine cqrdc(x,ldx,n,p,qraux,jpvt,work,job)
C        LINPACK routine, computes the QR factorization of x.
C     subroutine cqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
C        LINPACK routine, applies the QR factorization of x.
C     subroutine caxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     real cdotu(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C
C     Noel M. Nachtigal
C     May 31, 1993
C
C**********************************************************************
C
      INTRINSIC AMAX1, MAX0
      EXTERNAL CAXPBY, CQRDC, CQRSL
C
      INTEGER IERR, K, KSTAR, L, LSTAR, M, MK, MKSTAR, NL, NLSTAR, N
      INTEGER IDX(6,*), IWK(M,13), NDIM, NLEN, VF
      COMPLEX VECS(NDIM,*), ZWK(M,8*M+7)
      REAL ADJUST, DWK(M,11), NORM1, NORM2
C
C     Miscellaneous parameters.
C
      COMPLEX CONE, CZERO
      PARAMETER (CONE = (1.0E0,0.0E0),CZERO = (0.0E0,0.0E0))
      REAL SZERO
      PARAMETER (SZERO = 0.0E0)
C
C     Local variables.
C
      INTEGER I, IJ, J, KBLKSZ, NP1
      COMPLEX CTMP1, CTMP2
      REAL STMP1, STMP2
C
C     Find the index of the vector pair with the smallest pass value.
C
      IERR = 0
      IF (VF.NE.0) WRITE (VF,'(A23)') 'PQ block did not close:'
      J     = MK + 1
      STMP1 = DWK(IDX(1,J),11)
      DO 100 I = MK+2, N
         STMP2 = DWK(IDX(1,I),11)
         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
      NORM2 = AMAX1(NORM1,NORM2)
      IF (VF.NE.0) WRITE (VF,'(A40,I5,2E11.4)')
     $   '... updated norms, restarting from step:', IWK(IDX(1,J),7), 
     $NORM1, NORM2
      IF (IWK(IDX(1,J),7).EQ.N) RETURN
      N      = IWK(IDX(1,J),7)
      L      = IWK(IDX(1,N),9)
      NL     = IWK(IDX(1,L+1),2)
      K      = IWK(IDX(1,N),8)
      MK     = IWK(IDX(1,K+1),1)
      LSTAR  = IWK(IDX(1,N),11)
      NLSTAR = IWK(IDX(1,LSTAR+1),2)
      KSTAR  = IWK(IDX(1,N),10)
      MKSTAR = IWK(IDX(1,KSTAR+1),1)
C
C     Initialize local variables.
C
      ZWK(IDX(1,N),6*M+IDX(1,N))   = CONE
      NP1        = N + 1
      DWK(IDX(1,N),11) = SZERO
      KBLKSZ     = N - MK
C
C     Step 2:
C     Compute k^\star.
C
      I = KSTAR
      DO 110 J = I+1, K
         IF (IWK(IDX(1,J+1),1).LE.NL-1) KSTAR = J
 110  CONTINUE
      MKSTAR = IWK(IDX(1,KSTAR+1),1)
C
C     Compute U_{m_k:n-1,n}.  Save the old coefficients.
C
      IWK(IDX(1,K+1+1),1) = N
      IF (KBLKSZ.EQ.1) THEN
         ZWK(IDX(1,MK),6*M+IDX(1,NP1)) = ZWK(IDX(1,MK),6*M+IDX(1,N))
         ZWK(IDX(1,MK),6*M+IDX(1,N))   = ZWK(IDX(1,N),M+IDX(1,MK)) / 
     $ZWK(IDX(1,MK),5*M+IDX(1,MK)) * DWK(IDX(1,MK),3) / DWK(IDX(1,N),3)
      ELSE
         DO 130 J = MK, N-1
            ZWK(J-MK+1,8*M+1) = ZWK(IDX(1,N),M+IDX(1,J)) * DWK(IDX(1,J),
     $3) / DWK(IDX(1,N),3)
            DO 120 IJ = MK, N-1
               ZWK(J-MK+1,4*M+IJ-MK+1) = ZWK(IDX(1,J),5*M+IDX(1,IJ))
 120        CONTINUE
 130     CONTINUE
         CALL CQRDC (ZWK(1,4*M+1),M,KBLKSZ,KBLKSZ,ZWK(1,8*M+3),0,CZERO,0
     $)
         CALL CQRSL (ZWK(1,4*M+1),M,KBLKSZ,KBLKSZ,ZWK(1,8*M+3),ZWK(1,8*M
     $+1),
     $               CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),CZERO,
     $100,J)
         DO 140 J = MK, N-1
            ZWK(IDX(1,J),6*M+IDX(1,NP1)) = ZWK(IDX(1,J),6*M+IDX(1,N))
            ZWK(IDX(1,J),6*M+IDX(1,N))   = ZWK(J-MK+1,8*M+2)
 140     CONTINUE
      END IF
C
C     Convert inner vectors to regular vectors.
C
      DO 150 I = MK, N-1
         CTMP1 = ZWK(IDX(1,I),6*M+IDX(1,N)) - ZWK(IDX(1,I),6*M+IDX(1,NP1
     $))
         CTMP2 = CTMP1 * DWK(IDX(1,I),7) / DWK(IDX(1,N),7)
         CALL CAXPBY (NLEN,VECS(1,IDX(4,N)),CONE,VECS(1,IDX(4,N)),-CTMP2
     $,VECS(1,IDX(4,I)))
         CTMP2 = CTMP1 * DWK(IDX(1,I),8) / DWK(IDX(1,N),8) * DWK(IDX(1,N
     $),3) / DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,IDX(5,N)),CONE,VECS(1,IDX(5,N)),-CTMP2
     $,VECS(1,IDX(5,I)))
 150  CONTINUE
C
      RETURN
      END
C
C**********************************************************************
C
      SUBROUTINE CUCPL2 (NDIM,NLEN,M,N,K,KSTAR,L,LSTAR,MK,MKSTAR,NL,
     $                   NLSTAR,VF,IERR,ADJUST,NORM1,NORM2,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 coupled Lanczos code,  and is not meant
C     to be called directly by the user code.
C
C     Parameters:
C        See descriptions in the main routine CUCPL.
C
C     External routines used:
C     subroutine caxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     subroutine cqrdc(x,ldx,n,p,qraux,jpvt,work,job)
C        LINPACK routine, computes the QR factorization of x.
C     subroutine cqrsl(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
      INTRINSIC ABS, AMAX1, MAX0
      EXTERNAL CAXPBY, CQRDC, CQRSL
C
      INTEGER IERR, K, KSTAR, L, LSTAR, M, MK, MKSTAR, NL, NLSTAR, N
      INTEGER IDX(6,*), IWK(M,13), NDIM, NLEN, VF
      COMPLEX VECS(NDIM,*), ZWK(M,8*M+7)
      REAL ADJUST, DWK(M,11), NORM1, NORM2
C
C     Miscellaneous parameters.
C
      COMPLEX CONE, CZERO
      PARAMETER (CONE = (1.0E0,0.0E0),CZERO = (0.0E0,0.0E0))
      REAL SZERO
      PARAMETER (SZERO = 0.0E0)
C
C     Local variables.
C
      INTEGER I, IJ, J, LBLKSZ, NP1
      COMPLEX SCALV, SCALW, CTMP1, CTMP2
      REAL STMP1, STMP2
C
C     Find the index of the vector pair with the smallest pass value.
C
      IERR = 0
      IF (VF.NE.0) WRITE (VF,'(A23)') 'VW block did not close:'
      J     = NL
      STMP1 = DWK(IDX(1,J),6)
      DO 100 I = NL+1, N
         STMP2 = DWK(IDX(1,I),6)
         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
      NORM2 = ADJUST * STMP1
      NORM1 = AMAX1(NORM1,NORM2)
      IF (VF.NE.0) WRITE (VF,'(A40,I5,2E11.4)')
     $      '... updated norms, restarting from step:', IWK(IDX(1,J),7),
     $ NORM1, NORM2
      IF (IWK(IDX(1,J),7).EQ.N) RETURN
      N      = IWK(IDX(1,J),7)
      L      = IWK(IDX(1,N),9)
      NL     = IWK(IDX(1,L+1),2)
      K      = IWK(IDX(1,N),12)
      MK     = IWK(IDX(1,K+1),1)
      LSTAR  = IWK(IDX(1,N),11)
      NLSTAR = IWK(IDX(1,LSTAR+1),2)
      KSTAR  = IWK(IDX(1,N),13)
      MKSTAR = IWK(IDX(1,KSTAR+1),1)
C
C     Initialize local variables.
C
      NP1        = N + 1
      DWK(IDX(1,N),6) = SZERO
      LBLKSZ     = N - NL + 1
C
C     Step 15:
C     Compute l^\star.
C
      I = LSTAR
      DO 110 J = I+1, L
         IF (IWK(IDX(1,J+1),2).LE.MK) LSTAR = J
 110  CONTINUE
      NLSTAR = IWK(IDX(1,LSTAR+1),2)
C
C     Compute L_{n_l:n,n}.  Save the old coefficients.
C
      IWK(IDX(1,L+1+1),2) = 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 130 J = NL, N
            ZWK(J-NL+1,8*M+1) = ZWK(IDX(1,J),M+IDX(1,N))
            DO 120 IJ = NL, N
               ZWK(J-NL+1,4*M+IJ-NL+1) = ZWK(IDX(1,J),IDX(1,IJ))
 120        CONTINUE
 130     CONTINUE
         CALL CQRDC (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,8*M+3),0,CZERO,0
     $)
         CALL CQRSL (ZWK(1,4*M+1),M,LBLKSZ,LBLKSZ,ZWK(1,8*M+3),ZWK(1,8*M
     $+1),
     $               CZERO,ZWK(1,8*M+1),ZWK(1,8*M+2),ZWK(1,8*M+1),CZERO,
     $100,J)
         DO 140 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,8*M+2)
 140     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),4)
      SCALW = ZWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),5) * 
     $DWK(IDX(1,N),3) / DWK(IDX(1,NP1),3)
      ZWK(IDX(1,NP1),2*M+IDX(1,N)) = CZERO
      DO 150 I = NL, N
         CTMP1 = ZWK(IDX(1,I),2*M+IDX(1,N)) - ZWK(IDX(1,I),2*M+IDX(1,NP1
     $))
         CTMP2 = CTMP1 * DWK(IDX(1,I),4) / SCALV
         CALL CAXPBY (NLEN,VECS(1,IDX(2,NP1)),CONE,VECS(1,IDX(2,NP1)),-
     $CTMP2,VECS(1,IDX(2,I)))
         CTMP2 = CTMP1 * DWK(IDX(1,I),5) / SCALW * DWK(IDX(1,N),3) / 
     $DWK(IDX(1,I),3)
         CALL CAXPBY (NLEN,VECS(1,IDX(3,NP1)),CONE,VECS(1,IDX(3,NP1)),-
     $CTMP2,VECS(1,IDX(3,I)))
 150  CONTINUE
C
      RETURN
      END
C
C**********************************************************************
C
      COMPLEX FUNCTION CUCPLL(I,J)
C
C     Purpose:
C     Return the recurrence coefficients for inner vectors VW.
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     April 17, 1993
C
C**********************************************************************
C
C
C     Common block CUCPLX.
C
      REAL NORMA
      COMMON /CUCPLX/NORMA
C
      INTRINSIC CMPLX
C
      INTEGER I, J
C
      IF ((I.LT.1).OR.(J.LT.1)) THEN
         CUCPLL = (0.0E0,0.0E0)
      ELSE IF (I.EQ.J) THEN
         CUCPLL = CMPLX(NORMA,0.0E0)
      ELSE
         CUCPLL = (0.0E0,0.0E0)
      END IF
C
      RETURN
      END
C
C**********************************************************************
C
      REAL FUNCTION CUCPLO (I)
C
C     Purpose:
C     Return the scaling parameter OMEGA(I).
C
C     Parameters:
C     I = the index of the parameter OMEGA (input).
C
C     Noel M. Nachtigal
C     June 1, 1992
C
C**********************************************************************
C
      INTEGER I
C
      CUCPLO = 1.0E0
C
      RETURN
      END
C
C**********************************************************************
C
      COMPLEX FUNCTION CUCPLU(I,J)
C
C     Purpose:
C     Return the recurrence coefficients for inner vectors PQ.
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     April 17, 1993
C
C**********************************************************************
C
      INTEGER I, J
C
      IF ((I.LT.1).OR.(J.LT.1)) THEN
         CUCPLU = (0.0E0,0.0E0)
      ELSE
         CUCPLU = (0.0E0,0.0E0)
      END IF
C
      RETURN
      END
C
C**********************************************************************