SUBROUTINE BT(N,X,F,G,COMPFG,FM,EPS,MAXCOM,MAXIT,RESET,
     *              IPRINT,IFAIL,IWORK,WORK,LWORK)
C
C     B(undle)T(rust)-Algorithm (trying to minimize a convex function)
C
C     Double Precision version
C
C     N      : dimension of the problem
C     X      : vector, dimension N
C              Input : starting point x0
C              Output: solution point
C     F      : Input : f(x0)
C              Output: f(solution)
C     G      : Input : subgradient g(x0)
C     COMPFG : SUBROUTINE, computes f(x) and a subgradient g(x),
C              Parameter: N,X,F,G
C     FM     : estimate of the minimum of f
C     EPS    : stopping criterion,
C              eucl. norm of alpha-subgradient =< EPS and alpha =< eps
C     MAXCOM : maximal number of calls to COMPFG
C     MAXIT  : maximal number of iterations
C     RESET  : maximal number of stored subgradients, RESET > 3
C     IPRINT : to control printout
C              IPRINT=0: no printout
C              IPRINT=1: final printout
C              IPRINT=2: in iteration k
C                        niter= number of the iteration
C                        ncomp= number of calls to compfg
C                        f    = f(xk)
C                        gn   = norm of an alpha-subgradient at xk
C                      ( gives information about optimality:
C                        f(y) >= f(xk)-gn*norm(y-xk)-alpha for all y )
C              IPRINT>2: debug printout
C     IFAIL  : information about the result of BT
C              IFAIL=0: everything is okay
C              IFAIL=1: number of calls of compfg > MAXCOM
C              IFAIL=2: number of iterations > MAXIT
C              IFAIL=3: 'T too small'
C              IFAIL=4: not enough working space
C              IFAIL=5: failure in quadratic program
C     WORK   : working array of dimension LWORK
C     LWORK  : dimension of WORK,
C              LWORK  >= (RESET+N+9)*RESET+4*N+10
C
C***********************************************
C  Version 1.0 of Bundle Trust, March 2, 1989  *
C                                              *
C  Copyright by:                               *
C                                              *
C  Helga Schramm                               *
C  Mathematisches Institut                     *
C  Universitaet Bayreuth                       *
C  Postfach 101251                             *
C  D-8580 Bayreuth                             *
C***********************************************
C
      INTEGER          N,MAXIT,MAXCOM,IPRINT,IFAIL,LWORK,RESET
      DOUBLE PRECISION F,FM,EPS
      INTEGER          IWORK(RESET)
      DOUBLE PRECISION X(N),G(N),WORK(LWORK)
      EXTERNAL         COMPFG
C
C    local variables
C
      INTEGER          LALPHA,LXN,LA,LHESS,LY,LYA,LD,LDL,LDU,
     *                 LWORK1,LWORK2,LWORKT
C
C
      LALPHA= 1
      LXN   = LALPHA+RESET
      LA    = LXN+N
      LHESS = LA+RESET*N
      LY    = LHESS+RESET*(RESET+1)/2
      LYA   = LY+RESET
      LD    = LYA+RESET
      LDL   = LD+N
      LDU   = LDL+N
      LWORK1= LDU+N
      LWORK2= (RESET+11)*RESET/2+10
      LWORKT= LWORK1+LWORK2-1
      IFAIL = -1
      IF (LWORK .LT. LWORKT) IFAIL= 4
      IF (IPRINT .GT. 0) THEN
         WRITE(*,*)
         WRITE(*,*) 'BT-Algorithm'
         WRITE(*,*) '============'
         WRITE(*,*)
         WRITE(*,*) 'workspace provided is    WORK(',LWORK,')'
         WRITE(*,*) 'to solve problem we need WORK(',LWORKT,')'
         WRITE(*,*)
      ENDIF
      IF (IFAIL .GT. 0) THEN
         IF (IPRINT .GT. 5) THEN
            WRITE(*,*) 'Not enough workspace to solve problem'
            WRITE(*,*) 'IFAIL= ', IFAIL
         ENDIF
         RETURN
      ENDIF
      CALL BT1(N,X,F,G,COMPFG,WORK(LALPHA),FM,EPS,MAXCOM,MAXIT,
     *         RESET,IPRINT,IFAIL,WORK(LXN),WORK(LA),WORK(LHESS),
     *         IWORK,WORK(LWORK1),LWORK2,WORK(LY),WORK(LYA),
     *         WORK(LD),WORK(LDL),WORK(LDU))
      END
      SUBROUTINE BT1(N,X,F,G,COMPFG,ALPHA,FM,EPS,MAXCOM,MAXIT,RESET,
     *              IPRINT,IFAIL,XN,A,HESS,IWORK,WORK,LWORK,Y,YA,
     *              D,DL,DU)
      INTEGER          N,MAXIT,MAXCOM,IPRINT,IFAIL,LWORK,RESET
      DOUBLE PRECISION F,FM,EPS
      INTEGER          IWORK(RESET)
      DOUBLE PRECISION X(N),G(N),ALPHA(RESET),XN(N),A(RESET,N),
     *                 HESS(RESET*(RESET+1)/2),WORK(LWORK),
     *                 Y(RESET),YA(RESET),D(N),DL(N),DU(N)
      EXTERNAL         COMPFG,MLIS4
C
      DOUBLE PRECISION MU,ZERO
      PARAMETER (MU  = 1.0D-1,
     *           ZERO= 1.0D-10)
C
C    local variables
C
      INTEGER          I,J,NITER,NITER1,NN,NCOMP,INFORM,LALT,K,L,MODE,
     *                 LOGIC,IS,LIN,KK,K1,LINF,NAC
      DOUBLE PRECISION T,FN,V,VDA,VSS,ADA,VGXAL,W,VK,GD,TL,TU,SU,
     *                 ALN,ALM,ZM,ZP,ALP,GG,TP,TMIN,FPN,AM,D2,TPS,TNC,
     *                 TOL,NU,VV,TA,TLL,TUL,ALL,ALU
      LOGICAL          KRIT,LLIN
      INTRINSIC        DABS,DSQRT
c
c       Graph Partitioning additions:
c       temporary storage of eigenvalues
        real eigval(5), xbst
c       istop tells BT to abort if we get stuck in Lanczos
        common / stop / istop, mulply
c xbst is current best feasible solution
        common / eigvls / eigval, xbst
c       Only write out the nperm reliable eigenvalues
        common / iconst / matmul, nfig, nblock, nval, nperm
C
C    start
C
c     GPP change: DL should be initialised else may be used before assigned
      DO 338 I=1,N
         DL(I)= 0.0d0
  338 CONTINUE
      KRIT  = .TRUE.
      LLIN  = .FALSE.
      LINF  = 0
      NCOMP = 0
      NITER = 0
      NITER1= 0
      IFAIL = -1
      INFORM= 0
      IS= 0
      NN= 0
      ALPHA(1)= 0.0D0
      DO 1 I= 1,RESET
         YA(I)= -1.0D0
         Y(I) =  1.0D0
    1 CONTINUE
      SU = DMAX1(100.0D0*(F-FM),1.0D0)
      T  = 0.0D0
      TP = 0.0D0
      TLL= 0.0D0
      TUL= 0.0D0
      GG = 0.0D0
      DO 2 I=1,N
         GG= GG+G(I)*G(I)
    2 CONTINUE
      ZP = DSQRT(GG)
      ALP= 0.0D0
      ALM= 0.0D0
      ALN= 0.0D0
      ALL= 0.0D0
      ALU= 0.0D0
      IF (ZP .LE. EPS) THEN
         IFAIL= 0
         IF (IPRINT .GT. 1) WRITE(*,*) 'x already optimal'
         RETURN
      ENDIF
c gpp addition: print eigval and matmul
      IF (IPRINT .GT. 1) WRITE(*,'(3X,A,1X,A,7X,A,10X,A,9X,A)')
     *                        'niter','ncomp','ubd','bfs','gap'
        xgap = f / xbst - 1.0
      IF (IPRINT .GT. 1)
     *   WRITE(*,'(A,I4,A,I4,A,f15.2,A,f15.2,A,f7.3,a,f12.4,a,f12.4)')
     *           ' ',NITER1+1,' ',NCOMP+1,' ',F,' ',xbst,' ',xgap,
     *           ' ', zp,' ',alp

998     format(1x, 5f13.5)

C
C    iteration
C
  100 CONTINUE
      TA= T
      IF (KRIT) THEN
         IF (NITER .GE. RESET .OR. NN .EQ. RESET+100) THEN
C
C    reset
C
            IF (NN .EQ. RESET+100) THEN
               DO 3 I=1,N
                  A(1,I)= -D(I)/T
    3         CONTINUE
               ALPHA(1)= ALP
               DO 5 J=1,NITER
                  IF (ABS(ALPHA(J)) .LT. ZERO) THEN
                     ALPHA(2)= DMAX1(ALPHA(J),0.0D0)
                     DO 4 I=1,N
                        A(2,I)= A(J,I)
    4               CONTINUE
                  ENDIF
    5          CONTINUE
               NN= 2
               NITER= 2
               GOTO 111
            ENDIF
            NN= 0
            DO 7 J=1,NITER
            IF ((Y(J).GE.ZERO) .OR. (ABS(ALPHA(J)) .LT. ZERO)) THEN
               NN= NN+1
               ALPHA(NN)= DMAX1(ALPHA(J),0.0D0)
               DO 6 I=1,N
                  A(NN,I)= A(J,I)
    6          CONTINUE
            ENDIF
    7       CONTINUE
            IF (NN .GE. RESET) THEN
               NN= 2
               ALPHA(1)= ALP
               DO 8 I=1,N
                  A(1,I)= -D(I)/T
    8          CONTINUE
               DO 10 J=1,NITER
                  IF (ABS(ALPHA(J)) .LT. ZERO) THEN
                     ALPHA(2)= DMAX1(ALPHA(J),0.0D0)
                     DO 9 I=1,N
                        A(2,I)= A(J,I)
    9                CONTINUE
                  ENDIF
   10          CONTINUE
            ENDIF
  111       NITER= NN
            DO 11 I=1,RESET
               Y(I)= 0.0D0
               YA(I)= -1.0D0
   11       CONTINUE
            DO 12 I=NITER+1,RESET
               ALPHA(I)= 0.0D0
   12       CONTINUE
C
C    hess update
C
            DO 15 K=1,NITER
               LALT= K*(K-1)/2
               DO 14 J=1,K
                  K1= LALT+J
                  HESS(K1)= 0.0D0
                  DO 13 I=1,N
                     HESS(K1)= HESS(K1)+A(J,I)*A(K,I)
   13             CONTINUE
   14          CONTINUE
   15       CONTINUE
            KK= 1
            Y(1)= 1.0D0
            MODE= 1
            IF (IPRINT .GT. 1) THEN
               WRITE(*,*)
               WRITE(*,*) 'Reset, nn=',NN
               WRITE(*,*)
            ENDIF
         ENDIF
C
C    end of reset
C
         NITER1= NITER1+1
         IF (IFAIL .LT. 0 .AND. NITER1 .GE. MAXIT) THEN
            IFAIL= 2
            IF (IPRINT .GT. 1) WRITE(*,*) 'max. number of iterations'
            MAXCOM= NCOMP
         ENDIF
C
C    data for quad
C
         NITER= NITER+1
         ALPHA(NITER)= ALN
         DO 16 I=1,N
            A(NITER,I)= G(I)
   16    CONTINUE
         LALT= NITER*(NITER-1)/2
         DO 18 J=1,NITER
            K =LALT+J
            HESS(K)= 0.0D0
            DO 17 I=1,N
               HESS(K)= HESS(K)+A(J,I)*A(NITER,I)
   17       CONTINUE
   18    CONTINUE
         IF (NITER .EQ. 1) Y(1)= 1.0D0
         TL = 0.0D0
         TU = SU
         ZM = ZP
         ALM= ALP
      ENDIF
C
C    end krit
C
      IF (NITER .EQ. 1 .AND. KRIT) THEN
C
C    compute first t via line search
C
         DO 19 I=1,N
            D(I) = -G(I)
            XN(I)= X(I)
   19    CONTINUE
         FN= F
         T= 2.0D0*(F-FM)/ZP
         T= DMAX1(T,1.0D0/ZP)
         TMIN= 1.0D-20
         FPN= -GG
         D2= GG
         AM= 0.0D0
         TOL= DMAX1((F-FM),EPS)
         CALL MLIS4 (COMPFG,N,XN,FN,FPN,T,TMIN,SU,D,D2,G,0.2D0,0.1D0,
     *               LOGIC,NCOMP,MAXCOM,X,TOL,AM,TPS,TNC)
         IF (IPRINT. GE. 10) WRITE(*,'(5X,A,I4,3X,E16.8)')
     *                              'line search:',LOGIC,T
         V= -T*GG
         IF (IPRINT. GE. 5) THEN
            WRITE(*,'(48X,A,E16.8)') 't  =',T
            WRITE(*,'(48X,A,E16.8)') 'rho=',0.5*T*T*GG
         ENDIF
         IF (LOGIC .EQ. 3) THEN
C
C    null step
C
            ALN= DMAX1(TPS,0.0D0)
            DO 20 I=1,N
               X(I)= XN(I)
   20        CONTINUE
            KRIT= .TRUE.
            GOTO 100
         ENDIF
         IF (LOGIC .GE. 0 .AND. LOGIC. LE. 2) THEN
C
C    serious step
C
            DO 22 I=1,NITER
               GD= 0.0D0
               DO 21 J=1,N
                  GD= GD+T*D(J)*A(I,J)
   21          CONTINUE
               ALPHA(I)= ALPHA(I)+FN-F-GD
               ALPHA(I)= DMAX1(ALPHA(I),0.0D0)
   22       CONTINUE
            ALN= 0.0D0
            DO 23 I=1,N
               X(I)= XN(I)
   23       CONTINUE
            F= FN
            KRIT= .TRUE.
            GOTO 100
         ENDIF
         T= 0.1D0*(F-FM)/ZP
         T= DMAX1(T,1.0D0/ZP)
         DO 210 I=1,N
            X(I)= XN(I)
  210    CONTINUE
      ENDIF
C
C    compute d(t)
C
      MODE = 2
      LLIN= .FALSE.
      IF (KRIT .OR. NITER .LE. 2) THEN
          MODE = 1
          LIN  = 0
      ENDIF
      IF (MODE .EQ. 2) THEN
          Y(NITER)    = 1.0D0
          IWORK(NITER)= 1
      ENDIF
      IF (LIN .GT. 2 .AND. T .GT. TLL .AND. T .LT. TUL) THEN
C
C    piecewise linear
C
         LLIN= .TRUE.
         NU= (T-TLL)/(TUL-TLL)
         DO 24 I=1,N
            D(I)= DL(I)+NU*(DU(I)-DL(I))
   24    CONTINUE
         ALP= ALL+NU*(ALU-ALL)
         V= -ALPHA(1)+A(1,1)*D(1)
         DO 25 I=2,N
            V= V+A(1,I)*D(I)
   25    CONTINUE
         DO 27 J= 2,NITER
            VV= -ALPHA(J)+A(J,1)*D(1)
            DO 26 I=2,N
               VV= VV+A(J,I)*D(I)
   26       CONTINUE
            V= DMAX1(V,VV)
   27    CONTINUE
      ELSE
C
C    quadratic program
C
         CALL QUAD01(NITER,N,HESS,T,ALPHA,MODE,KK,IWORK,Y,ADA,VK,
     *               VDA,VSS,W,L,VGXAL,WORK,LWORK,INFORM)
C
         V= T*VDA
C
         DO 29 I=1,N
            D(I)= 0.0D0
            DO 28 J=1,NITER
               D(I)= D(I)-T*Y(J)*A(J,I)
   28       CONTINUE
   29    CONTINUE
C
         IF (IPRINT .GT. 10) THEN
            NAC= 0
            WRITE(*,*)
            WRITE(*,'(/A,I2,A,E23.15)') 'inform',INFORM,'   w',W
            WRITE(*,*) ' active subgradients:'
            DO 30 I=1,NITER
               IF (Y(I) .GE. ZERO) THEN
                 WRITE(*,*) I,Y(I),ALPHA(I)
                 NAC= NAC+1
               ENDIF
   30       CONTINUE
            WRITE(*,*) 'NAC=',NAC
         ENDIF
      ENDIF
C
      LIN= LIN+1
      DO 31 I=1,NITER
         IF (YA(I)*Y(I) .LT. ZERO .AND. (DABS(YA(I))+Y(I)) .GT. ZERO)
     *      LIN= 0
         YA(I)= Y(I)
   31 CONTINUE
C
C    next trialpoint
C
      DO 32 I=1,N
         XN(I)= X(I)+D(I)
   32 CONTINUE
C
C    compute zp,alpha
C
      ZP= 0.0D0
      DO 33 I=1,N
         ZP= ZP+D(I)*D(I)
   33 CONTINUE
      ZP= ZP/(T**2)
      IF (.NOT. LLIN) ALP= ADA
C
      IF (IPRINT .GT. 20) WRITE(*,*) 'zp2,alp  :',ZP,ALP
      IF (ZP .LE. -ZERO .OR. ALP .LE. -ZERO
     *    .OR. INFORM .GT. 0 .AND. LINF .LE. 3) THEN
         LINF= LINF+1
         KRIT= .TRUE.
         NN  = RESET+100
         IF (IPRINT .GT. 5) WRITE(*,*) 'linf= ',LINF
         IF (LINF .LE. 3) GOTO 100
      ENDIF
      IF (ZP .LE. -ZERO .OR. ALP .LE. -ZERO
     *   .OR. INFORM .GT. 0 .OR. LINF .GT. 3) THEN
         IF (IPRINT .GT. 1) WRITE(*,*) 'failure in qpdf2s'
         IF (IPRINT .GT. 1) WRITE(*,*) 'inform= ',INFORM
         IFAIL = 5
         MAXIT = NITER1
         MAXCOM= NCOMP
      ENDIF
C
      ZP= DSQRT(DABS(ZP))
C
c GPP addition: write out the eigenvalues
        xgap = f / xbst - 1.0
      IF ((IPRINT .GT. 1 .AND. KRIT) .OR. IPRINT .GE. 5) then
        WRITE(*,'(A,I4,A,I4,A,f15.2,A,f15.2,A,f7.3,a,f12.4,a,f12.4)')
     *           ' ',NITER1,' ',NCOMP+1,' ',F,' ',xbst,' ',xgap,
     *           ' ', zp,' ',alp

      endif

      IF (IPRINT. GE. 5) THEN
         WRITE(*,'(48X,A,E16.8)') 't  =',T
         WRITE(*,'(48X,A,E16.8)') 'rho=',0.5D0*(T*ZP)**2
      ENDIF
C
C    compute f(xn),g(xn)
C
      CALL COMPFG(N,XN,FN,G)
      NCOMP= NCOMP+1
C
      IF (IFAIL .LT. 0 .AND. NCOMP .GE. MAXCOM) THEN
         IFAIL= 1
         IF (IPRINT .GT. 1) WRITE(*,*) 'max. number of calls to compfg'
         MAXIT= NITER1
      ENDIF
c GPP change: abort if we have become stuck in Lanczos
      if (istop .eq. 1) then
         ifail = 1
         maxit = niter1
         IF (IPRINT .GT. 1)
     *       WRITE(*,*) ' Lanczos does not converge under given setup'
         return
      endif
C
C    stopping criterion
C
      IF ((DABS(ZP) .LE. EPS) .AND. (ALP .LE. EPS)) THEN
         IF (IPRINT .GE. 1) THEN
            WRITE(*,*)
            WRITE(*,*) 'convergence'
            WRITE(*,*)
         ENDIF
         IFAIL = 0
         MAXCOM= NCOMP
         MAXIT = NITER1
      ELSE
C
C    compute t
C
         KRIT= .FALSE.
C
C    test serious criterion
C
         IF (FN-F .LT. 0.1D0*V) THEN
            TL= T
            GD= 0.0D0
            DO 34 I=1,N
               GD= GD+G(I)*D(I)
   34       CONTINUE
            IF ((GD .GE. 0.2D0*V) .OR. (TU-T .LE. MU)) THEN
C
C    serious iteration
C
               KRIT= .TRUE.
               IS= MAX(1,IS+1)
               IF ((FN-F .LE. 0.7D0*V) .OR. (IS .GT. 2)) THEN
                  IF (ABS(V-FN+F).GE.ZERO)
     *               TP= -0.5D0*V*T/ABS(V-FN+F)
                  T= DMAX1(TP,5.0D0*T)
                  T= DMIN1(T,SU)
                  IS= 0
               ENDIF
C
C    update alpha
C
               DO 36 I=1,NITER
                  GD= 0.0D0
                  DO 35 J=1,N
                     GD= GD+A(I,J)*D(J)
   35             CONTINUE
                  ALPHA(I)= ALPHA(I)+FN-F-GD
                  ALPHA(I)= DMAX1(ALPHA(I),0.0D0)
   36          CONTINUE
               ALN= 0.0D0
               DO 37 I=1,N
                  X(I)= XN(I)
   37          CONTINUE
               F= FN
            ELSE
               IF (TU .EQ. SU) THEN
                  IF (ABS(V-FN+F).GE.ZERO) TP= -0.5D0*V*T/ABS(V-FN+F)
                  T= DMAX1(TP,5.0D0*T)
                  T= DMIN1(T,SU)
               ELSE
                  IF (ABS(V-FN+F).GE.ZERO) TP= -0.5D0*V*T/ABS(V-FN+F)
                  T= DMAX1(TP,TL+0.25D0*(TU-TL))
                  T= DMIN1(T,TL+0.5D0*(TU-TL))
               ENDIF
               IF (LIN .EQ. 1 .OR. LIN .EQ. 2) THEN
                  DO 38 I=1,N
                     DL(I)= D(I)
   38             CONTINUE
                  TLL= TL
                  ALL= ALP
               ENDIF
            ENDIF
         ELSE
            TU= T
            IF (TL .EQ. 0.0D0) THEN
C
C    test null criterion
C
               ALN= F-FN
               DO 39 I=1,N
                  ALN= ALN+G(I)*D(I)
  39           CONTINUE
C
               IF ((ALN. LE. DMAX1(ALP,EPS)) .OR.
     *            (DABS(F-FN). LE. 1.0D0*(ZM+ALM))) THEN
C
C    null iteration
C
                  KRIT= .TRUE.
                  ALN= DMAX1(ALN,0.0D0)
                  IS= MIN(-1,IS-1)
                  IF (IS .LT. -N) THEN
                     T= 0.5D0*T
                     IS= 0
                  ENDIF
               ELSE
                  T= TL+0.1D0*(TU-TL)
               ENDIF
            ELSE
               T= TL+0.5D0*(TU-TL)
               IF (LIN .EQ. 1 .OR. LIN .EQ. 2) THEN
                  DO 40 I=1,N
                     DU(I)= D(I)
   40             CONTINUE
               TUL= TU
               ALU= ALP
               ENDIF
            ENDIF
         ENDIF
      ENDIF
      IF (IFAIL .LT. 0 .AND. DABS(T-TA) .LT. ZERO*1.0D-3
     *    .AND. .NOT. KRIT) THEN
            IFAIL = 3
            IF (IPRINT .GT. 1) WRITE(*,*) 'T is wrong'
            MAXIT = NITER1
            MAXCOM= NCOMP
      ENDIF
      IF (IFAIL .LT. 0) GOTO 100
      END
      SUBROUTINE QUAD01(M,N,G,T,ALPHA,MODE,K,JS,Y,ADA,V,VDA,
     *                  VSS,W,L,VGXAL,WORK,LWORK,INFORM)
      INTEGER          M,N,MODE,K,L,LWORK,INFORM
      DOUBLE PRECISION T,ADA,V,VDA,VSS,W,VGXAL
      INTEGER          JS(M)
      DOUBLE PRECISION G(M*(M+1)/2),ALPHA(M),Y(M),WORK(LWORK)
C
C    local variables
C
      INTEGER          ITMAX,IDELMX,ITER
      DOUBLE PRECISION EPSTOP,EPSCON
      PARAMETER ( ITMAX  = 500,
     *            EPSTOP = 1.0D-12,
     *            EPSCON = 1.0D-12,
     *            IDELMX = 500      )
C
      CALL QUAD02(ITMAX,M,N,M,G,M*(M+1)/2,1.0D0/T,ALPHA,EPSTOP,
     *            EPSCON,IDELMX,MODE,K,JS,Y,ADA,V,VDA,VSS,
     *            W,L,VGXAL,WORK,LWORK,ITER,INFORM)
      END
      SUBROUTINE QUAD02(ITMAX,M,N,KMAX,G,LENG,U,A,EPSTOP,EPSCON,
     *                  IDELMX,MODE,K,JS,X,ADD,V,VDD,VSS,W,
     *                  L,VGXAL,WORK,LWORK,ITER,INFORM)
C
C
      INTEGER           ITMAX,M,N,KMAX,LENG,IDELMX,MODE,L,LWORK,
     *                  ITER,INFORM
      DOUBLE PRECISION  U,EPSTOP,EPSCON,ADD,V,VDD,VSS,W,VGXAL
      INTEGER           JS(KMAX)
      DOUBLE PRECISION  G(LENG),A(M),X(M),WORK(LWORK)
C
      INTEGER           LWORK1,K,KMAX1,K1,IDELET,KMAXA
      DOUBLE PRECISION  WOLD
C
C
      IF (MODE .LE. 1) THEN
         LWORK1=10+KMAX*(KMAX+11)/2
         IF (M .LT. 1 .OR. N .LT. 1 .OR. KMAX .LT. MIN0(M,N+1).OR.
     *       LENG .LT. M*(M+1)/2 .OR.
     *       EPSTOP .LT. 0.0D0 .OR. EPSCON .LE. 0.0D0 .OR.
     *       MODE .LT. 1 .OR. LWORK .LT. LWORK1)  THEN
            INFORM = 5
            RETURN
         ENDIF
         KMAX1=   KMAX
         WORK(1)= LWORK1
         WORK(2)= KMAX1
      ELSE
         LWORK1= WORK(1)+0.1D0
         IF (LWORK .LT. LWORK1) THEN
            INFORM= 5
            RETURN
         ENDIF
         KMAX1= WORK(2)+0.1D0
         K1   = WORK(3)+0.1D0
         IF (KMAX1 .LT. MIN0(M,N+1) .OR. EPSTOP .LT. 0.0D0 .OR.
     *       EPSCON .LE. 0.0D0 .OR. K .NE. K1) THEN
            INFORM= 5
            RETURN
         ENDIF
         WOLD  = WORK(4)
         IDELET= WORK(5)+0.1D0
      ENDIF
      KMAXA = MIN0(M,N+1)
      CALL QUAD03(ITMAX,M,N,KMAXA,G,LENG,U,A,EPSTOP,EPSCON,IDELMX,MODE,
     *            K,JS,X,ADD,V,VDD,VSS,W,L,VGXAL,ITER,INFORM,
     *            WORK(10),WORK(KMAX1+10),WORK(2*KMAX1+10),
     *            WORK(3*KMAX1+10),WORK(4*KMAX1+10),WORK(5*KMAX1+10),
     *            KMAXA*(KMAXA+1)/2,WOLD,IDELET)
      WORK(3)= K
      WORK(4)= WOLD
      WORK(5)= IDELET
      END
      SUBROUTINE QUAD03(ITMAX,M,N,KMAX,G,LENG,U,A,EPSTOP,EPSCON,IDELMX,
     *                  MODE,K,JS,X,ADD,V,VDD,VSS,W,L,VGXAL,ITER,
     *                  INFORM,Y,YHAT,YSS,YA,YE,R,MDIMR,WOLD,IDELET)
C
      INTEGER           ITMAX,M,N,KMAX,LENG,IDELMX,MODE,K,L,ITER,
     *                  INFORM,MDIMR,IDELET
      DOUBLE PRECISION  U,EPSTOP,EPSCON,ADD,V,VDD,VSS,W,VGXAL,
     *                  WOLD
      INTEGER           JS(KMAX)
      DOUBLE PRECISION  G(LENG),A(M),X(M),Y(KMAX),YHAT(KMAX),
     *                  YSS(KMAX),YA(KMAX),YE(KMAX),R(MDIMR)
C
      INTEGER           NINCRW,INCRW,INCRWM,KDIMR,ISETYA,ISETYE,NAUG,
     *                  NEXC,NDEL,NREF,IWOLD,JDEL,JADD,I,J,JA,LG,KG,
     *                  IA,JR,JG,IG,KR,JROT,JROTAT,ISTEP,KMIN1,KA,J1,
     *                  KREFAC
      DOUBLE PRECISION  T,ERRKT,ERRX,EPS1,EPS2,HALF,ONE,TWO,S,DELTA,
     *                  DELTAP,DELTAW,DMU,ERRORV,R1,R2,RHOSQ,VIOLAT,
     *                  WMIN,YBARL,YANEW,YENEW,DTAU,DGAMMA,DSIGMA,DNU,
     *                  DALPHA,S1,YAYE,YESQR,DABSS
C
      COMMON /QPLOG/ ERRKT,ERRX,NAUG,NEXC,NDEL
      DATA  EPS1, EPS2, HALF, ONE, TWO
     *    /0.5D0,1.0D-2,0.5D0,1.0D0,2.0D0/
C
C
      NINCRW= 3
      INCRW =-1
      INCRWM=-1
      IF (MODE .GT. 1) KDIMR= K*(K+1)/2
      ISETYA= 0
      IF (MODE .EQ. 2) ISETYA=1
      ISETYE= 0
      NAUG  = 0
      NEXC  = 0
      NDEL  = 0
      NREF  = 0
      ITER  =-1
      IWOLD = 0
      JDEL  = 0
      JADD  = 0
      T     = 0.0D0
      GOTO (10,500,100,151,500,500,500,513),MODE
   10 CONTINUE
   50 W= HALF*G(1)+U*A(1)
      L= 1
      KG=0
      DO 60 J=1,M
         KG= KG+J
         S= HALF*G(KG)+U*A(J)
         IF (W .GT. S) THEN
            W=S
            L=J
         ENDIF
   60 CONTINUE
   70 CONTINUE
      K      = 1
      KDIMR  = 1
      JS(1)  = L
      YHAT(1)= ONE
      LG     = L*(L+1)/2
      V      =-(G(LG)+U*A(L))
      IDELET = 0
      R(1)   = DSQRT(ONE+G(LG))
      YA(1)  = -A(L)/R(1)
      YE(1)  = ONE/R(1)
      JADD   = L
  100 CONTINUE
      DO 110 J=1,M
         X(J)=0.0D0
  110 CONTINUE
      ERRX= 0.0D0
      ADD= 0.0D0
      DO 120 J=1,K
         JA   = JS(J)
         X(JA)= YHAT(J)
         ERRX = ERRX+YHAT(J)
         ADD = ADD+A(JA)*YHAT(J)
  120 CONTINUE
      ERRX = ONE -ERRX
      VDD = 0.0D0
      VGXAL= 0.0D0
      ERRKT= 0.0D0
      DO 150 I=1,M
        S=0.0D0
        DO 140 J=1,K
           JA= JS(J)
           IG= MIN0(I,JA)
           JG= MAX0(I,JA)
           KG= (JG-1)*JG/2+IG
           S=  S+G(KG)*YHAT(J)
  140   CONTINUE
C       S=SUM OF G(I,J)*X(J) FOR J=1 THROUGH M.
        VDD= VDD+X(I)*S
        S   = S+U*A(I)
        IF (I .EQ. 1) VSS= -S
        VSS= DMAX1(VSS,-S)
        S= S+V
        DABSS= DABS(S)
        IF (X(I) .GT. 0.0D0) ERRKT=DMAX1(ERRKT,DABSS)
        IF (S .GE. VGXAL .AND. I .GT. 1) GOTO 150
        VGXAL=S
        L=I
  150 CONTINUE
      W   = HALF*VDD+U*ADD
      VDD=-(VDD+U*ADD)
  151 CONTINUE
      LG= L*(L+1)/2
      VIOLAT= VGXAL+EPSTOP*(ONE+G(LG))
      ERRORV= DMAX1(DABS(V-VDD),DABS(V-VSS),DABS(VDD-VSS))
     *        /(ONE+DABS(V))
      IF (IWOLD .EQ. 0) WOLD= W
      IWOLD = 1
      DELTAW= WOLD-W
      WOLD  = W
      IF (DELTAW .LE. 0.0D0) INCRW= INCRW+1
      IF(VIOLAT.GE.0.0D0) GOTO 800
      DO 160 J=1,K
         IF (L .EQ. JS(J)) GOTO 8160
  160 CONTINUE
      IF (ITER  .GE. ITMAX)  GOTO 8161
      IF (INCRW .EQ. NINCRW) GOTO 8162
      IF(INCRWM .EQ. -1) WMIN= W
      IF(W .EQ. WMIN) INCRWM= INCRWM+1
      IF(W .GE. WMIN) GOTO 170
         WMIN  = W
         INCRWM= 0
  170 CONTINUE
      IF (INCRWM .EQ. NINCRW) GOTO 8170
      DO 212 J=1,K
         JA  = JS(J)
         IG  = MIN0(L,JA)
         JG  = MAX0(L,JA)
         KG  = (JG-1)*JG/2+IG
         X(J)= G(KG)
         YSS(J)= ONE+G(KG)
  212 CONTINUE
      CALL QUAD04(R,YSS,KDIMR,K,1)
      DO 215 J=1,K
         Y(J)= YSS(J)
  215 CONTINUE
      RHOSQ= 0.0D0
      DO 220 J=1,K
         RHOSQ=RHOSQ+YSS(J)**2
  220 CONTINUE
      RHOSQ= (ONE+G(LG))-RHOSQ
      IF (RHOSQ .GT. EPSCON*(ONE+G(LG)) .AND. K .LT. KMAX) GOTO 270
      CALL QUAD04(R,YSS,KDIMR,K,0)
      YBARL= 0.0D0
      DO 230 J=1,K
         YBARL= YBARL+YSS(J)
  230 CONTINUE
      DELTA = YBARL-ONE
      DELTAP= 0.0D0
      DO 235 I=1,K
        IA= JS(I)
        S = 0.0D0
        DO 232 J=1,K
           JA= JS(J)
           IG= MIN0(IA,JA)
           JG= MAX0(IA,JA)
           KG= (JG-1)*JG/2+IG
           S=S+G(KG)*YSS(J)
  232   CONTINUE
  235 DELTAP= DELTAP+YSS(I)*S
      S= 0.0D0
      DO 240 J=1,K
         S=S+YSS(J)*X(J)
  240 CONTINUE
      DELTAW= (DELTAP+G(LG)*YBARL**2)-TWO*YBARL*S
      DELTAP= DABS((DELTAP+G(LG))-TWO*S)
      S     = RHOSQ-(DELTA**2+DELTAP)
      IF (DELTA .LT. -EPS1 .AND. K .LT. KMAX) GOTO 260
      T = 100./(ONE-EPS1)
      JR= 0
      DO 250 J=1,K
         IF (YSS(J) .LE. 0.0D0) GOTO 250
         IF (T*YSS(J) .LT. YHAT(J)) GOTO 250
         T = YHAT(J)/YSS(J)
         JR= J
  250 CONTINUE
      IF (HALF*T*DELTAW .LT. (EPS2-ONE-DELTA)*VGXAL) GOTO 400
      IF (K .EQ. KMAX .AND. JR .NE. 0) GOTO 400
      IF (K .EQ. KMAX) GOTO 8250
  260 CONTINUE
      RHOSQ= DMAX1(RHOSQ,DELTA**2+DELTAP)
      IF(RHOSQ .EQ. 0.0D0) GOTO 8250
  270 CONTINUE
      NAUG = NAUG+1
      ISTEP= 3
      YANEW=-A(L)
      YENEW= ONE
      KDIMR= KDIMR+1
      DO 310 J=1,K
         R(KDIMR)= Y(J)
         YANEW= YANEW-Y(J)*YA(J)
         YENEW= YENEW-Y(J)*YE(J)
         KDIMR= KDIMR+1
  310 CONTINUE
      R(KDIMR)= DSQRT(RHOSQ)
      K      = K+1
      JS(K)  = L
      YHAT(K)= 0.0D0
      YA(K)  = YANEW/R(KDIMR)
      YE(K)  = YENEW/R(KDIMR)
      JADD   = L
      T      = 0.0D0
      GOTO 500
  400 CONTINUE
      NEXC   = NEXC+1
      DO 410 J=1,K
         YHAT(J)= DMAX1(YHAT(J)-T*YSS(J),0.0D0)
  410 CONTINUE
      ISTEP= 4
  415 CONTINUE
      JDEL= JS(JR)
      IF (JR .EQ. K) GOTO 446
      IDELET= IDELET+1
      KMIN1 = K-1
      DO 420 J=JR,KMIN1
         YHAT(J)= YHAT(J+1)
         JS(J)=JS(J+1)
  420 CONTINUE
      IF (ISTEP .NE. 4) GOTO 422
      DO 421 J=JR,KMIN1
         X(J)=X(J+1)
  421 CONTINUE
  422 CONTINUE
      JROTAT= JR+1
      KR= JR*(JR+1)/2-1
      DO 431 JROT=JROTAT,K
         KR=     KR+JROT
         R1=     R(KR)
         R2=     R(KR+1)
         DMU=    DMAX1(DABS(R1),DABS(R2))
         DTAU=   DMU*DSQRT((R1/DMU)**2+(R2/DMU)**2)
         R(KR)=  DTAU
         DGAMMA= R1/DTAU
         DSIGMA= R2/DTAU
         DNU   = DSIGMA/(ONE+DGAMMA)
         DALPHA= DGAMMA*YA(JROT-1)+DSIGMA*YA(JROT)
         YA(JROT)  = DNU*(YA(JROT-1)+DALPHA)-YA(JROT)
         YA(JROT-1)= DALPHA
         DALPHA= DGAMMA*YE(JROT-1)+DSIGMA*YE(JROT)
         YE(JROT)= DNU*(YE(JROT-1)+DALPHA)-YE(JROT)
         YE(JROT-1)= DALPHA
         IF (JROT .EQ. K) GOTO 435
         KA= KR
         DO 430 J=JROT,KMIN1
            KA= KA+J
            DALPHA = DGAMMA*R(KA)+DSIGMA*R(KA+1)
            R(KA+1)= DNU*(R(KA)+DALPHA)-R(KA+1)
            R(KA)  = DALPHA
  430    CONTINUE
  431 CONTINUE
  435 KR= (JR-1)*JR/2+1
      KA= KR+JR
      DO 445 J=JR,KMIN1
        DO 440 J1=1,J
           R(KR)= R(KA)
           KR= KR+1
           KA= KA+1
  440   CONTINUE
        KA=KA+1
  445 CONTINUE
  446 CONTINUE
    6 KDIMR= KDIMR-K
      K    = K-1
      IF (DABS(A(JDEL)) .GT. ONE .AND. JR .LE. K) ISETYA= 1
      IF(ISTEP .EQ. 6) GOTO 500
  450 CONTINUE
      YANEW=-A(L)
      YENEW= ONE
      IF (K .GT. 0) GOTO 460
      KDIMR = 1
      IDELET= 0
      S1= 0.0D0
      GOTO 480
  460 CONTINUE
      IF(ISTEP.NE.4) GOTO 462
      DO 461 J=1,K
         YSS(J)=ONE+X(J)
  461 CONTINUE
      GOTO 472
  462 CONTINUE
      DO 471 J=1,K
         JA= JS(J)
         IG= MIN0(L,JA)
         JG= MAX0(L,JA)
         KG= (JG-1)*JG/2+IG
         YSS(J)=ONE+G(KG)
  471 CONTINUE
  472 CONTINUE
      CALL QUAD04(R,YSS,KDIMR,K,1)
      KDIMR= KDIMR+1
      S1   = 0.0D0
      DO 473 J=1,K
         S1= S1+YSS(J)**2
         R(KDIMR)= YSS(J)
         YANEW= YANEW-YSS(J)*YA(J)
         YENEW= YENEW-YSS(J)*YE(J)
         KDIMR= KDIMR+1
  473 CONTINUE
  480 R(KDIMR)= DSQRT(DMAX1((ONE+G(LG))-S1,0.0D0))
      K= K+1
      JS(K)= L
      IF (ISTEP .EQ. 4) YHAT(K)= T*(ONE+DELTA)
      IF (R(KDIMR) .EQ. 0.0D0) GOTO 8480
      YA(K)= YANEW/R(KDIMR)
      YE(K)= YENEW/R(KDIMR)
      IF(ISTEP .EQ. 7) GOTO 750
      JADD= L
  500 CONTINUE
  505 ITER= ITER+1
      IF (ISETYA .EQ. 0 .AND. ISETYE .EQ. 0) GOTO 513
      DO 510 J=1,K
         JA=JS(J)
         IF(ISETYA.EQ.1) YA(J)=-A(JA)
         IF(ISETYE.EQ.1) YE(J)=ONE
  510 CONTINUE
      IF(ISETYA .EQ. 1) CALL QUAD04(R,YA,KDIMR,K,1)
      IF(ISETYE .EQ. 1) CALL QUAD04(R,YE,KDIMR,K,1)
      ISETYA= 0
      ISETYE= 0
  513 CONTINUE
      YESQR= 0.0D0
      YAYE = 0.0D0
      DO 515 J=1,K
         YESQR= YESQR+YE(J)**2
         YAYE=YAYE+YE(J)*YA(J)
  515 CONTINUE
  516 CONTINUE
      V= (ONE-U*YAYE)/YESQR
      DO 517 J=1,K
         Y(J)=U*YA(J)+V*YE(J)
  517 CONTINUE
  540 CALL QUAD04(R,Y,KDIMR,K,0)
      V=  ONE-V
  550 CONTINUE
      JR= 0
      DO 551 J=1,K
         IF(Y(J).LE.0.0D0) JR=J
  551 CONTINUE
      JADD= 0
      IF (JR .GT. 0) GOTO 600
      DO 560 J=1,K
         YHAT(J)= Y(J)
  560 CONTINUE
      JDEL= 0
      T   = ONE
      GOTO 100
  600 CONTINUE
      NDEL= NDEL+1
      T   = ONE
      DO 610 J=1,K
        S1= YHAT(J)-Y(J)
        IF (S1 .LE. 0.0D0) GOTO 610
        S1= YHAT(J)/S1
        IF (T .LT. S1) GOTO 610
        T = S1
        JR= J
  610 CONTINUE
      DO 620 J=1,K
         YHAT(J)=DMAX1(YHAT(J)+T*(Y(J)-YHAT(J)),0.0D0)
  620 CONTINUE
      ISTEP= 6
      GOTO 415
  700 CONTINUE
      NREF  = NREF+1
      ISTEP = 7
      IDELET= 0
      KREFAC= K
      K     = 0
  750 CONTINUE
      IF (K .EQ. KREFAC) GOTO 760
      L = JS(K+1)
      LG= L*(L+1)/2
      GOTO 450
  760 CONTINUE
      ISETYA= 1
      ISETYE= 1
      JDEL  =-1000
      JADD  = 1000
      T     = 0.0D0
      GOTO 500
  800 INFORM= 0
  805 DO 810 J=1,M
         X(J)=0.0D0
  810 CONTINUE
      ERRX= 0.0D0
      DO 820 J=1,K
         JA   = JS(J)
         ERRX = ERRX+YHAT(J)
         X(JA)=YHAT(J)
  820 CONTINUE
      ERRX= ONE-ERRX
      RETURN
 8160 IF(IDELET .GT. IDELMX) GOTO 700
      INFORM= 1
      GOTO 805
 8161 INFORM= 4
      GOTO 805
 8162 INFORM= 2
      GOTO 805
 8170 INFORM= 2
      GOTO 805
 8250 IF(IDELET .GT. IDELMX) GOTO 700
      INFORM= 3
      GOTO 805
 8480 CONTINUE
      IF(IDELET .GT. 0) GOTO 700
      NDEL= NDEL+1
      IF (INCRW .EQ. NINCRW-1) INCRW=INCRW-1
      JDEL =-JS(K)
      JADD = 0
      T    = 0.0D0
      KDIMR= KDIMR-K
      K    = K-1
      S1   = 0.0D0
      DO 8484 J=1,K
         S1=S1+YHAT(J)
 8484 CONTINUE
      IF (S1 .EQ. 0.0D0) GOTO 8486
      DO 8485 J=1,K
         YHAT(J)=YHAT(J)/S1
 8485 CONTINUE
      ISETYA= 1
      ISETYE= 1
      GOTO 500
 8486 YHAT(1)= ONE
      ISETYA= 1
      ISETYE= 1
      GOTO 500
      END
      SUBROUTINE QUAD04(R,Y,MDIMR,M,ITRANS)
      INTEGER           MDIMR,M,ITRANS
      DOUBLE PRECISION  R(MDIMR),Y(M)
C
      INTEGER           I,J,K,KA,I1,IA,IOLD
      DOUBLE PRECISION  S
C
C
      IF (ITRANS.EQ.1) GOTO 3
      I= M
      K= MDIMR
      Y(I)= Y(I)/R(K)
    1 IF (I .EQ. 1) RETURN
      I1= I
      K = K-I
      KA= K
      I = I-1
      J = I
      S = Y(I)
      DO 2 IA=I1,M
        KA= KA+J
        S = S-R(KA)*Y(IA)
        J = J+1
    2 CONTINUE
      Y(I)= S/R(K)
      GOTO 1
    3 Y(1)= Y(1)/R(1)
      I= 1
      K= 1
    4 IF (I .EQ. M) RETURN
      IOLD= I
      I= I+1
      K= K+1
      S= Y(I)
      DO 5 IA=1,IOLD
        S= S-R(K)*Y(IA)
        K=K+1
    5 CONTINUE
      Y(I)= S/R(K)
      GOTO 4
      END
      SUBROUTINE MLIS4 (SIMUL,N,XN,FN,FPN,T,TMIN,TMAX,D,D2,G,
     *                  AMD,AMF,LOGIC,NAP,NAPMAX,X,TOL,A,TPS,TNC)
C
      INTEGER          N,LOGIC,NAP,NAPMAX
      DOUBLE PRECISION T,TMIN,TMAX,FN,FPN,A,D2,TOL,AMF,AMD,
     *                 TPS,TNC
      DOUBLE PRECISION XN(N),D(N),G(N),X(N)
      EXTERNAL SIMUL
C
C     Modification of mlis2 (C. Lemarechal)
C
      INTEGER          INDIC,INDICA,INDICD,I
      DOUBLE PRECISION TA,TG,TD,TC,TP,TESF,TESD,FPG,PS,Z1,FFN,FG,FA,FP,
     *                 FPA,FPD,SIGN,DEN,F,BARMIN,BARMUL,BARMAX,BARR,
     *                 FD,P,Z,TEST,ANUM,MC,MP
C
C
      INDIC=4
      TESF=AMF*FPN
      TESD=AMD*FPN
      TD=0.
      TG=0.
      FG=FN
      FPG=FPN
      TA=0.
      FA=FN
      FPA=FPN
      INDICA=1
      LOGIC=0
      BARMIN=0.01
      BARMUL=3.
      BARMAX=.3
      BARR=BARMIN
      IF (T.GT.TMIN) GO TO 20
      T=TMIN
      IF (T.LE.TMAX) GO TO 20
      TMIN=TMAX
   20 IF (FN+T*FPN.LT.FN+0.9*T*FPN) GO TO 30
      T=2.*T
      GO TO 20
   30  IF (T.GE.TMAX) THEN
          T=TMAX
          LOGIC=1
       ENDIF
      DO 50 I=1,N
   50 X(I)=XN(I)+T*D(I)
  100 NAP=NAP+1
      IF (NAP.GT.NAPMAX) THEN
          LOGIC=4
          IF (TG.EQ.0.) GO TO 999
          DO 120 I=1,N
  120     XN(I)=XN(I)+TG*D(I)
          INDIC=4
          CALL SIMUL (N,XN,FN,G)
          GO TO 999
      ENDIF
      INDIC=4
      CALL SIMUL(N,X,F,G)
      IF (INDIC.EQ.0) THEN
          LOGIC=5
          FN=F
          DO 170 I=1,N
  170     XN(I)=X(I)
          GO TO 999
      ENDIF
  200 IF(INDIC.GT.0) GO TO 210
      TD=T
      INDICD=INDIC
      LOGIC=0
      T=TG+0.1*(TD-TG)
      GO TO 905
  210 PS= 0.0
      DO 211 I=1,N
         PS= PS+G(I)*D(I)
  211 CONTINUE
      FP=PS
      FFN=F-FN
      IF(FFN.LT.T*TESF) GO TO 300
      TD=T
      FD=F
      FPD=FP
      INDICD=INDIC
      LOGIC=0
      IF(TG.NE.0.) GO TO 500
      IF(FPD.LT.TESD) GO TO 500
      TPS=(FN-F)+TD*FPD
      TNC=D2*TD*TD
      P=DMAX1(A*TNC,TPS)
      IF(P.GT.TOL) GO TO 500
      LOGIC=3
      FPN=FPD
      GO TO 999
  300 CONTINUE
      IF(FP.LT.TESD) GO TO 320
      LOGIC=0
      FN=F
      FPN=FP
      DO 310 I=1,N
  310 XN(I)=X(I)
      GO TO 999
  320 IF (LOGIC.EQ.0) GO TO 350
      FN=F
      FPN=FP
      DO 330 I=1,N
  330 XN(I)=X(I)
      GO TO 999
  350 TG=T
      FG=F
      FPG=FP
      TA=T
      T=9.*TG
      Z=FPN+3.*FP-4.*FFN/TG
      IF(Z.GT.0.) T=DMIN1(T,TG*DMAX1(1.0D0,-FP/Z))
      T=TG+T
      IF(T.LT.TMAX) GO TO 900
      LOGIC=1
      T=TMAX
      GO TO 900
  500 IF(INDICA.GT.0 .AND. INDICD.GT.0) GO TO 510
      TA=T
      T=0.9*TG+0.1*TD
      GO TO 900
  510 TEST=BARR*(TD-TG)
      PS=FP+FPA-3.D0*(FA-F)/(TA-T)
      Z1=PS*PS-FP*FPA
      IF (Z1.GE.0.D0) GO TO 520
      IF (FP.LT.0.) TC=TD
      IF (FP.GE.0.) TC=TG
      GO TO 600
  520 Z1=DSQRT(Z1)
      IF (T-TA.LT.0.) Z1=-Z1
      SIGN=(T-TA)/ABS(T-TA)
      IF ((PS+FP)*SIGN.GT.0.) GO TO 550
      DEN=2.D0*PS+FP+FPA
      ANUM=Z1-FP-PS
      IF (ABS((T-TA)*ANUM).GE.(TD-TG)*ABS(DEN)) GO TO 530
      TC=T+ANUM*(TA-T)/DEN
      GO TO 600
  530 TC=TD
      GO TO 600
  550 TC=T+FP*(TA-T)/(PS+FP+Z1)
  600 MC=0
      IF (TC.LT.TG) MC=-1
      IF (TC.GT.TD) MC=1
      TC=DMAX1(TC,TG+TEST)
      TC=DMIN1(TC,TD-TEST)
      PS=FPD-FPG
      IF (PS.NE.0.D0) GO TO 620
      TP=0.5*(TD+TG)
      GO TO 650
  620 TP=((FG-FPG*TG)-
     * (FD-FPD*TD))/PS
  650 MP=0
      IF (TP.LT.TG) MP=-1
      IF (TP.GT.TD) MP=1
      TP=DMAX1(TP,TG+TEST)
      TP=DMIN1(TP,TD-TEST)
      TA=T
      IF (MC.EQ.0  .AND.  MP.EQ.0) T=DMIN1(TC,TP)
      IF (MC.EQ.0  .AND.  MP.NE.0) T=TC
      IF (MC.NE.0  .AND.  MP.EQ.0) T=TP
      IF (MC.EQ.1  .AND.  MP.EQ.1) T=TD-TEST
      IF (MC.EQ.-1 .AND. MP.EQ.-1) T=TG+TEST
      IF (MC*MP.EQ.-1) T=0.5*(TG+TD)
      IF (MC.EQ.0 .AND. MP.EQ.0) THEN
          BARR=BARMIN
        ELSE
          BARR=DMIN1(BARMUL*BARR,BARMAX)
      ENDIF
  900 FA=F
      FPA=FP
  905 INDICA=INDIC
  920 IF (TD.EQ.0.) GO TO 990
      IF (TD-TG.LE.TMIN) GO TO 950
      DO 930 I=1,N
      Z=XN(I)+T*D(I)
      IF (Z.NE.X(I) .AND. Z.NE.XN(I)) GO TO 990
  930 CONTINUE
  950 LOGIC=6
      IF (INDICD.LT.0) LOGIC=INDICD
      IF (TG.EQ.0.) GO TO 970
      DO 960 I=1,N
  960 XN(I)=XN(I)+TG*D(I)
      INDIC=4
      CALL SIMUL(N,XN,FN,G)
  970 CONTINUE
      RETURN
  990 DO 995 I=1,N
  995 X(I)=XN(I)+T*D(I)
      GO TO 100
  999 RETURN
      END


C   VERSION 2    DOES NOT USE EISPACK
C
C ------------------------------------------------------------------
C
      SUBROUTINE DNLASO(OP, IOVECT, N, NVAL, NFIG, NPERM,
     *   NMVAL, VAL, NMVEC, VEC, NBLOCK, MAXOP, MAXJ, WORK,
     *   IND, IERR)
C
      INTEGER N, NVAL, NFIG, NPERM, NMVAL, NMVEC, NBLOCK,
     *   MAXOP, MAXJ, IND(1), IERR
      DOUBLE PRECISION VEC(NMVEC,1), VAL(NMVAL,1), WORK(1)
      EXTERNAL OP, IOVECT
C
C AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT
C
C COMPUTER SCIENCES DEPARTMENT
C UNIVERSITY OF TEXAS AT AUSTIN
C AUSTIN, TX 78712
C
C VERSION 2 ORIGINATED APRIL 1982
C
C CURRENT VERSION  JUNE 1983
C
C DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF
C THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX.  THE SUBROUTINE
C DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS
C THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND
C SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA.
C HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING
C PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS.
C DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES.  IF
C THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE
C NEGATIVE OF THE MATRIX.
C
C
C ON INPUT
C
C
C   OP   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
C     OP(N,M,P,Q).  P AND Q ARE N X M MATRICES AND Q IS
C     RETURNED AS THE MATRIX TIMES P.
C
C   IOVECT   A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE
C     IOVECT(N,M,Q,J,K).  Q IS AN N X M MATRIX.  IF K = 0
C     THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH
C     THE JTH LANCZOS VECTORS.  IF K = 1 THEN Q IS RETURNED
C     AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS.  SEE
C     DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES.
C
C   N   THE ORDER OF THE MATRIX.
C
C   NVAL   NVAL SPECIFIES THE EIGENVALUES TO BE FOUND.
C     DABS(NVAL)  IS THE NUMBER OF EIGENVALUES DESIRED.
C     IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST)
C     EIGENVALUES ARE FOUND.  IF NVAL > 0 THE ALGEBRAICALLY
C     LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND.  NVAL MUST NOT
C     BE ZERO.  DABS(NVAL) MUST BE LESS THAN  MAXJ/2.
C
C   NFIG   THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE
C     EIGENVALUES.  NFIG MUST BE GREATER THAN OR EQUAL TO 1.
C
C   NPERM   AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER
C     SUPPLIED EIGENPAIRS.  IN MOST CASES NPERM WILL BE ZERO.  SEE
C     DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER
C     THAN ZERO.  NPERM MUST NOT BE LESS THAN ZERO.
C
C   NMVAL   THE ROW DIMENSION OF THE ARRAY VAL.  NMVAL MUST BE GREATER
C     THAN OR EQUAL TO DABS(NVAL).
C
C   VAL   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
C     DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4.  IF NPERM
C     IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED
C     IN VAL.  SEE DOCUMENTATION FOR DETAILS.
C
C   NMVEC   THE ROW DIMENSION OF THE ARRAY VEC.  NMVEC MUST BE GREATER
C     THAN OR EQUAL TO N.
C
C   VEC   A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW
C     DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL).  IF
C     NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST
C     CONTAIN THE USER SUPPLIED EIGENVECTORS.
C
C   NBLOCK   THE BLOCK SIZE.  SEE DOCUMENTATION FOR CHOOSING
C     AN APPROPRIATE VALUE FOR NBLOCK.  NBLOCK MUST BE GREATER
C     THAN ZERO AND LESS THAN  MAXJ/6.
C
C   MAXOP   AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE
C     OP.  DNLASO TERMINATES WHEN MAXOP IS EXCEEDED.  SEE
C     DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP.
C
C   MAXJ   AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND
C     DOCUMENTATION ON IOVECT).  FOR THE FASTEST CONVERGENCE MAXJ
C     SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE
C     MAXJ LARGER THAN MAXOP*NBLOCK.
C
C   WORK   A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS
C     LARGE AS
C
C         2*N*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK*NBLOCK + 3*NV
C
C            + THE MAXIMUM OF
C                 N*NBLOCK
C                   AND
C         MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1)
C
C     WHERE NV = DABS(NVAL)
C
C     THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED
C     STARTING VECTORS.  SEE DOCUMENTATION FOR GUIDELINES IN
C     CHOOSING STARTING VECTORS.
C
C   IND   AN INTEGER ARRAY OF DIMENSION AT LEAST DABS(NVAL).
C
C   IERR   AN INTEGER VARIABLE.
C
C
C ON OUTPUT
C
C
C   NPERM   THE NUMBER OF EIGENPAIRS NOW KNOWN.
C
C   VEC   THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS.
C
C   VAL   THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING
C     EIGENVALUES.  THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF
C     THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN-
C     VALUES.  THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES
C     OF THE ACCURACY OF THE EIGENVALUES.  THE FOURTH COLUMN CONTAINS
C     ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS.  SEE
C     DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES.
C
C   WORK   IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2)
C     THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS
C     FOR RESTARTING THE ALGORITHM AND DNLASO CAN BE IMMEDIATELY
C     RECALLED TO CONTINUE WORKING ON THE PROBLEM.
C
C   IND   IND(1)  CONTAINS THE ACTUAL NUMBER OF CALLS TO OP.  ON SOME
C     OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER
C     THAN MAXOP.
C
C   IERR   AN ERROR COMPLETION CODE.  THE NORMAL COMPLETION CODE IS
C     ZERO.  SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO
C     COMPLETION CODES.
C
C
C INTERNAL VARIABLES.
C
C
      INTEGER I, I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11,
     *  I12, I13, M, NBAND, NOP, NV, IABS, MAX0
      LOGICAL RARITZ, SMALL
      DOUBLE PRECISION DELTA, EPS, DDOT, DNRM2, DABS, TARR(1)
C     EXTERNAL DNPPLA, DNWLA, DMVPC, DORTQR, DAXPY,
      EXTERNAL DNPPLA, DMVPC, DORTQR, DAXPY,
     *  DCOPY, DDOT, DNRM2, DSCAL, DSWAP, DLAEIG, DLABCM,
     *  DLABFC, DLAGER, DLARAN, DVSORT
C
C NOP   RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE
C   SUBROUTINE OP.
C
C NV   SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED,
C   AND PASSED TO DNWLA.
C
C SMALL   SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED.
C
C RARITZ   RETURNED FROM DNWLA AND PASSED TO DNPPLA.  RARITZ IS .TRUE.
C   IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED.
C
C DELTA   RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX
C   WHICH IS CLOSEST TO THE DESIRED EIGENVALUES.
C
C DNPPLA   A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED
C   BY DNWLA.
C
C DNWLA   A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM
C   WITH SELECTIVE ORTHOGONALIZATION.
C
C DMVPC   A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND
C   ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS.
C
C DORTQR   A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS
C   USING HOUSEHOLDER REFLECTIONS.
C
C DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP   A SUBSET OF THE BASIC LINEAR
C   ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION.
C
C DLARAN   A SUBROUTINE TO GENERATE RANDOM VECTORS
C
C DLAEIG, DLAGER, DLABCM, DLABFC   SUBROUTINES FOR BAND EIGENVALUE
C   CALCULATIONS.
C
C ------------------------------------------------------------------
C
C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS.
C
      NV = IABS(NVAL)
      IND(1) = 0
      IERR = 0
      IF (N.LT.6*NBLOCK) IERR = 1
      IF (NFIG.LE.0) IERR = IERR + 2
      IF (NMVEC.LT.N) IERR = IERR + 4
      IF (NPERM.LT.0) IERR = IERR + 8
      IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16
      IF (NV.LT.MAX0(1,NPERM)) IERR = IERR + 32
      IF (NV.GT.NMVAL) IERR = IERR + 64
      IF (NV.GT.MAXOP) IERR = IERR + 128
      IF (NV.GE.MAXJ/2) IERR = IERR + 256
      IF (NBLOCK.LT.1) IERR = IERR + 512
      IF (IERR.NE.0) RETURN
C
      SMALL = NVAL.LT.0
C
C ------------------------------------------------------------------
C
C THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS.
C IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION
C OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO  -1
C AND DNLASO TERMINATES.
C
      IF (NPERM.EQ.0) GO TO 110
C
C THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST
C EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE
C NEGATIVE OF THE MATRIX.
C
      IF (SMALL) GO TO 20
      DO 10 I=1,NPERM
         VAL(I,1) = -VAL(I,1)
   10 CONTINUE
C
C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS.
C
   20 CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TARR, NMVEC, N, VEC)
C
C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON.
C IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE.
C
      DO 60 I=1,NPERM
         VAL(I,2) = DABS(VAL(I,2))
         VAL(I,3) = DNRM2(N,VEC(1,I),1)
   60 CONTINUE
C
C THIS PERFORMS THE ORTHONORMALIZATION.
C
      M = N*NBLOCK + 1
      CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M))
      M = N*NBLOCK - NPERM
      DO 70 I = 1, NPERM
         M = M + NPERM + 1
         IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 70
         IERR = -1
         RETURN
C
   70 CONTINUE
C
C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN
C THE ARRAY WORK FOR LATER REFERENCE IN DNWLA.
C
      M = 2*N*NBLOCK + 1
      CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1)
C
C THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE
C PRECISION
C
C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT
C ***IN A PRODUCTION CODE
C
  110 EPS = 1.0D0
c      DO 120 I = 1,1000
c         EPS = 0.5D0*EPS
c         TEMP = 1.0D0 + EPS
c         IF(TEMP.EQ.1.0D0) GO TO 130
c  120 CONTINUE
        eps = 1.0d-14
C
C ------------------------------------------------------------------
C
C THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM
C WITH SELECTIVE ORTHOGONALIZATION.
C
  130 NBAND = NBLOCK + 1
      I1 = 1 + N*NBLOCK
      I2 = I1 + N*NBLOCK
      I3 = I2 + NV
      I4 = I3 + NV
      I5 = I4 + NV
      I6 = I5 + MAXJ*NBAND
      I7 = I6 + NBLOCK*NBLOCK
      I8 = I7 + NBLOCK*NBLOCK
      I9 = I8 + MAXJ*(NV+1)
      I10 = I9 + NBLOCK
      I11 = I10 + 2*NV + 6
      I12 = I11 + MAXJ*(2*NBLOCK+1)
      I13 = I12 + MAXJ
      CALL DNWLA(OP, IOVECT, N, NBAND, NV, NFIG, NPERM, VAL, NMVEC,
     *   VEC, NBLOCK, MAXOP, MAXJ, NOP, WORK(1), WORK(I1),
     *   WORK(I2), WORK(I3), WORK(I4), WORK(I5), WORK(I6),
     *   WORK(I7), WORK(I8), WORK(I9), WORK(I10), WORK(I11),
     *   WORK(I12), WORK(I13), IND, SMALL, RARITZ, DELTA, EPS, IERR)
C
C ------------------------------------------------------------------
C
C THIS SECTION CALLS DNPPLA (THE POST PROCESSOR).
C
      IF (NPERM.EQ.0) GO TO 140
      I1 = N*NBLOCK + 1
      I2 = I1 + NPERM*NPERM
      I3 = I2 + NPERM*NPERM
      I4 = I3 + MAX0(N*NBLOCK,2*NPERM*NPERM)
      I5 = I4 + N*NBLOCK
      I6 = I5 + 2*NPERM + 4
      CALL DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL, NMVEC,
     *  VEC, NBLOCK, WORK(I1), WORK(I2), WORK(I3), WORK(I4),
     *  WORK(I5), WORK(I6), DELTA, SMALL, RARITZ, EPS)
C
  140 IND(1) = NOP
      RETURN
      END
C
C ------------------------------------------------------------------
C
      SUBROUTINE DNWLA(OP, IOVECT, N, NBAND, NVAL, NFIG, NPERM, VAL,
     *   NMVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, P1, P0,
     *   RES, TAU, OTAU, T, ALP, BET, S, P2, BOUND, ATEMP, VTEMP, D,
     *   IND, SMALL, RARITZ, DELTA, EPS, IERR)
C
      INTEGER N, NBAND, NVAL, NFIG, NPERM, NMVEC, NBLOCK, MAXOP, MAXJ,
     *   NOP, IND(1), IERR
      LOGICAL RARITZ, SMALL
      DOUBLE PRECISION VAL(1), VEC(NMVEC,1), P0(N,1), P1(N,1),
     *   P2(N,1), RES(1), TAU(1), OTAU(1), T(NBAND,1),
     *   ALP(NBLOCK,1), BET(NBLOCK,1), BOUND(1), ATEMP(1),
     *   VTEMP(1), D(1), S(MAXJ,1), DELTA, EPS
      EXTERNAL OP, IOVECT
C
C DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE
C ORTHOGONALIZATION.
C
C   NBAND  NBLOCK + 1  THE BAND WIDTH OF T.
C
C   NVAL   THE NUMBER OF DESIRED EIGENVALUES.

C
C   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS
C     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE
C     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL
C     TO THE CURRENT KRYLOV SUBSPACE.
C
C   NOP   THE NUMBER OF CALLS TO OP.
C
C   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS.
C
C   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS.
C
C   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION.
C
C   T   THE BAND MATRIX.
C
C   ALP   THE CURRENT DIAGONAL BLOCK.
C
C   BET   THE CURRENT OFF DIAGONAL BLOCK.
C
C   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY THE BAND
C      EIGENVALUE SOLVER DLAEIG.
C
C   S   EIGENVECTORS OF T.
C
C   SMALL   .TRUE.  IF THE SMALL EIGENVALUES ARE DESIRED.
C
C   RARITZ  RETURNED AS  .TRUE.  IF A FINAL RAYLEIGH-RITZ PROCEDURE
C     IS TO BE DONE.
C
C   DELTA   RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE
C     OF THE MATRIX.  USED IN ESTIMATING THE ACCURACY OF THE
C     COMPUTED EIGENVALUES.
C
C
C INTERNAL VARIABLES USED
C
      INTEGER I, I1, II, J, K, L, M, NG, NGOOD,
     *   NLEFT, NSTART, NTHETA, NUMBER, MIN0, MTEMP
      LOGICAL ENOUGH, TEST
      DOUBLE PRECISION ALPMAX, ALPMIN, ANORM, BETMAX, BETMIN,
     *   EPSRT, PNORM, RNORM, TEMP,
     *   TMAX, TMIN, TOLA, TOLG, UTOL, DABS,
     *   DMAX1, DMIN1, DSQRT, DDOT, DNRM2, TARR(1), DZERO(1)
      EXTERNAL DMVPC, DORTQR, DAXPY, DCOPY, DDOT,
     *   DNRM2, DSCAL, DSWAP, DLAEIG, DLAGER, DLARAN, DVSORT
C
C J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT
C    KRYLOV SUBSPACE.
C
C NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS
C    LIE IN THE CURRENT KRYLOV SUBSPACE).
C
C NLEFT   THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED,
C    I.E.  NLEFT = NVAL - NPERM.
C
C NUMBER = NPERM + NGOOD.
C
C ANORM   AN ESTIMATE OF THE NORM OF THE MATRIX.
C
C EPS   THE RELATIVE MACHINE PRECISION.
C
C UTOL   THE USER TOLERANCE.
C
C TARR  AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO
C       DLAEIG
C
C DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE
C       CONSISTENCY IN CALLS TO DCOPY
C
      DZERO(1) = 0.0D0
      RNORM = 0.0D0
      IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM))
      PNORM = RNORM
      DELTA = 10.D28
      EPSRT = DSQRT(EPS)
      NLEFT = NVAL - NPERM
      NOP = 0
      NUMBER = NPERM
      RARITZ = .FALSE.
      UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE((-FLOAT(NFIG))))
      J = MAXJ
C
C ------------------------------------------------------------------
C
C ANY ITERATION OF THE ALGORITHM BEGINS HERE.
C
   30 DO 50 I=1,NBLOCK
         TEMP = DNRM2(N,P1(1,I),1)
         IF (TEMP.EQ.0D0) CALL DLARAN(N, P1(1,I))
   50 CONTINUE
      IF (NPERM.EQ.0) GO TO 70
      DO 60 I=1,NPERM
         TAU(I) = 1.0D0
         OTAU(I) = 0.0D0
   60 CONTINUE
   70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1)
      CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1)
      CALL DCOPY(J*NBAND, DZERO, 0, T, 1)
      MTEMP = NVAL + 1
      DO 75 I = 1, MTEMP
         CALL DCOPY(J, DZERO, 0, S(1,I), 1)
   75 CONTINUE
      NGOOD = 0
      TMIN = 1.0D28
      TMAX = -1.0D28
      TEST = .TRUE.
      ENOUGH = .FALSE.
      BETMAX = 0.0D0
      J = 0
C
C ------------------------------------------------------------------
C
C THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP.
C
   80 J = J + NBLOCK
C
C THIS IS THE SELECTIVE ORTHOGONALIZATION.
C
      IF (NUMBER.EQ.0) GO TO 110
      DO 100 I=1,NUMBER
         IF (TAU(I).LT.EPSRT) GO TO 100
         TEST = .TRUE.
         TAU(I) = 0.0D0
         IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0
         DO 90 K=1,NBLOCK
            TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1)
            CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1)
C
C THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A
C NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS
C TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST.
C
            IF (DABS(TEMP*BET(K,K)).GT.DBLE(FLOAT(N))*EPSRT*
     *        ANORM .AND. I.GT.NPERM) GO TO 380
   90    CONTINUE
  100 CONTINUE
C
C IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET.
C
  110 IF(.NOT. TEST) GO TO 160
      CALL DORTQR(N, N, NBLOCK, P1, ALP)
      TEST = .FALSE.
      IF(J .EQ. NBLOCK) GO TO 160
      DO 130 I = 1,NBLOCK
         IF(ALP(I,I) .GT. 0.0D0) GO TO 130
         M = J - 2*NBLOCK + I
         L = NBLOCK + 1
         DO 120 K = I,NBLOCK
            BET(I,K) = -BET(I,K)
            T(L,M) = -T(L,M)
            L = L - 1
            M = M + 1
  120    CONTINUE
  130 CONTINUE
C
C THIS IS THE LANCZOS STEP.
C
  160 CALL OP(N, NBLOCK, P1, P2)
      NOP = NOP + 1
      CALL IOVECT(N, NBLOCK, P1, J, 0)
C
C THIS COMPUTES P2=P2-P0*BET(TRANSPOSE)
C
      DO 180 I=1,NBLOCK
         DO 170 K=I,NBLOCK
            CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1)
  170    CONTINUE
  180 CONTINUE
C
C THIS COMPUTES ALP AND P2=P2-P1*ALP.
C
      DO 200 I=1,NBLOCK
         DO 190 K=1,I
            II = I - K + 1
            ALP(II,K) = DDOT(N,P1(1,I),1,P2(1,K),1)
            CALL DAXPY(N, -ALP(II,K), P1(1,I), 1, P2(1,K), 1)
            IF (K.NE.I) CALL DAXPY(N, -ALP(II,K), P1(1,K),
     *        1, P2(1,I), 1)
  190   CONTINUE
  200 CONTINUE
C
C  REORTHOGONALIZATION OF THE SECOND BLOCK
C
      IF(J .NE. NBLOCK) GO TO 220
      DO 215 I=1,NBLOCK
         DO 210 K=1,I
            TEMP = DDOT(N,P1(1,I),1,P2(1,K),1)
            CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1)
            IF (K.NE.I) CALL DAXPY(N, -TEMP, P1(1,K),
     *        1, P2(1,I), 1)
            II = I - K + 1
            ALP(II,K) = ALP(II,K) + TEMP
  210   CONTINUE
  215 CONTINUE
C
C THIS ORTHONORMALIZES THE NEXT BLOCK
C
  220 CALL DORTQR(N, N, NBLOCK, P2, BET)
C
C THIS STORES ALP AND BET IN T.
C
      DO 250 I=1,NBLOCK
         M = J - NBLOCK + I
         DO 230 K=I,NBLOCK
            L = K - I + 1
            T(L,M) = ALP(L,I)
  230    CONTINUE
         DO 240 K=1,I
            L = NBLOCK - I + K + 1
            T(L,M) = BET(K,I)
  240    CONTINUE
  250 CONTINUE
C
C THIS NEGATES T IF SMALL IS FALSE.
C
      IF (SMALL) GO TO 280
      M = J - NBLOCK + 1
      DO 270 I=M,J
         DO 260 K=1,L
            T(K,I) = -T(K,I)
  260    CONTINUE
  270 CONTINUE
C
C THIS SHIFTS THE LANCZOS VECTORS
C
  280 CALL DCOPY(NBLOCK*N, P1, 1, P0, 1)
      CALL DCOPY(NBLOCK*N, P2, 1, P1, 1)
      CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX)
      ANORM = DMAX1(RNORM, TMAX, -TMIN)
      IF (NUMBER.EQ.0) GO TO 305
C
C THIS COMPUTES THE EXTREME EIGENVALUES OF ALP.
C
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
     1   P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      ALPMIN = TARR(1)
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR,
     1  NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      ALPMAX = TARR(1)
C
C THIS COMPUTES ALP = BET(TRANSPOSE)*BET.
C
  305 DO 310 I = 1, NBLOCK
         DO 300 K = 1, I
            L = I - K + 1
            ALP(L,K) = DDOT(NBLOCK-I+1, BET(I,I), NBLOCK, BET(K,I),
     1        NBLOCK)
  300    CONTINUE
  310 CONTINUE
      IF(NUMBER .EQ. 0) GO TO 330
C
C THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET.
C
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
     1  P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, ANORM*ANORM)
      BETMIN = DSQRT(TARR(1))
C
C THIS UPDATES TAU AND OTAU.
C
      DO 320 I=1,NUMBER
         TEMP = (TAU(I)*DMAX1(ALPMAX-VAL(I),VAL(I)-ALPMIN)
     *     +OTAU(I)*BETMAX+EPS*ANORM)/BETMIN
         IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN
         OTAU(I) = TAU(I)
         TAU(I) = TEMP
  320 CONTINUE
C
C THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET.
C
  330 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR,
     1  NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0,
     2  ANORM*ANORM)
      BETMAX = DSQRT(TARR(1))
      IF (J.LE.2*NBLOCK) GO TO 80
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND
C LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK
C IS JUSTIFIED.
C
      TOLG = EPSRT*ANORM
      TOLA = UTOL*RNORM
      IF(MAXJ-J .LT. NBLOCK .OR. (NOP .GE. MAXOP .AND.
     1  NLEFT .NE. 0)) GO TO 390
      GO TO 400

C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO
C SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN
C ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY
C IS AT 400.
C
  380 J = J - NBLOCK
      IERR = -8
  390 IF (NLEFT.EQ.0) RETURN
      TEST = .TRUE.
  400 NTHETA = MIN0(J/2,NLEFT+1)
      CALL DLAEIG(J, NBAND, 1, NTHETA, T, VAL(NUMBER+1),
     1  MAXJ, S, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D)
C
C THIS CHECKS FOR TERMINATION OF A CHECK RUN
C
      IF(NLEFT .NE. 0 .OR. J .LT. 6*NBLOCK) GO TO 410
      IF(VAL(NUMBER+1)-ATEMP(1) .GT. VAL(NPERM) - TOLA) GO TO 790
C
C THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T
C TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED.
C
 410  IF (NTHETA.LE.NLEFT) GO TO 470
      IF (NPERM.EQ.0) GO TO 430
      M = NUMBER + NLEFT + 1
      IF (VAL(M).GE.VAL(NPERM)) GO TO 430
      NPERM = NPERM - 1
      NGOOD = 0
      NUMBER = NPERM
      NLEFT = NLEFT + 1
      GO TO 400
C
C THIS UPDATES DELTA.
C
  430 M = NUMBER + NLEFT + 1
      DELTA = DMIN1(DELTA,VAL(M))
      ENOUGH = .TRUE.
      IF(NLEFT .EQ. 0) GO TO 80
      NTHETA = NLEFT
      VTEMP(NTHETA+1) = 1
C
C ------------------------------------------------------------------
C
C THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL.
C
C THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES.
C
      IF (.NOT.(TEST .OR. ENOUGH)) GO TO 470
      DELTA = DMIN1(DELTA,ANORM)
      PNORM = DMAX1(RNORM,DMAX1(-VAL(NUMBER+1),DELTA))
      TOLA = UTOL*PNORM
      NSTART = 0
      DO 460 I=1,NTHETA
         M = NUMBER + I
         IF (DMIN1(ATEMP(I)*ATEMP(I)/(DELTA-VAL(M)),ATEMP(I))
     *     .GT.TOLA) GO TO 450
         IND(I) = -1
         GO TO 460
C
  450    ENOUGH = .FALSE.
         IF (.NOT.TEST) GO TO 470
         IND(I) = 1
         NSTART = NSTART + 1
  460 CONTINUE
C
C  COPY VALUES OF IND INTO VTEMP
C
      DO 465 I = 1,NTHETA
         VTEMP(I) = DBLE(FLOAT(IND(I)))
  465 CONTINUE
      GO TO 500
C
C THIS CHECKS FOR NEW GOOD VECTORS.
C
  470 NG = 0
      DO 490 I=1,NTHETA
         IF (VTEMP(I).GT.TOLG) GO TO 480
         NG = NG + 1
         VTEMP(I) = -1
         GO TO 490
C
  480    VTEMP(I) = 1
  490 CONTINUE
C
      IF (NG.LE.NGOOD) GO TO 80
      NSTART = NTHETA - NG
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS.
C IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED.
C
  500 TEST = TEST .AND. .NOT.ENOUGH
      NGOOD = NTHETA - NSTART
      NSTART = NSTART + 1
      NTHETA = NTHETA + 1
C
C THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND
C EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR
C FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE
C EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS
C PAUSE.
C
      CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1)
      IF (NSTART.EQ.0) GO TO 580
      IF (NSTART.EQ.NTHETA) GO TO 530
      CALL DVSORT(NTHETA, VTEMP, ATEMP, 1, VAL(NPERM+1), MAXJ,
     *  J, S)
C
C THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING
C VECTORS.
C
  530 IF (.NOT.TEST) NSTART = 0
      IF (.NOT.TEST) GO TO 580
C
C  FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW
C
      TEMP = ATEMP(1)
      DO 535 I = 1, NSTART
         TEMP = DMIN1(TEMP, ATEMP(I))
  535 CONTINUE
      M = NGOOD + 1
      L = NGOOD + MIN0(NSTART,NBLOCK)
      DO 540 I=M,L
         CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1)
  540 CONTINUE
      M = (NSTART-1)/NBLOCK
      IF (M.EQ.0) GO TO 570
      L = NGOOD + NBLOCK
      DO 560 I=1,M
         DO 550 K=1,NBLOCK
            L = L + 1
            IF (L.GT.NTHETA) GO TO 570
            I1 = NGOOD + K
            CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1)
  550    CONTINUE
  560 CONTINUE
  570 NSTART = MIN0(NSTART,NBLOCK)
C
C THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS.
C
  580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. ENOUGH)) GO TO 600
      DO 590 I=1,NGOOD
         M = NPERM + I
         RES(M) = ATEMP(I)
  590 CONTINUE
C
C THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE
C LANCZOS VECTORS.
C
  600 NUMBER = NPERM + NGOOD
      IF (TEST .OR. ENOUGH) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1)
      IF (NGOOD.EQ.0) GO TO 620
      M = NPERM + 1
      DO 610 I=M,NUMBER
         CALL DCOPY(N, DZERO, 0, VEC(1,I), 1)
  610 CONTINUE
  620 DO 670 I=NBLOCK,J,NBLOCK
         CALL IOVECT(N, NBLOCK, P2, I, 1)
         DO 660 K=1,NBLOCK
            M = I - NBLOCK + K
            IF (NSTART.EQ.0) GO TO 640
            DO 630 L=1,NSTART
               I1 = NGOOD + L
               CALL DAXPY(N, S(M,I1), P2(1,K), 1, P1(1,L), 1)
  630       CONTINUE
  640       IF (NGOOD.EQ.0) GO TO 660
            DO 650 L=1,NGOOD
               I1 = L + NPERM
               CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1)
  650       CONTINUE
  660    CONTINUE
  670 CONTINUE
      IF (TEST .OR. ENOUGH) GO TO 690
C
C THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE
C TAU RECURRENCE.
C
      M = NPERM + 1
      DO 680 I=M,NUMBER
         TEMP = 1.0D0/DNRM2(N,VEC(1,I),1)
         CALL DSCAL(N, TEMP, VEC(1,I), 1)
         TAU(I) = 1.0D0
         OTAU(I) = 1.0D0
  680 CONTINUE
C
C  SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG
C
      CALL DCOPY(NTHETA, VAL(NPERM+1), 1, VTEMP, 1)
      CALL DVSORT(NTHETA, VTEMP, ATEMP, 0, TARR, MAXJ, J, S)
      GO TO 80
C
C ------------------------------------------------------------------
C
C THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE
C PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING
C THE PERMANENT VECTORS.
C
  690 IF (NGOOD.EQ.0 .AND. NOP.GE.MAXOP) GO TO 810
      IF (NGOOD.EQ.0) GO TO 30
C
C THIS ORTHONORMALIZES THE VECTORS
C
      CALL DORTQR(NMVEC, N, NPERM+NGOOD, VEC, S)
C
C THIS SORTS THE VALUES AND VECTORS.
C
      IF(NPERM .NE. 0) CALL DVSORT(NPERM+NGOOD, VAL, RES, 0, TEMP,
     *   NMVEC, N, VEC)
      NPERM = NPERM + NGOOD
      NLEFT = NLEFT - NGOOD
      RNORM = DMAX1(-VAL(1),VAL(NPERM))
C
C THIS DECIDES WHERE TO GO NEXT.
C
      IF (NOP.GE.MAXOP .AND. NLEFT.NE.0) GO TO 810
      IF (NLEFT.NE.0) GO TO 30
      IF (VAL(NVAL)-VAL(1).LT.TOLA) GO TO 790
C
C THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED
C TO LOOK FOR UNDISCLOSED MULTIPLICITIES.
C
      M = NPERM - NBLOCK + 1
      IF (M.LE.0) RETURN
      DO 780 I=1,M
         L = I + NBLOCK - 1
         IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30
  780 CONTINUE
C
C THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ
C PROCEDURE IS NEEDED.
C
  790 M = NPERM - NBLOCK
      IF (M.LE.0) RETURN
      DO 800 I=1,M
         L = I + NBLOCK
         IF (VAL(L)-VAL(I).GE.TOLA) GO TO 800
         RARITZ = .TRUE.
         RETURN
  800 CONTINUE
C
      RETURN
C
C ------------------------------------------------------------------
C
C THIS REPORTS THAT MAXOP WAS EXCEEDED.
C
  810 IERR = -2
      GO TO 790
C
      END
C
C ------------------------------------------------------------------
C
      SUBROUTINE DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL,
     *   NMVEC, VEC, NBLOCK, H, HV, P, Q, BOUND, D, DELTA, SMALL,
     *   RARITZ, EPS)
C
      INTEGER N, NPERM, NOP, NMVAL, NMVEC, NBLOCK
      LOGICAL SMALL, RARITZ
      DOUBLE PRECISION VAL(NMVAL,1), VEC(NMVEC,1), H(NPERM,1),
     *   HV(NPERM,1), P(N,1), Q(N,1), BOUND(1), D(1), DELTA, EPS
      EXTERNAL OP, IOVECT
C
C THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS.  BLOCK MATRIX
C VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP.
C
      INTEGER I, J, JJ, K, KK, L, M, MOD
      DOUBLE PRECISION HMIN, HMAX, TEMP, DDOT, DNRM2, DZERO(1)
      EXTERNAL DAXPY, DCOPY, DDOT, DLAGER, DLAEIG
C
C IF RARITZ IS .TRUE.  A FINAL RAYLEIGH-RITZ PROCEDURE IS APPLIED
C TO THE EIGENVECTORS.
C
      DZERO(1) = 0.0D0
      IF (.NOT.RARITZ) GO TO 190
C
C ------------------------------------------------------------------
C
C THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE
C APPROXIMATE EIGENVECTORS.  TEMP = -1 IS USED WHEN SMALL IS
C FALSE TO AVOID HAVING TO RESORT THE EIGENVALUES AND EIGENVECTORS
C COMPUTED BY DLAEIG.
C
      CALL DCOPY(NPERM*NPERM, DZERO, 0, H, 1)
      TEMP = -1.0D0
      IF (SMALL) TEMP = 1.0D0
      M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 40
      DO 10 I=1,M
         CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1)
   10 CONTINUE
      CALL IOVECT(N, M, P, M, 0)
      CALL OP(N, M, P, Q)
      NOP = NOP + 1
      DO 30 I=1,M
         DO 20 J=I,NPERM
            JJ = J - I + 1
            H(JJ,I) = TEMP*DDOT(N,VEC(1,J),1,Q(1,I),1)
   20    CONTINUE
   30 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 90
   40 M = M + NBLOCK
      DO 80 I=M,NPERM,NBLOCK
         DO 50 J=1,NBLOCK
            L = I - NBLOCK + J
            CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1)
   50    CONTINUE
         CALL IOVECT(N, NBLOCK, P, I, 0)
         CALL OP(N, NBLOCK, P, Q)
         NOP = NOP + 1
         DO 70 J=1,NBLOCK
            L = I - NBLOCK + J
            DO 60 K=L,NPERM
               KK = K - L + 1
               H(KK,L) = TEMP*DDOT(N,VEC(1,K),1,Q(1,J),1)
   60       CONTINUE
   70      CONTINUE
   80 CONTINUE
C
C THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H.
C
   90 HMIN = H(1,1)
      HMAX = H(1,1)
      CALL DLAGER(NPERM, NPERM, 1, H, HMIN, HMAX)
      CALL DLAEIG(NPERM, NPERM, 1, NPERM, H, VAL, NPERM,
     1  HV, BOUND, P, D, Q, EPS, HMIN, HMAX)
C
C THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF
C Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H.
C
      DO 120 I=1,NPERM
         CALL DCOPY(N, DZERO, 0, VEC(1,I), 1)
  120 CONTINUE
      M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 150
      CALL IOVECT(N, M, P, M, 1)
      DO 140 I=1,M
         DO 130 J=1,NPERM
            CALL DAXPY(N, HV(I,J), P(1,I), 1, VEC(1,J), 1)
  130    CONTINUE
  140 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 190
  150 M = M + NBLOCK
      DO 180 I=M,NPERM,NBLOCK
         CALL IOVECT(N, NBLOCK, P, I, 1)
         DO 170 J=1,NBLOCK
            L = I - NBLOCK + J
            DO 160 K=1,NPERM
               CALL DAXPY(N, HV(L,K), P(1,J), 1, VEC(1,K), 1)
  160       CONTINUE
  170    CONTINUE
  180 CONTINUE
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1))
C AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS.
C
  190 IF (.NOT.SMALL) DELTA = -DELTA
      M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 220
      DO 200 I=1,M
         CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1)
  200 CONTINUE
      CALL OP(N, M, P, Q)
      NOP = NOP + 1
      DO 210 I=1,M
         VAL(I,1) = DDOT(N,P(1,I),1,Q(1,I),1)
         CALL DAXPY(N, -VAL(I,1), P(1,I), 1, Q(1,I), 1)
         VAL(I,2) = DNRM2(N,Q(1,I),1)
  210 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 260
  220 M = M + 1
      DO 250 I=M,NPERM,NBLOCK
         DO 230 J=1,NBLOCK
            L = I - 1 + J
            CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1)
  230    CONTINUE
         CALL OP(N, NBLOCK, P, Q)
         NOP = NOP + 1
         DO 240 J=1,NBLOCK
            L = I - 1 + J
            VAL(L,1) = DDOT(N,P(1,J),1,Q(1,J),1)
            CALL DAXPY(N, -VAL(L,1), P(1,J), 1, Q(1,J), 1)
            VAL(L,2) = DNRM2(N,Q(1,J),1)
  240    CONTINUE
  250 CONTINUE
C
C THIS COMPUTES THE ACCURACY ESTIMATES.  FOR CONSISTENCY WITH DILASO
C A DO LOOP IS NOT USED.
C
  260 I = 0
  270 I = I + 1
      IF (I.GT.NPERM) RETURN
      TEMP = DELTA - VAL(I,1)
      IF (.NOT.SMALL) TEMP = -TEMP
      VAL(I,4) = 0.0D0
      IF (TEMP.GT.0.0D0) VAL(I,4) = VAL(I,2)/TEMP
      VAL(I,3) = VAL(I,4)*VAL(I,2)
      GO TO 270
C
      END


      DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
      INTEGER I, INCX, J, N, NEXT, NN
      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
      DATA   ZERO, ONE /0.0D0, 1.0D0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C
      IF(N .GT. 0) GO TO 10
         DNRM2  = ZERO
         GO TO 300
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N * INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20    GO TO NEXT,(30, 50, 70, 110)
   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( DX(I) .EQ. ZERO) GO TO 200
      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM / DX(I)) / DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / DX(I))**2
         XMAX = DABS(DX(I))
         GO TO 200
C
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
   75 SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
   85 HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 95 J =I,NN,INCX
      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + DX(J)**2
      DNRM2 = DSQRT( SUM )
      GO TO 300
C
  200 CONTINUE
      I = I + INCX
      IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      DNRM2 = XMAX * DSQRT(SUM)
  300 CONTINUE
      RETURN
      END
C
C ------------------------------------------------------------------
C
      SUBROUTINE DORTQR(NZ, N, NBLOCK, Z, B)
C
      INTEGER NZ, N, NBLOCK
      DOUBLE PRECISION Z(NZ,1), B(NBLOCK,1)
C
C THIS SUBROUTINE COMPUTES THE QR FACTORIZATION OF THE N X NBLOCK
C MATRIX Z.  Q IS FORMED IN PLACE AND RETURNED IN Z.  R IS
C RETURNED IN B.
C
      INTEGER I, J, K, LENGTH, M
      DOUBLE PRECISION SIGMA, TAU, TEMP, DDOT, DNRM2, DSIGN
      EXTERNAL DAXPY, DDOT, DNRM2, DSCAL
C
C THIS SECTION REDUCES Z TO TRIANGULAR FORM.
C
      DO 30 I=1,NBLOCK
C
C THIS FORMS THE ITH REFLECTION.
C
         LENGTH = N - I + 1
         SIGMA = DSIGN(DNRM2(LENGTH,Z(I,I),1),Z(I,I))
         B(I,I) = -SIGMA
         Z(I,I) = Z(I,I) + SIGMA
         TAU = SIGMA*Z(I,I)
         IF (I.EQ.NBLOCK) GO TO 30
         J = I + 1
C
C THIS APPLIES THE ROTATION TO THE REST OF THE COLUMNS.
C
         DO 20 K=J,NBLOCK
            IF (TAU.EQ.0.0D0) GO TO 10
            TEMP = -DDOT(LENGTH,Z(I,I),1,Z(I,K),1)/TAU
            CALL DAXPY(LENGTH, TEMP, Z(I,I), 1, Z(I,K), 1)
   10       B(I,K) = Z(I,K)
            Z(I,K) = 0.0D0
   20    CONTINUE
   30 CONTINUE
C
C THIS ACCUMULATES THE REFLECTIONS IN REVERSE ORDER.
C
      DO 70 M=1,NBLOCK
C
C THIS RECREATES THE ITH = NBLOCK-M+1)TH REFLECTION.
C
         I = NBLOCK + 1 - M
         SIGMA = -B(I,I)
         TAU = Z(I,I)*SIGMA
         IF (TAU.EQ.0.0D0) GO TO 60
         LENGTH = N - NBLOCK + M
         IF (I.EQ.NBLOCK) GO TO 50
         J = I + 1
C
C THIS APPLIES IT TO THE LATER COLUMNS.
C
         DO 40 K=J,NBLOCK
            TEMP = -DDOT(LENGTH,Z(I,I),1,Z(I,K),1)/TAU
            CALL DAXPY(LENGTH, TEMP, Z(I,I), 1, Z(I,K), 1)
   40    CONTINUE
   50    CALL DSCAL(LENGTH, -1.0D0/SIGMA, Z(I,I), 1)
   60    Z(I,I) = 1.0D0 + Z(I,I)
   70 CONTINUE
      RETURN
      END
      SUBROUTINE  DSCAL(N,DA,DX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DA,DX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        DX(I) = DA*DX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DX(I) = DA*DX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DX(I) = DA*DX(I)
        DX(I + 1) = DA*DX(I + 1)
        DX(I + 2) = DA*DX(I + 2)
        DX(I + 3) = DA*DX(I + 3)
        DX(I + 4) = DA*DX(I + 4)
   50 CONTINUE
      RETURN
      END
C
C-------------------------------------------------------------------
C
      SUBROUTINE DVSORT(NUM, VAL, RES, IFLAG, V, NMVEC, N, VEC)
      INTEGER NUM, IFLAG, NMVEC, N
      DOUBLE PRECISION VAL(1), RES(1), V(1), VEC(NMVEC,1)
C
C  THIS SUBROUTINE SORTS THE EIGENVALUES (VAL) IN ASCENDING ORDER
C  WHILE CONCURRENTLY SWAPPING THE RESIDUALS AND VECTORS.
      INTEGER I, K, M
      DOUBLE PRECISION TEMP
      IF(NUM .LE. 1) RETURN
      DO 20 I = 2, NUM
         M = NUM - I + 1
         DO 10 K = 1, M
            IF(VAL(K) .LE. VAL(K+1)) GO TO 10
            TEMP = VAL(K)
            VAL(K) = VAL(K+1)
            VAL(K+1) = TEMP
            TEMP = RES(K)
            RES(K) = RES(K+1)
            RES(K+1) = TEMP
            CALL DSWAP(N, VEC(1,K), 1, VEC(1,K+1), 1)
            IF(IFLAG .EQ. 0) GO TO 10
            TEMP = V(K)
            V(K) = V(K+1)
            V(K+1) = TEMP
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DA
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF (DA .EQ. 0.0D0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DY(IY) + DA*DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        DY(I) = DY(I) + DA*DX(I)
        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
C
C     COPIES A VECTOR, X, TO A VECTOR, Y.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        DY(I) = DX(I)
        DY(I + 1) = DX(I + 1)
        DY(I + 2) = DX(I + 2)
        DY(I + 3) = DX(I + 3)
        DY(I + 4) = DX(I + 4)
        DY(I + 5) = DX(I + 5)
        DY(I + 6) = DX(I + 6)
   50 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      DDOT = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DTEMP + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      DDOT = DTEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
   60 DDOT = DTEMP
      RETURN
      END
C
C
      SUBROUTINE DLAEIG(N, NBAND, NL, NR, A, EIGVAL, LDE,
     1   EIGVEC, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
C
C  THIS IS A SPECIALIZED VERSION OF THE SUBROUTINE BNDEIG TAILORED
C  SPECIFICALLY FOR USE BY THE LASO PACKAGE.
C
      INTEGER N, NBAND, NL, NR, LDE
      DOUBLE PRECISION A(NBAND,1), EIGVAL(1),
     1   EIGVEC(LDE,1), BOUND(2,1), ATEMP(1), D(1), VTEMP(1),
     2   EPS, TMIN, TMAX
C
C  LOCAL VARIABLES
C
      INTEGER I, M, NVAL
      DOUBLE PRECISION ARTOL, ATOL
C
C  FUNCTIONS CALLED
C
      DOUBLE PRECISION DMAX1
C
C  SUBROUTINES CALLED
C
C     DLABCM, DLABFC, DLAGER, DCOPY
C
C  SET PARAMETERS
C
      ATOL = DBLE(FLOAT(N))*EPS*DMAX1(TMAX,-TMIN)
      ARTOL = ATOL/DSQRT(EPS)
      NVAL = NR - NL + 1
C
C   CHECK FOR SPECIAL CASE OF N = 1
C
      IF(N .NE. 1) GO TO 30
      EIGVAL(1) = A(1,1)
      EIGVEC(1,1) = 1.0D0
      RETURN
C
C   SET UP INITIAL EIGENVALUE BOUNDS
C
   30 M = NVAL + 1
      DO 50 I = 2, M
         BOUND(1,I) = TMIN
         BOUND(2,I) = TMAX
   50 CONTINUE
      BOUND(2,1) = TMAX
      BOUND(1,NVAL + 2) = TMIN
      IF(NL .EQ. 1) BOUND(2,1) = TMIN
      IF(NR .EQ. N) BOUND(1,NVAL + 2) = TMAX
C
   60 CALL DLABCM(N, NBAND, NL, NR, A, EIGVAL, LDE,
     1  EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP)
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE DLAGER(N, NBAND, NSTART, A, TMIN, TMAX)
C
C  THIS SUBROUTINE COMPUTES BOUNDS ON THE SPECTRUM OF A BY
C  EXAMINING THE GERSCHGORIN CIRCLES. ONLY THE NEWLY CREATED
C  CIRCLES ARE EXAMINED
C
C  FORMAL PARAMETERS
C
      INTEGER N, NBAND, NSTART
      DOUBLE PRECISION A(NBAND,1), TMIN, TMAX
C
C  LOCAL VARIABLES
C
      INTEGER I, K, L, M
      DOUBLE PRECISION TEMP
C
C  FUNCTIONS CALLED
C
      INTEGER MIN0
      DOUBLE PRECISION DMIN1, DMAX1
C
      DO 50 K = NSTART, N
         TEMP = 0.0D0
         DO 10 I = 2, NBAND
            TEMP = TEMP + DABS(A(I,K))
   10    CONTINUE
   20    L = MIN0(K,NBAND)
         IF(L .EQ. 1) GO TO 40
         DO 30 I = 2, L
            M = K - I + 1
            TEMP = TEMP + DABS(A(I,M))
   30    CONTINUE
   40    TMIN = DMIN1(TMIN, A(1,K)-TEMP)
         TMAX = DMAX1(TMAX, A(1,K)+TEMP)
   50 CONTINUE
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE DLARAN(N, X)
C
C  THIS SUBROUTINE SETS THE VECTOR X TO RANDOM NUMBERS
C
C  FORMAL PARAMETERS
C
      INTEGER N
      DOUBLE PRECISION X(N)
C
C  LOCAL VARIABLES
C
      INTEGER I, IURAND
C
C  FUNCTIONS CALLED
C
      REAL URAND
      DOUBLE PRECISION DBLE
C
C  SUBROUTINES CALLED
C
C     NONE
C
C  INITIALIZE SEED
C
c     julie 14.3.92- temp or maybe not
c      common /random/iurand
      DATA IURAND /0/
C
      DO 10 I = 1, N
         X(I) = DBLE(URAND(IURAND)) - 0.5D0
   10 CONTINUE
      RETURN
      END
C
C ------------------------------------------------------------------
C
      SUBROUTINE DMVPC(NBLOCK, BET, MAXJ, J, S, NUMBER, RESNRM,
     *     ORTHCF, RV)
C
      INTEGER NBLOCK, MAXJ, J, NUMBER
      DOUBLE PRECISION BET(NBLOCK,1), S(MAXJ,1), RESNRM(1),
     *     ORTHCF(1), RV(1)
C
C THIS SUBROUTINE COMPUTES THE NORM AND THE SMALLEST ELEMENT
C (IN ABSOLUTE VALUE) OF THE VECTOR BET*SJI, WHERE SJI
C IS AN NBLOCK VECTOR OF THE LAST NBLOCK ELEMENTS OF THE ITH
C EIGENVECTOR OF T.  THESE QUANTITIES ARE THE RESIDUAL NORM
C AND THE ORTHOGONALITY COEFFICIENT RESPECTIVELY FOR THE
C CORRESPONDING RITZ PAIR.  THE ORTHOGONALITY COEFFICIENT IS
C NORMALIZED TO ACCOUNT FOR THE LOCAL REORTHOGONALIZATION.
C
      INTEGER I, K, M
      DOUBLE PRECISION DDOT, DNRM2, DABS, DMIN1
C
      M = J - NBLOCK + 1
      DO 20 I=1,NUMBER
         DO 10 K=1,NBLOCK
            RV(K) = DDOT(NBLOCK,S(M,I),1,BET(K,1),NBLOCK)
            IF (K.EQ.1) ORTHCF(I) = DABS(RV(K))
            ORTHCF(I) = DMIN1(ORTHCF(I),DABS(RV(K)))
   10    CONTINUE
         RESNRM(I) = DNRM2(NBLOCK,RV,1)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE  DSWAP (N,DX,INCX,DY,INCY)
C
C     INTERCHANGES TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DX(IX)
        DX(IX) = DY(IY)
        DY(IY) = DTEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C       CLEAN-UP LOOP
C
   20 M = MOD(N,3)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
   30 CONTINUE
      IF( N .LT. 3 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,3
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
        DTEMP = DX(I + 1)
        DX(I + 1) = DY(I + 1)
        DY(I + 1) = DTEMP
        DTEMP = DX(I + 2)
        DX(I + 2) = DY(I + 2)
        DY(I + 2) = DTEMP
   50 CONTINUE
      RETURN
      END
      REAL FUNCTION URAND(IY)
      INTEGER  IY
C
C      URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED  ON  THEORY  AND
C  SUGGESTIONS  GIVEN  IN  D.E. KNUTH (1969),  VOL  2.   THE INTEGER  IY
C  SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL
C  TO URAND.  THE CALLING PROGRAM SHOULD  NOT  ALTER  THE  VALUE  OF  IY
C  BETWEEN  SUBSEQUENT CALLS TO URAND.  VALUES OF URAND WILL BE RETURNED
C  IN THE INTERVAL (0,1).
C
      INTEGER  IA,IC,ITWO,M2,M,MIC
      DOUBLE PRECISION  HALFM
      REAL  S
      DOUBLE PRECISION  DATAN,DSQRT
      DATA M2/0/,ITWO/2/
      IF (M2 .NE. 0) GO TO 20
C
C  IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH
C
      M = 1
   10 M2 = M
      M = ITWO*M2
      IF (M .GT. M2) GO TO 10
      HALFM = M2
C
C  COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
C
      IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5
      IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1
      MIC = (M2 - IC) + M2
C
C  S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
C
      S = 0.5/HALFM
C
C  COMPUTE NEXT RANDOM NUMBER
C
   20 IY = IY*IA
C
C  THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW
C  INTEGER OVERFLOW ON ADDITION
C
      IF (IY .GT. MIC) IY = (IY - M2) - M2
C
      IY = IY + IC
C
C  THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE
C  WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION
C
      IF (IY/2 .GT. M2) IY = (IY - M2) - M2
C
C  THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER
C  OVERFLOW AFFECTS THE SIGN BIT
C
      IF (IY .LT. 0) IY = (IY + M2) + M2
      URAND = FLOAT(IY)*S
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE DLABCM(N, NBAND, NL, NR, A, EIGVAL,
     1  LDE, EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP)
C
C  THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES
C  FOR THE BNDEIG PACKAGE.  EIGENVALUES ARE COMPUTED BY
C  A MODIFIED RAYLEIGH QUOTIENT ITERATION.  THE EIGENVALUE COUNT
C  OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE
C  THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO
C  INSURE CONVERGENCE TO THE DESIRED EIGENVALUES.
C
C  FORMAL PARAMETERS.
C
      INTEGER N, NBAND, NL, NR, LDE
      DOUBLE PRECISION A(NBAND,1), EIGVAL(1),
     1  EIGVEC(LDE,1), ATOL, ARTOL, BOUND(2,1), ATEMP(1),
     2  D(1), VTEMP(1)
C
C
C  LOCAL VARIABLES
C
      LOGICAL FLAG
      INTEGER I, J, L, M, NUML, NUMVEC, NVAL
      DOUBLE PRECISION ERRB, GAP, RESID, RQ, SIGMA, VNORM
C
C
C  FUNCTIONS CALLED
C
      INTEGER MIN0
      DOUBLE PRECISION DMAX1, DMIN1, DDOT, DNRM2
C
C  SUBROUTINES CALLED
C
C     DLABAX, DLABFC, DLARAN, DAXPY, DCOPY, DSCAL
C
C  REPLACE ZERO VECTORS BY RANDOM
C
      NVAL = NR - NL + 1
      FLAG = .FALSE.
      DO 5 I = 1, NVAL
         IF(DDOT(N, EIGVEC(1,I), 1, EIGVEC(1,I), 1) .EQ. 0.0D0)
     1      CALL DLARAN(N,EIGVEC(1,I))
    5 CONTINUE
C
C  LOOP OVER EIGENVALUES
C
      SIGMA = BOUND(2,NVAL+1)
      DO 400 J = 1, NVAL
         NUML = J
C
C  PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT
C
   10    CALL DLABAX(N, NBAND, A, EIGVEC(1,J), VTEMP)
         VNORM = DNRM2(N, VTEMP, 1)
         IF(VNORM .EQ. 0.0D0) GO TO 20
         CALL DSCAL(N, 1.0D0/VNORM, VTEMP, 1)
         CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
         CALL DAXPY(N, -SIGMA, EIGVEC(1,J), 1, VTEMP, 1)
C
C  LOOP OVER SHIFTS
C
C  COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE
C
   20       VNORM = DNRM2(N, EIGVEC(1,J), 1)
            IF(VNORM .NE. 0.0D0) GO TO 30
            CALL DLARAN(N, EIGVEC(1,J))
            GO TO 10
C
   30       RQ = SIGMA + DDOT(N, EIGVEC(1,J), 1, VTEMP, 1)
     1        /VNORM/VNORM
            CALL DAXPY(N, SIGMA-RQ, EIGVEC(1,J), 1, VTEMP, 1)
            RESID = DMAX1(ATOL, DNRM2(N, VTEMP, 1)/VNORM)
            CALL DSCAL(N, 1.0/VNORM, EIGVEC(1,J), 1)
C
C  ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH
C
            IF(BOUND(2,J+1) - BOUND(1,J+1) .LT. 3.0D0*ATOL) GO TO 300
C
C  COMPUTE MINIMAL ERROR BOUND
C
            ERRB = RESID
            GAP = DMIN1(BOUND(1,J+2) - RQ, RQ - BOUND(2,J))
            IF(GAP .GT. RESID) ERRB = DMAX1(ATOL, RESID*RESID/GAP)
C
C  TENTATIVE NEW SHIFT
C
            SIGMA = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
C
C  CHECK FOR TERMINALTION
C
            IF(RESID .GT. 2.0D0*ATOL) GO TO 40
            IF(RQ - ERRB .GT. BOUND(2,J) .AND.
     1        RQ + ERRB .LT. BOUND(1,J+2)) GO TO 310
C
C  RQ IS TO THE LEFT OF THE INTERVAL
C
   40       IF(RQ .GE. BOUND(1,J+1)) GO TO 50
            IF(RQ - ERRB .GT. BOUND(2,J)) GO TO 100
            IF(RQ + ERRB .LT. BOUND(1,J+1)) CALL DLARAN(N,EIGVEC(1,J))
            GO TO 200
C
C  RQ IS TO THE RIGHT OF THE INTERVAL
C
   50       IF(RQ .LE. BOUND(2,J+1)) GO TO 100
            IF(RQ + ERRB .LT. BOUND(1,J+2)) GO TO 100
C
C  SAVE THE REJECTED VECTOR IF INDICATED
C
            IF(RQ - ERRB .LE. BOUND(2,J+1)) GO TO 200
            DO 60 I = J, NVAL
               IF(BOUND(2,I+1) .GT. RQ) GO TO 70
   60       CONTINUE
            GO TO 80
C
   70       CALL DCOPY(N, EIGVEC(1,J), 1, EIGVEC(1,I), 1)
C
   80       CALL DLARAN(N, EIGVEC(1,J))
            GO TO 200
C
C  PERTURB RQ TOWARD THE MIDDLE
C
  100       IF(SIGMA .LT. RQ) SIGMA = DMAX1(SIGMA, RQ-ERRB)
            IF(SIGMA .GT. RQ) SIGMA = DMIN1(SIGMA, RQ+ERRB)
C
C  FACTOR AND SOLVE
C
  200       DO 210 I = J, NVAL
               IF(SIGMA .LT. BOUND(1,I+1)) GO TO 220
  210       CONTINUE
            I = NVAL + 1
  220       NUMVEC = I - J
            NUMVEC = MIN0(NUMVEC, NBAND + 2)
            IF(RESID .LT. ARTOL) NUMVEC = MIN0(1,NUMVEC)
            CALL DCOPY(N, EIGVEC(1,J), 1, VTEMP, 1)
            CALL DLABFC(N, NBAND, A, SIGMA, NUMVEC, LDE,
     1        EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
C
C  PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW
C
            IF(NUMVEC .EQ. 1) GO TO 227
            L = NUMVEC - 1
            DO 225 I = 1,L
               M = J + I
               CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,M), 1)
  225       CONTINUE
C
C  UPDATE INTERVALS
C
  227       NUML = NUML - NL + 1
            IF(NUML .GE. 0) BOUND(2,1) = DMIN1(BOUND(2,1), SIGMA)
            DO 230 I = J, NVAL
               IF(SIGMA .LT. BOUND(1,I+1)) GO TO 20
               IF(NUML .LT. I) BOUND(1,I+1) = SIGMA
               IF(NUML .GE. I) BOUND(2,I+1) = SIGMA
  230       CONTINUE
            IF(NUML .LT. NVAL + 1) BOUND(1,NVAL+2) = DMAX1(SIGMA,
     1        BOUND(1,NVAL+2))
            GO TO 20
C
C  ACCEPT AN EIGENPAIR
C
  300    CALL DLARAN(N, EIGVEC(1,J))
         FLAG = .TRUE.
         GO TO 310
C
  305    FLAG = .FALSE.
         RQ = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
         CALL DLABFC(N, NBAND, A, RQ, NUMVEC, LDE,
     1    EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
         VNORM = DNRM2(N, EIGVEC(1,J), 1)
         IF(VNORM .NE. 0.0) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C  ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES
C
  310    EIGVAL(J) = RQ
         IF(J .EQ. 1) GO TO 330
         M = J - 1
         DO 320 I = 1, M
            CALL DAXPY(N, -DDOT(N,EIGVEC(1,I),1,EIGVEC(1,J),1),
     1        EIGVEC(1,I), 1, EIGVEC(1,J), 1)
  320    CONTINUE
  330    VNORM = DNRM2(N, EIGVEC(1,J), 1)
         IF(VNORM .EQ. 0.0D0) GO TO 305
         CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE
C
         IF(FLAG) GO TO 305
         IF(J .EQ. NVAL) RETURN
         M = J + 1
         DO 340 I = M, NVAL
            CALL DAXPY(N, -DDOT(N,EIGVEC(1,J),1,EIGVEC(1,I),1),
     1        EIGVEC(1,J), 1, EIGVEC(1,I), 1)
  340    CONTINUE
  400 CONTINUE
      RETURN
C
  500 CONTINUE
      END
C
C***********************************************************************
C
      SUBROUTINE DLABFC(N, NBAND, A, SIGMA, NUMBER, LDE,
     1  EIGVEC, NUML, LDAD, ATEMP, D, ATOL)
C
C  THIS SUBROUTINE FACTORS (A-SIGMA*I) WHERE A IS A GIVEN BAND
C  MATRIX AND SIGMA IS AN INPUT PARAMETER.  IT ALSO SOLVES ZERO
C  OR MORE SYSTEMS OF LINEAR EQUATIONS.  IT RETURNS THE NUMBER
C  OF EIGENVALUES OF A LESS THAN SIGMA BY COUNTING THE STURM
C  SEQUENCE DURING THE FACTORIZATION.  TO OBTAIN THE STURM
C  SEQUENCE COUNT WHILE ALLOWING NON-SYMMETRIC PIVOTING FOR
C  STABILITY, THE CODE USES A GUPTA'S MULTIPLE PIVOTING
C  ALGORITHM.
C
C  FORMAL PARAMETERS
C
      INTEGER N, NBAND, NUMBER, LDE, NUML, LDAD
      DOUBLE PRECISION A(NBAND,1), SIGMA, EIGVEC(LDE,1),
     1  ATEMP(LDAD,1), D(LDAD,1), ATOL
C
C  LOCAL VARIABLES
C
      INTEGER I, J, K, KK, L, LA, LD, LPM, M, NB1
      DOUBLE PRECISION ZERO(1)
C
C  FUNCTIONS CALLED
C
      INTEGER MIN0
      DOUBLE PRECISION DABS
C
C  SUBROUTINES CALLED
C
C     DAXPY, DCOPY, DSWAP
C
C
C  INITIALIZE
C
      ZERO(1) = 0.0D0
      NB1 = NBAND - 1
      NUML = 0
      CALL DCOPY(LDAD*NBAND, ZERO, 0, D, 1)
C
C   LOOP OVER COLUMNS OF A
C
      DO 100 K = 1, N
C
C   ADD A COLUMN OF A TO D
C
         D(NBAND, NBAND) = A(1,K) - SIGMA
         M = MIN0(K, NBAND) - 1
         IF(M .EQ. 0) GO TO 20
         DO 10 I = 1, M
            LA = K - I
            LD = NBAND - I
            D(LD,NBAND) = A(I+1, LA)
   10    CONTINUE
C
   20    M = MIN0(N-K, NB1)
         IF(M .EQ. 0) GO TO 40
         DO 30 I = 1, M
            LD = NBAND + I
            D(LD, NBAND) = A(I+1, K)
   30    CONTINUE
C
C   TERMINATE
C
   40    LPM = 1
         IF(NB1 .EQ. 0) GO TO 70
         DO 60 I = 1, NB1
            L = K - NBAND + I
            IF(D(I,NBAND) .EQ. 0.0D0) GO TO 60
            IF(DABS(D(I,I)) .GE. DABS(D(I,NBAND))) GO TO 50
            IF((D(I,NBAND) .LT. 0.0D0 .AND. D(I,I) .LT. 0.0D0)
     1        .OR. (D(I,NBAND) .GT. 0.0D0 .AND. D(I,I) .GE. 0.0D0))
     2        LPM = -LPM
            CALL DSWAP(LDAD-I+1, D(I,I), 1, D(I,NBAND), 1)
            CALL DSWAP(NUMBER, EIGVEC(L,1), LDE, EIGVEC(K,1), LDE)
   50       CALL DAXPY(LDAD-I, -D(I,NBAND)/D(I,I), D(I+1,I), 1,
     1        D(I+1,NBAND), 1)
            CALL DAXPY(NUMBER, -D(I,NBAND)/D(I,I), EIGVEC(L,1),
     1        LDE, EIGVEC(K,1), LDE)
   60    CONTINUE
C
C  UPDATE STURM SEQUENCE COUNT
C
   70    IF(D(NBAND,NBAND) .LT. 0.0D0) LPM = -LPM
         IF(LPM .LT. 0) NUML = NUML + 1
         IF(K .EQ. N) GO TO 110
C
C   COPY FIRST COLUMN OF D INTO ATEMP
         IF(K .LT. NBAND) GO TO 80
         L = K - NB1
         CALL DCOPY(LDAD, D, 1, ATEMP(1,L), 1)
C
C   SHIFT THE COLUMNS OF D OVER AND UP
C
         IF(NB1 .EQ. 0) GO TO 100
   80    DO 90 I = 1, NB1
            CALL DCOPY(LDAD-I, D(I+1,I+1), 1, D(I,I), 1)
            D(LDAD,I) = 0.0D0
   90    CONTINUE
  100 CONTINUE
C
C  TRANSFER D TO ATEMP
C
  110 DO 120 I = 1, NBAND
         L = N - NBAND + I
         CALL DCOPY(NBAND-I+1, D(I,I), 1, ATEMP(1,L), 1)
  120 CONTINUE
C
C   BACK SUBSTITUTION
C
      IF(NUMBER .EQ. 0) RETURN
      DO 160 KK = 1, N
         K = N - KK + 1
         IF(DABS(ATEMP(1,K)) .LE. ATOL)
     1     ATEMP(1,K) = DSIGN(ATOL,ATEMP(1,K))
C
  130    DO 150 I = 1, NUMBER
            EIGVEC(K,I) = EIGVEC(K,I)/ATEMP(1,K)
            M = MIN0(LDAD, K) - 1
            IF(M .EQ. 0) GO TO 150
            DO 140 J = 1, M
                L = K - J
                EIGVEC(L,I) = EIGVEC(L,I) - ATEMP(J+1,L)*EIGVEC(K,I)
  140       CONTINUE
  150    CONTINUE
  160 CONTINUE
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE DLABAX(N, NBAND, A, X, Y)
C
C  THIS SUBROUTINE SETS Y = A*X
C  WHERE X AND Y ARE VECTORS OF LENGTH N
C  AND A IS AN  N X NBAND  SYMMETRIC BAND MATRIX
C
C  FORMAL PARAMETERS
C
      INTEGER N, NBAND
      DOUBLE PRECISION A(NBAND,1), X(1), Y(1)
C
C  LOCAL VARIABLES
C
      INTEGER I, K, L, M
      DOUBLE PRECISION ZERO(1)
C
C  FUNCTIONS CALLED
C
      INTEGER MIN0
C
C  SUBROUTINES CALLED
C
C     DCOPY
C
      ZERO(1) = 0.0D0
      CALL DCOPY(N, ZERO, 0, Y, 1)
      DO 20 K = 1, N
         Y(K) = Y(K) + A(1,K)*X(K)
         M = MIN0(N-K+1, NBAND)
         IF(M .LT. 2) GO TO 20
         DO 10 I = 2, M
            L = K + I - 1
            Y(L) = Y(L) + A(I,K)*X(K)
            Y(K) = Y(K) + A(I,K)*X(L)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END