C-----HLEVAL  (EIGENVALUES OF HERMITIAN MATRICES)-----------------------
C
C     CONTAINS MAIN PROGRAM FOR COMPUTING DISTINCT EIGENVALUES OF
C     A HERMITIAN MATRIX USING LANCZOS TRIDIAGONALIZATION WITHOUT
C     REORTHOGONALIZATION
C
C     PORTABILITY:
C     THIS PROGRAM IS NOT PORTABLE DUE TO THE USE OF COMPLEX*16
C     VARIABLES.  MOREOVER, THE PFORT VERIFIER IDENTIFIED THE
C     FOLLOWING ADDITIONAL NONPORTABLE CONSTRUCTIONS:
C
C     1.  DATA/MACHEP/ STATEMENT
C     2.  ALL READ(5,*) STATEMENTS (FREE FORMAT)
C     3.  FORMAT(20A4) USED WITH EXPLANATORY HEADER EXPLAN.
C     4.  HEXADECIMAL FORMAT (4Z20) USED IN ALPHA/BETA FILES 1 AND 2.
C
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION  ALPHA(5000),BETA(5001),VS(5000)
      COMPLEX*16  V1(5000),V2(5000)
      DOUBLE PRECISION  GR(1500),GC(1500),LB(20),UB(20)
      DOUBLE PRECISION  BTOL,GAPTOL,TTOL,MACHEP,EPSM,RELTOL
      DOUBLE PRECISION  SCALE1,SCALE2,SCALE3,SCALE4,BISTOL,CONTOL,MULTOL
      DOUBLE PRECISION  ONE,ZERO,TEMP,TKMAX,BETAM,BKMIN,T0,T1
      REAL  G(5000),EXPLAN(20)
      INTEGER  MP(5000),NMEV(20)
      INTEGER  SVSEED,RHSEED,SVSOLD
      INTEGER IABS
      REAL  ABS
      DOUBLE PRECISION  DABS, DSQRT, DFLOAT
      EXTERNAL CMATV
C
C-----------------------------------------------------------------------
      DATA MACHEP/Z3410000000000000/
      EPSM = 2.0D0*MACHEP
C-----------------------------------------------------------------------
C     THE ARRAYS V1 AND V2 ARE DEFINED AS COMPLEX*16 IN THE MAIN PROGRAM
C     AND IN THE SUBROUTINE LANCZS.  HOWEVER, IN THE OTHER SUBROUTINES
C     THEY ARE DECLARED AS DOUBLE PRECISION ARRAYS. NOTE THAT THE
C     DIMENSION OF V2 ASSUMES THAT NO MORE THAN KMAX/2 EIGENVALUES OF
C     THE T-MATRICES ARE BEING COMPUTED IN ANY ONE OF THE SUB-INTERVALS
C     BEING CONSIDERED.  V2 MUST CONTAIN UPPER AND LOWER BOUNDS
C     ON EACH T-EIGENVALUE COMPUTED BY BISEC IN ANY ONE GIVEN INTERVAL.
C
C     ARRAYS MUST BE DIMENSIONED AS FOLLOWS:
C          1. ALPHA:   >= KMAX.  BETA:  >= (KMAX+1)
C          2. V1:  >= MAX(N,(KMAX+1)/2).   V2:  >= MAX(N,KMAX/2)
C          3. VS:  >= MAX(2*N,KMAX).
C          4. GR,GC:  >=  N
C          5. G:  >= MAX(2*KMAX,N)
C          6. MP:   >= KMAX
C          7. LB,UB:  >= NUMBER OF SUB-INTERVALS SPECIFIED
C          8. NMEV:  >= NUMBER OF T-MATRICES SPECIFIED
C          9. EXPLAN:  DIMENSION IS 20.
C
C
C     IMPORTANT TOLERANCES OR SCALES THAT ARE USED REPEATEDLY
C     THROUGHOUT THE PROGRAM ARE THE FOLLOWING:
C     SCALED MACHINE EPSILON:  TTOL = TKMAX*EPSM WHERE
C     EPSM = 2*MACHINE EPSILON AND
C     TKMAX = MAX(|ALPHA(J)|,BETA(J), J = 1,MEV)
C     BISEC CONVERGENCE TOLERANCE:  BISTOL = DSQRT(1000+MEV)*TTOL
C     BISEC MULTIPLICITY TOLERANCE:  MULTOL = (1000+MEV)*TTOL
C     LANCZOS CONVERGENCE TOLERANCE:   CONTOL = BETA(MEV+1)*1.D-10
C-----------------------------------------------------------------------
C     OUTPUT HEADER
      WRITE(6,10)
   10 FORMAT(/' LANCZOS PROCEDURE FOR HERMITIAN MATRICES'/)
C
C     SET PROGRAM PARAMETERS
C     SCALEK ARE USED IN TOLERANCES NEEDED IN SUBROUTINES LUMP,
C     ISOEV AND PRTEST.  USER MUST NOT MODIFY THESE SCALES.
      SCALE1 = 5.0D2
      SCALE2 = 5.0D0
      SCALE3 = 5.0D0
      SCALE4 = 1.0D4
      ONE  = 1.0D0
      ZERO = 0.0D0
C     BTOL = EPSM
      BTOL = 1.0D-8
      GAPTOL = 1.0D-8
      ICONV = 0
      MOLD = 0
      MOLD1 = 1
      ICT = 0
      MMB = 0
      IPROJ = 0
C
C     READ USER-SPECIFIED PARAMETERS FROM INPUT FILE 5 (FREE FORMAT)
C
C     READ USER-PROVIDED HEADER FOR RUN
      READ(5,20) EXPLAN
      WRITE(6,20) EXPLAN
      READ(5,20) EXPLAN
      WRITE(6,20) EXPLAN
   20 FORMAT(20A4)
C
C     READ ORDER OF MATRICES (N) , MAXIMUM ORDER OF T-MATRIX (KMAX),
C     NUMBER OF T-MATRICES ALLOWED (NMEVS), AND MATRIX IDENTIFICATION
C     NUMBERS (MATNO)
      READ(5,20) EXPLAN
      READ(5,*) N,KMAX,NMEVS,MATNO
C
C     READ SEEDS FOR LANCZS AND INVERR SUBROUTINES (SVSEED AND RHSEED)
C     READ MAXIMUM NUMBER OF ITERATIONS ALLOWED FOR EACH INVERSE
C     ITERATION (MXINIT) AND MAXIMUM NUMBER OF STURM SEQUENCES
C     ALLOWED (MXSTUR)
      READ(5,20) EXPLAN
      READ(5,*) SVSEED,RHSEED,MXINIT,MXSTUR
C
C     ISTART = (0,1):  ISTART = 0 MEANS ALPHA/BETA FILE IS NOT
C     AVAILABLE.  ISTART = 1 MEANS ALPHA/BETA FILE IS AVAILABLE ON
C     FILE 2.
C     ISTOP = (0,1):  ISTOP = 0 MEANS PROCEDURE GENERATES ALPHA/BETA
C     FILE AND THEN TERMINATES.  ISTOP = 1 MEANS PROCEDURE GENERATES
C     ALPHAS/BETAS IF NEEDED AND THEN COMPUTES EIGENVALUES AND ERROR
C     ESTIMATES AND THEN TERMINATES.
      READ(5,20) EXPLAN
      READ(5,*) ISTART,ISTOP
C
C     IHIS = (0,1):  IHIS = 0 MEANS ALPHA/BETA FILE IS NOT WRITTEN
C     TO FILE 1.  IHIS = 1 MEANS ALPHA/BETA FILE IS WRITTEN TO FILE 1.
C     IDIST = (0,1):  IDIST = 0 MEANS DISTINCT T-EIGENVALUES
C     ARE NOT WRITTEN TO FILE 11.  IDIST = 1 MEANS DISTINCT
C     T-EIGENVALUES ARE WRITTEN TO FILE 11.
C     IWRITE = (0,1):  IWRITE = 0 MEANS NO INTERMEDIATE OUTPUT
C     FROM THE COMPUTATIONS IS WRITTEN TO FILE 6.  IWRITE = 1 MEANS
C     T-EIGENVALUES AND ERROR ESTIMATES ARE WRITTEN TO FILE 6
C     AS THEY ARE COMPUTED.
      READ(5,20) EXPLAN
      READ(5,*) IHIS,IDIST,IWRITE
C
C     READ IN THE RELATIVE TOLERANCE (RELTOL) FOR USE IN THE
C     SPURIOUS, T-MULTIPLICITY, AND PRTESTS.
      READ(5,20) EXPLAN
      READ(5,*) RELTOL
C
C     READ IN THE SIZES OF THE T-MATRICES TO BE CONSIDERED.
      READ(5,20) EXPLAN
      READ(5,*) (NMEV(J), J=1,NMEVS)
C
C     READ IN THE NUMBER OF SUBINTERVALS TO BE CONSIDERED.
      READ(5,20) EXPLAN
      READ(5,*) NINT
C
C     READ IN THE LEFT-END POINTS OF THE SUBINTERVALS TO BE CONSIDERED.
C     THESE MUST BE IN ALGEBRAICALLY-INCREASING ORDER
      READ(5,20) EXPLAN
      READ(5,*) (LB(J), J=1,NINT)
C
C     READ IN THE RIGHT-END POINTS OF THE SUBINTERVALS TO BE CONSIDERED.
C     THESE MUST BE IN ALGEBRAICALLY-INCREASING ORDER
      READ(5,20) EXPLAN
      READ(5,*) (UB(J), J=1,NINT)
C
C-----------------------------------------------------------------------
C
C     INITIALIZE THE ARRAYS FOR THE USER-SPECIFIED MATRIX
C     AND PASS THE STORAGE LOCATIONS OF THESE ARRAYS TO THE
C     MATRIX-VECTOR MULTIPLY SUBROUTINE CMATV.
C
      CALL USPEC(N,MATNO)
C
C-----------------------------------------------------------------------
C
C     MASK UNDERFLOW AND OVERFLOW
C
      CALL MASK
C
C-----------------------------------------------------------------------
C
C     WRITE TO FILE 6, A SUMMARY OF THE PARAMETERS FOR THIS RUN
C
      WRITE(6,30) MATNO,N,KMAX
   30 FORMAT(/3X,'MATRIX ID',4X,'ORDER OF A',4X,'MAX ORDER OF T'/
     1 I12,I14,I18/)
C
      WRITE(6,40) ISTART,ISTOP
   40 FORMAT(/2X,'ISTART',3X,'ISTOP'/2I8/)
C
      WRITE(6,50) IHIS,IDIST,IWRITE
   50 FORMAT(/4X,'IHIS',3X,'IDIST',2X,'IWRITE'/3I8/)
C
      WRITE(6,60) SVSEED,RHSEED
   60 FORMAT(/' SEEDS FOR RANDOM NUMBER GENERATOR'//
     1 4X,'LANCZS SEED',4X,'INVERR SEED'/2I15/)
C
      WRITE(6,70) (NMEV(J), J=1,NMEVS)
   70 FORMAT(/' SIZES OF T-MATRICES TO BE CONSIDERED'/(6I12))
C
      WRITE(6,80) RELTOL,GAPTOL,BTOL
   80 FORMAT(/' RELATIVE TOLERANCE USED TO COMBINE COMPUTED T-EIGENVALUE
     1S'/E15.3/' RELATIVE GAP TOLERANCES USED IN INVERSE ITERATION'/
     1E15.3/' RELATIVE TOLERANCE FOR CHECK ON SIZE OF BETAS'/E15.3/)
C
      WRITE(6,90) (J,LB(J),UB(J), J=1,NINT)
   90 FORMAT(/' BISEC WILL BE USED ON THE FOLLOWING INTERVALS'/
     1 (I6,2E20.6))
C
      IF (ISTART.EQ.0) GO TO 140
C
C     READ IN ALPHA BETA HISTORY
C
      READ(2,100)MOLD,NOLD,SVSOLD,MATOLD
  100 FORMAT(2I6,I12,I8)
C
      IF (KMAX.LT.MOLD) KMAX = MOLD
      KMAX1 = KMAX + 1
C
C     CHECK THAT ORDER N, MATRIX ID MATNO, AND RANDOM SEED SVSEED
C     AGREE WITH THOSE IN THE HISTORY FILE.  IF NOT PROCEDURE STOPS.
C
      ITEMP = (NOLD-N)**2+(MATNO-MATOLD)**2+(SVSEED-SVSOLD)**2
C
      IF (ITEMP.EQ.0) GO TO 120
C
      WRITE(6,110)
  110 FORMAT(' PROGRAM TERMINATES'/    ' READ FROM FILE 2 CORRESPONDS TO
     1 DIFFERENT MATRIX THAN MATRIX SPECIFIED'/)
      GO TO 640
C
  120 CONTINUE
      MOLD1 = MOLD+1
C
      READ(2,130)(ALPHA(J), J=1,MOLD)
      READ(2,130)(BETA(J), J=1,MOLD1)
  130 FORMAT(4Z20)
C
      IF (KMAX.EQ.MOLD) GO TO 160
C
      READ(2,130)(V1(J), J=1,N)
      READ(2,130)(V2(J), J=1,N)
C
  140 CONTINUE
      IIX = SVSEED
C
C-----------------------------------------------------------------------
C
      CALL LANCZS(CMATV,V1,V2,VS,ALPHA,BETA,GR,GC,G,KMAX,MOLD1,N,IIX)
C
C-----------------------------------------------------------------------
C
      KMAX1 = KMAX + 1
C
      IF (IHIS.EQ.0.AND.ISTOP.GT.0) GO TO 160
C
      WRITE(1,150) KMAX,N,SVSEED,MATNO
  150 FORMAT(2I6,I12,I8,' = KMAX,N,SVSEED,MATNO')
C
C     TO AVOID PERTURBATIONS CAUSED BY HEX TO DECIMAL AND DECIMAL TO HEX
C     CONVERSIONS, THE ALPHA AND BETA MUST BE WRITTEN OUT IN HEX.
      WRITE(1,130)(ALPHA(I), I=1,KMAX)
      WRITE(1,130)(BETA(I), I=1,KMAX1)
C
      WRITE(1,130)(V1(I), I=1,N)
      WRITE(1,130)(V2(I), I=1,N)
C
      IF (ISTOP.EQ.0) GO TO 540
C
  160 CONTINUE
      BKMIN = BTOL
      WRITE(6,170)
  170 FORMAT(/' T-MATRICES (ALPHA AND BETA) ARE NOW AVAILABLE'/)
C
C-----------------------------------------------------------------------
C     SUBROUTINE TNORM CHECKS MIN(BETA)/(ESTIMATED NORM(A)) > BTOL .
C     IF THIS IS VIOLATED IB IS SET EQUAL TO THE NEGATIVE OF THE INDEX
C     OF THE MINIMAL BETA.  IF(IB < 0) THEN SUBROUTINE TNORM IS
C     CALLED FOR EACH VALUE OF MEV TO DETERMINE WHETHER OR NOT THERE
C     IS A BETA IN THE T-MATRIX SPECIFIED THAT VIOLATES THIS TEST.
C     IF THERE IS SUCH A BETA THE PROGRAM TERMINATES FOR THE USER
C     TO DECIDE WHAT TO DO.  THIS TEST CAN BE OVER-RIDDEN BY
C     SIMPLY MAKING BTOL SMALLER, BUT THEN THERE IS THE POSSIBILITY
C     THAT LOSSES IN THE LOCAL ORTHOGONALITY MAY HURT THE COMPUTATIONS.
C     BTOL = 1.D-8 IS HOWEVER A CONSERVATIVE CHOICE FOR BTOL.
C
C     TNORM ALSO COMPUTES TKMAX = MAX(|ALPHA(K)|,BETA(K), K=1,KMAX).
C     TKMAX IS USED TO SCALE THE TOLERANCES USED IN THE
C     T-MULTIPLICITY AND SPURIOUS TESTS IN BISEC. TKMAX IS ALSO USED IN
C     THE PROJECTION TEST FOR HIDDEN EIGENVALUES THAT HAD 'TOO SMALL'
C     A PROJECTION ON THE STARTING VECTOR.
C
      CALL TNORM(ALPHA,BETA,BKMIN,TKMAX,KMAX,IB)
C
C-----------------------------------------------------------------------
C
      TTOL = EPSM*TKMAX
C
C     LOOP ON THE SIZE OF THE T-MATRIX
C
  180 CONTINUE
      MMB = MMB + 1
      MEV = NMEV(MMB)
C     IS MEV TOO LARGE ?
      IF(MEV.LE.KMAX) GO TO 200
      WRITE(6,190) MMB, MEV, KMAX
  190 FORMAT(/' TERMINATE PRIOR TO CONSIDERING THE',I6,'TH T-MATRIX'/
     1' BECAUSE THE SIZE REQUESTED',I6,' IS GREATER THAN THE MAXIMUM SIZ
     1E ALLOWED',I6/)
      GO TO 540
C
  200 MP1 = MEV + 1
      BETAM = BETA(MP1)
C
      IF (IB.GE.0) GO TO 210
C
      T0 = BTOL
C
C-----------------------------------------------------------------------
C
      CALL TNORM(ALPHA,BETA,T0,T1,MEV,IBMEV)
C
C-----------------------------------------------------------------------
C
      TEMP = T0/TKMAX
      IBMEV = IABS(IBMEV)
      IF (TEMP.GE.BTOL) GO TO 210
      IBMEV = -IBMEV
      GO TO 600
C
  210 CONTINUE
      IC = MXSTUR-ICT
C
C-----------------------------------------------------------------------
C     BISEC LOOP. THE SUBROUTINE BISEC INCORPORATES DIRECTLY THE
C     T-MULTIPLICITY AND SPURIOUS TESTS. T-EIGENVALUES WILL BE
C     CALCULATED BY BISEC SEQUENTIALLY ON INTERVALS
C     (LB(J),UB(J)), J = 1,NINT).
C
C     ON RETURN FROM BISEC
C     NDIS = NUMBER OF DISTINCT EIGENVALUES OF T(1,MEV) ON UNION
C            OF THE (LB,UB) INTERVALS
C     VS = DISTINCT T-EIGENVALUES IN ALGEBRAICALLY INCREASING ORDER
C     MP = MULTIPLICITIES OF THE T-EIGENVALUES IN VS
C     MP(I) = (0,1,MI), MI>1, I=1,NDIS  MEANS:
C        (0)  VS(I) IS SPURIOUS
C        (1)  VS(I) IS T-SIMPLE AND GOOD
C        (MI) VS(I) IS MULTIPLE AND IS THEREFORE NOT ONLY GOOD BUT
C             ALSO A CONVERGED GOOD T-EIGENVALUE.
C     WITHIN BISEC V1 AND V2 ARE DEFINED AS DOUBLE PRECISION ARRAYS
C
C
      CALL BISEC(ALPHA,BETA,V1,V2,VS,LB,UB,EPSM,TTOL,MP,NINT,
     1 MEV,NDIS,IC,IWRITE)
C
C-----------------------------------------------------------------------
C
      IF (NDIS.EQ.0) GO TO 620
C
C     COMPUTE THE TOTAL NUMBER OF STURM SEQUENCES USED TO DATE
C     COMPUTE THE BISEC CONVERGENCE AND T-MULTIPLICITY TOLERANCES USED.
C     COMPUTE THE CONVERGENCE TOLERANCE FOR EIGENVALUES OF A.
      ICT = ICT + IC
      TEMP = DFLOAT(MEV+1000)
      MULTOL = TEMP*TTOL
      TEMP = DSQRT(TEMP)
      BISTOL = TTOL*TEMP
      CONTOL = BETAM*1.D-10
C
C-----------------------------------------------------------------------
C    SUBROUTINE LUMP 'COMBINES' T-EIGENVALUES THAT ARE 'TOO CLOSE'.
C    NOTE HOWEVER THAT CLOSE SPURIOUS T-EIGENVALUES ARE NOT AVERAGED
C    WITH GOOD ONES. HOWEVER, THEY MAY BE USED TO INCREASE THE
C    MULTIPLICITY OF A GOOD T-EIGENVALUE.
C
      LOOP = NDIS
      CALL LUMP(VS,RELTOL,MULTOL,SCALE2,MP,LOOP)
C
C-----------------------------------------------------------------------
C
      IF(NDIS.EQ.LOOP) GO TO 230
C
      WRITE(6,220) NDIS, MEV, LOOP
  220 FORMAT(/I6,' DISTINCT T-EIGENVALUES WERE COMPUTED IN BISEC AT MEV
     1',I6/ 2X,' LUMP SUBROUTINE REDUCES NUMBER OF DISTINCT EIGENVALUES
     1TO',I6)
C
  230 CONTINUE
      NDIS = LOOP
      BETA(MP1) = BETAM
C
C-----------------------------------------------------------------------
C     THE SUBROUTINE ISOEV LABELS THOSE SIMPLE EIGENVALUES OF T(1,MEV)
C     WITH VERY SMALL GAPS BETWEEN NEIGHBORING EIGENVALUES OF T(1,MEV)
C     TO AVOID COMPUTING ERROR ESTIMATES FOR ANY SIMPLE GOOD
C     T-EIGENVALUE THAT IS TOO CLOSE TO A SPURIOUS T-EIGENVALUE.
C     ON RETURN FROM ISOEV, G CONTAINS CODED MINIMAL GAPS
C     BETWEEN THE DISTINCT EIGENVALUES OF T(1,MEV). (G IS REAL).
C     G(I) < 0 MEANS MINGAP IS DUE TO LEFT GAP G(I) > 0 MEANS DUE TO
C     RIGHT GAP. MP(I) = -1 MEANS THAT THE GOOD T-EIGENVALUE IS SIMPLE
C     AND HAS A VERY SMALL MINGAP IN T(1,MEV) DUE TO A SPURIOUS
C     T-EIGENVALUE.   NG = NUMBER OF GOOD EIGENVALUES.
C     NISO = NUMBER OF ISOLATED GOOD T-EIGENVALUES.
C
      CALL ISOEV(VS,GAPTOL,MULTOL,SCALE1,G,MP,NDIS,NG,NISO)
C
C-----------------------------------------------------------------------
C
      WRITE(6,240)NG,NISO,NDIS
  240 FORMAT(/I6,' GOOD T-EIGENVALUES HAVE BEEN COMPUTED'/
     1 I6,' OF THESE ARE T-ISOLATED'/
     2 I6,' = NUMBER OF DISTINCT T-EIGENVALUES COMPUTED'/)
C
C     DO WE WRITE DISTINCT EIGENVALUES OF T-MATRIX TO FILE 4?
      IF (IDIST.EQ.0) GO TO 280
C
      WRITE(11,250) NDIS,NISO,MEV,N,SVSEED,MATNO
  250 FORMAT(/4I6,I12,I8,' = NDIS,NISO,MEV,N,SVSEED,MATNO'/)
C
      WRITE(11,260) (MP(I),VS(I),G(I), I=1,NDIS)
  260 FORMAT(2(I3,E25.16,E12.3))
C
      WRITE(11,270) NDIS, (MP(I), I=1,NDIS)
  270 FORMAT(/I6,' = NDIS, T-MULTIPLICITIES (0 MEANS  SPURIOUS)'/(20I4))
C
  280 CONTINUE
C
      IF (NISO.NE.0) GO TO 310
C
      WRITE(4,290) MEV
  290 FORMAT(/' AT MEV = ',I6,' THERE ARE NO ISOLATED T-EIGENVALUES'/
     1' SO NO ERROR ESTIMATES WERE COMPUTED/')
C
      WRITE(6,300)
  300 FORMAT(/' ALL COMPUTED GOOD T-EIGENVALUES ARE MULTIPLE'/
     1 ' THEREFORE ALL SUCH EIGENVALUES ARE ASSUMED TO HAVE CONVERGED')
C
      ICONV = 1
      GO TO 350
C
  310 CONTINUE
C
C-----------------------------------------------------------------------
C     SUBROUTINE INVERR COMPUTES ERROR ESTIMATES FOR ISOLATED GOOD
C     T-EIGENVALUES USING INVERSE ITERATION ON T(1,MEV). ON RETURN
C     G(J) = MINIMUM GAP IN T(1,MEV) FOR EACH VS(J), J=1,NDIS
C     G(MEV+I) = BETAM*|U(MEV)| = ERROR ESTIMATE FOR ISOLATED GOOD
C              T-EIGENVALUES, WHERE I = 1, NISO AND  BETAM = BETA(MEV+1)
C              U(MEV) IS MEVTH COMPONENT OF THE UNIT EIGENVECTOR OF T
C              CORRESPONDING TO THE ITH ISOLATED GOOD T-EIGENVALUE.
C     A NEGATIVE ERROR ESTIMATE MEANS THAT FOR THAT PARTICULAR
C     EIGENVALUE THE INVERSE ITERATION DID NOT CONVERGE IN <= MXINIT
C     STEPS AND THAT THE CORRESPONDING ERROR ESTIMATE IS QUESTIONABLE.
C
C     V2 CONTAINS THE ISOLATED GOOD T-EIGENVALUES
C     V1 CONTAINS THE MINGAPS TO THE NEAREST DISTINCT  EIGENVALUE
C        OF T(1,MEV) FOR EACH ISOLATED GOOD EIGENVALUE IN V2.
C     VS CONTAINS THE NDIS DISTINCT EIGENVALUES OF T(1,MEV)
C     MP CONTAINS THE CORRESPONDING CODED T-MULTIPLICITIES
C     WITHIN INVERR V1 AND V2 ARE DOUBLE PRECISION ARRAYS
C
      IT = MXINIT
      CALL INVERR(ALPHA,BETA,V1,V2,VS,EPSM,G,MP,MEV,MMB,NDIS,NISO,N,
     1 RHSEED,IT,IWRITE)
C
C-----------------------------------------------------------------------
C
C     SIMPLE CHECK FOR CONVERGENCE. CHECKS TO SEE IF ALL OF THE ERROR
C     ESTIMATES ARE SMALLER THAN CONTOL.
C     IF THIS TEST IS SATISFIED, THEN CONVERGENCE FLAG, ICONV IS SET
C     TO 1.  TYPICALLY ERROR ESTIMATES ARE VERY CONSERVATIVE.
C
      WRITE(6,320) CONTOL
  320 FORMAT(/' CONVERGENCE IS TESTED USING THE CONVERGENCE TOLERANCE',
     1E13.4/)
C
      II = MEV +1
      IF = MEV+NISO
      DO 330 I = II,IF
      IF (ABS(G(I)).GT.CONTOL) GO TO 350
  330 CONTINUE
      ICONV = 1
      MMB = NMEVS
C
      WRITE(6,340) CONTOL
  340 FORMAT(' ALL COMPUTED ERROR ESTIMATES WERE LESS THAN',E15.4/
     1 ' THEREFORE PROCEDURE TERMINATES'/)
C
  350 CONTINUE
C
C     IF CONVERGENCE IS INDICATED, THAT IS ICONV = 1 ,THEN
C     THE SUBROUTINE PRTEST IS CALLED TO CHECK FOR ANY CONVERGED
C     EIGENVALUES THAT HAVE BEEN MISLABELLED AS SPURIOUS BECAUSE
C     THE PROJECTION OF THEIR EIGENVECTOR(S) ON THE STARTING
C     VECTOR WERE TOO SMALL.
C     NUMERICAL TESTS INDICATE THAT SUCH EIGENVALUES ARE RARE.
C     IF FOR SOME REASON MANY OF THESE HIDDEN EIGENVALUES APPEAR
C     ON SOME RUN, YOU CAN BE CERTAIN THAT SOMETHING IS FOULED UP.
C
      IF (ICONV.EQ.0) GO TO 480
C
C-----------------------------------------------------------------------
C
      CALL PRTEST (ALPHA,BETA,VS,TKMAX,EPSM,RELTOL,SCALE3,SCALE4,
     1 MP,NDIS,MEV,IPROJ)
C
C-----------------------------------------------------------------------
C
      IF(IPROJ.EQ.0) GO TO 470
C
      IF(IDIST.EQ.1)  WRITE(11,360) IPROJ
  360 FORMAT(' SUBROUTINE PRTEST WANTS TO RELABEL',I6,' SPURIOUS EIGENVA
     1LUES'/' WE ACCEPT RELABELLING ONLY IF LAST COMPONENT OF T-EIGENVEC
     1TOR IS L.T. 1.D-10'/)
C
      IIX = RHSEED
C
C-----------------------------------------------------------------------
C
      CALL GENRAN(IIX,G,MEV)
C
C-----------------------------------------------------------------------
C
      ITEN = -10
      NISOM = NISO + MEV
      IWRITO = IWRITE
      IWRITE = 0
C
      DO 390 J = 1,NDIS
      IF(MP(J).NE.ITEN) GO TO 390
      T0 = VS(J)
C
C-----------------------------------------------------------------------
C
      IT = MXINIT
      CALL INVERM(ALPHA,BETA,V1,V2,T0,TEMP,T1,EPSM,G,MEV,IT,IWRITE)
C
C-----------------------------------------------------------------------
C
      IF(TEMP.LE.1.D-10) GO TO 380
C     ERROR ESTIMATE WAS NOT SMALL REJECT RELABELLING OF THIS EIGENVALUE
      IF(IDIST.EQ.1)  WRITE(11,370) J,T0,TEMP
  370 FORMAT(/' LAST COMPONENT FOR',I6,'TH EIGENVALUE',E20.12/' IS TOO L
     1ARGE = ',E15.6,' SO DO NOT ACCEPT PRTEST RELABELLING'/)
      MP(J) = 0
      IPROJ = IPROJ - 1
      GO TO 390
C     RELABELLING ACCEPTED
  380 NISOM = NISOM + 1
      G(NISOM) = BETAM*TEMP
  390 CONTINUE
      IWRITE = IWRITO
C
      IF(IPROJ.EQ.0) GO TO 430
      WRITE(6,400) IPROJ
  400 FORMAT(/I6,' T-EIGENVALUES WERE RECLASSIFIED AS GOOD.'/
     1' THESE ARE IDENTIFIED IN FILE 3 BY A T-MULTIPLICITY OF -10'/' USE
     2R SHOULD INSPECT EACH TO MAKE SURE NEIGHBORS HAVE CONVERGED'/)
C
      IF(IDIST.EQ.1)  WRITE(11,410) IPROJ
  410 FORMAT(/I6,' T-EIGENVALUES WERE RELABELLED AS GOOD'/
     1' BELOW IS CORRECTED T-MULTIPLICITY PATTERN'/)
C
      WRITE(6,420) NDIS, (MP(I), I=1,NDIS)
      IF(IDIST.EQ.1)  WRITE(11,420) NDIS, (MP(I), I=1,NDIS)
  420 FORMAT(/I6,' = NDIS, T-MULTIPLICITIES (0 MEANS  SPURIOUS)'/
     1 6X, ' (-10) MEANS SPURIOUS T-EIGENVALUE RELABELLED AS GOOD'/(20I4
     1))
C
C     RECALCULATE MINGAPS FOR DISTINCT T(1,MEV) EIGENVALUES.
  430 NM1 = NDIS - 1
      G(NDIS) = VS(NM1)-VS(NDIS)
      G(1) = VS(2)-VS(1)
C
      DO 440 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
  440 CONTINUE
      IF(IPROJ.EQ.0) GO TO 470
C     WRITE TO FILE 4 ERROR ESTIMATES FOR THOSE T-EIGENVALUES RELABELLED
      NGOOD = 0
      DO 450 J = 1,NDIS
      IF(MP(J).EQ.0) GO TO 450
      NGOOD = NGOOD + 1
      IF(MP(J).NE.ITEN) GO TO 450
      T0 = VS(J)
      NISO = NISO + 1
      NISOM = MEV + NISO
      WRITE(4,460) NGOOD,T0,G(NISOM),G(J)
  450 CONTINUE
  460 FORMAT(I10,E25.16,2E14.3)
C
  470 CONTINUE
C
C     WRITE THE GOOD T-EIGENVALUES TO FILE 3.  FIRST TRANSFER THEM
C     TO V2 AND THEIR T-MULTIPLICITIES TO THE CORRESPONDING POSITIONS
C     IN MP AND COMPUTE THE A-MINGAPS, THE MINIMAL GAPS BETWEEN THE
C     GOOD T-EIGENVALUES.  THESE GAPS WILL BE PUT IN THE ARRAY G.
C     SINCE G CURRENTLY CONTAINS THE MINIMAL GAPS BETWEEN THE DISTINCT
C     EIGENVALUES OF THE T-MATRIX, THESE GAPS WILL FIRST BE
C     TRANSFERRED TO GC.  NOTE THAT GC<0 MEANS THAT THAT MINIMAL GAP
C     IN THE T-MATRIX IS DUE TO A SPURIOUS T-EIGENVALUE.
C     ALL THIS INFORMATION IS PRINTED TO FILE 3
C
  480 CONTINUE
C
      NG = 0
      DO 490 I = 1,NDIS
      IF (MP(I).EQ.0) GO TO 490
      NG = NG+1
      MP(NG) = MP(I)
      GR(NG) = VS(I)
      TEMP = G(I)
      TEMP = DABS(TEMP)
      J = I+1
      IF (G(I).LT.ZERO) J = I-1
      IF (MP(J).EQ.0) TEMP = -TEMP
      GC(NG) = TEMP
  490 CONTINUE
C
      WRITE(6,500)MEV
  500 FORMAT(//' T-EIGENVALUE CALCULATION AT MEV = ',I6,'    IS COMPLETE
     1')
C
C     NG = NUMBER OF COMPUTED DISTINCT GOOD T-EIGENVALUES.  NEXT
C     GENERATE GAPS BETWEEN GOOD T-EIGENVALUES (AMINGAPS) AND PUT THEM
C     IN G.  G(J) < 0 MEANS THE AMINGAP IS DUE TO THE LEFT-HAND GAP.
C
      NGM1 = NG - 1
      G(NG) = GR(NGM1)-GR(NG)
      G(1) = GR(2)-GR(1)
C
      DO 510 J = 2,NGM1
      T0 = GR(J)-GR(J-1)
      T1 = GR(J+1)-GR(J)
      G(J) = T1
      IF (T0.LT.T1) G(J) = -T0
  510 CONTINUE
C
C     WRITE GOOD T-EIGENVALUES OUT TO FILE 3.
C
      WRITE(3,520)NG,NDIS,MEV,N,SVSEED,MATNO,MULTOL,IB,BTOL
  520 FORMAT(4I6,I12,I8,' = NG,NDIS,MEV,N,SVEED,MATNO'/
     1 E20.12,I6,E13.4,' = MUTOL,INDEX MINIMAL BETA,BTOL'/
     1' EV NO',2X,'TMULT',7X,'GOOD T-EIGENVALUE',7X,'TMINGAP',7X,'AMINGA
     1P')
C
      WRITE(3,530)(I,MP(I),GR(I),GC(I),G(I), I=1,NG)
  530 FORMAT(2I6,E25.16,2E14.3)
C
C     IF CONVERGENCE FLAG ICONV.NE.1 AND NUMBER OF T-MATRICES
C     CONSIDERED TO DATE IS LESS THAN NUMBER ALLOWED, INCREMENT MEV.
C     AND LOOP BACK TO 210 TO REPEAT COMPUTATIONS.  RESTORE BETA(MEV+1).
C
      BETA(MP1) = BETAM
C
      IF (MMB.LT.NMEVS.AND.ICONV.NE.1) GO TO 180
C
C     END OF LOOP ON DIFFERENT SIZE T-MATRICES ALLOWED.
C
  540 CONTINUE
C
      IF(ISTOP.EQ.0)  WRITE(6,550)
  550 FORMAT(/' T-MATRICES (ALPHA AND BETA) ARE NOW AVAILABLE, TERMINATE
     1')
      IF (IHIS.EQ.1.AND.KMAX.NE.MOLD) WRITE(1,560)
  560 FORMAT(/' ABOVE ARE THE FOLLOWING VECTORS '/
     1 '  ALPHA(I), I = 1,KMAX'/
     2 '  BETA(I), I = 1,KMAX+1'/
     3 ' FINAL TWO LANCZOS VECTORS OF ORDER N FOR I = KMAX,KMAX+1'/
     4 ' ALL ENTRIES IN THIS FILE HAVE FORMAT 4Z20 '/
     5 ' ----- END OF FILE 1 NEW ALPHA, BETA HISTORY---------------'///)
C
      IF (ISTOP.EQ.0) GO TO 640
C
      WRITE(3,570)
  570 FORMAT(/' ABOVE ARE COMPUTED GOOD T-EIGENVALUES'/
     1 ' NG = NUMBER OF GOOD T-EIGENVALUES COMPUTED'/
     2 ' NDIS = NUMBER OF COMPUTED DISTINCT EIGENVALUES OF T(1,MEV)'/
     3 ' N = ORDER OF A,  MATNO = MATRIX IDENT'/
     4 ' MULTOL = MULTIPLICITY TOLERANCE FOR T-EIGENVALUES IN BISEC'/
     4 ' TMULT IS THE T-MULTIPLICITY OF GOOD T-EIGENVALUE'/
     5 ' TMULT = -1 MEANS SPURIOUS T-EIGENVALUE TOO CLOSE'/
     6 ' DO NOT COMPUTE ERROR ESTIMATES FOR SUCH EIGENVALUES'/
     7 ' AMINGAP = MINIMAL GAP BETWEEN THE COMPUTED A-EIGENVALUES'/
     8 ' AMINGAP .LT. 0. MEANS MINIMAL GAP IS DUE TO LEFT-HAND GAP'/
     9 ' TMINGAP= MINIMAL GAP W.R.T.  DISTINCT EIGENVALUES IN T(1,MEV)'/
     1 ' TMINGAP .LT. 0. MEANS MINGAP IS DUE TO SPURIOUS EIGENVALUE'/
     2 ' ----- END OF FILE 3 GOODEIGENVALUES-----------------------'///)
C
      IF (IDIST.EQ.1) WRITE(11,580)
  580 FORMAT(/' ABOVE ARE THE DISTINCT EIGENVALUES OF T(1,MEV).'/
     2 ' THE FORMAT IS      T-MULTIPLICITY    EIGENVALUE   TMINGAP'/
     3 '        THIS FORMAT IS REPEATED TWICE ON EACH LINE.'/
     4 ' T-MULTIPLICITY = -1 MEANS THAT THE SUBROUTINE ISOEV HAS TAGGED'
     5 /'   THIS SIMPLE T-EIGENVALUE AS HAVING A VERY CLOSE SPURIOUS'/
     6 '   T-EIGENVALUE SO THAT NO ERROR ESTIMATE WILL BE COMPUTED'/
     7 '   FOR THAT EIGENVALUE IN SUBROUTINE INVERR.'/
     8 ' TMINGAP .LT. 0, TMINGAP IS DUE TO LEFT GAP .GT. 0, RIGHT GAP.'/
     9 ' EACH OF THE DISTINCT T-EIGENVALUE TABLES IS FOLLOWED'/
     9 ' BY THE T-MULTIPLICITY PATTERN.'/
     1 ' NDIS = NUMBER OF COMPUTED DISTINCT EIGENVALUES OF T(1,MEV).'/
     2 ' NG = NUMBER OF GOOD T-EIGENVALUES. '/
     3 ' NISO = NUMBER OF ISOLATED GOOD T-EIGENVALUES. '/
     4 ' NISO ALSO IS THE COUNT OF +1 ENTRIES IN T-MULTIPLICITY PATTERN.
     5 '/' ----- END OF FILE 4 DISTINCT T-EIGENVALUES----------------'//
     6 /)
C
      IF(NISO.NE.0)  WRITE(4,590)
  590 FORMAT(/' ABOVE ARE THE ERROR ESTIMATES OBTAINED FOR THE ISOLATED
     1GOOD EIGENVALUES'/
     1' OBTAINED VIA INVERSE ITERATION IN THE SUBROUTINE INVERR.'/
     1' ALL OTHER GOOD EIGENVALUES HAVE CONVERGED.'/
     2' ERROR ESTIMATE = BETAM*ABS(UM)'/
     2' WHERE BETAM = BETA(MEV+1) AND UM = U(MEV).'/
     3' U = UNIT EIGENVECTOR OF T WHERE T*U = EV*U AND EV = ISOLATED GOO
     3D EIGENVALUE.'/
     4' TMINGAP = GAP TO NEAREST DISTINCT EIGENVALUE OF T(1,MEV).'/
     5' TMINGAP .LT. 0. MEANS MINGAP IS DUE TO LEFT NEIGHBOR.'/
     6' ERROR ESTIMATE L.T. 0 MEANS INVERSE ITERATION DID NOT CONVERGE'/
     7' ------ END OF FILE 7 ERRINV -------------------------------'//)
      GO TO 640
C
  600 CONTINUE
C
      IBB = IABS(IBMEV)
      IF (IBMEV.LT.0) WRITE(6,610) MEV,IBB,BETA(IBB)
  610 FORMAT(/' PROGRAM TERMINATES BECAUSE MEV REQUESTED = ',I6,' IS .GT
     1',I6/' AT WHICH AN ABNORMALLY SMALL BETA = ' , E13.4,' OCCURRED'/)
      GO TO 640
C
  620 IF (NDIS.EQ.0.AND.ISTOP.GT.0) WRITE(6,630)
  630 FORMAT(/' INTERVALS SPECIFIED FOR BISECT DID NOT CONTAIN ANY EIGEN
     1VALUES'/' PROGRAM TERMINATES')
C
  640 CONTINUE
C
      STOP
C-----END OF MAIN PROGRAM FOR LANCZOS HERMITIAN EIGENVALUE COMPUTATIONS-
      END

.