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