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