C-----LSVAL (SINGULAR VALUES OF REAL, RECTANGULAR MATRICES------------- C C CONTAINS MAIN PROGRAM FOR COMPUTING DISTINCT SINGULAR VALUES OF C A REAL M X N MATRIX USING LANCZOS TRIDIAGONALIZATION WITHOUT C REORTHOGONALIZATION AND WITH SPECIAL STARTING VECTORS. C C FOR A GIVEN REAL MATRIX A OF ORDER M X N THE LANCZOS RECURSION C IS APPLIED TO THE ASSOCIATED REAL SYMMETRIC MATRIX B OF ORDER C MN = M + N C C ---- ---- C | 0 A | C B = | | C | A-TRANSPOSE 0 | C ---- ---- C C USING SPECIAL STARTING VECTORS. PLEASE NOTE: ONLY EVEN ORDER C LANCZOS TRIDIAGONAL MATRICES AND ONLY NONNEGATIVE SUBINTERVALS C ARE PERMISSIBLE. 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 BETA FILES. C C----------------------------------------------------------------------- DOUBLE PRECISION BETA(5001),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 SVMAT,STRAN C----------------------------------------------------------------------- DATA MACHEP/Z3410000000000000/ EPSM = 2.0D0*MACHEP C----------------------------------------------------------------------- C C ARRAYS MUST BE DIMENSIONED AS FOLLOWS: C 1. BETA: >= (KMAX+1) WHERE KMAX IS READ IN AND IS C THE SIZE OF THE LARGEST T-MATRIX THAT CAN BE CONSIDERED. C 2. V1: >= MAX(M,KMAX+1) C 3. V2: >= MAX(N,KMAX) C 4. VS: >= KMAX C 5. G: >= MAX(2*KMAX,M,N) 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(BETA(J), J = 1,MEV) C BISEC CONVERGENCE TOLERANCE: BISTOL = DSQRT(1000+MEV)*TTOL C BISEC 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, RECTANGULAR 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. 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 HEADERS FOR RUN READ(5,20) EXPLAN WRITE(6,20) EXPLAN READ(5,20) EXPLAN WRITE(6,20) EXPLAN 20 FORMAT(20A4) C C READ THE ROW ORDER M OF THE MATRIX AND THE COLUMN ORDER N. C READ THE MAXIMUM ORDER OF THE T-MATRICES ALLOWED (KMAX), C THE NUMBER OF T-MATRICES ALLOWED (NMEVS), AND A C MATRIX IDENTIFICATION NUMBER (MATNO). READ(5,20) EXPLAN READ(5,*) M,N,KMAX,NMEVS,MATNO NM = M + N 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 BETA FILE IS NOT C AVAILABLE. ISTART = 1 MEANS BETA FILE IS AVAILABLE ON C FILE 2. C ISTOP = (0,1): ISTOP = 0 MEANS PROCEDURE GENERATES BETA C FILE AND THEN TERMINATES. ISTOP = 1 MEANS PROCEDURE GENERATES C BETAS IF NEEDED AND THEN COMPUTES SINGULAR VALUES AND C ERROR ESTIMATES AND THEN TERMINATES. READ(5,20) EXPLAN READ(5,*) ISTART,ISTOP C C IHIS = (0,1): IHIS = 0 MEANS BETA FILE IS NOT WRITTEN C TO FILE 1. IHIS = 1 MEANS 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. SPECIFY THE PARITY (IPAR) OF THE C LANCZOS STARTING VECTOR. IF M > N, THEN IPAR = 1, C IF M < N, THEN IPAR = 2. READ(5,20) EXPLAN READ(5,*) IHIS,IDIST,IWRITE,IPAR IF(M.GT.N) IPAR = 1 IF(M.LT.N) IPAR = 2 IPARO = IPAR C C READ IN THE RELATIVE TOLERANCE (RELTOL) FOR USE IN THE C SPURIOUS, T-MULTIPLICITY, AND PRTEST TESTS. READ(5,20) EXPLAN READ(5,*) RELTOL C C READ IN THE SIZES OF THE T-MATRICES TO BE CONSIDERED. C NOTE THAT ONLY EVEN ORDER T-SIZES ARE PERMISSIBLE. READ(5,20) EXPLAN READ(5,*) (NMEV(J), J=1,NMEVS) C C CHECK TO SEE THAT ALL T-SIZES PROVIDED ARE EVEN ORDERED. C TERMINATE IF THAT IS NOT THE CASE. DO 30 I = 1,NMEVS NMEV2 = NMEV(I)/2 IF(2*NMEV2.NE.NMEV(I)) GO TO 670 30 CONTINUE 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 SUBROUTINES SVMAT AND STRAN. C CALL USPEC(M,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,40) MATNO,M,N,KMAX 40 FORMAT(/3X,'MATRIX ID',5X,'M',5X,'N',4X,'MAX ORDER OF T'/ 1 I12,2I6,I18/) C WRITE(6,50) ISTART,ISTOP 50 FORMAT(/2X,'ISTART',3X,'ISTOP'/2I8/) C WRITE(6,60) IHIS,IDIST,IWRITE,IPAR 60 FORMAT(/4X,'IHIS',3X,'IDIST',2X,'IWRITE',4X,'IPAR'/4I8/) C WRITE(6,70) SVSEED,RHSEED 70 FORMAT(/' SEEDS FOR RANDOM NUMBER GENERATOR'// 1 4X,'LANCZS SEED',4X,'INVERR SEED'/2I15/) C WRITE(6,80) (NMEV(J), J=1,NMEVS) 80 FORMAT(/' SIZES OF T-MATRICES TO BE CONSIDERED'/(6I12)) C WRITE(6,90) RELTOL,GAPTOL,BTOL 90 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,100) (J,LB(J),UB(J), J=1,NINT) 100 FORMAT(/' BISEC WILL BE USED ON THE FOLLOWING INTERVALS'/ 1 (I6,2E20.6)/) C IF (ISTART.EQ.0.AND.IPAR.EQ.1) WRITE(6,110) IF (ISTART.EQ.0.AND.IPAR.EQ.2) WRITE(6,120) 110 FORMAT(/' STARTING VECTOR IS OF FORM (0,V2)'/) 120 FORMAT(/' STARTING VECTOR IS OF FORM (V1,0)'/) C IF (ISTART.EQ.0) GO TO 170 C C READ IN BETA HISTORY FROM FILE 2 C READ(2,130)MOLD,MO,NO,IPARO,IPAR,SVSOLD,MATOLD 130 FORMAT(3I6,2I3,I12,I8) C IF (KMAX.LT.MOLD) KMAX = MOLD KMAX1 = KMAX + 1 C C CHECK THAT M, N, MATRIX ID MATNO, AND RANDOM SEED SVSEED C AGREE WITH THOSE IN THE HISTORY FILE. IF NOT PROCEDURE STOPS. C ITEMP = (MO-M)**2+(NO-N)**2+(MATNO-MATOLD)**2+(SVSEED-SVSOLD)**2 C IF (ITEMP.EQ.0) GO TO 150 C WRITE(6,140) 140 FORMAT(' PROGRAM TERMINATES'/ ' READ FROM FILE 2 CORRESPONDS TO 1 DIFFERENT MATRIX THAN MATRIX SPECIFIED'/) GO TO 690 C 150 CONTINUE MOLD1 = MOLD+1 C READ(2,160)(BETA(J), J=1,MOLD1) 160 FORMAT(4Z20) C IF (KMAX.EQ.MOLD) GO TO 190 C READ(2,160)(V1(J), J=1,M) READ(2,160)(V2(J), J=1,N) C 170 CONTINUE IIX = SVSEED C C----------------------------------------------------------------------- C CALL LANCZS(SVMAT,STRAN,BETA,V1,V2,G,KMAX,MOLD1,M,N,IPAR,IIX) C C----------------------------------------------------------------------- C KMAX1 = KMAX + 1 C IF (IHIS.EQ.0.AND.ISTOP.GT.0) GO TO 190 C WRITE(1,180) KMAX,M,N,IPARO,IPAR,SVSEED,MATNO 180 FORMAT(3I6,2I3,I12,I8,' = KMAX,M,N,IPARO,IPAR,SVSEED,MATNO') C WRITE(1,160)(BETA(I), I=1,KMAX1) C WRITE(1,160)(V1(I), I=1,M) WRITE(1,160)(V2(I), I=1,N) C IF (ISTOP.EQ.0) GO TO 570 C 190 CONTINUE BKMIN = BTOL WRITE(6,200) 200 FORMAT(/' T-MATRICES (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(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 T-EIGENVALUES THAT HAD 'TOO SMALL' C A PROJECTION ON THE STARTING VECTOR. C CALL TNORM(BETA,BKMIN,TKMAX,KMAX,IB) C C----------------------------------------------------------------------- C TTOL = EPSM*TKMAX C C LOOP ON THE SIZE OF THE T-MATRIX C 210 CONTINUE MMB = MMB + 1 C NOTE THAT ONLY EVEN ORDER T-SIZES ARE PERMISSIBLE. MEV = NMEV(MMB) C IS MEV TOO LARGE ? IF(MEV.LE.KMAX) GO TO 230 WRITE(6,220) MMB, MEV, KMAX 220 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 570 C 230 MP1 = MEV + 1 BETAM = BETA(MP1) C IF (IB.GE.0) GO TO 240 C T0 = BTOL C C----------------------------------------------------------------------- C CALL TNORM(BETA,T0,T1,MEV,IBMEV) C C----------------------------------------------------------------------- C TEMP = T0/TKMAX IBMEV = IABS(IBMEV) IF (TEMP.GE.BTOL) GO TO 240 IBMEV = -IBMEV GO TO 630 C 240 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 = T-MULTIPLICITIES OF THE T-EIGENVALUES STORED 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 T-MULTIPLE AND IS THEREFORE NOT ONLY GOOD BUT C ALSO A CONVERGED GOOD T-EIGENVALUE. C C CALL BISEC(BETA,V1,V2,VS,LB,UB,EPSM,TTOL,MP,NINT, 1 MEV,NDIS,IC,IWRITE) C C----------------------------------------------------------------------- C IF (NDIS.EQ.0) GO TO 650 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 T-EIGENVALUES. 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 T-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 260 C WRITE(6,250) NDIS, MEV, LOOP 250 FORMAT(/I6,' DISTINCT T-EIGENVALUES WERE COMPUTED IN BISEC AT MEV 1=',I6/ 2X,' LUMP SUBROUTINE REDUCES NUMBER OF DISTINCT T-EIGENVALU 1ES TO',I6) C 260 CONTINUE NDIS = LOOP BETA(MP1) = BETAM C C----------------------------------------------------------------------- C THE SUBROUTINE ISOEV LABELS THOSE SIMPLE T-EIGENVALUES OF T(1,MEV) C WITH VERY SMALL GAPS BETWEEN NEIGHBORING T-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 T-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 T-EIGENVALUE. C 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,270)NG,NISO,NDIS 270 FORMAT(/I6,' SINGULAR VALUES HAVE BEEN COMPUTED'/ 1 I6,' OF THESE ARE ISOLATED'/ 2 I6,' = NUMBER OF DISTINCT T-EIGENVALUES COMPUTED'/) C C DO WE WRITE DISTINCT T-EIGENVALUES TO FILE 11? IF (IDIST.EQ.0) GO TO 310 C WRITE(11,280) NDIS,NISO,MEV,M,N,SVSEED,MATNO 280 FORMAT(5I5,I12,I8,' = NDIS,NISO,MEV,M,N,SVSEED,MATNO'/) C WRITE(11,290) (MP(I),VS(I),G(I), I=1,NDIS) 290 FORMAT(2(I3,E25.16,E12.3)) C WRITE(11,300) NDIS, (MP(I), I=1,NDIS) 300 FORMAT(/I6,' = NDIS, T-MULTIPLICITIES (0 MEANS SPURIOUS)'/(20I4)) C 310 CONTINUE C IF (NISO.NE.0) GO TO 340 C WRITE(4,320) MEV 320 FORMAT(/' AT MEV = ',I6,' THERE ARE NO ISOLATED T-EIGENVALUES'/ 1' SO NO ERROR ESTIMATES WERE COMPUTED/') C WRITE(6,330) 330 FORMAT(/' ALL COMPUTED SINGULAR VALUES ARE T-MULTIPLE'/ 1 ' THEREFORE ALL COMPUTED SINGULAR VALUES ARE ASSUMED TO HAVE CONV 1ERGED'/) C ICONV = 1 GO TO 380 C 340 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 T-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 IT = MXINIT CALL INVERR(BETA,V1,V2,VS,EPSM,G,MP,MEV,MMB,NDIS,NISO,NM, 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. C IF THIS TEST IS SATISFIED, THEN CONVERGENCE FLAG, ICONV IS SET C TO 1. TYPICALLY ERROR ESTIMATES ARE VERY CONSERVATIVE. C WRITE(6,350) CONTOL 350 FORMAT(/' CONVERGENCE IS TESTED USING THE CONVERGENCE TOLERANCE', 1E13.4/) C II = MEV +1 IF = MEV+NISO DO 360 I = II,IF IF (ABS(G(I)).GT.CONTOL) GO TO 380 360 CONTINUE ICONV = 1 MMB = NMEVS C WRITE(6,370) CONTOL 370 FORMAT(' ALL COMPUTED ERROR ESTIMATES WERE LESS THAN',E15.4/ 1 ' THEREFORE PROCEDURE TERMINATES'/) C 380 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 SINGULAR VECTOR ON THE STARTING C VECTOR WAS TOO SMALL. NUMERICAL TESTS INDICATE THAT C SUCH SINGULAR VALUES ARE RARE. THEREFORE, IF MANY OF C THESE HIDDEN SINGULAR VALUES APPEAR ON SOME RUN, THE USER C CAN BE CERTAIN THAT SOMETHING IS FOULED UP. C IF (ICONV.EQ.0) GO TO 510 C C----------------------------------------------------------------------- C CALL PRTEST (BETA,VS,TKMAX,EPSM,RELTOL,SCALE3,SCALE4, 1 MP,NDIS,MEV,IPROJ) C C----------------------------------------------------------------------- C IF(IPROJ.EQ.0) GO TO 500 C IF(IDIST.EQ.1) WRITE(11,390) IPROJ 390 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 420 J = 1,NDIS IF(MP(J).NE.ITEN) GO TO 420 T0 = VS(J) C C----------------------------------------------------------------------- C IT = MXINIT CALL INVERM(BETA,V1,V2,T0,TEMP,T1,EPSM,G,MEV,IT,IWRITE) C C----------------------------------------------------------------------- C IF(TEMP.LE.1.D-10) GO TO 410 C ERROR ESTIMATE WAS NOT SMALL REJECT RELABELLING OF THIS C T-EIGENVALUE. IF(IDIST.EQ.1) WRITE(11,400) J,T0,TEMP 400 FORMAT(/' LAST COMPONENT FOR',I6,'TH T-EIGENVALUE',E20.12/' IS TOO 1 LARGE = ',E15.6,' SO DO NOT ACCEPT PRTEST RELABELLING'/) MP(J) = 0 IPROJ = IPROJ - 1 GO TO 420 C RELABELLING ACCEPTED 410 NISOM = NISOM + 1 G(NISOM) = BETAM*TEMP 420 CONTINUE IWRITE = IWRITO C IF(IPROJ.EQ.0) GO TO 460 WRITE(6,430) IPROJ 430 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,440) IPROJ 440 FORMAT(/I6,' T-EIGENVALUES WERE RELABELLED AS GOOD'/ 1' BELOW IS CORRECTED T-MULTIPLICITY PATTERN'/) C WRITE(6,450) NDIS, (MP(I), I=1,NDIS) IF(IDIST.EQ.1) WRITE(11,450) NDIS, (MP(I), I=1,NDIS) 450 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. 460 NDIS1 = NDIS - 1 G(NDIS) = VS(NDIS1)-VS(NDIS) G(1) = VS(2)-VS(1) C DO 470 J = 2,NDIS1 T0 = VS(J)-VS(J-1) T1 = VS(J+1)-VS(J) G(J) = T1 IF (T0.LT.T1) G(J) = -T0 470 CONTINUE IF(IPROJ.EQ.0) GO TO 500 C WRITE TO FILE 4 ERROR ESTIMATES FOR THOSE T-EIGENVALUES RELABELLED NGOOD = 0 DO 480 J = 1,NDIS IF(MP(J).EQ.0) GO TO 480 NGOOD = NGOOD + 1 IF(MP(J).NE.ITEN) GO TO 480 T0 = VS(J) NISO = NISO + 1 NISOM = MEV + NISO WRITE(4,490) NGOOD,T0,G(NISOM),G(J) 480 CONTINUE 490 FORMAT(I10,E25.16,2E14.3) C 500 CONTINUE C C WRITE THE COMPUTED SINGULAR VALUES TO FILE 3. FIRST TRANSFER THEM C TO V2 AND THEIR T-MULTIPLICITIES TO THE CORRESPONDING POSITIONS C IN MP AND COMPUTE THE B-MINGAPS, THE MINIMAL GAPS BETWEEN THE C SINGULAR VALUES CONSIDERED AS EIGENVALUES OF THE B-MATRIX. C 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 510 CONTINUE C NG = 0 DO 520 I = 1,NDIS IF (MP(I).EQ.0) GO TO 520 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 520 CONTINUE C WRITE(6,530)MEV 530 FORMAT(//' SINGULAR VALUE CALCULATION AT MEV = ',I6,' IS COMPLE 1TE'//) C C NG = NUMBER OF COMPUTED DISTINCT GOOD T-EIGENVALUES. NEXT C GENERATE GAPS BETWEEN GOOD T-EIGENVALUES (BMINGAPS) AND PUT THEM C IN G. G(J) < 0 MEANS THE BMINGAP 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 540 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 540 CONTINUE C C WRITE GOOD T-EIGENVALUES (COMPUTED SINGULAR VALUES) OUT TO FILE 3. C WRITE(3,550)NG,NDIS,MEV,M,N,SVSEED,MATNO,IPARO,MULTOL,IB,BTOL 550 FORMAT(5I6,I12,I8,I2,'=NG,ND,MEV,M,N,SEED,MN,IPARO'/ 1 E20.12,I6,E13.4,' = MUTOL,INDEX MINIMAL BETA,BTOL'/ 1' SV NO',2X,'T-MULT',10X,'SINGULAR VALUE',7X,'BMINGAP',7X,'TMINGAP 1') C WRITE(3,560)(I,MP(I),V2(I),G(I),V1(I), I=1,NG) 560 FORMAT(I6,I8,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 210 C C END OF LOOP ON DIFFERENT SIZE T-MATRICES ALLOWED. C 570 CONTINUE C IF(ISTOP.EQ.0) WRITE(6,580) 580 FORMAT(/' T-MATRICES (BETA) ARE NOW AVAILABLE, TERMINATE'/) IF (IHIS.EQ.1.AND.KMAX.NE.MOLD) WRITE(1,590) 590 FORMAT(/' ABOVE ARE THE FOLLOWING VECTORS '/ 2 ' BETA(I), I = 1,KMAX+1'/ 3 ' FINAL TWO LANCZOS VECTORS OF ORDERS M,N FOR I = KMAX,KMAX+1'/ 4 ' ALL VECTORS IN THIS FILE HAVE FORMAT 4Z20'/ 5 ' ----- END OF FILE 1 NEW BETA HISTORY---------------'///) C IF (ISTOP.EQ.0) GO TO 690 C WRITE(3,600) 600 FORMAT(/' ABOVE ARE COMPUTED SINGULAR VALUES'/ 1 ' NG = NUMBER OF SINGULAR VALUES COMPUTED'/ 2 ' NDIS = NUMBER OF COMPUTED DISTINCT EIGENVALUES OF T(1,MEV)'/ 3 ' M = ROW ORDER OF A N = COLUMN ORDER, MATNO = MATRIX IDENT'/ 4 ' MULTOL = T-MULTIPLICITY TOLERANCE FOR T-EIGENVALUES IN BISEC'/ 4 ' T-MULT IS THE T-MULTIPLICITY OF SINGULAR VALUE'/ 5 ' T-MULT = -1 MEANS SPURIOUS T-EIGENVALUE TOO CLOSE'/ 6 ' DO NOT COMPUTE ERROR ESTIMATES FOR SUCH T-EIGENVALUES'/ 7 ' BMINGAP = MINIMAL GAP BETWEEN THE COMPUTED SINGULAR VALUES'/ 8 ' BMINGAP .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 SINGULAR VALUES------------------------'//) C IF (IDIST.EQ.1) WRITE(11,610) 610 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 COMPUTED SINGULAR VALUE AS HAVING A VERY CLOSE SPURIOUS 6 '/' T-EIGENVALUE SO THAT NO ERROR ESTIMATE WILL BE COMPUTED'/ 7 ' FOR THAT SINGULAR VALUE 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 COMPUTED SINGULAR VALUES. '/ 3 ' NISO = NUMBER OF ISOLATED (IN T-MATRIX) SINGULAR VALUES. '/ 4 ' NISO ALSO IS THE COUNT OF +1 ENTRIES IN T-MULTIPLICITY PATTERN. 5 '/' ----- END OF FILE 11 DISTINCT T-EIGENVALUES--------------'//) C IF(NISO.NE.0) WRITE(4,620) 620 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 = SV*U AND SV = 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 SPURIOUS T-EIGENVALUE.' 6/' ------ END OF FILE 4 ERRINV -------------------------------'//) GO TO 690 C 630 CONTINUE C IBB = IABS(IBMEV) IF (IBMEV.LT.0) WRITE(6,640) MEV,IBB,BETA(IBB) 640 FORMAT(/' PROGRAM TERMINATES BECAUSE MEV REQUESTED = ',I6,' IS .GT 1',I6/' AT WHICH AN ABNORMALLY SMALL BETA = ' , E13.4,' OCCURRED'/) GO TO 690 C 650 IF (NDIS.EQ.0.AND.ISTOP.GT.0) WRITE(6,660) 660 FORMAT(/' INTERVALS SPECIFIED FOR BISECT DID NOT CONTAIN ANY T-EIG 1ENVALUES'/' PROGRAM TERMINATES') GO TO 690 C 670 WRITE(6,680) I, NMEV(I) 680 FORMAT(//I6,'TH T-SIZE REQUESTED ',I6,' IS ODD'/ 1' BUT ONLY EVEN T-SIZES ARE PERMISSIBLE. PROGRAM TERMINATES FOR U 1SER TO FIX'//) GO TO 690 C 690 CONTINUE C STOP C-----END OF MAIN PROGRAM FOR LANCZOS SINGULAR VALUE COMPUTATIONS------- END .