C-----LSSUB-------(SINGULAR VALUES AND VECTORS)-------------------------
C     ACCORDING TO PFORT THESE SUBROUTINES ARE PORTABLE
C
C
C     SUBROUTINES    BISEC, INVERR, TNORM, LUMP, ISOEV, PRTEST, AND
C                    INVERM ARE USED WITH LANCZOS SINGULAR VALUE
C                    PROGRAM LSVAL.  STURMI, INVERM, LBISEC, TNORM
C                    ARE USED WITH THE LANCZOS SINGULAR VECTOR
C                    PROGRAM LSVEC.
C
C
C-----COMPUTE T-EIGENVALUES BY BISECTION--------------------------------
C
      SUBROUTINE BISEC(BETA,BETA2,VB,VS,LBD,UBD,EPS,TTOL,MP,
     1 NINT,MEV,NDIS,IC,IWRITE)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  BETA(1),BETA2(1),VB(1),VS(1)
      DOUBLE PRECISION  LBD(1),UBD(1),EPS,EPT,EP0,EP1,TEMP,TTOL
      DOUBLE PRECISION  ZERO,ONE,HALF,YU,YV,LB,UB,XL,XU,X1,X0,XS,BETAM
      INTEGER  MP(1),IDEF(10)
      DOUBLE PRECISION  DABS, DSQRT, DMAX1, DMIN1, DFLOAT
C-----------------------------------------------------------------------
C     COMPUTES EIGENVALUES OF T(1,MEV) BY LOOPING INTERNALLY ON THE
C     USER-SPECIFIED INTERVALS,  (LB(J),UB(J)), J = 1,NINT.  INTERVALS
C     ARE TREATED AS OPEN ON THE LEFT AND CLOSED ON THE RIGHT.
C     THE BISEC SUBROUTINE SIMULTANEOUSLY LABELS SPURIOUS T-EIGENVALUES
C     AND DETERMINES THE T-MULTIPLICITIES OF EACH GOOD T-EIGENVALUE.
C     SPURIOUS T-EIGENVALUES ARE LABELLED BY A T-MULTIPLICITY = 0.
C     ANY T-EIGENVALUE WITH A T-MULTIPLICITY >= 1 IS 'GOOD'.
C
C     IF IWRITE = 0 THEN MOST OF THE WRITES TO FILE 6 ARE NOT
C     ACTIVATED.
C
C     NOTE THAT PROGRAM ASSUMES THAT NO MORE THAN MMAX/2 T-EIGENVALUES
C     OF T(1,MEV) ARE TO BE COMPUTED IN ANY ONE OF THE SUBINTERVALS
C     CONSIDERED, WHERE MMAX = DIMENSION OF VB SPECIFIED BY THE USER
C     IN THE MAIN PROGRAM LEVAL.
C
C     ON ENTRY
C     BETA2(J) IS SET = BETA(J)*BETA(J).  THE STORAGE FOR BETA2 COULD
C     BE ELIMINATED BY RECOMPUTING THE BETA(J)**2 FOR EACH STURM
C     SEQUENCE.
C
C     EPS = 2*MACHEP =  4.4 * 10**-16 ON IBM 3081.
C     TTOL = EPS*TKMAX WHERE
C     TKMAX = MAX(BETA(K), K=1,KMAX)
C
C     ON EXIT
C     NDIS = TOTAL NUMBER OF COMPUTED DISTINCT T-EIGENVALUES OF
C            T(1,MEV) ON THE UNION OF THE (LB,UB) INTERVALS.
C     VS = COMPUTED DISTINCT T-EIGENVALUES OF T(1,MEV) IN ALGEBRAICALLY-
C          INCREASING ORDER
C     MP = CORRESPONDING T-MULTIPLICITIES OF THESE T-EIGENVALUES
C     MP(I) = (0,1,MI), MI>1, I=1,NDIS  MEANS:
C        (0)  V(I) IS SPURIOUS
C        (1)  V(I) IS T-ISOLATED AND GOOD
C        (MI) V(I) IS T-MULTIPLE AND HENCE A CONVERGED GOOD T-EIGENVALUE
C     IC = TOTAL NUMBER OF STURMS USED
C
C     DEFAULTS
C     ISKIP = 0 INITIALLY. IF DEFAULT OCCURS ON J-TH SUB-INTERVAL, SET
C             ISKIP=ISKIP+1 AND IDEF(ISKIP) = J
C             DEFAULTS OCCUR IF THERE ARE NO T-EIGENVALUES  IN THE
C             SUBINTERVAL SPECIFIED OR IF THE NUMBER
C             OF STURMS SEQUENCES REQUIRED EXCEEDS MXSTUR.
C             WHEN A DEFAULT OCCURS THE PROGRAM
C             SKIPS THE INTERVAL INVOLVED AND GOES ON TO THE NEXT
C             INTERVAL.
C
C-----------------------------------------------------------------------
C     SPECIFY PARAMETERS
      ZERO = 0.0D0
      ONE  = 1.0D0
      HALF = 0.5D0
      MXSTUR = IC
      NDIS = 0
      IC = 0
      ISKIP = 0
      MP1 = MEV+1
C     SAVE THEN SET BETA(MEV+1) = 0. GENERATE BETA**2
      BETAM = BETA(MP1)
      BETA(MP1) = ZERO
C
      DO 10 I = 1,MP1
   10 BETA2(I) = BETA(I)*BETA(I)
C
C     EP0 IS USED IN T-MULTIPLICITY AND SPURIOUS TESTS
C     EP1 AND EPS ARE USED IN THE BISEC CONVERGENCE TEST
C
      TEMP = DFLOAT(MEV+1000)
      EP0  = TEMP*TTOL
      EP1  = DSQRT(TEMP)*TTOL
C
      WRITE(6,20)MEV,NINT
   20 FORMAT(/' BISEC CALCULATION'/' ORDER OF T IS',I6/
     1' NUMBER OF INTERVALS IS',I6/)
C
      WRITE(6,30) EP0,EP1
   30 FORMAT(/' MULTOL, TOLERANCE USED IN T-MULTIPLICITY AND SPURIOUS TE
     1STS = ',E10.3/' BISTOL, TOLERANCE USED IN BISEC CONVERGENCE TEST =
     1',E10.3/)
C
C     LOOP ON THE NINT INTERVALS  (LB(J),UB(J)), J=1,NINT
      DO 430 JIND = 1,NINT
      LB = LBD(JIND)
      UB = UBD(JIND)
C
      WRITE(6,40)JIND,LB,UB
   40 FORMAT(//1X,'BISEC INTERVAL NO',2X,'LOWER BOUND',2X,'UPPER BOUND'/
     1I18,2E13.5/)
C
C     INITIALIZATION AND PARAMETER SPECIFICATION
C     ICT IS TOTAL STURM COUNT ON (LB,UB)
C
      NA = 0
      MD = 0
      NG = 0
      ICT = 0
C
C     START OF T-EIGENVALUE CALCULATIONS
      X1 = UB
      ISTURM = 1
      GO TO 330
C     FORWARD STURM CALCULATION TO DETERMINE NA = NO. T-EIGENVALUES > UB
   50 NA = NEV
C
      X1 = LB
      ISTURM = 2
      GO TO 330
C     FORWARD STURM CALC TO DETERMINE MT = NO. T-EIGENVALUES ON (LB,UB)
   60 CONTINUE
      MT=NEV
      ICT = ICT +2
C
      WRITE(6,70)MT,NA
   70 FORMAT(/2I6,' = NO. TMEV ON (LB,UB) AND NO. .GT. UB'/)
C
C     DEFAULT TEST: IS ESTIMATED NUMBER OF STURMS > MXSTUR?
      IEST = 30*MT
      IF (IEST.LT.MXSTUR) GO TO 90
C
      WRITE(6,80)
   80 FORMAT(//' ESTIMATED NUMBER OF STURMS REQUIRED EXCEEDS USER LIMIT'
     1/' SKIP THIS SUBINTERVAL')
      GO TO 110
C
   90 CONTINUE
C
      IF (MT.GE.1) GO TO 120
C
      WRITE(6,100)
  100 FORMAT(//' THERE ARE NO T-EIGENVALUES ON THIS INTERVAL)'/)
C
  110 ISKIP = ISKIP+1
      IDEF(ISKIP) = JIND
      GO TO 430
C
C     REGULAR CASE.
  120 CONTINUE
C
      IF (IWRITE.NE.0) WRITE(6,130)
  130 FORMAT(/' DISTINCT T-EIGENVALUES COMPUTED USING BISEC'/
     1 13X,'T-EIGENVALUE',2X,'TMULT',3X,'MD',4X,'NG')
C
C     SET UP INITIAL UPPER AND LOWER BOUNDS FOR T-EIGENVALUES
      DO 140 I=1,MT
      VB(I) = LB
      MTI = MT + I
  140 VB(MTI) = UB
C
C     CALCULATE T-EIGENVALUES FROM LB UP TO UB  K = MT,...,1
C     MAIN LOOP FOR FINDING KTH T-EIGENVALUE
C
      K = MT
  150 CONTINUE
      IC0 = 0
      XL = VB(K)
      MTK = MT+K
      XU = VB(MTK)
C
      ISTURM = 3
      X1 = XU
      IC0 = IC0 + 1
      GO TO 330
C     FORWARD STURM CALCULATION AT XU
  160 NU=NEV
C
C     BISECTION LOOP FOR KTH T-EIGENVALUE. TEST X1=MIDPOINT OF (XL,XU)
      ISTURM = 4
  170 CONTINUE
      X1 = (XL+XU)*HALF
      XS = DABS(XL)+DABS(XU)
      X0 = XU-XL
      EPT = EPS*XS+EP1
C
C     EPT IS CONVERGENCE TOLERANCE FOR KTH T-EIGENVALUE
C
      IF (X0.LE.EPT) GO TO 230
C
C     T-EIGENVALUE HAS NOT YET CONVERGED
C
      IC0 = IC0 + 1
      GO TO 330
C     FORWARD STURM CALCULATION AT CURRENT T-EIGENVALUE APPROXIMATION.
  180 CONTINUE
C
C     UPDATE T-EIGENVALUE INTERVAL (XL,XU)
C
      IF (NEV.LT.K) GO TO 190
C
C     NUMBER OF T-EIGENVALUES NEV = K
      XL = X1
      GO TO 170
  190 CONTINUE
C     NUMBER OF T-EIGENVALUES NEV<K
      XU = X1
      NU = NEV
C
C     UPDATE OF T-EIGENVALUE BOUNDS
C
      IF (NEV.EQ.0) GO TO 210
C
      DO 200 I = 1,NEV
  200 VB(I) = DMAX1(X1,VB(I))
C
  210 NEV1 = NEV+1
C
      DO 220 II = NEV1,K
      I = MT+II
  220 VB(I) = DMIN1(X1,VB(I))
C
      GO TO 170
C
C     END (XL,XU) BISECTION LOOP FOR KTH T-EIGENVALUE ON (LB,UB)
C     TEST FOR T-MULTIPLICITY AND IF SIMPLE THEN TEST FOR SPURIOUSNESS
C
  230 CONTINUE
      NDIS = NDIS+1
      MD = MD+1
      VS(NDIS) = X1
C
      JSTURM = 1
      X1 = XL-EP0
      GO TO 370
C     BACKWARD STURM CALCULATION
  240 KL = KEV
      JL = JEV
C
      JSTURM = 2
      IC0 = IC0 + 2
      X1 = XU+EP0
      GO TO 370
C     BACKWARD STURM CALCULATION
  250 JU = JEV
      KU = KEV
C
C     FOR T(1,MEV)
C     NU - KU = NO. T-EIGENVALUES ON (XU, XU + EP0)
C     KL - KU = NO. T-EIGENVALUES ON (XL - EP0, XU + EP0)
C
C     FOR T(2,MEV)
C     JL -JU = NO. T-EIGENVALUES ON (XL - EP0, XU + EP0)
C
C     IS THIS A SIMPLE T-EIGENVALUE?
C
      IF (KL-KU-1.EQ.0) GO TO 290
C
C     VS(NDIS) = KTH-T-EIGENVALUE OF (LB,UB) IS T-MULTIPLE AND HENCE
C     GOOD
      IF (KU.EQ.NU) GO TO 280
C     CONTINUE TO CHECK FOR T-MULTIPLICITY
  260 CONTINUE
      ISTURM = 5
      X1 = X1+EP0
      IC0 = IC0 + 1
      GO TO 330
C     FORWARD STURM CALCULATION
  270 KNE = KU-NEV
      KU = NEV
      IF (KNE.NE.0) GO TO 260
C     SPECIFY T-MULTIPLICITY = MP(NDIS)
  280 MPEV = KL-KU
      KNEW = KU
      GO TO 300
C     END T-MULTIPLE CASE
C
C     T-EIGENVALUE IS SIMPLE   CHECK IF IT IS SPURIOUS
  290 CONTINUE
      MPEV = 1
      IF (JU.LT.JL) MPEV=0
      KNEW = K-1
C
C     X1 >= XU+EP0
C     SPURIOUS TEST AND SIMPLE CASE COMPLETED
C     START OF NEXT T-EIGENVALUE COMPUTATION
C
  300 K = KNEW
      MP(NDIS) = MPEV
      IF (MPEV.GE.1) NG = NG + 1
C
      IF (IWRITE.NE.0) WRITE(6,310) VS(NDIS),MPEV,MD,NG
  310 FORMAT(E25.16,3I6)
C
C     UPDATE STURM COUNT. IC0 = STURM COUNT FOR KTH T-EIGENVALUE
      ICT = ICT + IC0
C
C     EXIT TEST FOR K DO LOOP
C
      IF (K.LE.0) GO TO 410
C
C     UPDATE LOWER BOUNDS
      DO 320 I=1,KNEW
  320 VB(I) = DMAX1(X1,VB(I))
C
      GO TO 150
C     END OF BISECTION LOOP FOR KTH EIGENVALUE
C
C     FORWARD STURM CALCULATION
  330 NEV = -NA
      YU  = ONE
C
      DO 360 I = 1,MEV
      IF (YU.NE.ZERO) GO TO 340
      YV = BETA(I)/EPS
      GO TO 350
  340 YV = BETA2(I)/YU
  350 YU = X1 - YV
      IF (YU.GE.ZERO) GO TO 360
      NEV = NEV + 1
  360 CONTINUE
C     NEV = NUMBER OF T-EIGENVALUES ON (X1,UB)
C
      GO TO (50,60,160,180,270), ISTURM
C
C     BACKWARD STURM CALCULATION FOR T(1,MEV) AND T(2,MEV)
  370 KEV = -NA
      YU  = ONE
C
      DO 400 II = 1,MEV
      I = MP1-II
      IF (YU.NE.ZERO) GO TO 380
      YV = BETA(I+1)/EPS
      GO TO 390
  380 YV = BETA2(I+1)/YU
  390 YU = X1-YV
      JEV = 0
      IF (YU.GE.ZERO) GO TO 400
      KEV = KEV+1
      JEV = 1
  400 CONTINUE
      JEV = KEV-JEV
C
      GO TO (240,250), JSTURM
C
C     KEV = -NA + (NUMBER OF T(1,MEV) T-EIGENVALUES) > X1
C     JEV = -NA + (NUMBER OF T(2,MEV) T-EIGENVALUES) > X1
C
C     SET PARAMETERS FOR NEXT INTERVAL
  410 CONTINUE
      IC = ICT+IC
      MXSTUR = MXSTUR-ICT
C
      WRITE(6,420) JIND,NG,MD,ICT
  420 FORMAT(/' T-EIGENVALUE CALCULATION ON INTERVAL',I6,'  IS COMPLETE'
     1 /3X,'NO. GOOD',3X,'NO. DISTINCT',4X,'STURMS'/I10,I13,I10)
C
  430 CONTINUE
C
C     END LOOP ON THE SUBINTERVALS (LB(J),UB(J)), J=1,NINT
C     ISKIP OUTPUT
C
      IF (ISKIP.GT.0) WRITE(6,440)ISKIP
  440 FORMAT(' BISEC DEFAULTED ON',I3,3X,'INTERVALS'/
     1 ' DEFAULTS OCCUR IF AN INTERVAL HAS NO T-EIGENVALUES'/
     2 ' OR THE STURM ESTIMATE EXCEEDS THE USER-SPECIFIED LIMIT'/)
C
      IF (ISKIP.GT.0) WRITE(6,450)(IDEF(I), I=1,ISKIP)
  450 FORMAT(' BISEC DEFAULTED ON INTERVALS'/(10I8))
C
C     RESET BETA AT I = MP1
      BETA(MP1) = BETAM
C-----END OF BISEC------------------------------------------------------
      RETURN
      END
C
C-----INVERSE ITERATION ON T(1,MEV)-------------------------------------
C
      SUBROUTINE INVERR(BETA,V1,V2,VS,EPS,G,MP,MEV,MMB,NDIS,NISO,
     1 NM,IKL,IT,IWRITE)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION BETA(1),V1(1),V2(1),VS(1)
      DOUBLE PRECISION X1,U,Z,EST,TEMP,T0,T1,RATIO,SUM,XU,NORM,TSUM
      DOUBLE PRECISION  BETAM,EPS,EPS3,EPS4,ZERO,ONE
      REAL G(1)
      INTEGER  MP(1)
C-----------------------------------------------------------------------
      DOUBLE PRECISION FINPRO
      REAL  ABS
      DOUBLE PRECISION  DABS, DMIN1, DSQRT, DFLOAT
C-----------------------------------------------------------------------
C     COMPUTES ERROR ESTIMATES FOR COMPUTED ISOLATED GOOD T-EIGENVALUES
C     IN VS AND WRITES THESE T-EIGENVALUES AND ESTIMATES TO FILE 4.
C     BY DEFINITION A GOOD T-EIGENVALUE IS ISOLATED IF ITS
C     CLOSEST NEIGHBOR IS ALSO GOOD, OR IF ONE OF ITS NEIGHBORS IS
C     SPURIOUS BUT THAT NEIGHBOR IS FAR ENOUGH AWAY.  SO
C     IN PARTICULAR, WE COMPUTE ESTIMATES FOR GOOD T-EIGENVALUES
C     THAT ARE IN CLUSTERS OF GOOD T-EIGENVALUES.
C
C     USES INVERSE ITERATION ON T(1,MEV) SOLVING THE EQUATION
C     (T - X1*I)V2 = RIGHT-HAND SIDE (RANDOMLY-GENERATED)
C     FOR EACH SUCH GOOD T-EIGENVALUE X1.
C
C     PROGRAM REFACTORS T-X1*I ON EACH ITERATION OF INVERSE ITERATION.
C     TYPICALLY ONLY ONE ITERATION IS NEEDED PER T-EIGENVALUE X1.
C
C     POSSIBLE STORAGE COMPRESSION
C     G STORAGE COULD BE ELIMINATED BY REGENERATING THE RANDOM
C     RIGHT-HAND SIDE ON EACH ITERATION AND PRINTING OUT THE
C     ERROR ESTIMATES AS THEY ARE GENERATED.
C
C     ON ENTRY AND EXIT
C     MEV = ORDER OF T
C     BETA CONTAINS THE NONZERO ENTRIES OF THE T-MATRIX
C     VS = COMPUTED DISTINCT EIGENVALUES OF T(1,MEV)
C     MP = T-MULTIPLICITY OF EACH EIGENVALUE IN VS. MP(I) = -1 MEANS
C          VS(I) IS A GOOD T-EIGENVALUE BUT THAT IT IS SITTING CLOSE TO
C          A SPURIOUS T-EIGENVALUE.  MP(I) = 0 MEANS VS(I) IS SPURIOUS.
C          ESTIMATES ARE COMPUTED ONLY FOR THOSE T-EIGENVALUES
C          WITH MP(I) = 1. FLAGGING WAS DONE IN SUBROUTINE ISOEV
C          PRIOR TO ENTERING INVERR.
C     NISO = NUMBER OF ISOLATED GOOD T-EIGENVALUES CONTAINED IN VS
C     NDIS =  NUMBER OF DISTINCT T-EIGENVALUES IN VS
C     IKL = SEED FOR RANDOM NUMBER GENERATOR
C     EPS = 2. * MACHINE EPSILON
C
C     IN PROGRAM:
C     ITER = MAXIMUM NUMBER OF INVERSE ITERATION STEPS ALLOWED FOR EACH
C            X1.  ITER = IT ON ENTRY.
C     G = ARRAY OF DIMENSION AT LEAST MEV + NISO.  USED TO STORE
C         RANDOMLY-GENERATED RIGHT-HAND SIDE.  THIS IS NOT
C         REGENERATED FOR EACH X1. G IS ALSO USED TO STORE ERROR
C         ESTIMATES AS THEY ARE COMPUTED FOR LATER PRINTOUT.
C     V1,V2 = WORK SPACES USED IN THE FACTORIZATION OF T(1,MEV).
C     AT THE END OF THE INVERSE ITERATION COMPUTATION FOR X1, V2
C     CONTAINS THE UNIT EIGENVECTOR OF T(1,MEV) CORRESPONDING TO X1.
C     V1 AND V2 MUST BE OF DIMENSION AT LEAST MEV.
C
C     ON EXIT
C     G(J) = MINIMUM GAP IN T(1,MEV) FOR EACH VS(J), J=1,NDIS
C     G(MEV+I) = BETAM*|V2(MEV)| = ERROR ESTIMATE FOR ISOLATED GOOD
C              T-EIGENVALUES, WHERE I = 1,NISO  AND  BETAM = BETA(MEV+1)
C              V2(MEV) IS LAST COMPONENT OF THE UNIT T-EIGENVECTOR OF
C              T(1,MEV) CORRESPONDING TO ITH ISOLATED GOOD T-EIGENVALUE.
C
C     IF FOR SOME X1 IT.GT.ITER THEN THE ERROR ESTIMATE IN G IS MARKED
C     WITH A - SIGN.
C
C     V2 = ISOLATED GOOD T-EIGENVALUES
C     V1 = MINIMAL T-GAPS FOR THE EIGENVALUES IN V2.
C     THESE ARE CONSTRUCTED FOR WRITE-OUT PURPOSES ONLY AND NOT
C     NEEDED ELSEWHERE IN THE PROGRAM.
C-----------------------------------------------------------------------
C
C     LABEL OUTPUT FILE  4
      IF (MMB.EQ.1) WRITE(4,10)
   10 FORMAT(' INVERSE ITERATION ERROR ESTIMATES'/)
C
C     FILE 6 (TERMINAL) OUTPUT OF ERROR ESTIMATES
      IF (IWRITE.NE.0.AND.NISO.NE.0) WRITE(6,20)
   20 FORMAT(/' INVERSE ITERATION ERROR ESTIMATES'/'  JISO',' JDIST',8X
     1,'GOOD T-EIGENVALUE',4X,'BETAM*UM',5X,'TMINGAP')
C
C     INITIALIZATION AND PARAMETER SPECIFICATION
      ZERO = 0.0D0
      ONE = 1.0D0
      NG = 0
      NISO = 0
      ITER = IT
      MP1 = MEV+1
      MM1 = MEV-1
      BETAM = BETA(MP1)
      BETA(MP1) = ZERO
C
C     CALCULATE SCALE AND TOLERANCES
      TSUM = ZERO
      DO 30 I = 2,MEV
   30 TSUM = TSUM + BETA(I)
C
      EPS3 = EPS*TSUM
      EPS4 = DFLOAT(MEV)*EPS3
C
C     GENERATE SCALED RANDOM RIGHT-HAND SIDE
      ILL = IKL
C
C-----------------------------------------------------------------------
      CALL GENRAN(ILL,G,MEV)
C-----------------------------------------------------------------------
C
      GSUM = ZERO
      DO 40 I = 1,MEV
   40 GSUM = GSUM+ABS(G(I))
      GSUM = EPS4/GSUM
C
      DO 50 I = 1,MEV
   50 G(I) = GSUM*G(I)
C
C     LOOP ON ISOLATED GOOD T-EIGENVALUES IN VS (MP(I) = 1) TO
C     CALCULATE CORRESPONDING UNIT T-EIGENVECTOR OF T(1,MEV)
C
      DO 180 JEV = 1,NDIS
      IF (MP(JEV).EQ.0) GO TO 180
      NG = NG + 1
      IF (MP(JEV).NE.1) GO TO 180
      IT = 1
      NISO = NISO + 1
      X1 = VS(JEV)
C
C     INITIALIZE RIGHT HAND SIDE FOR INVERSE ITERATION
      DO 60 I = 1,MEV
   60 V2(I) = G(I)
C
C     TRIANGULAR FACTORIZATION WITH NEAREST NEIGHBOR PIVOT
C     STRATEGY. INTERCHANGES ARE LABELLED BY SETTING BETA < 0.
C
   70 CONTINUE
      U = -X1
      Z = BETA(2)
C
      DO 90 I = 2,MEV
      IF (BETA(I).GT.DABS(U)) GO TO 80
C     NO INTERCHANGE
      V1(I-1) = Z/U
      V2(I-1) = V2(I-1)/U
      V2(I) = V2(I)-BETA(I)*V2(I-1)
      RATIO = BETA(I)/U
      U = -X1-Z*RATIO
      Z = BETA(I+1)
      GO TO 90
   80 CONTINUE
C     INTERCHANGE CASE
      RATIO = U/BETA(I)
      BETA(I) = -BETA(I)
      V1(I-1) = -X1
      U = Z-RATIO*V1(I-1)
      Z = -RATIO*BETA(I+1)
      TEMP = V2(I-1)
      V2(I-1) = V2(I)
      V2(I) = TEMP-RATIO*V2(I)
   90 CONTINUE
      IF (U.EQ.ZERO) U = EPS3
C
C     SMALLNESS TEST AND DEFAULT VALUE FOR LAST COMPONENT
C     PIVOT(I-1) = |BETA(I)| FOR INTERCHANGE CASE
C     (I-1,I+1) ELEMENT IN RIGHT FACTOR = BETA(I+1)
C     END OF FACTORIZATION AND FORWARD SUBSTITUTION
C
C     BACK SUBSTITUTION
      V2(MEV) = V2(MEV)/U
      DO 110 II = 1,MM1
      I = MEV-II
      IF (BETA(I+1).LT.ZERO) GO TO 100
C     NO INTERCHANGE
      V2(I) = V2(I)-V1(I)*V2(I+1)
      GO TO 110
C     INTERCHANGE CASE
  100 BETA(I+1) = -BETA(I+1)
      V2(I) = (V2(I)-V1(I)*V2(I+1)-BETA(I+2)*V2(I+2))/BETA(I+1)
  110 CONTINUE
C
C     TESTS FOR CONVERGENCE OF INVERSE ITERATION
C     IF SUM |V2| COMPS. LE. 1 AND IT. LE. ITER DO ANOTHER INVIT STEP
C
      NORM = DABS(V2(MEV))
      DO 120 II = 1,MM1
      I = MEV-II
  120 NORM = NORM+DABS(V2(I))
C
      IF (NORM.GE.ONE) GO TO 140
      IT = IT+1
      IF (IT.GT.ITER) GO TO 140
      XU = EPS4/NORM
C
      DO 130 I = 1,MEV
  130 V2(I) = V2(I)*XU
C
      GO TO 70
C     ANOTHER INVERSE ITERATION STEP
C
C     INVERSE ITERATION FINISHED
C     NORMALIZE COMPUTED T-EIGENVECTOR : V2 = V2/||V2||
  140 CONTINUE
      SUM = FINPRO(MEV,V2(1),1,V2(1),1)
      SUM = ONE/DSQRT(SUM)
C
      DO 150 II = 1,MEV
  150 V2(II) = SUM*V2(II)
C
C     SAVE ERROR ESTIMATE FOR LATER OUTPUT
      EST = BETAM*DABS(V2(MEV))
      IF (IT.GT.ITER) EST = -EST
      MEVPNI = MEV + NISO
      G(MEVPNI) = EST
      IF (IWRITE.EQ.0) GO TO 180
C
C     FILE 6 (TERMINAL) OUTPUT OF ERROR ESTIMATES.
      IF (JEV.EQ.1) GAP = VS(2) - VS(1)
      IF (JEV.EQ.MEV) GAP = VS(MEV) - VS(MEV-1)
      IF (JEV.EQ.MEV.OR.JEV.EQ.1) GO TO 160
      TEMP = DMIN1(VS(JEV+1)-VS(JEV),VS(JEV)-VS(JEV-1))
      GAP = TEMP
  160 CONTINUE
C
      WRITE(6,170) NISO,JEV,X1,EST,GAP
  170 FORMAT(2I6,E25.16,2E12.3)
C
  180 CONTINUE
C
C     END ERROR ESTIMATE LOOP ON ISOLATED GOOD T-EIGENVALUES.
C     GENERATE DISTINCT MINGAPS FOR T(1,MEV).  THIS IS USEFUL AS AN
C     INDICATOR OF THE GOODNESS OF THE INVERSE ITERATION ESTIMATES.
C     TRANSFER ISOLATED GOOD T-EIGENVALUES AND CORRESPONDING TMINGAPS
C     TO V2 AND V1 FOR OUTPUT PURPOSES ONLY.
C
      NM1 = NDIS - 1
      G(NDIS) = VS(NM1)-VS(NDIS)
      G(1) = VS(2)-VS(1)
C
      DO 190 J = 2,NM1
      T0 = VS(J)-VS(J-1)
      T1 = VS(J+1)-VS(J)
      G(J) = T1
      IF (T0.LT.T1) G(J)=-T0
  190 CONTINUE
      ISO = 0
      DO 200 J = 1,NDIS
      IF (MP(J).NE.1) GO TO 200
      ISO = ISO+1
      V1(ISO) = G(J)
      V2(ISO) = VS(J)
  200 CONTINUE
C
      IF(NISO.EQ.0) GO TO 250
C
C     ERROR ESTIMATES ARE WRITTEN TO FILE 4
      WRITE(4,210)MEV,NDIS,NG,NISO,NM,IKL,ITER,BETAM
  210 FORMAT(1X,'TSIZE',2X,'NDIS',1X,'NGOOD',2X,'NISO',3X,'M+N'/5I6/
     1 4X,'RHSEED',2X,'MXINIT',5X,'BETAM'/I10,I8,E10.3/
     2 2X,'GOODEVNO',8X,'GOOD T-EIGENVALUE',6X,'BETAM*UM',7X,'TMINGAP')
C
      ISPUR = 0
      I = 0
      DO 240 J = 1,NDIS
      IF(MP(J).NE.0) GO TO 220
      ISPUR = ISPUR + 1
      GO TO 240
  220 IF(MP(J).NE.1) GO TO 240
      I = I + 1
      MEVI = MEV + I
      IGOOD = J - ISPUR
      WRITE(4,230) IGOOD,V2(I),G(MEVI),V1(I)
  230 FORMAT(I10,E25.16,2E14.3)
  240 CONTINUE
      GO TO 270
C
  250 WRITE(4,260)
  260 FORMAT(/' THERE ARE NO ISOLATED T-EIGENVALUES SO NO ERROR ESTIMATE
     1S WERE COMPUTED')
C     RESTORE BETA(MEV+1) = BETAM
  270 BETA(MP1) = BETAM
C-----END OF INVERR-----------------------------------------------------
      RETURN
      END
C
C-----START OF TNORM----------------------------------------------------
C
      SUBROUTINE TNORM(BETA,BMIN,TMAX,MEV,IB)
C
C-----------------------------------------------------------------------
      DOUBLE PRECISION  BETA(1)
      DOUBLE PRECISION  TMAX,BMIN,BSIZE,BTOL
      DOUBLE PRECISION  DABS, DMAX1
C-----------------------------------------------------------------------
C     COMPUTE SCALING FACTOR USED IN THE T-MULTIPLICITY, SPURIOUS AND
C     PRTESTS.   CHECK RELATIVE SIZE OF THE BETA(K), K=1,MEV
C     AS A TEST ON THE LOCAL ORTHOGONALITY OF THE LANCZOS VECTORS.
C
C             TMAX = MAX (BETA(I),  I=1,MEV)
C             BMIN = MIN (BETA(I) I=2,MEV)
C             BSIZE = BMIN/TMAX
C             |IB| = INDEX OF MINIMAL(BETA)
C             IB < 0 IF BMIN/TMAX < BTOL
C-----------------------------------------------------------------------
C     SPECIFY PARAMETERS
      IB = 2
      BTOL = BMIN
      BMIN = BETA(2)
      TMAX = BETA(2)
C
      DO 20 I = 2,MEV
      IF (BETA(I).GE.BMIN) GO TO 10
      IB = I
      BMIN = BETA(I)
   10 TMAX = DMAX1(TMAX,BETA(I))
   20 CONTINUE
C
C     TEST OF LOCAL ORTHOGONALITY USING SCALED BETAS
      BSIZE = BMIN/TMAX
      IF (BSIZE.GE.BTOL) GO TO 40
C
C     DEFAULT.  BSIZE IS SMALLER THAN TOLERANCE BTOL SPECIFIED IN MAIN
C     PROGRAM.  PROGRAM TERMINATES FOR USER TO DECIDE WHAT TO DO
C     BECAUSE LOCAL ORTHOGONALITY OF THE LANCZOS VECTORS COULD BE
C     LOST.
C
      IB = -IB
      WRITE(6,30) MEV
   30 FORMAT(/' BETA TEST INDICATES POSSIBLE LOSS OF LOCAL ORTHOGONALITY
     1OVER 1ST',I6,' LANCZOS VECTORS'/)
C
   40 CONTINUE
C
      WRITE(6,50) IB
   50 FORMAT(/' MINIMUM BETA RATIO OCCURS AT',I6,' TH BETA'/)
C
      WRITE(6,60) MEV,BMIN,TMAX,BSIZE
   60 FORMAT(/1X,'TSIZE',6X,'MIN BETA',5X,'TKMAX',6X,'MIN RATIO'/
     1 I6,E14.3,E10.3,E15.3/)
C
C-----END OF TNORM------------------------------------------------------
      RETURN
      END

.