C
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 FOR 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 T-EIGENVALUES IS SET
C EQUAL TO THE VALUE OF THAT T-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 T-EIGENVALUES ARE
C NEVER COMBINED 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 T-EIGENVALUES OF T(1,MEV)
C MP CONTAINS THE CORRESPONDING T-MULTIPLICITIES
C NDIS = NUMBER OF DISTINCT T-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 TO 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(BETA,TEIG,TKMAX,EPSM,RELTOL,SCALE3,SCALE4,
1 TMULT,NDIST,MEV,IPROJ)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION BETA(1),TEIG(1),SIGMA(4)
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(4)
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 SINGULAR VALUES OF A HAVE BEEN
C MISSED BY LANCZOS PROCEDURE. A SINGULAR VALUE WHOSE
C SINGULAR VECTOR(S) HAS A VERY SMALL PROJECTION ON THE
C STARTING VECTOR (< SINGLE PRECISION) CAN BE MISSED BECAUSE
C IT WILL THEN ALSO BE AN EIGENVALUE OF T(2,MEV) TO WITHIN
C THE SQUARE OF THIS ORIGINAL PROJECTION. HOWEVER,
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 SINGULAR VALUES, THEN ONE CAN EXPECT TO ALSO HAVE
C CONVERGENCE ON ANY SUCH 'HIDDEN' SINGULAR VALUES. (IF THERE
C ARE ANY). PROCEDURE CONSIDERS ONLY SPURIOUS T-EIGENVALUES AND
C ONLY THOSE SPURIOUS T-EIGENVALUES THAT ARE ISOLATED FROM GOOD
C T-EIGENVALUES. FOR EACH SUCH T-EIGENVALUE IT DOES 2 STURM
C SEQUENCES AND A FEW SCALAR MULTIPLICATIONS. UPON RETURN TO MAIN
C PROGRAM ERROR ESTIMATES WILL BE COMPUTED FOR ANY T-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 = -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 = -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(BETA,X1,TOLN,EPSM,MMAX,MK1,MK2,IC,IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION 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 GOOD T-EIGENVALUE 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 BETA ARRAY 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 SINGULAR VECTOR PROGRAM
C USES THESE VALUES TO DETERMINE A 1ST GUESS AT AN APPROPRIATE
C SIZE T-MATRIX FOR THE SINGULAR VALUE 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 GOOD T-EIGENVALUES BEING CONSIDERED WERE MULTIPLE
C AS SINGULAR VALUES 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 SINGULAR VECTOR PROGRAM TO COMPUTE ADDITIONAL
C SINGULAR VECTORS CORRESPONDING TO THESE MULTIPLE SINGULAR
C VALUES. THE MAIN PROGRAM LSVEC PROVIDED DOES NOT INCLUDE
C 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 - 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 - 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(BETA,V1,V2,X1,ERROR,ERRORV,EPS,G,MEV,IT,
1 IWRITE)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION 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 BETA CONTAINS THE 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 = ZERO
DO 10 I = 2,MEV
10 TSUM = TSUM + 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 = -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 = -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) = -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) = BETA(MEV)*V2(MEV-1)-X1*V2(MEV)
DO 130 J = 2,MM1
JM = MP1 - J
V1(JM) = BETA(JM)*V2(JM-1) + BETA(JM+1)*V2(JM+1
1) - X1*V2(JM)
130 CONTINUE
C
V1(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(' INVERSE ITERATION OUTPUT'/
1 2X,'TSIZE',13X,'T-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(BETA,EPSM,EVAL,EVALN,LB,UB,TTOL,M,NEVT)
C
C-----------------------------------------------------------------------
DOUBLE PRECISION 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(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.EPT) 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 EIGENVALUES OF T(1,M) ON (X1,XU)
C THERE IS EXACTLY ONE 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 - 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