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= 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