C-----DOCUMENTATION FOR THE SINGLE-VECTOR------------------------------- C LANCZOS SINGULAR VALUE/VECTOR PROGRAMS C FOR REAL, RECTANGULAR MATRICES C C----------------------------------------------------------------------- C REFERENCE: C LANCZOS ALGORITHMS FOR LARGE SYMMETRIC EIGENVALUE COMPUTATIONS C VOL. 1 THEORY VOL. 2 PROGRAMS. JANUARY 1985. BIRKHAUSER, C BASEL. C C AUTHORS: JANE CULLUM AND RALPH A. WILLOUGHBY, IBM RESEARCH, C YORKTOWN HEIGHTS, NY 10598. PHONE: 914-945-2227 C----------------------------------------------------------------------- C C GIVEN A REAL RECTANGULAR MATRIX A OF ORDER M X N THE THREE C SETS OF FORTRAN FILES LABELLED LSVAL, LSSUB, AND LSMULT C CAN BE USED TO COMPUTE DISTINCT SINGULAR VALUES OF A IN C USER-SPECIFIED INTERVALS. C C CORRESPONDING SINGULAR VECTORS FOR SELECTED, COMPUTED C SINGULAR VALUES CAN BE COMPUTED USING THE SETS OF FILES C LABELLED LSVEC, LSSUB AND LSMULT. C C THESE PROGRAMS USE LANCZOS TRIDIAGONALIZATION WITHOUT C REORTHOGONALIZATION ON THE ASSOCIATED REAL SYMMETRIC MATRIX C C ---- ---- C | 0 A | C B = | | C | A-TRANSPOSE 0 | C ---- ---- C C OF ORDER M + N TO GENERATE REAL SYMMETRIC TRIDIAGONAL C MATRICES, T(1,MEV), OF ORDER MEV. SUBSETS OF THE EIGENVALUES OF C THESE T-MATRICES, LABELLED AS THE 'GOOD EIGENVALUES' OF T(1,MEV), C ARE APPROXIMATIONS TO THE DESIRED SINGULAR VALUES OF A. C CORRESPONDING RITZ VECTORS FOR B ARE APPROXIMATIONS TO C EIGENVECTORS OF B WHICH IN TURN CONTAIN APPROXIMATIONS TO C THE DESIRED LEFT AND RIGHT SINGULAR VECTORS OF A. THIS C PROCEDURE USES A SPECIAL STARTING VECTOR SUGGESTED BY GOLUB C AND KAHAN. THUS, THE STARTING LANCZOS VECTOR IS EITHER OF C THE FORM (V1,0) OR (0,V2) WHERE V1 IS MX1 AND V2 IS NX1 AND C ALL SUCCEEDING LANCZOS VECTORS GENERATED ALTERNATE BETWEEN C THESE 2 FORMS. THIS SPECIAL CHOICE OF STARTING VECTOR RESULTS C IN SIGNIFICANT GAINS IN STORAGE AND OPERATION COUNTS AND C ALSO IN CONVERGENCE RELATIVE TO A 'BRUTE FORCE' APPLICATION C OF THE REAL SYMMETRIC LANCZOS PROCEDURE DIRECTLY TO THE C MATRIX B ABOVE. FOR MORE DETAILS SEE REFERENCE 1 BELOW. C IN THE DISCUSSIONS T(1,MEV) DENOTES THE LANCZOS T-MATRIX C OF SIZE MEV. C C THE IDEAS USED IN THESE PROGRAMS ARE DISCUSSED IN THE FOLLOWING C REFERENCES. C C 1. JANE CULLUM, RALPH A. WILLOUGHBY AND MARK LAKE, A LANCZOS C ALGORITHM FOR COMPUTING SINGULAR VALUES AND VECTORS OF LARGE C MATRICES, SIAM J. SCIENTIFIC AND STATISTICAL COMPUTING, C VOL. 4, JUNE 1983, PP. 197-215. C C 2. JANE CULLUM AND RALPH A. WILLOUGHBY, LANCZOS ALGORITHMS C FOR LARGE SYMMETRIC MATRICES, PROGRESS IN C SCIENTIFIC COMPUTING, EDITORS, G. GOLUB, H.O. KREISS, C S. ARBARBANEL, AND R. GLOWINSKI, BIRKHAUSER BOSTON INC., C CAMBRIDGE, MASSACHUSETTS, 1984. C C 3. JANE CULLUM AND RALPH A. WILLOUGHBY, COMPUTING EIGENVECTORS C (AND EIGENVALUES) OF LARGE, SYMMETRIC MATRICES USING C LANCZOS TRIDIAGONALIZATION, LECTURE NOTES IN MATHEMATICS, C 773, NUMERICAL ANALYSIS PROCEEDINGS, DUNDEE 1979, EDITED BY C G. A. WATSON, SPRINGER-VERLAG, (1980), BERLIN, PP.46-63. C C 4. IBID, LANCZOS AND THE COMPUTATION IN SPECIFIED INTERVALS OF C THE SPECTRUM OF LARGE SPARSE, REAL SYMMETRIC MATRICES, SPARSE C MATRIX PROCEEDINGS 1978, ED. I.S. DUFF AND G. W. STEWART, C SIAM, PHILADELPHIA, PP.220-255, 1979. C C 5. IBID, COMPUTING EIGENVALUES OF VERY LARGE SYMMETRIC MATRICES- C AN IMPLEMENTATION OF A LANCZOS ALGORITHM WITHOUT C REORTHOGONALIZATION, J. COMPUT. PHYS. 44(1981), 329-358. C C C-----PORTABILITY------------------------------------------------------- C C C PROGRAMS WERE TESTED FOR PORTABILITY USING THE PFORT VERIFIER. C FOR DETAILS OF THE VERIFIER SEE FOR EXAMPLE, B. G. RYDER AND C A. D. HALL, "THE PFORT VERIFIER", COMPUTING SCIENCE TECHNICAL C REPORT 12, BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974, C (REVISED), JANUARY 1981. C C EXCEPT FOR THE FOLLOWING CONSTRUCTIONS WHICH CAN BE EASILY C MODIFIED BY THE USER TO MATCH THE PARTICULAR COMPUTER BEING C USED, THE PROGRAM STATEMENTS ARE PORTABLE. C C NONPORTABLE CONSTRUCTIONS. C C IN LSVAL AND IN LSVEC C 1. DATA/MACHEP STATEMENT C 2. ALL READ(5,*) STATEMENTS (FREE FORMAT) C 3. FORMAT(20A4) USED FOR THE EXPLANATORY HEADER ARRAY, EXPLAN C 4. FORMAT(4Z20) USED TO READ AND WRITE BETA FILES 1 AND 2. C IN LSMULT C 1. IN SVMAT, STRAN, AND USPEC THE ENTRY THAT PASSES THE C STORAGE LOCATIONS OF THE ARRAYS DEFINING THE C USER-SPECIFIED MATRIX. C 2. IN SAMPLE USPEC FOR 'DIAGONAL' MATRICES: THE FREE C FORMAT (8,*) AND THE FORMAT (20A4). C IN LSSUB C 1. ALL STATEMENTS ARE PORTABLE. C C C IN THE COMMENTS BELOW: C COMPLEX*16 = COMPLEX VARIABLE, 16 BYTES OF STORAGE C REAL*8 = REAL VARIABLE, 8 BYTES OF STORAGE C REAL*4 = REAL VARIABLE, 4 BYTES OF STORAGE C INTEGER*4 = INTEGER VARIABLE, 4 BYTES C C C-----A-MATRIX SPECIFICATION-------------------------------------------- C C C SUBROUTINE USPEC IS USED TO SPECIFY THE USER-SUPPLIED A-MATRIX. C SUBROUTINES SVMAT AND STRAN ARE, RESPECTIVELY, CORRESPONDING C MATRIX-VECTOR MULTIPLE SUBROUTINES FOR A AND FOR A-TRANSPOSE. C THESE SUBROUTINES SHOULD BE DESIGNED TO TAKE ADVANTAGE OF C ANY SPECIAL PROPERTIES OF THE USER-SUPPLIED MATRIX. THE C MATRIX-VECTOR MULTIPLIES REQUIRED BY THE LANCZOS PROCEDURES C MUST BE COMPUTED RAPIDLY AND ACCURATELY. C C SUBROUTINE USPEC HAS THE CALLING SEQUENCE C C CALL USPEC(M,N,MATNO) C C WHERE M IS THE NUMBER OF ROWS IN THE USER-SPECIFIED C A-MATRIX AND N IS THE NUMBER OF COLUMNS. MATNO IS A C <= 8 DIGIT INTEGER USED AS A MATRIX AND TEST IDENTIFICATION C NUMBER. THIS SUBROUTINE DEFINES (DIMENSIONS) THE ARRAYS C REQUIRED TO SPECIFY THE A-MATRIX. THIS SUBROUTINE ALSO C INITIALIZES THESE ARRAYS AND ANY OTHER PARAMETERS NEEDED TO C DEFINE THE MATRIX. THE STORAGE LOCATIONS OF THESE PARAMETERS C AND ARRAYS ARE THEN PASSED TO THE MATRIX-VECTOR MULTIPLY C SUBROUTINES SVMAT AND STRAN VIA ENTRIES. SAMPLE SUBROUTINES C ARE INCLUDED IN THE FORTRAN FILE LSMULT. C C IMPORTANT NOTE: C THE SAMPLE MATRIX-VECTOR MULTIPLY SUBROUTINES IN LSMULT C ASSUME THAT M >= N. THEY ALSO ASSUME THAT THE USER-SUPPLIED C INFORMATION ABOUT THE GIVEN MATRIX IS STORED ON FILE 8. C THE USER SHOULD SEE THE LSMULT PROGRAMS FOR MORE DETAILS. C C SUBROUTINE SVMAT HAS THE CALLING SEQUENCE C C CALL SVMAT(W,U,SUM) C C WHERE U AND W ARE REAL*8 VECTORS AND SUM IS A REAL*8 C SCALAR. SVMAT CALCULATES U = A*W - SUM*U FOR THE C USER-SPECIFIED A-MATRIX. SUBROUTINE STRAN HAS THE C CALLING SEQUENCE C C CALL STRAN(W,U,SUM) C C STRAN CALCULATES U = (A-TRANSPOSE)*W - SUM*U FOR THE C TRANSPOSE OF THE USER-SUPPLIED A-MATRIX. THE ARRAY AND PARAMETER C INFORMATION NEEDED TO PERFORM THE MATRIX-VECTOR MULTIPLIES C IS PASSED TO THE SVMAT AND THE STRAN SUBROUTINES FROM THE C USPEC SUBROUTINE VIA ENTRIES. ONE SET OF THE SAMPLE SVMAT C AND STRAN SUBROUTINES INCLUDED IN LSMULT COMPUTES C MATRIX-VECTOR MULTIPLIES FOR AN ARBITRARY SPARSE, C RECTANGULAR MATRIX STORED IN THE SPARSE FORMAT SPECIFIED C IN THE CORRESPONDING SAMPLE USPEC SUBROUTINE. THE LANCZS C SUBROUTINE CALLS SVMAT AND STRAN IN THE GENERATION OF C THE LANCZOS T-MATRICES FOR THE B MATRIX. C C THE DATA FOR THE A-MATRIX IS ASSUMED TO BE ON FILE 8 AND C IN THE FOLLOWING SPARSE FORMAT: C NZ = NUMBER OF NONZERO ELEMENTS OF A C ICOL(K), K = 1,N, NUMBER OF NONZEROS OF A IN COLUMN K. C IROW(K), K = 1,NZ, ROW INDEX OF A(K). C A(K), K=1,NZ CONTAINS THE ELEMENTS OF A BY COLUMN. C C C SVMATV AND STRAN ARE CALLED FROM THE SUBROUTINE LANCZS C WHICH GENERATES THE LANCZOS TRIDIAGONAL MATRICES, THE C BETA HISTORY. SIMILARLY, THSE SUBROUTINES ARE CALLED FROM C THE CORRESPONDING SINGULAR VECTOR PROGRAM, LSVEC. C SVMAT AND STRAN ARE DECLARED AS EXTERNAL VARIABLES. C EACH IS AN ARGUMENT FOR THE LANCZS SUBROUTINE. C C USPEC, SVMAT, AND STRAN SUBROUTINES SUITABLE FOR THE C USER-SPECIFIED MATRIX MUST BE SUPPLIED BY THE USER. C C THE MAIN PROGRAMS FOR THE SINGULAR VALUE AND SINGULAR VECTOR C CALCULATIONS ASSUME THAT INPUT FILE 5 CONTAINS THE ROW ORDER C M AND THE COLUMN ORDER N OF THE GIVEN A-MATRIX AND MATNO, C AN IDENTIFICATION NUMBER OF <= 8 DIGITS FOR THE GIVEN MATRIX. C C C-----MACHEP------------------------------------------------------------ C C C MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING THE RELATIVE C PRECISION OF THE FLOATING POINT ARITHMETIC USED. C MACHEP = 2.2 * 10**-16 FOR DOUBLE PRECISION ARITHMETIC ON C IBM 370-3081. C C THE USER WILL HAVE TO RESET THIS PARAMETER TO C THE CORRESPONDING VALUE FOR THE MACHINE BEING USED. NOTE THAT C IF A MACHINE WITH A MACHINE EPSILON THAT IS MUCH LARGER THAN THE C VALUE GIVEN HERE IS BEING USED, THEN THERE COULD BE C PROBLEMS WITH THE TOLERANCES. C C C-----SUBROUTINES AND FUNCTIONS USER MUST SUPPLY------------------------ C C C GENRAN, FINPRO, MASK, USPEC, SVMAT AND STRAN C C GENRAN = COMPUTES K PSEUDO-RANDOM NUMBERS AND STORES THEM IN C THE REAL ARRAY, G. THIS SUBROUTINE IS USED TO C GENERATE A STARTING VECTOR FOR THE LANCZOS PROCEDURE C IN THE SUBROUTINE LANCZS AND A STARTING RIGHT-HAND SIDE C FOR INVERSE ITERATION IN THE SUBROUTINE INVERR. C C TESTS REPORTED IN THE REFERENCES USED EITHER GGL1 OR C GGL2 FROM THE IBM LIBRARY SLMATH. C THE EXISTING CALLING SEQUENCE IS: C C CALL GENRAN(IIX,G,K). C C WHERE IIX =INTEGER SEED, G = REAL*4 ARRAY WHOSE C DIMENSION MUST BE >= K. K RANDOM NUMBERS ARE GENERATED C AND PLACED IN G. C C FINPRO = DOUBLE PRECISION FUNCTION WHICH COMPUTES THE INNER C PRODUCT OF 2 DOUBLE PRECISION VECTORS OF DIMENSION K. C TESTS REPORTED IN THE REFERENCES USED THE HARWELL C LIBRARY SUBROUTINE FM02AD. C EXISTING CALLING SEQUENCE IS C C CALL FINPRO(N,V,J,W,K). C C COMPUTES THE INNER PRODUCT OF DIMENSION N OF THE VECTORS C V AND W. SUCCESSIVE COMPONENTS OF V AND OF W ARE STORED C AT LOCATIONS THAT ARE ,RESPECTIVELY, J AND K UNITS APART. C C MASK = MASKS OVERFLOW AND UNDERFLOW. C USER MUST SUPPLY OR COMMENT OUT CALL. C C USPEC = DIMENSIONS AND INITIALIZES ARRAYS NEEDED TO SPECIFY C USER-SUPPLIED A-MATRIX. SEE A-MATRIX SPECIFICATION SECTION C C SVMAT = MATRIX-VECTOR MULTIPLY FOR USER-SUPPLIED A-MATRIX. C SEE A-MATRIX SPECIFICATION SECTION. C C STRAN = MATRIX-VECTOR MULTIPLY FOR TRANSPOSE OF USER-SUPPLIED C A-MATRIX. SEE A-MATRIX SPECIFICATION SECTION. C C C----------------------------------------------------------------------- C C COMMENTS FOR SINGULAR VALUE COMPUTATIONS C C----------------------------------------------------------------------- C C C-----PARAMETER CONTROLS FOR SINGULAR VALUE PROGRAMS-------------------- C C C PARAMETER CONTROLS ARE INTRODUCED TO ALLOW SEGMENTATION OF THE C SINGULAR VALUE COMPUTATIONS AND TO ALLOW VARIOUS COMBINATIONS C OF READ/WRITES. C C THE FLAG ISTART CONTROLS THE T-MATRIX (BETA HISTORY) C GENERATION. C C ISTART = (0,1) MEANS C C (0) THERE IS NO EXISTING BETA HISTORY AND ONE C MUST BE GENERATED. C C (1) THERE IS AN EXISTING BETA HISTORY AND IT IS C TO BE READ IN FROM FILE 2 AND EXTENDED IF NECESSARY. C C THE FLAG ISTOP CAN BE USED IN CONJUNCTION WITH THE FLAG ISTART TO C ALLOW SEGMENTATION OF THE SINGULAR VALUE COMPUTATIONS. C C ISTOP = (0,1) MEANS C C (0) PROGRAM COMPUTES ONLY THE REQUESTED BETAS, C STORES THEM AND THE LAST 2 LANCZOS VECTORS GENERATED C IN FILE 1 AND THEN TERMINATES. C C (1) PROGRAM COMPUTES REQUESTED BETAS AND THEN C USES THE BISEC SUBROUTINE TO CALCULATE EIGENVALUES C OF THE TRIDIAGONAL MATRICES GENERATED FOR THE ORDERS C SPECIFIED BY THE USER AND ON THE USER-SPECIFIED C INTERVALS. PROGRAM THEN USES THE SUBROUTINE INVERR C TO COMPUTE ERROR ESTIMATES FOR THE ISOLATED GOOD C T-EIGENVALUES WHICH ARE USED TO CHECK THE C CONVERGENCE OF THESE T-EIGENVALUES. C C CONTROL PARAMETERS FOR WRITES C C IHIS = (0,1) MEANS C C (0) IF ISTOP .GT. 0 THEN BETAS ARE NOT SAVED ON FILE 1. C C (1) PROGRAM WRITES BETAS AND LAST 2 LANCZOS C VECTORS TO FILE 1 SO THAT THE T-MATRIX GENERATION C MAY BE REUSED OR CONTINUED LATER IF NECESSARY. C TYPICALLY ONE WOULD ALWAYS DO THIS ON ANY RUN WHERE C A HISTORY FILE IS BEING GENERATED. HISTORY MUST BE C SAVED IN MACHINE FORMAT ((4Z20) FOR IBM/3081) SO C THAT NO ERRORS DUE TO FORMAT CONVERSIONS OCCUR. C C IDIST = (0,1) MEANS C C (0) DISTINCT EIGENVALUES OF T-MATRICES ARE NOT SAVED. C C (1) PROGRAM WRITES COMPUTED DISTINCT EIGENVALUES OF C T-MATRICES ALONG WITH THEIR T-MULTIPLICITIES C TO FILE 11. C C IWRITE = (0,1) MEANS C C (0) NO EXTENDED OUTPUT FROM SUBROUTINES BISEC AND INVERR C IS SENT TO FILE 6. C C (1) INDIVIDUAL COMPUTED T-EIGENVALUES AND CORRESPONDING C ERROR ESTIMATES FROM THE SUBROUTINES BISEC AND INVERR C ARE PRINTED OUT TO FILE 6 AS THEY ARE COMPUTED. C C THE PROGRAM ALWAYS MAKES A SEPARATE LIST OF THE COMPUTED GOOD C EIGENVALUES OF THE LANCZOS MATRICES T(1,MEV) CONSIDERED, C THESE ARE THE APPROXIMATIONS TO THE DESIRED SINGULAR VALUES, C ALONG WITH THEIR MINIMAL GAPS AS SINGULAR VALUES OF A AND C WRITES THEM TO FILE 3. CORRESPONDING ERROR ESTIMATES FOR ANY C ISOLATED COMPUTED GOOD T-EIGENVALUES (SINGULAR VALUES OF A) C ARE ALWAYS WRITTEN TO FILE 4. C C C-----INPUT/OUTPUT FILES FOR SINGULAR VALUE PROGRAMS-------------------- C C ANY INPUT DATA OTHER THAN THE BETA HISTORY SHOULD BE STORED C ON FILE 5. SEE SAMPLE INPUT/OUTPUT FROM TYPICAL RUN. C THE READ STATEMENTS IN THE GIVEN FORTRAN PROGRAM ASSUME THAT C THE DATA STORED ON FILE 5 IS IN FREE FORMAT. USER SHOULD NOTE C THAT 'FREE FORMAT' IS NOT CLASSIFIED AS PORTABLE BY PFORT SO THAT C THE USER MAY HAVE TO MODIFY THE READ STATEMENTS FROM FILE 5 TO C CONFORM TO WHAT IS PERMISSIBLE ON THE MACHINE BEING USED. C C FILE 6 WAS USED AS THE INTERACTIVE TERMINAL OUTPUT FILE. C THIS FILE PROVIDES A RUNNING ACCOUNT OF THE PROGRESS OF THE C COMPUTATIONS. THE AMOUNT OF INFORMATION PRINTED OUT IS C CONTROLLED BY THE PARAMETER IWRITE. C C DESCRIPTION OF OTHER I/O FILES C C FILE (K) CONTAINS: C C (1) OUTPUT FILE: C HISTORY FILE OF NEWLY-GENERATED T-MATRIX C (BETA VECTOR) AND LAST 2 LANCZOS VECTORS USED C IN THE T-MATRIX GENERATION. C IF IHIS = 0 AND ISTOP = 1, FILE 1 IS NOT WRITTEN. C C (2) INPUT FILE: C SAME AS FILE 1 EXCEPT THAT IT CONTAINS A C PREVIOUSLY-GENERATED T-MATRIX (IF ANY). IF ISTART = 1, C PROGRAM ASSUMES THAT THERE IS A HISTORY FILE OF C BETAS ON FILE 2. THESE BETAS AND THE LAST TWO LANCZOS C VECTORS USED IN THE T-MATRIX GENERATION ARE READ IN. C C (3) OUTPUT FILE: C COMPUTED GOOD EIGENVALUES OF THE T-MATRICES CONSIDERED. C ALSO CONTAINS T-MULTIPLICITIES OF THESE T-EIGENVALUES AS C EIGENVALUES OF THE T-MATRIX, AND THEIR GAPS AS C EIGENVALUES IN THE B MATRIX AND IN THE T-MATRIX. C NOTE THAT THESE GOOD T-EIGENVALUES ARE THE COMPUTED C SINGULAR VALUES OF THE A-MATRIX AND THAT THE GAPS C OF THESE EIGENVALUES AS EIGENVALUES OF THE B-MATRIX C ARE EQUAL TO THEIR GAPS AS SINGULAR VALUES OF A. FILE C 3 IS ALWAYS WRITTEN. C C (4) OUTPUT FILE: C ERROR ESTIMATES FOR THE ISOLATED COMPUTED SINGULAR C SINGULAR VALUES (ISOLATED GOOD EIGENVALUES OF T(1,MEV)) C THESE ARE OBTAINED USING THE SUBROUTINE INVERR. THESE C ESTIMATES USE THE LAST COMPONENTS OF THE ASSOCIATED C T-EIGENVECTORS WHICH ARE COMPUTED USING INVERSE C ITERATION. FILE 4 IS ALWAYS WRITTEN. C C C (8) INPUT FILE: C SAMPLE USPEC SUBROUTINE ASSUMES THAT THE ARRAYS C REQUIRED TO SPECIFY THE USER'S MATRIX ARE STORED ON C FILE 8. USERS MUST MAKE WHATEVER DEFINITIONS ARE C APPROPRIATE FOR THEIR MATRICES. C C (9) OUTPUT FILE: OPTIONAL C CAN BE USED TO STORE THE TRUE SINGULAR VALUES OF C A GIVEN TEST MATRIX, WHEN THE SINGULAR VALUE PROCEDURE C IS BEING EXERCISED ON A TEST MATRIX. C C (11) OUTPUT FILE: C COMPUTED DISTINCT EIGENVALUES OF T-MATRICES USED. C ALSO CONTAINS THEIR T-MULTIPLICITIES AND T-GAPS TO C NEAREST DISTINCT T-EIGENVALUES, AND THE T-MULTIPLICITY C PATTERN OF THE GOOD AND THE SPURIOUS T-EIGENVALUES. C FILE 11 IS WRITTEN ONLY IF IDIST = 1. C C C-----PARAMETERS SET BY THE SINGULAR VALUE PROGRAMS--------------------- C C C THESE PARAMETERS ARE SET INTERNALLY IN THE PROGRAM C C SCALEK K = 1,2,3,4 C C THE SCALING FACTORS SCALEK HAVE BEEN INTRODUCED IN AN C ATTEMPT TO MAKE THE TOLERANCES USED IN THE C T-MULTIPLICITY, SPURIOUS, ISOLATION AND PRTESTS ADJUST C TO THE SCALE OF THE GIVEN MATRIX. THESE FACTORS MUST C NOT BE MODIFIED. C C NOTE: THE USER SHOULD NOTE THAT IF THE MATRIX BEING C PROCESSED IS VERY STIFF, THAT IS THE RATIO OF THE LARGEST C SINGULAR VALUE TO THE SMALLEST SINGULAR VALUE IS VERY C LARGE, THEN THE TOLERANCES BEING USED IN BISEC, LUMP, ISOEV C AND PRTEST MAY NOT TREAT THE SMALLEST SINGULAR VALUES C VERY WELL. IN SOME SUCH CASES A USER-INTRODUCED REDUCTION C IN THE SIZE OF TKMAX AND THE SUBSEQUENT RECOMPUTATION OF C THE T-MATRIX EIGENVALUES CORRESPONDING TO THE SMALLEST C SINGULAR VALUES USING THIS TKMAX MAY RESULT IN IMPROVED C COMPUTATIONS AT THE LOW END. C C THE LUMP, ISOEV, AND PRTEST TOLERANCES THAT WERE USED C MOST IN THE TESTING OF THIS ALGORITHM WERE NOT C SCALE INVARIANT BUT SEEMED TO WORK WELL ON MATRICES THAT C HAD SINGULAR VALUES BOTH GREATER THAN AND LESS C THAN 1. THESE TOLERANCES ARE ALSO INCLUDED IN THESE THREE C SUBROUTINES BUT AS COMMENTED OUT STATEMENTS. THEY CAN BE C REVIVED BY COMMENTING OUT THE CORRESPONDING TOLERANCES C SPECIFIED IN THE STATEMENT ABOVE EACH OF THESE. 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 T-MULTIPLICITY TOLERANCE: MULTOL = (1000+MEV)*TTOL C LANCZOS CONVERGENCE TOLERANCE: CONTOL = BETA(MEV+1)*1.D-10 C C C BTOL = RELATIVE TOLERANCE USED TO ESTIMATE ANY LOSS OF LOCAL C ORTHOGONALITY OF THE LANCZOS VECTORS AFTER THE T-MATRIX C HAS BEEN GENERATED. THE LANCZOS PROCEDURE WORKS WELL C ONLY IF LOCAL ORTHOGONALITY BETWEEN SUCCESSIVE LANCZOS C VECTORS IS MAINTAINED. THE TNORM SUBROUTINE TESTS C WHETHER OR NOT C C MINIMUM |BETA(I)|/||A|| > BTOL. C I=2,KMAX C C IF THIS TEST IS VIOLATED BY SOME BETA AND A T-MATRIX THAT C WOULD INCLUDE SUCH A BETA IS REQUESTED, THEN THE LANCZOS C PROCEDURE WILL TERMINATE FOR THE USER TO DECIDE WHAT TO C DO. THE USER CAN OVER-RIDE THIS TEST BY SIMPLY DECREASING C THE SIZE OF BTOL, BUT THEN CONVERGENCE IS NOT AS CERTAIN. C THE PROGRAM SETS BTOL = 1.D-8 WHICH IS A VERY CONSERVATIVE C CHOICE. THE || A || IS ESTIMATED BY USING AN ESTIMATE C OF THE NORM OF THE T-MATRIX, T(1,KMAX). C C GAPTOL = RELATIVE TOLERANCE USED IN THE SUBROUTINE ISOEV C TO DETERMINE FOR WHICH OF THE GOOD T-EIGENVALUES, C THE COMPUTED SINGULAR VALUES, ERROR ESTIMATES SHOULD C BE COMPUTED. THE PROGRAM SETS GAPTOL = 1.D-8. C IF FOR A GIVEN 'GOOD' T-EIGENVALUE OF THE GIVEN C T-MATRIX THE COMPUTED GAP IN THE T-MATRIX IS TOO C SMALL AND IS DUE TO A 'SPURIOUS' EIGENVALUE OF C THE T-MATRIX, THEN THE 'GOOD' T-EIGENVALUE IS ASSUMED C TO HAVE CONVERGED AND AN ERROR ESTIMATE IS NOT C COMPUTED. C C C-----USER-SPECIFIED PARAMETERS FOR SINGULAR VALUE PROGRAMS------------- C C C RELTOL = RELATIVE TOLERANCE USED IN 'COMBINING' COMPUTED C EIGENVALUES OF T(1,MEV) PRIOR TO COMPUTING ERROR C ESTIMATES. C C THE LUMPING OF T-EIGENVALUES OCCURS IN SUBROUTINE LUMP. C LUMPING IS NECESSARY BECAUSE IT IS IMPOSSIBLE TO ACCURATELY C PREDICT THE ACCURACY OF THE BISEC SUBROUTINE. LUMP 'COMBINES' C T-EIGENVALUES THAT HAVE SLIPPED BY THE TOLERANCE THAT WAS USED C IN THE T-MULTIPLICITY TESTS. IN PARTICULAR IF FOR SOME J, C C |EVALUE(J)-EVALUE(J-1)| < DMAX1(RELTOL*|EVALUE(J)|,SCALE2*MULTOL) C C THEN THESE T-EIGENVALUES ARE 'COMBINED'. MULTOL IS THE TOLERANCE C THAT WAS USED IN THE T-MULTIPLICITY TEST IN BISEC. SEE THE HEADER C ON THE LUMP SUBROUTINE FOR MORE DETAILS. C C RELTOL IS SET TO 1.D-10. C C MXINIT = MAXIMUM NUMBER OF INVERSE ITERATIONS ALLOWED IN C SUBROUTINE INVERR FOR EACH ISOLATED GOOD T-EIGENVALUE C CONSIDERED. TYPICALLY ONLY ONE IS REQUIRED. C C SEEDS FOR RANDOM NUMBER GENERATORS = INTEGER*4 SCALARS. C C (1) SVSEED = SEED FOR STARTING VECTOR USED IN C T-MATRIX GENERATION IN LANCZS SUBROUTINE C C (2) RHSEED = SEED FOR RIGHT-HAND SIDE USED IN C INVERSE ITERATION COMPUTATIONS IN INVERR. C C BISEC DATA C C (1) NINT = NUMBER OF SUBINTERVALS ON WHICH SINGULAR VALUES C ARE TO BE COMPUTED. C C (2) LB(J) = (J = 1,NINT) = LEFT END POINTS OF THESE INTERVALS. C MUST BE PROVIDED IN INCREASING ORDER. THAT IS, C LB(J) < LB(J+1) FOR J = 1,NINT. C C (3) UB(J) = (J = 1,NINT) = RIGHT END POINTS OF THESE INTERVALS. C MUST BE PROVIDED IN INCREASING ORDER. THAT IS, C UB(J) < UB(J+1) FOR J = 1,NINT. C C (4) MXSTUR = MAXIMUM NUMBER OF STURM ITERATIONS ALLOWED FOR C ENTIRE SET OF SINGULAR VALUE CALCULATIONS OVER C ALL SPECIFIED SIZE T-MATRICES. PROGRAM WILL C TERMINATE IF THIS LIMIT IS EXCEEDED. C C T-MATRICES C C SIZES OF T-MATRICES C C (1) KMAX= MAXIMUM ORDER FOR T-MATRIX THAT USER IS WILLING C TO CONSIDER. C C (2) NMEVS = MAXIMUM NUMBER OF T-MATRICES THAT WILL BE C CONSIDERED. C C (3) NMEV(J) (J=1,NMEVS) = SIZES OF T-MATRIX TO BE C CONSIDERED SEQUENTIALLY. C C T-MATRIX-GENERATION C C IPAR = (1,2) MEANS C C (1) STARTING VECTOR IS OF FORM (0,V2) WHERE V2 IS C NX1. USE WHEN M > N . C C (2) STARTING VECTOR IF OF FORM (V1,0) WHERE V1 IS C MX1. USE WHEN M < N . C C USER SHOULD NOTE THAT THIS PROGRAM FIRST COMPUTES A T-MATRIX C OF ORDER KMAX AND THEN CYCLES THROUGH THE T-MATRICES SPECIFIED C A PRIORI BY THE USER, USING THE SUBROUTINE BISEC TO COMPUTE C EIGENVALUES OF THE T-MATRICES ON THE INTERVALS SPECIFIED BY C THE USER. SUBSETS OF THESE T-EIGENVALUES ARE THEN SELECTED C AS APPROXIMATIONS TO THE DESIRED SINGULAR VALUES. C C IDEALLY, ONE WOULD COMPUTE THE SINGULAR VALUE APPROXIMATIONS C AT A REASONABLE SIZE T-MATRIX, LOOK AT THE ACCURACY OF THE C COMPUTED RESULTS AND USE THAT TO DETERMINE AN APPROPRIATE C INCREMENT FOR THE SIZE OF THE T-MATRIX BASED UPON WHAT C HAS ALREADY CONVERGED AND UPON THE SIZES OF THE ERROR ESTIMATES C ON THOSE SINGULAR VALUES THAT ARE DESIRED BUT THAT HAVE NOT C YET CONVERGED. HOWEVER, IN THE INTERESTS OF GENERALITY AND C SIMPLICITY WE CHOSE NOT TO DO THAT HERE. C C .