C
C-----LESUB (1) REAL SYMMETRIC-----------------------------------------
C (2) HERMITIAN MATRICES
C (3) FACTORED INVERSES OF REAL SYMMETRIC MATRICES AND
C (4) REAL SYMMETRIC GENERALIZED, A*X = EVAL*B*X WHERE
C B IS POSITIVE DEFINITE, CHOLESKY FACTOR AVAILABLE
C
C ACCORDING TO PFORT THESE SUBROUTINES ARE PORTABLE EXCEPT FOR:
C (1) THE COMPLEX*16 VARIABLES AND THE CORRESPONDING FUNCTIONS
C FOR COMPLEX VARIABLES, DCMPLX, DREAL AND DCONJG USED IN
C THE SUBROUTINE CINPRD (USED ONLY IN CASE (2), HERMITIAN)
C (2) THE ENTRY IN THE SUBROUTINE LPERM USED TO PASS THE
C PERMUTATION FROM THE UPSEC SUBROUTINE TO LPERM. (USED
C ONLY IN CASES (3) AND (4), INVERSE AND GENERALIZED).
C
C SUBROUTINES BISEC, INVERR, TNORM, LUMP, ISOEV, PRTEST, AND
C INVERM ARE USED WITH THE LANCZOS EIGENVALUE
C PROGRAMS LEVAL, HLEVAL, LIVAL AND LGVAL. STURMI,
C INVERM, LBISEC, AND TNORM ARE USED WITH THE
C EIGENVECTOR PROGRAMS LEVEC, HLEVEC, LIVEC AND
C LGVEC. LPERM IS USED WITH LIVEC AND LGVEC.
C IN THE HERMITIAN CASE, THE SUBROUTINE CINPRD
C IS ALSO USED.
C
C-----COMPUTE T-EIGENVALUES BY BISECTION--------------------------------
C
SUBROUTINE BISEC(ALPHA,BETA,BETA2,VB,VS,LBD,UBD,EPS,TTOL,MP,
1 NINT,MEV,NDIS,IC,IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),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 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(|ALPHA(K)|,BETA(K), K=1,KMAX)
C
C ON EXIT
C NDIS = TOTAL NUMBER OF COMPUTED DISTINCT EIGENVALUES OF
C T(1,MEV) ON THE UNION OF THE (LB,UB) INTERVALS.
C VS = COMPUTED DISTINCT EIGENVALUES OF T(1,MEV) IN ALGEBRAICALLY-
C INCREASING ORDER
C MP = CORRESPONDING T-MULTIPLICITIES OF THESE EIGENVALUES
C MP(I) = (0,1,MI), MI>1, I=1,NDIS MEANS:
C (0) V(I) IS SPURIOUS
C (1) V(I) IS ISOLATED AND GOOD
C (MI) V(I) IS 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 MULTIPLE AND HENCE 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 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 T-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 T-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 - ALPHA(I) - 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-ALPHA(I)-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) EIGENVALUES) > X1
C JEV = -NA + (NUMBER OF T(2,MEV) EIGENVALUES) > X1
C SET PARAMETERS FOR NEXT INTERVAL
410 CONTINUE
IC = ICT+IC
MXSTUR = MXSTUR-ICT
C
WRITE(6,420) JIND,NG,MD
420 FORMAT(/' T-EIGENVALUE CALCULATION ON INTERVAL',I6,' IS COMPLETE'
1 /3X,'NO. GOOD',3X,'NO. DISTINCT'/I10,I13)
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(ALPHA,BETA,V1,V2,VS,EPS,G,MP,MEV,MMB,NDIS,NISO,
1 N,IKL,IT,IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),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)
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 T-NEIGHBOR IS ALSO GOOD, OR ITS CLOSEST NEIGHBOR 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 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 ALPHA, BETA CONTAIN THE NONZERO ENTRIES OF THE T-MATRIX
C VS = COMPUTED DISTINCT EIGENVALUES OF T(1,MEV)
C MP = T-MULTIPLICITY OF EACH T-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 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 T-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 = DABS(ALPHA(1))
DO 30 I = 2,MEV
30 TSUM = TSUM + DABS(ALPHA(I)) + 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 EIGENVECTOR OF T(1,MEV)
C
DO 180 JEV = 1,NDIS
C
IF (MP(JEV).EQ.0) GO TO 180
NG = NG + 1
IF (MP(JEV).NE.1) GO TO 180
C
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 = ALPHA(1)-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 = ALPHA(I)-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) = ALPHA(I)-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,N,IKL,ITER,BETAM
210 FORMAT(1X,'TSIZE',2X,'NDIS',1X,'NGOOD',2X,'NISO',1X,'ASIZE'/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
C RESTORE BETA(MEV+1) = BETAM
270 BETA(MP1) = BETAM
C-----END OF INVERR-----------------------------------------------------
RETURN
END
C
C-----START OF TNORM----------------------------------------------------
C
SUBROUTINE TNORM(ALPHA,BETA,BMIN,TMAX,MEV,IB)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),BETA(1)
DOUBLE PRECISION TMAX,BMIN,BMAX,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 (|ALPHA(I)|, 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)
BMAX = BETA(2)
TMAX = DABS(ALPHA(1))
C
DO 20 I = 2,MEV
IF (BETA(I).GE.BMIN) GO TO 10
IB = I
BMIN = BETA(I)
10 TMAX = DMAX1(TMAX,DABS(ALPHA(I)))
BMAX = DMAX1(BETA(I),BMAX)
20 CONTINUE
TMAX = DMAX1(BMAX,TMAX)
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
C
C-----START OF LUMP-----------------------------------------------------
C
SUBROUTINE LUMP(V1,RELTOL,MULTOL,SCALE2,LINDEX,LOOP)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION V1(1),SUM,RELTOL,MULTOL,THOLD,ZERO,SCALE2
INTEGER LINDEX(1)
DOUBLE PRECISION DABS, DFLOAT, DMAX1
C-----------------------------------------------------------------------
C LINDEX(J) = T-MULTIPLICITY OF JTH DISTINCT T-EIGENVALUE
C LOOP = NUMBER OF DISTINCT T-EIGENVALUES
C LUMP 'COMBINES' COMPUTED 'GOOD' T-EIGENVALUES THAT ARE
C 'TOO CLOSE'.
C VALUE OF RELTOL IS 1.D-10.
C
C IF IN A SET OF T-EIGENVALUES TO BE COMBINED THERE IS AN EIGENVALUE
C WITH LINDEX=1, THEN THE VALUE OF THE COMBINED EIGENVALUES IS SET
C EQUAL TO THE VALUE OF THAT EIGENVALUE. NOTE THAT IF A SPURIOUS
C T-EIGENVALUE IS TO BE 'COMBINED' WITH A GOOD T-EIGENVALUE, THEN
C THIS IS DONE ONLY BY INCREASING THE INDEX, LINDEX, FOR THAT
C T-EIGENVALUE. NUMERICAL VALUES OF SPURIOUS EIGENVALUES ARE NEVER
C COMBINE WITH THOSE OF GOOD T-EIGENVALUES.
C-----------------------------------------------------------------------
ZERO = 0.0D0
NLOOP = 0
J = 0
ICOUNT = 1
JI = 1
THOLD = DMAX1(RELTOL*DABS(V1(1)),SCALE2*MULTOL)
C THOLD = DMAX1(RELTOL*DABS(V1(1)),RELTOL)
C
10 J = J+1
IF (J.EQ.LOOP) GO TO 20
SUM = DABS(V1(J)-V1(J+1))
IF (SUM.LT.THOLD) GO TO 60
20 JF = JI + ICOUNT - 1
INDSUM = 0
ISPUR = 0
C
DO 30 KK = JI,JF
IF (LINDEX(KK).NE.0) GO TO 30
ISPUR = ISPUR + 1
INDSUM = INDSUM + 1
30 INDSUM = INDSUM + LINDEX(KK)
C
C IF (JF-JI.GE.1) WRITE(6,40) (V1(KKK), KKK=JI,JF)
40 FORMAT(/' LUMP LUMPS THE T-EIGENVALUES'/(4E20.13))
C
C COMPUTE THE 'COMBINED' T-EIGENVALUE AND THE RESULTING
C T-MULTIPLICITY
K = JI - 1
50 K = K+1
IF (K.GT.JF) GO TO 70
IF (LINDEX(K) .NE.1) GO TO 50
NLOOP = NLOOP + 1
V1(NLOOP) = V1(K)
GO TO 100
60 ICOUNT = ICOUNT + 1
GO TO 10
C
C ALL INDICES WERE 0 OR >1
70 NLOOP = NLOOP + 1
IDIF = INDSUM - ISPUR
IF (IDIF.EQ.0) GO TO 90
C
SUM = ZERO
DO 80 KK = JI,JF
80 SUM = SUM + V1(KK) * DFLOAT(LINDEX(KK))
C
V1(NLOOP) = SUM/DFLOAT(IDIF)
GO TO 100
90 V1(NLOOP) = V1(JI)
100 LINDEX(NLOOP) = INDSUM
IDIF = INDSUM - ISPUR
IF (IDIF.EQ.0.AND.ISPUR.EQ.1) LINDEX(NLOOP) = 0
IF (J.EQ.LOOP) GO TO 110
ICOUNT = 1
JI= J+1
THOLD = DMAX1(RELTOL*DABS(V1(JI)),SCALE2*MULTOL)
C THOLD = DMAX1(RELTOL*DABS(V1(JI)),RELTOL)
IF (JI.LT.LOOP) GO TO 10
NLOOP = NLOOP + 1
V1(NLOOP)= V1(JI)
LINDEX(NLOOP) = LINDEX(JI)
110 CONTINUE
C
C ON RETURN V1 CONTAINS THE DISTINCT T-EIGENVALUES
C LINDEX CONTAINS THE CORRESPONDING T-MULTIPLICITIES
C
LOOP = NLOOP
RETURN
C-----END OF LUMP-------------------------------------------------------
END
C
C
C-----START OF ISOEV----------------------------------------------------
C
SUBROUTINE ISOEV(VS,GAPTOL,MULTOL,SCALE1,G,MP,NDIS,NG,NISO)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION VS(1),T0,T1,MULTOL,GAPTOL,SCALE1,TEMP
REAL G(1),GAP
INTEGER MP(1)
REAL ABS
DOUBLE PRECISION DABS, DMAX1
C-----------------------------------------------------------------------
C GENERATE DISTINCT TMINGAPS AND USE THEM TO LABEL THE ISOLATED
C GOOD T-EIGENVALUES THAT ARE VERY CLOSE TO SPURIOUS ONES.
C ERROR ESTIMATES WILL NOT BE COMPUTED FOR THESE T-EIGENVALUES.
C
C ON ENTRY AND EXIT
C VS CONTAINS THE COMPUTED DISTINCT EIGENVALUES OF T(1,MEV)
C MP CONTAINS THE CORRESPONDING T-MULTIPLICITIES
C NDIS = NUMBER OF DISTINCT EIGENVALUES
C GAPTOL = RELATIVE GAP TOLERANCE SET IN MAIN
C
C ON EXIT
C G CONTAINS THE TMINGAPS.
C G(I) < 0 MEANS MINGAP IS DUE TO LEFT GAP
C MP(I) IS NOT CHANGED EXCEPT THAT MP(I)=-1, IF MP(I)=1,
C TMINGAP WAS TOO SMALL AND DUE TO A SPURIOUS T-EIGENVALUE.
C
C IF MP(I)=-1 THAT SIMPLE GOOD T-EIGENVALUE WILL BE SKIPPED
C IN THE SUBSEQUENT ERROR ESTIMATE COMPUTATIONS IN INVERR
C THAT IS, WE COMPUTE ERROR ESTIMATES ONLY FOR THOSE GOOD
C T-EIGENVALUES WITH MP(I)=1.
C-----------------------------------------------------------------------
C CALCULATE MINGAPS FOR DISTINCT T(1,MEV) EIGENVALUES.
NM1 = NDIS - 1
G(NDIS) = VS(NM1)-VS(NDIS)
G(1) = VS(2)-VS(1)
C
DO 10 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
10 CONTINUE
C
C SET MP(I)=-1 FOR SIMPLE GOOD T-EIGENVALUES WHOSE MINGAPS ARE
C 'TOO SMALL' AND DUE TO SPURIOUS T-EIGENVALUES.
C
NISO = 0
NG = 0
DO 20 J = 1,NDIS
IF (MP(J).EQ.0) GO TO 20
NG = NG+1
IF (MP(J).NE.1) GO TO 20
C VS(J) IS NEXT SIMPLE GOOD T-EIGENVALUE
NISO = NISO + 1
I = J+1
IF (G(J).LT.0.0) I = J-1
IF (MP(I).NE.0) GO TO 20
GAP = ABS(G(J))
T0 = DMAX1(SCALE1*MULTOL,GAPTOL*DABS(VS(J)))
C T0 = DMAX1(GAPTOL,GAPTOL*DABS(VS(J)))
TEMP = T0
IF (GAP.GT.TEMP) GO TO 20
MP(J) = -MP(J)
NISO = NISO-1
20 CONTINUE
C
C-----END OF ISOEV------------------------------------------------------
RETURN
END
C
C-----START OF PRTEST---------------------------------------------------
C
SUBROUTINE PRTEST(ALPHA,BETA,TEIG,TKMAX,EPSM,RELTOL,SCALE3,SCALE4,
1 TMULT,NDIST,MEV,IPROJ)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1), BETA(1),TEIG(1),SIGMA(10)
DOUBLE PRECISION EPSM,RELTOL,PRTOL,TKMAX,LRATIO,URATIO
DOUBLE PRECISION EPS,EPS1,BETAM,LBD,UBD,SIG,YU,YV,LRATS,URATS
DOUBLE PRECISION ZERO,ONE,TEN,BISTOL,SCALE3,SCALE4,AEV,TEMP
INTEGER TMULT(1),ISIGMA(10)
DOUBLE PRECISION DABS, DMAX1, DSQRT, DFLOAT
C-----------------------------------------------------------------------
C AFTER CONVERGENCE HAS BEEN ESTABLISHED, SUBROUTINE PRTEST
C TESTS COMPUTED EIGENVALUES OF T(1,MEV) THAT HAVE BEEN LABELLED
C SPURIOUS TO DETERMINE IF ANY EIGENVALUES OF A HAVE BEEN
C MISSED BY LANCZOS PROCEDURE. AN EIGENVALUE WITH A VERY SMALL
C PROJECTION ON THE STARTING VECTOR (< SINGLE PRECISION)
C CAN BE MISSED BECAUSE IT IS ALSO AN EIGENVALUE OF T(2,MEV) TO
C WITHIN THE SQUARE OF THIS ORIGINAL PROJECTION.
C OUR EXPERIENCE IS THAT SUCH SMALL PROJECTIONS OCCUR ONLY
C VERY INFREQUENTLY.
C
C THIS SUBROUTINE IS CALLED ONLY AFTER CONVERGENCE HAS BEEN
C ESTABLISHED. ONCE CONVERGENCE HAS BEEN OBSERVED ON THE
C OTHER EIGENVALUES THEN ONE CAN EXPECT TO ALSO HAVE CONVERGENCE
C ON ANY SUCH HIDDEN EIGENVALUES.(IF THERE ARE ANY). THIS
C PROCEDURE CONSIDERS ONLY SPURIOUS T-EIGENVALUES AND ONLY THOSE
C SPURIOUS T-EIGENVALUES THAT ARE ISOLATED FROM GOOD T-EIGENVALUES.
C FOR EACH SUCH T-EIGENVALUE IT DOES 2 STURM SEQUENCES
C AND A FEW SCALAR MULTIPLICATIONS. UPON RETURN TO MAIN
C PROGRAM ERROR ESTIMATES WILL BE COMPUTED FOR ANY EIGENVALUES
C THAT HAVE BEEN LABELLED AS 'HIDDEN'. SUCH T-EIGENVALUES
C WILL BE RELABELLED AS 'GOOD' ONLY IF THESE ERROR ESTIMATES
C ARE SUFFICIENTLY SMALL.
C-----------------------------------------------------------------------
ZERO = 0.0D0
ONE = 1.0D0
TEN = 10.0D0
PRTOL = 1.D-6
TEMP = DFLOAT(MEV+1000)
TEMP = DSQRT(TEMP)
BISTOL = TKMAX*EPSM*TEMP
NSIGMA = 4
SIGMA(1) = TEN*TKMAX
C
DO 10 J = 2,NSIGMA
10 SIGMA(J) = TEN*SIGMA(J-1)
C
IFIN = 0
MF = 1
ML = MEV
BETAM = BETA(MF)
BETA(MF) = ZERO
IPROJ = 0
J = 1
C
IF (TMULT(1).NE.0) GO TO 110
C
AEV = DABS(TEIG(1))
TEMP = PRTOL*AEV
EPS1 = DMAX1(TEMP,SCALE4*BISTOL)
C EPS1 = DMAX1(TEMP,PRTOL)
TEMP = RELTOL*AEV
EPS = DMAX1(TEMP,SCALE3*BISTOL)
C EPS = DMAX1(TEMP,RELTOL)
C
IF (TEIG(2)-TEIG(1).LT.EPS1.AND.TMULT(2).NE.0) GO TO 110
C
20 LBD = TEIG(J) - EPS
UBD = TEIG(J) + EPS
MEVL = 0
IL = 0
YU = ONE
C
DO 50 I=MF,ML
IF (YU.NE.ZERO) GO TO 30
YV = BETA(I)/EPSM
GO TO 40
30 YV = BETA(I)*BETA(I)/YU
40 YU = ALPHA(I)-LBD-YV
IF (YU.GE.ZERO) GO TO 50
C MEVL INCREMENTED
MEVL = MEVL + 1
IL = I
50 CONTINUE
C
LRATIO = YU
MEV1L = MEVL
IF (IL.EQ.ML) MEV1L=MEVL-1
C
C MEVL = NUMBER OF EVS OF T(1,MEV) WHICH ARE < LBD
C MEV1L = NUMBER OF EVS OF T(1,MEV-1) WHICH ARE < LBD
C LRATIO = DET(T(1,MEV)-LBD)/DET(T(1,MEV-1)-LBD):
C
MEVU = 0
IL = 0
YU = ONE
C
DO 80 I=MF,ML
IF (YU.NE.ZERO) GO TO 60
YV = BETA(I)/EPSM
GO TO 70
60 YV = BETA(I)*BETA(I)/YU
70 YU = ALPHA(I)-UBD-YV
IF (YU.GE.ZERO) GO TO 80
C MEVU INCREMENTED
MEVU = MEVU + 1
IL = I
80 CONTINUE
C
URATIO = YU
MEV1U = MEVU
IF (IL.EQ.ML) MEV1U=MEVU-1
C
C MEVU = NUMBER OF EVS OF T(MEV) WHICH ARE < UBD
C MEV1U = NUMBER OF EVS OF T(MEV-1) WHICH ARE < UBD
C URATIO = DET(TM-UBD)/DET(T(M-1)-UBD): TM=T(MF,ML)
C
NEV1 = MEV1U-MEV1L
C
DO 90 K=1,NSIGMA
SIG = SIGMA(K)
LRATS = LRATIO-SIG
URATS = URATIO-SIG
C NOTE THE INCREMENT IS ON NUMBER OF EVALUES OF T(M-1)
MEVLS = MEV1L
IF (LRATS.LT.0.) MEVLS=MEV1L+1
MEVUS = MEV1U
IF (URATS.LT.0.) MEVUS=MEV1U+1
ISIGMA(K) = MEVUS - MEVLS
90 CONTINUE
C
ICOUNT = 0
DO 100 K=1,NSIGMA
100 IF (ISIGMA(K).EQ.1) ICOUNT=ICOUNT + 1
C
IF (ICOUNT.LT.2.OR.NEV1.EQ.0) GO TO 110
TMULT(J) = -10
IPROJ=IPROJ+1
C
110 J=J+1
C
IF (J.GE.NDIST) GO TO 120
IF (TMULT(J).NE.0) GO TO 110
C
AEV = DABS(TEIG(J))
TEMP = PRTOL*AEV
EPS1 = DMAX1(TEMP,SCALE4*BISTOL)
C EPS1 = DMAX1(TEMP,PRTOL)
TEMP = RELTOL*AEV
EPS = DMAX1(TEMP,SCALE3*BISTOL)
C EPS = DMAX1(TEMP,RELTOL)
C
IF (TEIG(J)-TEIG(J-1).LT.EPS1.AND.TMULT(J-1).NE.0) GO TO 110
IF (TEIG(J+1)-TEIG(J).LT.EPS1.AND.TMULT(J+1).NE.0) GO TO 110
C
GO TO 20
C
120 IF (IFIN.EQ.1) GO TO 130
IF (TMULT(NDIST).NE.0) GO TO 130
C
AEV = DABS(TEIG(NDIST))
TEMP = PRTOL*AEV
EPS1 = DMAX1(TEMP,SCALE4*BISTOL)
C EPS1 = DMAX1(TEMP,PRTOL)
TEMP = RELTOL*AEV
EPS = DMAX1(TEMP,SCALE3*BISTOL)
C EPS = DMAX1(TEMP,RELTOL)
C
NDIST1=NDIST -1
TEMP = TEIG(NDIST)-TEIG(NDIST1)
IF (TEMP.LT.EPS1.AND.TMULT(NDIST1).NE.0) GO TO 130
IFIN = 1
C
GO TO 20
C
130 BETA(MF) = BETAM
C
C-----END OF PRTEST-----------------------------------------------------
RETURN
END
C
C------START OF STURMI--------------------------------------------------
C
SUBROUTINE STURMI(ALPHA,BETA,X1,TOLN,EPSM,MMAX,MK1,MK2,IC,IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),BETA(1)
DOUBLE PRECISION EPSM,X1,TOLN,EVL,EVU,BETA2
DOUBLE PRECISION U1,U2,V1,V2,ZERO,ONE
INTEGER I,IC,ICD,IC0,IC1,IC2,MK1,MK2,MMAX
C-----------------------------------------------------------------------
C
C FOR ANY EIGENVALUE OF A THAT HAS CONVERGED AS AN EIGENVALUE
C OF THE T-MATRICES THIS SUBROUTINE CALCULATES
C THE SMALLEST SIZE OF THE T-MATRIX, T(1,MK1) DEFINED
C BY THE ALPHA AND BETA ARRAYS SUCH THAT MK1.LE.MMAX
C AND THE INTERVAL (X1-TOLN,X1+TOLN) CONTAINS AT LEAST ONE
C EIGENVALUE OF T(1,MK1). IT ALSO CALCULATES MK2 <= MMAX
C AS THE SMALLEST SIZE T-MATRIX (IF ANY) SUCH THAT THIS INTERVAL
C CONTAINS AT LEAST TWO EIGENVALUES OF T(1,MK2).
C IF NO T-MATRIX OF ORDER < MMAX SATISFIES THIS REQUIREMENT
C THEN MK2 IS SET EQUAL TO MMAX. THE EIGENVECTOR PROGRAM
C USES THESE VALUES TO DETERMINE AN APPROPRIATE 1ST GUESS AT
C AN APPROPRIATE SIZE T-MATRIX FOR THE EIGENVALUE X1.
C
C ON EXIT IC = NUMBER OF EIGENVALUES OF T(1,MK2) IN THIS INTERVAL
C
C STURMI REGENERATES THE QUANTITIES BETA(I)**2 EACH TIME IT IS
C CALLED, OBVIOUSLY FOR THE PRICE OF ANOTHER VECTOR OF LENGTH
C MMAX THIS GENERATION COULD BE DONE ONCE IN THE MAIN
C PROGRAM BEFORE THE LOOP ON THE CALLS TO SUBROUTINE STURMI.
C
C IF ANY OF THE EIGENVALUES BEING CONSIDERED WERE MULTIPLE
C AS EIGENVALUES OF THE USER-SPECIFIED MATRIX, THEN
C THIS SUBROUTINE COULD BE MODIFIED TO COMPUTE ADDITIONAL
C SIZES MKJ, J = 3, ... WHICH COULD THEN BE USED IN THE
C MAIN LANCZOS EIGENVECTOR PROGRAM TO COMPUTE ADDITIONAL
C EIGENVECTORS CORRESPONDING TO THESE MULTIPLE EIGENVALUES.
C THE MAIN PROGRAM PROVIDED DOES NOT INCLUDE THIS OPTION.
C
C-----------------------------------------------------------------------
C INITIALIZATION OF PARAMETERS
MK1 = 0
MK2 = 0
ZERO = 0.0D0
ONE = 1.0D0
BETA(1) = ZERO
EVL = X1-TOLN
EVU = X1+TOLN
U1 = ONE
U2 = ONE
IC0 = 0
IC1 = 0
IC2 = 0
C
C MAIN LOOP FOR CALCULATING THE SIZES MK1,MK2
DO 60 I = 1,MMAX
BETA2 = BETA(I)*BETA(I)
IF (U1.NE.ZERO) GO TO 10
V1 = BETA(I)/EPSM
GO TO 20
10 V1 = BETA2/U1
20 U1 = EVL - ALPHA(I) - V1
IF (U1.LT.ZERO) IC1 = IC1+1
IF (U2.NE.ZERO) GO TO 30
V2 = BETA(I)/EPSM
GO TO 40
30 V2 = BETA2/U2
40 U2 = EVU - ALPHA(I) - V2
IF (U2.LT.ZERO) IC2 = IC2+1
C TEST FOR CHANGE IN NUMBER OF T-EIGENVALUES ON (EVL,EVU)
ICD = IC1-IC2
IC = ICD-IC0
IF (IC.GE.1) GO TO 50
GO TO 60
50 CONTINUE
IF (IC0.EQ.0) MK1 = I
IC0 = IC0+1
IF (IC0.GT.1) GO TO 70
60 CONTINUE
C
I = I-1
IF (IC0.EQ.0) MK1 = MMAX
70 MK2 = I
IC = ICD
C
IF (IWRITE.EQ.1) WRITE(6,80) X1,MK1,MK2,IC
80 FORMAT(' EVAL =',E20.12,' MK1 =',I6,' MK2 =',I6,' IC =',I3/)
C
RETURN
C-----END OF STURMI-----------------------------------------------------
END
C
C
C-----START OF INVERM---------------------------------------------------
C
SUBROUTINE INVERM(ALPHA,BETA,V1,V2,X1,ERROR,ERRORV,EPS,G,MEV,IT,
1 IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),BETA(1),V1(1),V2(1)
DOUBLE PRECISION X1,U,Z,TEMP,RATIO,SUM,XU,NORM,TSUM,BETAM
DOUBLE PRECISION EPS,EPS3,EPS4,ERROR,ERRORV,ZERO,ONE
REAL G(1)
DOUBLE PRECISION DABS, DSQRT, DFLOAT
DOUBLE PRECISION FINPRO
REAL ABS
C-----------------------------------------------------------------------
C
C COMPUTES T-EIGENVECTORS FOR ISOLATED GOOD T-EIGENVALUES X1
C USING INVERSE ITERATION ON T(1,MEV(X1)) SOLVING EQUATION
C (T - X1*I)V2 = RIGHT-HAND SIDE (RANDOMLY-GENERATED) .
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 IF IWRITE = 1 THEN THERE ARE EXTENDED WRITES TO FILE 6 (TERMINAL)
C
C ON ENTRY G CONTAINS A REAL*4 RANDOM VECTOR WHICH WAS GENERATED
C IN MAIN PROGRAM.
C
C ON ENTRY AND EXIT
C MEV = ORDER OF T
C ALPHA, BETA CONTAIN THE DIAGONAL AND OFFDIAGONAL ENTRIES OF T.
C EPS = 2. * MACHINE EPSILON
C
C IN PROGRAM:
C ITER = MAXIMUM NUMBER STEPS ALLOWED FOR INVERSE ITERATION
C ITER = IT ON ENTRY.
C V1,V2 = WORK SPACES USED IN THE FACTORIZATION OF T(1,MEV).
C V1 AND V2 MUST BE OF DIMENSION AT LEAST MEV.
C
C ON EXIT
C V2 = THE UNIT EIGENVECTOR OF T(1,MEV) CORRESPONDING TO X1.
C ERROR = |V2(MEV)| = ERROR ESTIMATE FOR CORRESPONDING
C RITZ VECTOR FOR X1.
C
C ERRORV = || T*V2 - X1*V2 || = ERROR ESTIMATE ON T-EIGENVECTOR.
C IF IT.GT.ITER THEN ERRORV = -ERRORV
C IT = NUMBER OF ITERATIONS ACTUALLY REQUIRED
C-----------------------------------------------------------------------
C INITIALIZATION AND PARAMETER SPECIFICATION
ONE = 1.0D0
ZERO = 0.0D0
ITER = IT
MP1 = MEV+1
MM1 = MEV-1
BETAM = BETA(MP1)
BETA(MP1) = ZERO
C
C CALCULATE SCALE AND TOLERANCES
TSUM = DABS(ALPHA(1))
DO 10 I = 2,MEV
10 TSUM = TSUM + DABS(ALPHA(I)) + BETA(I)
C
EPS3 = EPS*TSUM
EPS4 = DFLOAT(MEV)*EPS3
C
C GENERATE SCALED RANDOM RIGHT-HAND SIDE
GSUM = ZERO
DO 20 I = 1,MEV
20 GSUM = GSUM+ABS(G(I))
GSUM = EPS4/GSUM
C
C INITIALIZE RIGHT HAND SIDE FOR INVERSE ITERATION
DO 30 I = 1,MEV
30 V2(I) = GSUM*G(I)
IT = 1
C
C CALCULATE UNIT EIGENVECTOR OF T(1,MEV) FOR ISOLATED GOOD
C T-EIGENVALUE X1.
C
C TRIANGULAR FACTORIZATION WITH NEAREST NEIGHBOR PIVOT
C STRATEGY. INTERCHANGES ARE LABELLED BY SETTING BETA < 0.
C
40 CONTINUE
U = ALPHA(1)-X1
Z = BETA(2)
C
DO 60 I=2,MEV
IF (BETA(I).GT.DABS(U)) GO TO 50
C NO PIVOT 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 = ALPHA(I)-X1-Z*RATIO
Z = BETA(I+1)
GO TO 60
C PIVOT INTERCHANGE
50 CONTINUE
RATIO = U/BETA(I)
BETA(I) = -BETA(I)
V1(I-1) = ALPHA(I)-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)
60 CONTINUE
C
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 80 II = 1,MM1
I = MEV-II
IF (BETA(I+1).LT.ZERO) GO TO 70
C NO PIVOT INTERCHANGE
V2(I) = V2(I)-V1(I)*V2(I+1)
GO TO 80
C PIVOT INTERCHANGE
70 BETA(I+1) = -BETA(I+1)
V2(I) = (V2(I)-V1(I)*V2(I+1)-BETA(I+2)*V2(I+2))/BETA(I+1)
80 CONTINUE
C
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 90 II = 1,MM1
I = MEV-II
90 NORM = NORM+DABS(V2(I))
C
C IS DESIRED GROWTH IN VECTOR ACHIEVED ?
C IF NOT, DO ANOTHER INVERSE ITERATION STEP UNLESS NUMBER ALLOWED IS
C EXCEEDED.
IF (NORM.GE.ONE) GO TO 110
C
IT=IT+1
IF (IT.GT.ITER) GO TO 110
C
XU = EPS4/NORM
DO 100 I=1,MEV
100 V2(I) = V2(I)*XU
C
GO TO 40
C
C NORMALIZE COMPUTED T-EIGENVECTOR : V2 = V2/||V2||
C
110 CONTINUE
C
SUM = FINPRO(MEV,V2(1),1,V2(1),1)
SUM = ONE/DSQRT(SUM)
DO 120 II = 1,MEV
120 V2(II) = SUM*V2(II)
C
C SAVE ERROR ESTIMATE FOR LATER OUTPUT
ERROR = DABS(V2(MEV))
C
C GENERATE ERRORV = ||T*V2 - X1*V2||.
V1(MEV) = ALPHA(MEV)*V2(MEV)+BETA(MEV)*V2(MEV-1)-X1*V2(MEV)
DO 130 J = 2,MM1
JM = MP1 - J
V1(JM) = ALPHA(JM)*V2(JM) + BETA(JM)*V2(JM-1) + BETA(JM+1)*V2(JM+1
1) - X1*V2(JM)
130 CONTINUE
C
V1(1) = ALPHA(1)*V2(1) + BETA(2)*V2(2) - X1*V2(1)
ERRORV = FINPRO(MEV,V1(1),1,V1(1),1)
ERRORV = DSQRT(ERRORV)
IF (IT.GT.ITER) ERRORV = -ERRORV
IF (IWRITE.EQ.0) GO TO 150
C
C FILE 6 (TERMINAL) OUTPUT OF ERROR ESTIMATES.
WRITE(6,140) MEV,X1,ERROR,ERRORV
140 FORMAT(2X,'TSIZE',15X,'EIGENVALUE',11X,'U(M)',9X,'ERRORV'/
1 I6,E25.16,2E15.5)
C
C RESTORE BETA(MEV+1) = BETAM
150 CONTINUE
BETA(MP1) = BETAM
C-----END OF INVERM-----------------------------------------------------
RETURN
END
C
C-----START OF LBISEC---------------------------------------------------
C
SUBROUTINE LBISEC(ALPHA,BETA,EPSM,EVAL,EVALN,LB,UB,TTOL,M,NEVT)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION ALPHA(1),BETA(1),X0,X1,XL,XU,YU,YV,LB,UB
DOUBLE PRECISION EPSM,EP1,EVAL,EVALN,EVD,EPT
DOUBLE PRECISION ZERO,ONE,HALF,TTOL,TEMP
DOUBLE PRECISION DABS,DSQRT,DFLOAT
C-----------------------------------------------------------------------
C SPECIFY PARAMETERS
ZERO = 0.0D0
HALF = 0.5D0
ONE = 1.0D0
XL = LB
XU = UB
C
C EP1 = DSQRT(1000+M)*TTOL TTOL = EPSM*TKMAX
C TKMAX = MAX(|ALPHA(K)|,BETA(K), K= 1,KMAX)
C
TEMP = DFLOAT(1000+M)
EP1 = DSQRT(TEMP)*TTOL
C
NA = 0
X1 = XU
JSTURM = 1
GO TO 60
C FORWARD STURM CALCULATION
10 NA = NEV
X1 = XL
JSTURM = 2
GO TO 60
C FORWARD STURM CALCULATION
20 NEVT = NEV
C
C WRITE(6,30) M,EVAL,NEVT,EP1
30 FORMAT(/3X,'TSIZE',23X,'EV',9X/I8,E25.16/
1 I6,' = NUMBER OF T(1,M) EIGENVALUES ON TEST INTERVAL'/
1 E12.3,' = CONVERGENCE TOLERANCE'/)
C
IF (NEVT.NE.1) GO TO 120
C
C BISECTION LOOP
JSTURM = 3
40 X1 = HALF*(XL+XU)
X0 = XU-XL
EPT = EPSM*(DABS(XL) + DABS(XU)) + EP1
C CONVERGENCE TEST
IF (X0.LE.EP1) GO TO 100
GO TO 60
C FORWARD STURM CALCULATION
50 CONTINUE
IF(NEV.EQ.0) XU = X1
IF(NEV.EQ.1) XL = X1
GO TO 40
C NEV = NUMBER OF T-EIGENVALUES OF T(1,M) ON (X1,XU)
C THERE IS EXACTLY ONE T-EIGENVALUE OF T(1,M) ON (XL,XU)
C
C FORWARD STURM CALCULATION
60 NEV = -NA
YU = ONE
DO 90 I = 1,M
IF (YU.NE.ZERO) GO TO 70
YV = BETA(I)/EPSM
GO TO 80
70 YV = BETA(I)*BETA(I)/YU
80 YU = X1 - ALPHA(I) - YV
IF (YU.GE.ZERO) GO TO 90
NEV = NEV+1
90 CONTINUE
GO TO (10,20,50), JSTURM
C
100 CONTINUE
C
EVALN = X1
EVD = DABS(EVALN-EVAL)
C WRITE(6,110) EVALN,EVAL,EVD
110 FORMAT(/20X,'EVALN',21X,'EVAL',6X,'CHANGE'/2E25.16,E12.3/)
C
120 CONTINUE
RETURN
C-----END OF LBISEC-----------------------------------------------------
END
C-----START OF COMPLEX INNER PRODUCT------------------------------------
C
C COMPLEX INNER PRODUCT
C
SUBROUTINE CINPRD(V2,V1,SUM,N)
C-----------------------------------------------------------------------
DOUBLE PRECISION ZERO,SUM
COMPLEX*16 V2(1),V1(1),SUMC
C-----------------------------------------------------------------------
C
C NOTE THAT THE ORDER MATTERS HERE
C COMPUTES THE INNER PRODUCT OF THE CONJUGATE OF V2 WITH V1.
ZERO = 0.D0
SUMC = DCMPLX(ZERO,ZERO)
DO 10 J=1,N
10 SUMC = SUMC + DCONJG(V2(J))*V1(J)
SUM = DREAL(SUMC)
C
RETURN
C-----END OF COMPLEX INNER PRODUCT SUBROUTINE---------------------------
END
C
C-----LPERM-PERMUTES VECTORS--------------------------------------------
C
SUBROUTINE LPERM(W,U,IPERM)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION U(1),W(1)
INTEGER IPR(1),IPT(1)
C-----------------------------------------------------------------------
C SUBROUTINE HAS 2 BRANCHES: IPERM = 1, CALCULATES
C U = P*W
C LET J = IPR(K) THEN U(K) = W(J), K = 1,N.
C IPERM = 2, CALCULATES
C U = P'*W LET J = IPT(K) THEN U(K) = W(J), K=1,N.
C-----------------------------------------------------------------------
C
IF(IPERM.EQ.2) GO TO 30
C
C-----------------------------------------------------------------------
C IPERM = 1
DO 10 K = 1,N
J = IPR(K)
10 U(K) = W(J)
GO TO 60
C----------------------------------------------------------------------
C IPERM = 2
30 DO 40 K = 1,N
J = IPT(K)
40 U(K) = W(J)
C----------------------------------------------------------------------
60 CONTINUE
DO 50 K = 1,N
50 W(K) = U(K)
C
RETURN
C
C-----------------------------------------------------------------------
ENTRY LPERME(IPR,IPT,N)
C-----------------------------------------------------------------------
C
C-----END OF LPERM-------------------------------------------------
RETURN
END