C**********************************************************************
C
C     Copyright (C) 1991-1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file "cpyrit.doc" in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
C     This  file  contains  the  routines  for  the  QMR algorithm  for
C     symmetric  matrices, using  the three-term recurrence  variant of
C     the look-ahead Lanczos algorithm.
C
C**********************************************************************
C
      SUBROUTINE DSQMR (NDIM,NLEN,NLIM,MAXVW,M,NORM,DWK,IDX,IWK,
     $                  VECS,TOL,INFO)
C
C     Purpose:
C     This subroutine  uses the QMR algorithm based on  the three-term
C     variant  of the  look-ahead  Lanczos  process  to  solve  linear
C     systems.  It runs  the QMR algorithm  to convergence or  until a
C     user-specified iteration limit is reached. It is set up to solve
C     symmetric systems starting with identical starting vectors.
C
C     The  code is  set up  to solve  the system  A x = b  with initial
C     guess x_0 = 0.  Here  A x = b  denotes the preconditioned system,
C     and  it  is  connected with the  original system as follows.  Let
C     B y = c be the original unpreconditioned system to be solved, and
C     let y_0 be an arbitrary initial guess for its solution.  Then:
C          A x = b, where  A = M_1^{-1} B M_2^{-1},
C                          x = M_2 (y - y_0), b = M_1^{-1} (c - B y_0).
C     Here M = M_1 M_2 is the preconditioner.
C
C     To recover the final iterate y_n for  the original system B y = c
C     from the final iterate x_n for the preconditioned system A x = b,
C     set
C               y_n = y_0 + M_2^{-1} x_n.
C
C     The algorithm  was first described  in the RIACS Technical Report
C     90.46, `An Implementation of the Look-Ahead Lanczos Algorithm for
C     Non-Hermitian   Matrices,  Part  II`,  by  R.W. Freund  and  N.M.
C     Nachtigal, November 1990.
C
C     Parameters:
C     For a description of  the parameters, see the file "dsqmr.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 ddot(n,dx,incx,dy,incy)
C        BLAS-1 routine, computes y^H * x.
C     double precision dlamch(ch)
C        LAPACK routine, computes machine-related constants.
C     double precision dnrm2(n,dx,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine drandn(n,x,seed)
C        Library routine, fills x with random numbers.
C     subroutine drotg(da,db,dcos,dsin)
C        BLAS-1 routine, computes the Givens rotation which rotates the
C        vector [a; b] into [ sqrt(a**2 + b**2); 0 ].
C     subroutine dsqmr1(ndim,nlen,m,maxvw,numchk,n,l,lstar,nl,nlstar,vf,
C        ierr,adjust,norm1,norma,tmax,tnrm,dwk,iwk,idx,vecs)
C        Low-level routine, builds the vectors v_{n+1} and w_{n+1}.
C     double precision dsqmro(n)
C        User-supplied routine, specifies the QMR scaling factors.
C
C     Noel M. Nachtigal
C     October 24, 1990
C
C**********************************************************************
C
      INTRINSIC DABS, DBLE, DMAX1, DSQRT, MAX0, MIN0, MOD
      EXTERNAL DAXPBY, DDOT, DLAMCH, DNRM2, DRANDN, DROTG, DSQMR1
      EXTERNAL DSQMRO
      DOUBLE PRECISION DDOT, DLAMCH, DNRM2, DSQMRO
C
      INTEGER INFO(4), M, MAXVW, NDIM, NLEN, NLIM
      INTEGER IDX(3,NLIM+2), IWK(M,4)
      DOUBLE PRECISION DWK(M,5*M+13), NORM, TOL
      DOUBLE PRECISION VECS(NDIM,3*MAXVW+2+3)
C
C     Common block variables.
C
C
C     Common block DSQMRX.
C
      DOUBLE PRECISION NORMA
      COMMON /DSQMRX/NORMA
C
C     Miscellaneous parameters.
C
      DOUBLE PRECISION DHUN, DONE, DTEN, DZERO
      PARAMETER (DHUN = 1.0D2,DONE = 1.0D0,DTEN = 1.0D1,DZERO = 0.0D0)
C
C     Local variables, permanent.
C
      INTEGER IERR, IEND, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NQMR
      SAVE    IERR, IEND, L, LSTAR, MVWBLT, N, NL, NLSTAR, NMAX, NQMR
      INTEGER NUMCHK, RETLBL, TF, TRES, VF
      SAVE    NUMCHK, RETLBL, TF, TRES, VF
      DOUBLE PRECISION ADJUST, MAXOMG, NORM1, R0, RESN, TMAX, TMIN
      SAVE             ADJUST, MAXOMG, NORM1, R0, RESN, TMAX, TMIN
      DOUBLE PRECISION TNRM, UCHK, UNRM
      SAVE             TNRM, UCHK, UNRM
C
C     Local variables, transient.
C
      INTEGER I, IBASE, INIT, ISJ, ISN, J, REVCOM
      DOUBLE PRECISION DTMP, DTMP1, DTMP2
C
C     Initialize some of the permanent variables.
C
      DATA RETLBL /0/
C
C     Check the reverse communication flag to see where to branch.
C        REVCOM   RETLBL      Comment
C           0        0    first call, go to label 10
C           1       40    returning from AXB, go to label 40
C           1      100    returning from AXB, go to label 100
C           2       50    returning from ATXB, go to label 50
C
      REVCOM  = INFO(2)
      INFO(2) = 0
      IF (REVCOM.EQ.0) THEN
         N      = 0
         NQMR   = 0
         MVWBLT = 0
         IF (RETLBL.EQ.0) GO TO 10
      ELSE IF (REVCOM.EQ.1) THEN
         IF (RETLBL.EQ.40) THEN
            GO TO 40
         ELSE IF (RETLBL.EQ.100) THEN
            GO TO 100
         END IF
      ELSE IF (REVCOM.EQ.2) THEN
         IF (RETLBL.EQ.50) GO TO 50
      END IF
      IERR = 1
      GO TO 130
C
C     Check whether the inputs are valid.
C
 10   IERR = 0
      IF (NDIM.LT.1)       IERR = 2
      IF (NLEN.LT.1)       IERR = 2
      IF (NLIM.LT.1)       IERR = 2
      IF (MAXVW.LT.1)      IERR = 2
      IF (NLEN.GT.NDIM)    IERR = 2
      IF (M.LT.2*MAXVW+2)  IERR = 2
      IF (IERR.NE.0) GO TO 130
C
C     Extract from INFO the output units TF and VF, the true residual
C     flag TRES, and the left starting vector flag INIT.
C
      VF   = MAX0(INFO(1),0)
      INIT = VF / 100000
      VF   = VF - INIT * 100000
      TRES = VF / 10000
      VF   = VF - TRES * 10000
      TF   = VF / 100
      VF   = VF - TF * 100
C
C     Extract the norms.
C
      NORMA = DZERO
      NORM1 = DABS(NORM)
C
C     Set the adjustment parameters.
C
      NUMCHK = 25
      ADJUST = DTEN
C
C     Extract and check the various tolerances and norm estimates.
C
      TNRM = DLAMCH('E') * DTEN
      TMIN = DSQRT(DSQRT(DLAMCH('S')))
      TMAX = DONE / TMIN
      IF (TOL.LE.DZERO) TOL = DSQRT(DLAMCH('E'))
C
C     Start the trace messages and convergence history.
C
      IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') 0, DONE, DONE
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') 0, DONE, DONE
C
C     Set up wrapped indices.  The following indices are used:
C        IDX(1,I) = indices used for the work arrays row-dimensioned M;
C        IDX(2,I) = indices used for v_i;
C        IDX(3,I) = indices used for w_i.
C
      DO 20 I = 1, NLIM+2
         IDX(1,I) = MOD(I-1,M) + 1
         IDX(2,I) = 3 + 1 + 0*(MAXVW+2) + MOD(I-1,MAXVW+2)
         IDX(3,I) = IDX(2,I)
 20   CONTINUE
C
C     Set x_0 = 0 and compute the norm of the initial residual.
C
      CALL DAXPBY (NLEN,VECS(1,IDX(2,1)),DONE,VECS(1,2),DZERO,VECS(1,
     $IDX(2,1)))
      CALL DAXPBY (NLEN,VECS(1,1),DZERO,VECS(1,1),DZERO,VECS(1,1))
      R0 = DNRM2(NLEN,VECS(1,IDX(2,1)),1)
      IF ((TOL.GE.DONE).OR.(R0.EQ.DZERO)) GO TO 130
C
C     Check whether the auxiliary vector must be supplied.
C
CSYM  IF (INIT.EQ.0) CALL DRANDN (NLEN,VECS(1,3),1)
CSYM  CALL DAXPBY (NLEN,VECS(1,IDX(3,1)),DONE,VECS(1,3),DZERO,VECS(1,IDX(3,1)))
C
C     Scale the first pair of Lanczos vectors and check for invariant
C     subspaces.
C
      DTMP1 = R0
CSYM  DTMP2 = DNRM2(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 130
      DWK(IDX(1,1),5*M+8)   = DONE
      DWK(IDX(1,1),IDX(1,1)) = DDOT(NLEN,VECS(1,IDX(2,1)),1,VECS(1,IDX(3
     $,1)),1) / ( DTMP1 * DTMP2 )
      IF ((DTMP1.GE.TMAX).OR.(DTMP1.LE.TMIN)) THEN
         DTMP = DONE / DTMP1
         CALL DAXPBY (NLEN,VECS(1,IDX(2,1)),DTMP,VECS(1,IDX(2,1)),DZERO,
     $VECS(1,IDX(2,1)))
         DTMP1 = DONE
      END IF
      IF ((DTMP2.GE.TMAX).OR.(DTMP2.LE.TMIN)) THEN
         DTMP = DONE / DTMP2
CSYM     CALL DAXPBY (NLEN,VECS(1,IDX(3,1)),DTMP,VECS(1,IDX(3,1)),DZERO,VECS(1,IDX(3,1)))
         DTMP2 = DONE
      END IF
      DWK(IDX(1,1),5*M+9) = DONE / DTMP1
      DWK(IDX(1,1),5*M+10) = DONE / DTMP2
C
C     Initialize the counters.
C
      L      = 1
      N      = 1
      NMAX   = 0
      IWK(IDX(1,0+1),1) = 1
      IWK(IDX(1,1+1),1) = 1
      MVWBLT = 1
      NL     = IWK(IDX(1,L+1),1)
C
C     Set up the QMR iteration.
C
      RESN   = DONE
      DWK(IDX(1,1),5*M+12) = DSQMRO(1)
      DWK(IDX(1,1),5*M+4) = DWK(IDX(1,1),5*M+12) * R0
      MAXOMG = DONE / DWK(IDX(1,1),5*M+12)
C
C     This is one step of the three-term Lanczos algorithm.
C     V_n, W_n, D_{n-1}, F_{n-1}, H_{n-1}, and D_{nn} are given.
C
C     Note that the previous step was  not necessarily step N-1,  as it
C     could have been a restart.
C     Except at the first step, the following hold:
C       NL.LE.N, N-NL.LE.MAXVW-1.
C
 30   IWK(IDX(1,N),4) = L
      IWK(IDX(1,N),3) = N
C
C     Build v_{n+1} and w_{n+1}.
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,IDX(2,N)),VECS(1,3))
C
      INFO(2) = 1
      INFO(3) = IDX(2,N)
      INFO(4) = 3
      RETLBL  = 40
      RETURN
C
C     Have the caller carry out ATXB, then return here.
C        CALL ATXB (VECS(1,IDX(3,N)),VECS(1,3))
C
 40   INFO(2) = 2
      INFO(3) = IDX(3,N)
      INFO(4) = 4
      RETLBL  = 50
CSYM  RETURN
 50   IF ((VF.NE.0).AND.(N.LE.NMAX)) WRITE (VF,'(A11,I8)')
     $    'Rebuilding:', N+1
      CALL DSQMR1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR,VF,IERR,
     $            ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,DWK,IWK,IDX,VECS)
      IF ((IERR.NE.0).AND.(IERR.NE.16).AND.(IERR.NE.32)
     $    .AND.(IERR.NE.48)) GO TO 130
      MVWBLT = MAX0(MVWBLT,NL-IWK(IDX(1,L-1+1),1))
      IWK(IDX(1,N),2) = NLSTAR
C
C     Update the counter for steps taken.
C
      NMAX = MAX0(NMAX,N)
      IF (N.LT.NL-1) GO TO 120
C
C     The QMR code starts here.
C     At this point, steps up to NL-1 are guaranteed not to be rebuilt.
C     No errors are allowed in IERR,  other than possibly  having found
C     one  or both  invariant subspaces,  in which  case  all remaining
C     iterates are computed.
C
      IEND = NL-1
      IF (IERR.NE.0) IEND = N
 60   IF (NQMR.GT.IEND-1) GO TO 120
C
C     Update the QMR iteration counter.
C
      NQMR = NQMR + 1
C
C     Get the next scaling factor omega(i) and update MAXOMG.
C
      DWK(IDX(1,NQMR+1),5*M+12) = DSQMRO(NQMR+1)
      MAXOMG      = DMAX1(MAXOMG,DONE/DWK(IDX(1,NQMR+1),5*M+12))
C
C     Compute the starting index IBASE for the column of \hat{R}.
C
      IBASE       = MAX0(1,IWK(IDX(1,NQMR),2)-1)
      DWK(IDX(1,IBASE),5*M+5) = DZERO
C
C     Multiply the new column by the previous omegas.
C
      DO 70 J = IWK(IDX(1,NQMR),2), NQMR+1
         DWK(IDX(1,J),5*M+5) = DWK(IDX(1,J),5*M+12) * DWK(IDX(1,J),2*M+
     $IDX(1,NQMR))
 70   CONTINUE
C
C     Apply the previous rotations.
C     The loop below explicitly implements a call to DROT:
C     CALL DROT (1,DWK(IDX(1,J-1),5*M+5),1,DWK(IDX(1,J),5*M+5),1,DWK(IDX(1,J),5*M+13),DWK(IDX(1,J),5*M+6))
C
      DO 80 J = IBASE+1, NQMR
         DTMP1     = DWK(IDX(1,J),5*M+5)
         DTMP2     = DWK(IDX(1,J-1),5*M+5)
         DWK(IDX(1,J-1),5*M+5) = DWK(IDX(1,J),5*M+13) * DTMP2 + 
     $DWK(IDX(1,J),5*M+6) * DTMP1
         DWK(IDX(1,J),5*M+5)   = DWK(IDX(1,J),5*M+13) * DTMP1 - 
     $DWK(IDX(1,J),5*M+6) * DTMP2
 80   CONTINUE
C
C     Compute the rotation for the last element (this also applies it).
C
      CALL DROTG (DWK(IDX(1,NQMR),5*M+5),DWK(IDX(1,NQMR+1),5*M+5),
     $DWK(IDX(1,NQMR+1),5*M+13),DWK(IDX(1,NQMR+1),5*M+6))
C
C     Apply the new rotation to the right-hand side vector.
C     Could be replaced with:
C     DWK(IDX(1,NQMR+1),5*M+4) = DZERO
C     CALL DROT (1,DWK(IDX(1,NQMR),5*M+4),1,DWK(IDX(1,NQMR+1),5*M+4),1,DWK(IDX(1,NQMR+1),5*M+13),DWK(IDX(1,NQMR+1),5*M+6))
C
      DWK(IDX(1,NQMR+1),5*M+4) = -DWK(IDX(1,NQMR+1),5*M+6) * DWK(IDX(1,
     $NQMR),5*M+4)
      DWK(IDX(1,NQMR),5*M+4)   =  DWK(IDX(1,NQMR+1),5*M+13) * DWK(IDX(1,
     $NQMR),5*M+4)
C
C     Compute the next search direction s_i.
C     This is more complicated than it might have to be because storage
C     for the vectors VECS(1,ISN) is minimized.
C
      DTMP2 = DZERO
      DTMP  = DWK(IDX(1,NQMR),5*M+9)
      ISN   = 3 + 1 + MAXVW+2 + MOD(NQMR-1,2*MAXVW)
      DO 90 J = IBASE, NQMR-1
         DTMP1 = DWK(IDX(1,J),5*M+5) * DWK(IDX(1,J),5*M+7) / DTMP
         IF (DTMP1.EQ.DZERO) GO TO 90
         ISJ   = 3 + 1 + MAXVW+2 + MOD(J-1,2*MAXVW)
         CALL DAXPBY (NLEN,VECS(1,ISN),DTMP2,VECS(1,ISN),-DTMP1,VECS(1,
     $ISJ))
         DTMP2 = DONE
 90   CONTINUE
      CALL DAXPBY (NLEN,VECS(1,ISN),DTMP2,VECS(1,ISN),DONE,VECS(1,IDX(2,
     $NQMR)))
      DTMP = DTMP / DWK(IDX(1,NQMR),5*M+5)
C
C     Compute the new QMR iterate, then scale the search direction.
C
      DTMP1 = DTMP * DWK(IDX(1,NQMR),5*M+4)
      CALL DAXPBY (NLEN,VECS(1,1),DONE,VECS(1,1),DTMP1,VECS(1,ISN))
      DWK(IDX(1,NQMR),5*M+7) = DTMP
      DTMP      = DABS(DTMP)
      IF ((DTMP.GE.TMAX).OR.(DTMP.LE.TMIN)) THEN
         DWK(IDX(1,NQMR),5*M+7) = DONE
         CALL DAXPBY (NLEN,VECS(1,ISN),DTMP,VECS(1,ISN),DZERO,VECS(1,ISN
     $))
      END IF
C
C     Compute the residual norm upper bound.
C     If the scaled upper bound is within one order of magnitude of the
C     target convergence norm, compute the true residual norm.
C
      UNRM = DSQRT(DBLE(NQMR+1)) * MAXOMG * DABS(DWK(IDX(1,NQMR+1),5*M+4
     $)) / R0
      UCHK = UNRM
      IF ((TRES.EQ.0).AND.(UNRM/TOL.GT.DTEN).AND.(N.LT.NLIM)) GO TO 110
C
C     Have the caller carry out AXB, then return here.
C        CALL AXB (VECS(1,1),VECS(1,3))
C
      INFO(2) = 1
      INFO(3) = 1
      INFO(4) = 3
      RETLBL  = 100
      RETURN
 100  CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,2),-DONE,VECS(1,3))
      RESN = DNRM2(NLEN,VECS(1,3),1) / R0
      UCHK = RESN
C
C     Output the convergence history.
C
 110  IF (VF.NE.0) WRITE (VF,'(I8,2E11.4)') NQMR, UNRM, RESN
      IF (TF.NE.0) WRITE (TF,'(I8,2E11.4)') NQMR, UNRM, RESN
C
C     Check for convergence or termination.  Stop if:
C         1. algorithm converged;
C         2. there is an error condition;
C         3. the residual norm upper bound is smaller than the computed
C     residual norm by a factor of at least 100;
C         4. algorithm exceeded the iterations limit.
C
      IF (RESN.LE.TOL) THEN
         IERR = 0
         GO TO 130
      ELSE IF (IERR.NE.0) THEN
         GO TO 130
      ELSE IF (UNRM.LT.UCHK/DHUN) THEN
         IERR = 4
         GO TO 130
      ELSE IF (NQMR.GE.NLIM) THEN
         IERR = 4
         GO TO 130
      END IF
      GO TO 60
C
C     Update the running counter.
C
 120  IF (IERR.NE.0) GO TO 130
      IERR = 4
      IF (N.GE.NLIM) GO TO 130
      N = N + 1
      GO TO 30
C
C     Done.
C
 130  RETLBL  = 0
      NLIM    = NQMR
      INFO(1) = IERR
      NORM    = NORM1
      MAXVW   = MVWBLT
C
      RETURN
      END
C
C**********************************************************************
C
      DOUBLE PRECISION FUNCTION DSQMRH(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 DSQMRX.
C
      DOUBLE PRECISION NORMA
      COMMON /DSQMRX/NORMA
C
C
      INTEGER I, J
C
      IF ((I.LT.1).OR.(J.LT.1)) THEN
         DSQMRH = 0.0D0
      ELSE IF (I.EQ.J) THEN
         DSQMRH = NORMA
      ELSE
         DSQMRH = 0.0D0
      END IF
C
      RETURN
      END
C
C**********************************************************************
C
      DOUBLE PRECISION FUNCTION DSQMRO (I)
C
C     Purpose:
C     Returns the scaling parameter OMEGA(I).
C
C     Parameters:
C     I = the index of the parameter OMEGA (input).
C
C     Noel M. Nachtigal
C     October 7, 1990
C
C**********************************************************************
C
      INTEGER I
C
      DSQMRO = 1.0D0
C
      RETURN
      END
C
C**********************************************************************
C
      SUBROUTINE DSQMR1 (NDIM,NLEN,M,MAXVW,NUMCHK,N,L,LSTAR,NL,NLSTAR,
     $                  VF,IERR,ADJUST,NORM1,NORMA,TMAX,TMIN,TNRM,DWK,
     $                  IWK,IDX,VECS)
C
C     Purpose:
C     This subroutine builds a new pair of vectors v_{n+1} and w_{n+1}.
C     It is called internally by  the Lanczos code, and is not meant to
C     to be called directly by the user code.
C
C     Parameters:
C        See descriptions in "dulal.doc" or "dsqmr.doc".
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 ddot(n,x,incx,y,incy)
C        BLAS-1 routine, computes y^H * x.
C     double precision dlamch(ch)
C        LAPACK routine, computes machine-related constants.
C     double precision dnrm2(n,x,incx)
C        BLAS-1 routine, computes the 2-norm of x.
C     subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job)
C        LINPACK routine, computes the QR factorization of x.
C     subroutine dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
C        LINPACK routine, applies the QR factorization of x.
C     subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
C        LINPACK routine, computes the SVD of x.
C     subroutine dsqmr2 (ndim,nlen,m,dwk,iwk,idx,vecs)
C        Forces closure of an inner block.
C     double precision dsqmrh(i,n)
C        User-supplied routine, computes inner recurrence coefficients.
C
C     Noel M. Nachtigal
C     July 9, 1993
C
C**********************************************************************
C
      INTRINSIC DABS, DMAX1, MAX0
      EXTERNAL DAXPBY, DDOT, DLAMCH, DNRM2, DQRDC, DQRSL, DSVDC
      EXTERNAL DSQMR2, DSQMRH
      DOUBLE PRECISION DDOT, DLAMCH, DNRM2, DSQMRH
C
      INTEGER IERR, L, LSTAR, M, MAXVW, N, NDIM, NL, NLEN, NLSTAR
      INTEGER IDX(3,3), IWK(M,4), NUMCHK, VF
      DOUBLE PRECISION ADJUST, NORM1, NORMA, TMAX, TMIN, TNRM
      DOUBLE PRECISION DWK(M,5*M+13), VECS(NDIM,3*MAXVW+2+3)
C
C     Miscellaneous parameters.
C
      DOUBLE PRECISION DONE, DZERO
      PARAMETER (DONE = 1.0D0,DZERO = 0.0D0)
C
C     Local variables.
C
      INTEGER I, IJ, J, LBLKSZ, NP1
      DOUBLE PRECISION DTMP, DTMP1, DTMP2, DTMP3, DTMP4, DTMP5
      DOUBLE PRECISION IDTMP1, IDTMP2, IDTMP3, IDTMP4
      LOGICAL IBUILT, INNER, RERUN
C
C     Initialize local variables.
C
      NP1        = N + 1
      DWK(IDX(1,N),5*M+11) = DZERO
      RERUN      = .FALSE.
      IBUILT     = .FALSE.
      IERR       = 0
      LBLKSZ     = N - NL + 1
C
C     Update the norm estimates.
C
      IF (N.LE.NUMCHK) THEN
         DTMP1 = DNRM2(NLEN,VECS(1,3),1) * DWK(IDX(1,N),5*M+9)
CSYM     DTMP2 = DNRM2(NLEN,VECS(1,3),1) * DWK(IDX(1,N),5*M+10)
         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
         DWK(IDX(1,I),2*M+IDX(1,N)) = DZERO
 100  CONTINUE
C
C     Update D^{(n-1)} to D^{(n)}.
C
      DO 120 I = NL, N-1
         DTMP = DZERO
         DO 110 J = NL, N-1
            DTMP = DTMP + DWK(IDX(1,I),IDX(1,J)) * DWK(IDX(1,J),2*M+
     $IDX(1,N-1))
 110     CONTINUE
         DTMP     = ( DWK(IDX(1,I),M+IDX(1,N-1)) - DTMP ) / DWK(IDX(1,N)
     $,2*M+IDX(1,N-1))
         DWK(IDX(1,I),IDX(1,N)) = DTMP
         DWK(IDX(1,N),IDX(1,I)) = DTMP * DWK(IDX(1,N),5*M+8) / DWK(IDX(1
     $,I),5*M+8)
 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
         DTMP = DZERO
         DO 130 J = NL, I+1
            DTMP = DTMP + DWK(IDX(1,N),IDX(1,J)) * DWK(IDX(1,J),2*M+
     $IDX(1,I))
 130     CONTINUE
         DWK(IDX(1,N),M+IDX(1,I)) = DTMP
         DWK(IDX(1,I),M+IDX(1,N)) = DTMP * DWK(IDX(1,I),5*M+8) / 
     $DWK(IDX(1,N),5*M+8)
         NORMA    = DMAX1(NORMA,DABS(DWK(IDX(1,N),M+IDX(1,I))),
     $DABS(DWK(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 = DABS(DWK(IDX(1,NL),IDX(1,NL))).EQ.DZERO
      ELSE
         DO 160 I = NL, N
            DO 150 J = NL, N
               DWK(I-NL+1,4*M+J-NL+1) = DWK(IDX(1,I),IDX(1,J))
 150        CONTINUE
 160     CONTINUE
         CALL DSVDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+1),DWK(1,5*M
     $+3),DZERO,0,DZERO,0,DWK(1,5*M+2),0,I)
         IF (I.NE.0) THEN
            IERR = -I
            GOTO 370
         END IF
         DTMP = DZERO
         IF (DWK(1,5*M+1).NE.DZERO) DTMP = DWK(LBLKSZ,5*M+1) / DWK(1,5*M
     $+1)
         INNER = DTMP.LT.(DBLE(NLEN) * DLAMCH('E'))
      END IF
      IF ((VF.NE.0).AND.INNER)
     $    WRITE (VF,'(A31)') '... moment matrix D is singular'
C
C     Compute H_{n_{l-1}:n_l-1,n}, for n_l > 1.
C     Also, build the part common to both inner and regular vectors.
C     Assume that VECS(1,3) = A v_n and VECS(1,3) = A^T w_n.
C
      DWK(IDX(1,NP1),2*M+IDX(1,N)) = DZERO
      IF (NL.EQ.1) GO TO 180
      IF (NL.EQ.NLSTAR+1) THEN
         DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) = DWK(IDX(1,NLSTAR),M+IDX(1,N))
     $ / DWK(IDX(1,NLSTAR),IDX(1,NLSTAR))
         DTMP = DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),5*M+
     $9) / DWK(IDX(1,N),5*M+9)
         CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP,VECS(1,IDX(2,
     $NLSTAR)))
         DTMP = DWK(IDX(1,NLSTAR),2*M+IDX(1,N)) * DWK(IDX(1,NLSTAR),5*M+
     $10) / DWK(IDX(1,N),5*M+10) * DWK(IDX(1,N),5*M+8) / DWK(IDX(1,
     $NLSTAR),5*M+8)
CSYM     CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP,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
         DTMP2 = DZERO
         DO 170 I = NL-1, NLSTAR, -1
            DWK(IDX(1,I),2*M+IDX(1,N)) = DWK(IDX(1,NL-1),M+IDX(1,N)) * 
     $DWK(IDX(1,I),3*M+IDX(1,L-1))
            IF (N.NE.NL) GO TO 170
            DTMP1 = DWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),5*M+9)
            CALL DAXPBY (NLEN,VECS(1,IDX(2,NL-1)),DTMP2,VECS(1,IDX(2,NL-
     $1)),DTMP1,VECS(1,IDX(2,I)))
            DTMP1 = DWK(IDX(1,I),3*M+IDX(1,L-1)) * DWK(IDX(1,I),5*M+10) 
     $/ DWK(IDX(1,I),5*M+8)
CSYM        CALL DAXPBY (NLEN,VECS(1,IDX(3,NL-1)),DTMP2,VECS(1,IDX(3,NL-1)),DTMP1,VECS(1,IDX(3,I)))
            DTMP2 = DONE
 170     CONTINUE
         DTMP1 = DWK(IDX(1,NL-1),M+IDX(1,N)) / DWK(IDX(1,N),5*M+9)
         CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP1,VECS(1,IDX(2,
     $NL-1)))
         DTMP1 = DWK(IDX(1,NL-1),M+IDX(1,N)) * DWK(IDX(1,N),5*M+8) / 
     $DWK(IDX(1,N),5*M+10)
CSYM     CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP1,VECS(1,IDX(3,NL-1)))
      END IF
C
C     Compute F_{nn}.
C
 180  DWK(IDX(1,N),M+IDX(1,N)) = DDOT(NLEN,VECS(1,IDX(3,N)),1,VECS(1,3),
     $1) * DWK(IDX(1,N),5*M+9) * DWK(IDX(1,N),5*M+10)
      NORMA    = DMAX1(NORMA,DABS(DWK(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
         DWK(IDX(1,NL),2*M+IDX(1,N)) = DWK(IDX(1,NL),M+IDX(1,N)) / 
     $DWK(IDX(1,NL),IDX(1,NL))
      ELSE
         DO 200 J = NL, N
            DWK(J-NL+1,5*M+1) = DWK(IDX(1,J),M+IDX(1,N))
            DO 190 IJ = NL, N
               DWK(J-NL+1,4*M+IJ-NL+1) = DWK(IDX(1,J),IDX(1,IJ))
 190        CONTINUE
 200     CONTINUE
         CALL DQRDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),0,DZERO,0
     $)
         CALL DQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M
     $+1),
     $               DZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),DZERO,
     $100,J)
         DO 210 J = NL, N
            DWK(IDX(1,J),2*M+IDX(1,N)) = DWK(J-NL+1,5*M+2)
 210     CONTINUE
      END IF
C
C     Either D_l is nonsingular, or VW is being rerun.  For the latter,
C     just finish building the vectors.  For the former,  the check for
C     H_{n_{l-1}:n,n} could be done.  However, the look-ahead  strategy
C     requires  the smaller of  the norm checks for H_{n_{l-1}:n,n} and
C     H_{n_l:n,n+1}, so both norms are computed first and then checked.
C
C     Build regular vectors.
C
      DTMP = DWK(IDX(1,N),2*M+IDX(1,N))
      CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,3),-DTMP,VECS(1,
     $IDX(2,N)))
CSYM  CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,3),-DTMP,VECS(1,IDX(3,N)))
      DO 220 I = NL, N-1
         DTMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+9) / 
     $DWK(IDX(1,N),5*M+9)
         CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1)),-
     $DTMP,VECS(1,IDX(2,I)))
         DTMP = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+10) / 
     $DWK(IDX(1,N),5*M+10) * DWK(IDX(1,N),5*M+8) / DWK(IDX(1,I),5*M+8)
CSYM     CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)),-DTMP,VECS(1,IDX(3,I)))
 220  CONTINUE
C
C     Compute scale factors for the new vectors.
C
 230  DTMP3      = DNRM2(NLEN,VECS(1,IDX(2,NP1)),1)
CSYM  DTMP4      = DNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
      DTMP4      = DTMP3
      DTMP1      = DWK(IDX(1,N),5*M+9) * DTMP3
      DTMP2      = DWK(IDX(1,N),5*M+10) * DTMP4
      DWK(IDX(1,NP1),2*M+IDX(1,N)) = DTMP1
      IF (DTMP1.LT.TNRM) IERR = IERR + 16
      IF (DTMP2.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GOTO 360
      DWK(IDX(1,NP1),5*M+8)     = DWK(IDX(1,N),5*M+8) * DTMP1 / DTMP2
      DWK(IDX(1,NP1),IDX(1,NP1)) = DDOT(NLEN,VECS(1,IDX(3,NP1)),1,VECS(1
     $,IDX(2,NP1)),1) / ( DTMP3 * DTMP4 )
      IF ((DTMP3.GE.TMAX).OR.(DTMP3.LE.TMIN)) THEN
         DTMP = DONE / DTMP3
         CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DTMP,VECS(1,IDX(2,NP1)),
     $DZERO,VECS(1,IDX(2,NP1)))
         DTMP3 = DONE
      END IF
      IF ((DTMP4.GE.TMAX).OR.(DTMP4.LE.TMIN)) THEN
         DTMP = DONE / DTMP4
CSYM     CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DTMP,VECS(1,IDX(3,NP1)),DZERO,VECS(1,IDX(3,NP1)))
         DTMP4 = DONE
      END IF
      DWK(IDX(1,NP1),5*M+9) = DONE / DTMP3
      DWK(IDX(1,NP1),5*M+10) = DONE / DTMP4
C
C     Compute the last column of D_l^{-1}.
C
      IF (LBLKSZ.EQ.1) THEN
         DWK(IDX(1,NL),3*M+IDX(1,L)) = DONE / DWK(IDX(1,NL),IDX(1,NL))
      ELSE
         DO 240 I = 1, LBLKSZ-1
            DWK(I,5*M+1) = DZERO
 240     CONTINUE
         DWK(LBLKSZ,5*M+1) = DONE
         CALL DQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M
     $+1),
     $               DZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),DZERO,
     $100,J)
         DO 250 I = NL, N
            DWK(IDX(1,I),3*M+IDX(1,L)) = DWK(I-NL+1,5*M+2)
 250     CONTINUE
      END IF
C
C     If VW is being rerun, then skip to the end.
C
      IF (RERUN) GO TO 360
C
C     Compute the norm of H_{n_{l-1}:n,n}.
C
      DTMP1 = DZERO
      DTMP2 = DZERO
      DO 260 I = IWK(IDX(1,L-1+1),1), N
         DTMP  = DABS(DWK(IDX(1,I),2*M+IDX(1,N)))
         DTMP1 = DTMP1 + DTMP
         DTMP2 = DTMP2 + DTMP / DWK(IDX(1,I),5*M+8)
 260  CONTINUE
      DTMP2 = DTMP2 * DWK(IDX(1,N),5*M+8)
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
      DTMP5 = DWK(IDX(1,NP1),IDX(1,NP1)) * DWK(IDX(1,NP1),2*M+IDX(1,N)) 
     $* DWK(IDX(1,N),5*M+8) / DWK(IDX(1,NP1),5*M+8)
      DO 270 I = NL, N
         DTMP  = DABS(DTMP5 * DWK(IDX(1,I),3*M+IDX(1,L)))
         DTMP3 = DTMP3 + DTMP
         DTMP4 = DTMP4 + DTMP / DWK(IDX(1,I),5*M+8)
 270  CONTINUE
      DTMP4 = DTMP4 * DWK(IDX(1,NP1),5*M+8)
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),5*M+11) = DTMP
C
C     If H_{n_{l-1}:n,n} is bad, build inner vectors.
C
      IF (DMAX1(DTMP1,DTMP2).GT.NORM1) GO TO 320
C
C     If H_{n_l:n-1,n+1} is bad, check the inner vectors.
C     This only applies if n_l > 1.
C
      IF (NL.LE.1) GO TO 320
C
C     Build the inner vectors to compute the 2nd term at the next step.
C     Get the coefficients H_{n_l:n+1,n} in a temporary location.
C     Build inner vectors in VECS(1,3) and VECS(1,3).
C
      IBUILT = .TRUE.
      DO 280 I = NL, N
         DWK(IDX(1,I),2*M+IDX(1,NP1)) = DSQMRH(I,N)
         DTMP       = DWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),5*M+9)
     $ / DWK(IDX(1,N),5*M+9)
         CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP,VECS(1,IDX(2,I
     $)))
         DTMP       = DWK(IDX(1,I),2*M+IDX(1,NP1)) * DWK(IDX(1,I),5*M+10
     $) / DWK(IDX(1,N),5*M+10) * DWK(IDX(1,N),5*M+8) / DWK(IDX(1,I),5*M+
     $8)
CSYM     CALL DAXPBY (NLEN,VECS(1,3),DONE,VECS(1,3),-DTMP,VECS(1,IDX(3,I)))
 280  CONTINUE
      IDTMP3       = DNRM2(NLEN,VECS(1,3),1)
CSYM  IDTMP4       = DNRM2(NLEN,VECS(1,3),1)
      IDTMP4       = IDTMP3
      IDTMP1       = DWK(IDX(1,N),5*M+9) * IDTMP3
      IDTMP2       = DWK(IDX(1,N),5*M+10) * IDTMP4
      DWK(IDX(1,NP1),2*M+IDX(1,NP1)) = IDTMP1
      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
      DTMP = DZERO
      DO 290 J = NLSTAR, NL-1
         DTMP = DTMP + DWK(IDX(1,NL),IDX(1,J)) * DWK(IDX(1,I),2*M+IDX(1,
     $N))
 290  CONTINUE
      DO 300 J = NL, N
         DTMP = DTMP + DWK(IDX(1,NL),IDX(1,J)) * DWK(IDX(1,I),2*M+IDX(1,
     $NP1))
 300  CONTINUE
      DTMP1 = DZERO
      DTMP2 = DZERO
      DTMP3 = DMAX1(DTMP3,DTMP4)
      DTMP  = ( DWK(IDX(1,NL),M+IDX(1,N)) - DTMP ) / DWK(IDX(1,NP1),2*M+
     $IDX(1,NP1))
      DTMP5 = DTMP * DWK(IDX(1,NL),2*M+IDX(1,NL-1)) * DWK(IDX(1,NL-1),
     $5*M+8) / DWK(IDX(1,NL),5*M+8)
      DO 310 I = IWK(IDX(1,L-1+1),1), NL-1
         DTMP  = DABS(DWK(IDX(1,I),3*M+IDX(1,L-1)) * DTMP5)
         DTMP1 = DTMP1 + DTMP
         DTMP2 = DTMP2 + DTMP / DWK(IDX(1,I),5*M+8)
 310  CONTINUE
      DTMP2 = DTMP2 * DWK(IDX(1,N),5*M+8)
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 DSQMR2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR,ADJUST,
     $               NORM1,DWK,IWK,IDX,VECS)
         IF (IERR.NE.0) GO TO 370
         LBLKSZ = N - NL + 1
         INNER  = .FALSE.
         RERUN  = .TRUE.
         NP1    = N + 1
         GO TO 230
      END IF
C
C     The temporary vectors  contain either just partial inner vectors,
C     or the  completed ones,  depending on whether  IBUILT is  TRUE or
C     FALSE.  In either case, replace the regular vectors.
C
 330  CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,3),DZERO,VECS(1,
     $IDX(2,NP1)))
CSYM  CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,3),DZERO,VECS(1,IDX(3,NP1)))
C
C     Get the coefficients H_{n_l:n+1,n} and build inner vectors.
C
      IF (IBUILT) THEN
         DO 340 I = NL, NP1
            DWK(IDX(1,I),2*M+IDX(1,N)) = DWK(IDX(1,I),2*M+IDX(1,NP1))
 340     CONTINUE
         DTMP1 = IDTMP1
         DTMP2 = IDTMP2
         DTMP3 = IDTMP3
         DTMP4 = IDTMP4
      ELSE
         DO 350 I = NL, N
            DWK(IDX(1,I),2*M+IDX(1,N)) = DSQMRH(I,N)
            DTMP     = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+9) 
     $/ DWK(IDX(1,N),5*M+9)
            CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1))
     $,-DTMP,VECS(1,IDX(2,I)))
            DTMP     = DWK(IDX(1,I),2*M+IDX(1,N)) * DWK(IDX(1,I),5*M+10)
     $ / DWK(IDX(1,N),5*M+10) * DWK(IDX(1,N),5*M+8) / DWK(IDX(1,I),5*M+8
     $)
CSYM        CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)),-DTMP,VECS(1,IDX(3,I)))
 350     CONTINUE
         DTMP3 = DNRM2(NLEN,VECS(1,IDX(2,NP1)),1)
CSYM     DTMP4 = DNRM2(NLEN,VECS(1,IDX(3,NP1)),1)
         DTMP4 = DTMP3
         DTMP1 = DWK(IDX(1,N),5*M+9) * DTMP3
         DTMP2 = DWK(IDX(1,N),5*M+10) * DTMP4
      END IF
      DWK(IDX(1,NP1),2*M+IDX(1,N)) = DTMP1
      IF (DTMP1.LT.TNRM) IERR = IERR + 16
      IF (DTMP2.LT.TNRM) IERR = IERR + 32
      IF (IERR.NE.0) GOTO 360
      DWK(IDX(1,NP1),5*M+8)     = DWK(IDX(1,N),5*M+8) * DTMP1 / DTMP2
      DWK(IDX(1,NP1),IDX(1,NP1)) = DDOT(NLEN,VECS(1,IDX(3,NP1)),1,VECS(1
     $,IDX(2,NP1)),1) / ( DTMP3 * DTMP4 )
      IF ((DTMP3.GE.TMAX).OR.(DTMP3.LE.TMIN)) THEN
         DTMP = DONE / DTMP3
         CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DTMP,VECS(1,IDX(2,NP1)),
     $DZERO,VECS(1,IDX(2,NP1)))
         DTMP3 = DONE
      END IF
      IF ((DTMP4.GE.TMAX).OR.(DTMP4.LE.TMIN)) THEN
         DTMP = DONE / DTMP4
CSYM     CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DTMP,VECS(1,IDX(3,NP1)),DZERO,VECS(1,IDX(3,NP1)))
         DTMP4 = DONE
      END IF
      DWK(IDX(1,NP1),5*M+9) = DONE / DTMP3
      DWK(IDX(1,NP1),5*M+10) = 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 DSQMR2 (NDIM,NLEN,M,N,L,LSTAR,NL,NLSTAR,VF,IERR,
     $                  ADJUST,NORM1,DWK,IWK,IDX,VECS)
C
C     Purpose:
C     This subroutine rebuilds the data for vectors v_{n+1} and w_{n+1}
C     at the point where a block  was forced  to close.  The routine is
C     called internally by the look-ahead Lanczos code and is not meant
C     to be called directly by the user code.
C
C     Parameters:
C        See descriptions in "dulal.doc" or "dsqmr.doc".
C
C     External routines used:
C     subroutine daxpby(n,z,a,x,b,y)
C        Library routine, computes z = a * x + b * y.
C     subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job)
C        LINPACK routine, computes the QR factorization of x.
C     subroutine dqrsl(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 DAXPBY, DQRDC, DQRSL
C
      INTEGER IERR, L, LSTAR, M, N, NDIM, NL, NLEN, NLSTAR
      INTEGER IDX(3,3), IWK(M,4), VF
      DOUBLE PRECISION ADJUST, DWK(M,5*M+13), NORM1, VECS(NDIM,*)
C
C     Miscellaneous parameters.
C
      DOUBLE PRECISION DONE, DZERO
      PARAMETER (DONE = 1.0D0,DZERO = 0.0D0)
C
C     Local variables.
C
      INTEGER I, IJ, J, LBLKSZ, NP1
      DOUBLE PRECISION DTMP1, DTMP2, SCALV, SCALW
C
C     Find the index of the vector pair with the smallest pass value.
C
      IERR = 0
      IF (VF.NE.0) WRITE (VF,'(A20)') 'block did not close:'
      J     = NL
      DTMP1 = DWK(IDX(1,J),5*M+11)
      DO 100 I = NL+1, N
         DTMP2 = DWK(IDX(1,I),5*M+11)
         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),5*M+11) = 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
         DWK(IDX(1,NL),2*M+IDX(1,NP1)) = DWK(IDX(1,NL),2*M+IDX(1,N))
         DWK(IDX(1,NL),2*M+IDX(1,N))   = DWK(IDX(1,NL),M+IDX(1,N)) / 
     $DWK(IDX(1,NL),IDX(1,NL))
      ELSE
         DO 120 J = NL, N
            DWK(J-NL+1,5*M+1) = DWK(IDX(1,J),M+IDX(1,N))
            DO 110 IJ = NL, N
               DWK(J-NL+1,4*M+IJ-NL+1) = DWK(IDX(1,J),IDX(1,IJ))
 110        CONTINUE
 120     CONTINUE
         CALL DQRDC (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),0,DZERO,0
     $)
         CALL DQRSL (DWK(1,4*M+1),M,LBLKSZ,LBLKSZ,DWK(1,5*M+3),DWK(1,5*M
     $+1),
     $               DZERO,DWK(1,5*M+1),DWK(1,5*M+2),DWK(1,5*M+1),DZERO,
     $100,J)
         DO 130 J = NL, N
            DWK(IDX(1,J),2*M+IDX(1,NP1)) = DWK(IDX(1,J),2*M+IDX(1,N))
            DWK(IDX(1,J),2*M+IDX(1,N)) = DWK(J-NL+1,5*M+2)
 130     CONTINUE
      END IF
C
C     Convert inner vectors to regular vectors.
C
      SCALV = DWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),5*M+9)
      SCALW = DWK(IDX(1,NP1),2*M+IDX(1,N)) * DWK(IDX(1,NP1),5*M+10) * 
     $DWK(IDX(1,N),5*M+8) / DWK(IDX(1,NP1),5*M+8)
      DWK(IDX(1,NP1),2*M+IDX(1,N)) = DZERO
      DO 140 I = NL, N
         DTMP1 = DWK(IDX(1,I),2*M+IDX(1,N)) - DWK(IDX(1,I),2*M+IDX(1,NP1
     $))
         DTMP2 = DTMP1 * DWK(IDX(1,I),5*M+9) / SCALV
         CALL DAXPBY (NLEN,VECS(1,IDX(2,NP1)),DONE,VECS(1,IDX(2,NP1)),-
     $DTMP2,VECS(1,IDX(2,I)))
         DTMP2 = DTMP1 * DWK(IDX(1,I),5*M+10) / SCALW * DWK(IDX(1,N),5*M
     $+8) / DWK(IDX(1,I),5*M+8)
CSYM     CALL DAXPBY (NLEN,VECS(1,IDX(3,NP1)),DONE,VECS(1,IDX(3,NP1)),-DTMP2,VECS(1,IDX(3,I)))
 140  CONTINUE
C
      RETURN
      END
C
C**********************************************************************