C C-----DOCUMENTATION FOR SINGLE-VECTOR----------------------------------- C LANCZOS EIGENVALUE/EIGENVECTOR PROGRAMS FOR C (1) REAL SYMMETRIC MATRICES C (2) HERMITIAN MATRICES C (3) FACTORED INVERSES OF REAL SYMMETRIC MATRICES C (4) REAL SYMMETRIC, GENERALIZED PROBLEMS WHERE ONE OF THE C MATRICES IS POSITIVE DEFINITE AND ITS CHOLESKY FACTORS ARE C AVAILABLE 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 REAL SYMMETRIC MATRICES: C C GIVEN A REAL SYMMETRIC MATRIX A OF ORDER N THE THREE SETS OF C FORTRAN FILES LABELLED LEVAL, LESUB, AND LEMULT CAN BE USED TO C COMPUTE DISTINCT EIGENVALUES OF THE USER-SPECIFIED MATRIX C IN USER-SPECIFIED INTERVALS. C C CORRESPONDING EIGENVECTORS FOR SELECTED, COMPUTED EIGENVALUES CAN C BE COMPUTED USING THE SETS OF FILES LABELLED LEVEC, LESUB, AND C LEMULT. C C C HERMITIAN MATRICES: C C GIVEN A HERMITIAN MATRIX A OF ORDER N THE THREE SETS OF C FORTRAN FILES LABELLED HLEVAL, LESUB, AND HLEMULT CAN BE USED C TO COMPUTE DISTINCT EIGENVALUES IN USER-SPECIFIED INTERVALS. C C CORRESPONDING EIGENVECTORS FOR SELECTED, COMPUTED EIGENVALUES C CAN BE COMPUTED USING THE SETS OF PROGRAMS LABELLED HLEVEC, C LESUB, AND HLEMULT. C C C FACTORED INVERSES OF REAL SYMMETRIC MATRICES: C C GIVEN A REAL SYMMETRIC MATRIX A, THE LANCZOS RECURSION IS C APPLIED TO THE INVERSE OF A, USING A FACTORIZATION C OF A. THE SETS OF FILES LIVAL, LESUB, AND LIMULT C CAN BE USED TO COMPUTE THE DISTINCT EIGENVALUES OF THE C INVERSE OF THE A-MATRIX AND OF A IN USER-SPECIFIED C INTERVALS. THE PROGRAMS ACTUALLY ALLOW ONE TO WORK WITH C ANY MATRIX B = PCP' WHERE C = S0*A + SHIFT*I, WHERE C S0 AND SHIFT ARE SCALARS CHOSEN BY THE USER AND P IS A C PERMUTATION MATRIX CHOSEN SUCH THAT THE FACTORIZATION C OF THE B-MATRIX RETAINS SPARSITY. IN THE C SAMPLE LIMULT SUBROUTINES PROVIDED, S0 AND SHIFT MUST BE C CHOSEN SO THAT THE RESULTING B-MATRIX IS POSITIVE DEFINITE, C AND THE CHOLESKY FACTORS ARE USED TO SOLVE B*U = V. C HOWEVER, THE USER CAN EASILY REPLACE THE SAMPLE USPEC AND C BSOLV SUBROUTINES PROVIDED BY SUBROUTINES THAT ALLOW THE C GENERAL FACTORIZATION L*D*(L-TRANSPOSE). THESE LANCZOS C PROGRAMS APPLY THE LANCZOS RECURSION TO B-INVERSE, USING C THE FACTORIZATION PROVIDED. OPTIONAL PREPROCESSING PROGRAMS C PERMUT, LORDER, LFACT, AND LTEST ARE PROVIDED FOR SET-UP PURPOSES. C PERMUT USES THE SPARSPAK PACKAGE OF A. GEORGE, J. LIU AND C E. NG TO OBTAIN A REORDERING OF THE GIVEN MATRIX THAT C PRESERVES SPARSENESS ON SUBEQUENT FACTORIZATION. LORDER C CAN BE USED TO REORDER A GIVEN MATRIX, USING A GIVEN C PERMUTATION. LFACT CAN BE USED TO COMPUTE THE CHOLESKY C FACTORS OF A GIVEN POSITIVE DEFINITE B-MATRIX. LTEST CAN C BE USED TO ESTIMATE THE NUMERICAL CONDITION OF THE C B-MATRIX. C C CORRESPONDING EIGENVECTORS FOR SELECTED, COMPUTED C EIGENVALUES CAN BE COMPUTED USING THE SETS OF FILES C LABELLED LIVEC, LESUB, AND LIMULT. C C GENERALIZED REAL SYMMETRIC PROBLEMS: C C GIVEN 2 REAL SYMMETRIC MATRICES A AND B WHERE IN ADDITION B IS C POSITIVE DEFINITE AND ITS CHOLESKY FACTORS ARE AVAILABLE, C THE SETS OF FILES LGVAL, LGMULT, AND LESUB CAN BE USED C TO COMPUTE THE DISTINCT EIGENVALUES OF THE GENERALIZED C PROBLEM A*X = EVAL*B*X. C C CORRESPONDING EIGENVECTORS CAN BE COMPUTED USING THE PROGRAMS C LGVEC, LGMULT, AND LESUB. NOTE THAT THE PREPROCESSING PROGRAMS C AVAILABLE FOR USE IN CASE (3) (PERMUT, LORDER, LFACT, AND LTEST) C CAN ALSO BE USED IN THIS CASE TO OBTAIN A SUITABLE PERMUTATION, C AND A FACTORIZATION OF THE RESULTING B-MATRIX. THE A-MATRIX C CAN THEN BE PERMUTED USING LORDER. C C C THESE PROGRAMS ALL USE LANCZOS TRIDIAGONALIZATION WITHOUT C REORTHOGONALIZATION 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', YIELD C APPROXIMATIONS TO THE DESIRED EIGENVALUES. CORRESPONDING C RITZ VECTORS ARE APPROXIMATIONS TO THE DESIRED EIGENVECTORS. C NOTE THAT FOR CASE (4) THE GENERALIZED LANCZOS RECURSION C B*V(I+1)*BETA(I+1) = A*V(I) - B*V(I)*ALPHA(I) - B*V(I-1)*BETA(I) C IS USED, ALONG WITH THE B-NORM. C C THE IDEAS USED IN THESE PROGRAMS ARE DISCUSSED IN THE FOLLOWING C REFERENCES. C C 1. JANE CULLUM AND RALPH A. WILLOUGHBY, LANCZOS ALGORITHMS C FOR LARGE SYMMETRIC MATRICES, VOLUME ?, PROGRESS IN C SCIENTIFIC COMPUTING, EDITORS, G. GOLUB, H.O. KREISS, C S. ARBARBANEL, AND R. GLOWINSKI, BIRKHAUSER BOSTON INC., C CAMBRIDGE, MASSACHUSETTS, 1983. C C 2. 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 3. 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 4. 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 WITH THE EXCEPTION OF THE PROGRAMS FOR HERMITIAN MATRICES WHICH C ARE NOT PORTABLE BECAUSE OF THEIR USE OF COMPLEX*16 VARIABLES, C THE OTHER PROGRAMS INCLUDED ARE PORTABLE EXCEPT FOR A FEW C CONSTRUCTIONS WHICH, IF NECESSARY, WILL HAVE TO BE MODIFIED C BY THE USER FOR THE PARTICULAR COMPUTER BEING USED. C C NONPORTABLE CONSTRUCTIONS: C C REAL SYMMETRIC MATRICES: C IN LEVAL AND IN LEVEC 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 ALPHA/BETA FILES. C IN LEMULT C 1. IN CMATV AND USPEC THE ENTRY THAT PASSES THE STORAGE C LOCATIONS OF THE ARRAYS DEFINING THE USER-SPECIFIED C MATRIX. C 2. IN THE SAMPLE USPEC PROVIDED: FREE FORMAT (8,*), C THE FORMAT (20A4), AND DATA/MACHEP STATEMENT. C C HERMITIAN MATRICES: C IN HLEVAL AND IN HLEVEC 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. COMPLEX*16 VARIABLES AND FUNCTIONS SUCH AS DCMPLX. C 5. FORMAT (4Z20) USED TO READ AND WRITE ALPHA/BETA FILES. C IN HLEMULT C 1. IN CMATV AND USPEC THE ENTRY THAT PASSES THE STORAGE C LOCATIONS OF THE ARRAYS DEFINING THE USER-SPECIFIED C MATRIX. C 2. COMPLEX*16 VARIABLES AND FUNCTIONS SUCH AS DCMPLX. C 3. IN THE SAMPLE USPEC PROVIDED: FREE FORMAT (8,*), C THE FORMAT (20A4), AND DATA/MACHEP STATEMENT. C C FACTORED INVERSES OF REAL SYMMETRIC MATRICES: C IN LIVAL AND IN LIVEC 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 ALPHA/BETA FILES. C IN LIMULT C 1. IN USPEC AND BSOLV, THE ENTRIES THAT PASS C THE STORAGE LOCATIONS OF THE ARRAYS DEFINING THE C USER-SPECIFIED MATRIX. C 2. IN THE SAMPLE USPEC SUBROUTINES PROVIDED: C FORMATS (20A4) AND (4Z20), FREE FORMAT (8,*), AND C DATA/MACHEP STATEMENTS. C C C GENERALIZED SYMMETRIC PROBLEM, B-MATRIX POSITIVE C DEFINITE AND CHOLESKY FACTORS AVAILABLE: C IN LGVAL AND IN LGVEC 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 ALPHA/BETA FILES. C IN LGMULT C 1. IN USPECA, USPECB, AMATV AND LSOLV THE ENTRIES C THAT PASS THE STORAGE LOCATIONS OF THE ARRAYS DEFINING C THE USER-SPECIFIED MATRICES. C 2. IN THE SAMPLE USPECA AND USPECB SUBROUTINES PROVIDED: C FORMATS (20A4) AND (4Z20), FREE FORMAT (8,*), AND C DATA/MACHEP STATEMENTS. C C ALL 4 CASES USE THE FORTRAN FILE LESUB: C IN LESUB ALL STATEMENTS ARE PORTABLE EXCEPT FOR: C (1) THE ENTRY IN SUBROUTINE LPERM THAT PASSES THE C PERMUTATION FROM THE USPEC SUBROUTINE TO LPERM. C (THIS IS USED ONLY IN CASES (3) AND (4)). C (2) THE COMPLEX*16 VARIABLES AND FUNCTIONS USED IN C SUBROUTINE CINPRD. (THIS IS USED ONLY IN CASE (2)). C C IN THE COMMENTS BELOW: C 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-----MATRIX SPECIFICATION---------------------------------------------- C C C IN CASES (1) AND (2), SUBROUTINE USPEC IS USED TO SPECIFY THE C USER-SUPPLIED A-MATRIX. SIMILARLY, IN CASE (4) SUBROUTINES C USPECA AND USPECB DEFINE THE USER-SUPPLIED A-MATRIX AND B-MATRIX. C IN CASE (3) ((4)), SUBROUTINE USPECB DEFINES THE FACTORIZATION C OF THE MATRIX (B-MATRIX) USED BY THE LANCZOS PROCEDURE. C (IN CASE (3) THE A-MATRIX IS NOT USED DIRECTLY.) C C IN CASES (1) AND (2), SUBROUTINE CMATV IS A CORRESPONDING C MATRIX-VECTOR MULTIPLY SUBROUTINE WHICH SHOULD BE DESIGNED C TO TAKE ADVANTAGE OF ANY SPECIAL PROPERTIES OF THE GIVEN C MATRIX. IN CASE (4) THIS SUBROUTINE IS NEEDED FOR THE C A-MATRIX AND THUS IS CALLED AMATV. IN CASES (3) AND (4) C SUBROUTINES THAT CAN SOLVE B*U = V, USING A SPARSE C FACTORIZATION OF B ARE NEEDED. THESE SUBROUTINES ARE C CALLED RESPECTIVELY, BSOLV AND LSOLV. IN ALL CASES, C ANY MATRIX-VECTOR MULTIPLY AND SOLVE SUBROUTINES USED C MUST BE DESIGNED TO COMPUTE RAPIDLY AND ACCURATELY. C C IN ALL CASES: C SUBROUTINE USPEC(A OR B) HAS THE CALLING SEQUENCE C C CALL USPEC(N,MATNO) C C WHERE N IS THE ORDER OF THE USER-SUPPLIED MATRIX A, AND C MATNO IS A <= 8 DIGIT INTEGER USED AS A MATRIX AND C TEST IDENTIFICATION NUMBER. IN ALL CASES THIS (THESE) C SUBROUTINE(S) DEFINES (DIMENSIONS) THE ARRAYS REQUIRED C TO SPECIFY THE MATRIX (MATRICES IN CASE (4)) THAT WILL BE C USED BY THE LANCZS SUBROUTINE. IN CASES (1) AND (2) C THIS IS THE A-MATRIX; IN CASE (3) THIS IS THE FACTORIZATION C OF A SCALED, SHIFTED AND PERMUTED VERSION OF THE C USER-SPECIFIED A-MATRIX. IN CASE (4) THE A-MATRIX C IS SPECIFIED AS WELL AS THE FACTORIZATION OF THE C B-MATRIX. THIS SUBROUTINE ALSO INITIALIZES THE ARRAYS C AND ANY OTHER PARAMETERS NEEDED TO DEFINE THE MATRIX C (MATRICES). THE STORAGE LOCATIONS OF THESE PARAMETERS C AND ARRAYS ARE THEN PASSED TO THE MATRIX-VECTOR MULTIPLY C SUBROUTINE CMATV IN CASES (1) AND (2), TO THE SUBROUTINE C BSOLV IN CASE (3), AND TO THE SUBROUTINES AMATV C AND LSOLV IN CASE (4) VIA ENTRY CALLS. IN CASES (3) AND (4) C WHENEVER A MATRIX HAS BEEN PERMUTED, THERE IS ALSO AN C ENTRY INTO THE SUBROUTINE LPERM TO PASS THE LOCATIONS OF C THE PERMUTATIONS IPR AND IPRT USED. SAMPLE USPECS, CMATV, C AMATV, BSOLV AND LSOLV SUBROUTINES ARE INCLUDED C IN THE RELEVANT FILES. THESE SAMPLE PROGRAMS ASSUME THAT C THE USER-SUPPLIED A-MATRIX IS STORED ON FILE 8 IN CASES (1), C (2), AND (4), AND THAT THE FACTORIZATION OF THE B-MATRIX C IS ON FILE 7 IN CASES (3) AND (4). THE USER SHOULD SEE C THE INDIVIDUAL SAMPLE SUBROUTINES FOR MORE DETAILS. C C IN CASES (1) AND (2): C SUBROUTINE CMATV HAS THE CALLING SEQUENCE C C CALL CMATV(W,U,SUM) C C IN THE REAL SYMMETRIC CASE, U AND W ARE REAL*8 VECTORS C AND SUM IS A REAL*8 SCALAR. IN THE HERMITIAN CASE, U C AND W ARE COMPLEX*16 VECTORS AND SUM IS A REAL*8 SCALAR. C CMATV CALCULATES U = A*W - SUM*U FOR THE USER-SPECIFIED C MATRIX A. ONE OF THE SAMPLE CMATV SUBROUTINES INCLUDED C COMPUTES MATRIX-VECTOR MULTIPLIES FOR AN ARBITRARY SPARSE, C SYMMETRIC MATRIX STORED IN THE SPARSE FORMAT SPECIFIED IN THE C CORRESPONDING SAMPLE USPEC SUBROUTINE. FOR CASES (1) AND C (2) CMATV IS THE SUBROUTINE USED BY THE LANCZS SUBROUTINE C THAT GENERATES THE T-MATRICES. IN CASE (4) SUBROUTINE C AMATV HAS THE SAME CALLING SEQUENCE AS CMATV IN CASE (1). C C IN CASES (3) AND (4): C ALPHA/BETA HISTORY IS GENERATED USING SPARSE MATRIX INVERSION. C IN CASE (3), AT EACH ITERATION OF THE LANCZOS RECURSION C GIVEN A FACTORIZATION OF THE MATRIX BEING USED, THE C SUBROUTINE BSOLV FOR A GIVEN V, COMPUTES U SUCH THAT B*U = V. C THE CALLING SEQUENCE OF BSOLV IS C C CALL BSOLV(V,U,IBSOLV) C C WHEN IBSOLV = 2, U = (B-INVERSE)*V IS RETURNED. IN CASE (4), C AT EACH ITERATION OF THE GENERALIZED LANCZOS RECURSION BOTH THE C SUBROUTINE AMATV AND THE SUBROUTINE LSOLV ARE USED. THE C CALLING SEQUENCE OF LSOLV IS C C CALL LSOLV(V,U,ISOLV) C C WHERE U AND V ARE REAL*8 VECTORS. LSOLV PERFORMS 4 FUNCTIONS. C LET L DENOTE THE CHOLESKY FACTOR OF THE B-MATRIX USED IN LANCZS. C WHEN ISOLV = 1, LSOLV COMPUTES U = L*V. WHEN ISOLV = 2, C LSOLV COMPUTES U = (L-TRANSPOSE)*V. WHEN ISOLV = 3, LSOLV C COMPUTES U = (L-INVERSE)*V. WHEN ISOLV = 4, LSOLV C COMPUTES U = ((L-TRANSPOSE)-INVERSE)*V. C C SAMPLE PROGRAMS ASSUME THAT THE A-MATRIX (CASES (1),(2),(4)) C IS ON FILE 8 AND STORED IN THE FOLLOWING SPARSE FORMAT: C ICOL(K), K = 1,NZL, NUMBER OF SUBDIAGONAL NONZEROS IN COLUMN K. C IROW(K), K = 1,NZS, ROW INDEX OF ASD(K). C AD(K), K=1,N, CONTAINS THE DIAGONAL ELEMENTS OF THE A-MATRIX. C ASD(K), K=1,NZS CONTAINS THE SUBDIAGONAL ELEMENTS OF A BY COLUMN. C NZS = NUMBER OF NONZERO ELEMENTS BELOW THE DIAGONAL OF A C NZL = INDEX OF LAST COLUMN WITH NONZERO SUBDIAGONAL ENTRIES C N = ORDER OF THE A-MATRIX. C C NOTE THAT THE OPTIONAL PREPROCESSING PROGRAMS PERMUT AND C LORDER ASSUME THAT THE GIVEN MATRIX IS ON FILE 8. CASES (3) C AND (4) ASSUME THAT THE SPARSE FACTORIZATION OF B IS STORED ON C FILE 7. THE SAMPLE BSOLV SUBROUTINE SUPPLIED ASSUMES C THAT THE B-MATRIX IS POSITIVE DEFINITE AND THAT ITS CHOLESKY C FACTOR IS PROVIDED ON FILE 7, STORED IN SPARSE FORMAT IN C ARRAYS BD AND BSD. THE USER CAN EASILY REPLACE THIS SAMPLE C BSOLV SUBROUTINE AND THE CORRESPONDING SAMPLE USPEC C SUBROUTINE BY SUBROUTINES THAT DEFINE AND USE A GENERAL C FACTORIZATION L*D*(L-TRANSPOSE). C C THE SAMPLE USPEC, CMATV (CASES (1) AND (2)), AMATV (CASE (4)), C BSOLV (CASE (3)), AND LSOLV (CASE(4)) MUST BE MODIFIED BY C THE USER TO ACCOMODATE THE USER-SPECIFIED MATRIX OR MATRICES. 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, AND C CASES (1) AND (2), CMATV: CASE (3), BSOLV: C CASE (4), AMATV AND LSOLV. C C GENRAN = COMPUTES K PSEUDO-RANDOM NUMBERS AND STORES THEM IN C THE REAL*4 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 N. 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 MATRIX THAT WILL BE USED BY LANCZS SUBROUTINE. C IN CASE (4) A-MATRIX AND B-MATRIX MUST BOTH BE SPECIFIED. C C CMATV = MATRIX-VECTOR MULTIPLY FOR USER-SUPPLIED MATRIX. C CASES (1) AND (2). SEE MATRIX SPECIFICATION SECTION. C C AMATV = MATRIX-VECTOR MULTIPLY FOR USER-SUPPLIED A-MATRIX. C CASES (4) ONLY. SEE MATRIX SPECIFICATION SECTION. C C BSOLV = GIVEN A VECTOR V COMPUTES U SUCH THAT B*U = V, C USING THE FACTORIZATION OF B. USED IN CASE (3) ONLY. C SEE MATRIX SPECIFICATION SECTION. C C LSOLV = PERFORMS 4 FUNCTIONS. GIVEN A VECTOR V COMPUTES C U = L*V, U = (L-TRANSPOSE)*V, U = (L-INVERSE)*V OR C U = (L-TRANSPOSE)-INVERSE*V, USING THE CHOLESKY C FACTORS OF B. USED ONLY IN CASE (4). SEE MATRIX C SPECIFICATION SECTION. C C C----------------------------------------------------------------------- C C COMMENTS FOR EIGENVALUE COMPUTATIONS C C----------------------------------------------------------------------- C C C-----PARAMETER CONTROLS FOR EIGENVALUE PROGRAMS------------------------ C C C PARAMETER CONTROLS ARE INTRODUCED TO ALLOW SEGMENTATION OF THE C EIGENVALUE COMPUTATIONS AND TO ALLOW VARIOUS COMBINATIONS OF C READ/WRITES. C C THE FLAG ISTART CONTROLS THE T-MATRIX (ALPHA/BETA HISTORY) C GENERATION. C C ISTART = (0,1) MEANS C C (0) THERE IS NO EXISTING ALPHA/BETA HISTORY AND ONE C MUST BE GENERATED. C C (1) THERE IS AN EXISTING ALPHA/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 EIGENVALUE COMPUTATIONS. C C ISTOP = (0,1) MEANS C C (0) PROGRAM COMPUTES ONLY THE REQUESTED ALPHAS/BETAS, C STORES THEM AND THE LAST 2 LANCZOS VECTORS GENERATED C IN FILE 1 AND THEN TERMINATES. IN CASE (4) THERE C ARE ACTUALLY 3 VECTORS TO BE SAVED. C C (1) PROGRAM COMPUTES REQUESTED ALPHAS/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 ALPHA/BETAS ARE NOT SAVED ON C FILE 1. C C (1) PROGRAM WRITES ALPHAS/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 ARE INTRODUCED BY FORMAT CONVERSIONS. 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 T-EIGENVALUES ALONG WITH THEIR MINIMAL GAPS AND WRITES THEM OUT C TO FILE 3. CORRESPONDING ERROR ESTIMATES FOR ANY ISOLATED C GOOD T-EIGENVALUES ARE ALWAYS WRITTEN TO FILE 4. C C C-----INPUT/OUTPUT FILES FOR EIGENVALUE PROGRAMS------------------------ C C ANY INPUT DATA OTHER THAN THE ALPHA/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 (ALPHA AND C BETA VECTORS) AND LAST 2 LANCZOS VECTORS USED C IN THE T-MATRIX GENERATION. NOTE THAT IN CASE (4) C THREE 'LANCZOS' VECTORS ARE WRITTEN TO FILE 1. 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 ALPHAS C AND BETAS ON FILE 2. THESE ALPHAS AND BETAS ARE C READ IN ALONG WITH THE LAST 2 LANCZOS VECTORS THAT C WERE GENERATED. IN CASE (4) THREE 'LANCZOS' VECTORS C ARE READ IN FROM FILE 2. C C (3) OUTPUT FILE: C COMPUTED GOOD EIGENVALUES OF THE T-MATRICES USED. ALSO C CONTAINS T-MULTIPLICITIES OF THESE EIGENVALUES AS C EIGENVALUES OF THE T-MATRIX, AND THEIR GAPS AS C EIGENVALUES IN THE A-MATRIX AND IN THE T-MATRIX. C FILE 3 IS ALWAYS WRITTEN. IN CASE (3) THIS OUTPUT C CONTAINS THE EIGENVALUES OF THE B-INVERSE MATRIX C SINCE IN THIS CASE THE T-MATRICES CORRESPOND TO C THE B-INVERSE MATRIX AND NOT TO THE A-MATRIX. IN C THIS CASE, 3 SETS OF GAPS ARE GIVEN, THOSE IN C THE T-MATRIX, IN THE B-INVERSE MATRIX AND THOSE C FOR THE CORRESPONDING EIGENVALUES IN THE A-MATRIX. C C (4) OUTPUT FILE: C ERROR ESTIMATES FOR THE ISOLATED GOOD T-EIGENVALUES C WHICH ARE OBTAINED USING THE SUBROUTINE INVERR. THESE C ESITMATES 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 (7) INPUT FILE: C USED ONLY IN CASES (3) AND (4), FACTORED INVERSES C OF REAL SYMMETRIC MATRICES AND GENERALIZED EIGENVALUE C PROBLEM. CONTAINS THE REQUIRED FACTORIZATION OF THE C B-MATRIX. 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. NOTE THAT IN CASE C (3) THE LANCZS SUBROUTINE DOES NOT USE THE MATRIX C ON FILE 8 IN THE T-MATRIX GENERATION, RATHER IT C USES THE FACTORIZATION OF AN ASSOCIATED C B-MATRIX WHICH IS STORED ON FILE 7. IN CASE (4), C THE INFORMATION STORED ON BOTH FILES 7 AND 8 IS USED. C C (9) INPUT AND OUTPUT FILE: C CAN BE USED TO STORE THE TRUE EIGENVALUES OF THE C GIVEN PROBLEM, WHEN THE 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 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 EIGENVALUE 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. THESE TOLERANCES OCCUR IN LUMP, C ISOEV, AND PRTEST. 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 EIGENVALUE IN MAGNITUDE TO THE SMALLEST IN MAGNITUDE IS VERY C LARGE, THEN THE TOLERANCES BEING USED IN BISEC, LUMP, ISOEV C AND PRTEST MAY NOT TREAT THE SMALL END (SMALL IN MAGNITUDE) 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 IN (ONLY) THE LOWER END OF THE C SPECTRUM WITH THIS TKMAX MAY RESULT IN IMPROVED COMPUTATIONS C AT THE LOW END. C C THE LUMP, ISOEV, AND PRTEST TOLERANCES THAT WERE USED C MOST IN THE TESTING OF THESE ALGORITHMS WERE NOT C SCALE INVARIANT BUT SEEMED TO WORK WELL ON MATRICES THAT C HAD EIGENVALUES WITH MAGNITUDES 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 THE LANCZOS EIGENVALUE PROGRAMS 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 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 WHICH OF THE GOOD T-EIGENVALUES NEED C ERROR ESTIMATES. THE PROGRAM SETS GAPTOL = 1.D-8. C IF FOR A GIVEN 'GOOD' T-EIGENVALUE THE COMPUTED GAP C IS TOO SMALL AND IS DUE TO A 'SPURIOUS' T-EIGENVALUE C THEN THE 'GOOD' T-EIGENVALUE IS ASSUMED TO HAVE CONVERGED C AND NO ERROR ESTIMATES ARE COMPUTED. C C C-----USER-SPECIFIED PARAMETERS FOR EIGENVALUE 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 TYPICALLY ONLY ONE ITERATION 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 EIGENVALUES ARE C 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 EIGENVALUE CALCULATIONS OVER ALL C 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 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 THE C EIGENVALUES OF THE T-MATRICES ON THE INTERVALS SPECIFIED BY C THE USER. C C IDEALLY, ONE WOULD COMPUTE THE EIGENVALUE APPROXIMATIONS AT A C REASONABLE SIZE T-MATRIX, LOOK AT THE ACCURACY OF THE COMPUTED C 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 EIGENVALUES THAT ARE DESIRED BUT THAT HAVE NOT YET C CONVERGED. HOWEVER, IN THE INTERESTS OF GENERALITY AND C SIMPLICITY WE DID NOT DO THAT HERE. C C C-----CONVERGENCE TESTS FOR THE EIGENVALUE PROGRAMS-------------------- C C C THE CONVERGENCE TEST INCORPORATED IN THIS PROGRAM IS C BASED UPON THE ASSUMPTION THAT THOSE T-EIGENVALUES AND THEIR C ASSOCIATED T-EIGENVECTORS WHICH CORRESPOND TO THE C EIGENVALUES AND RITZVECTORS WHICH WE WISH TO COMPUTE C CONVERGE AS THE T-SIZE IS INCREASED. C C AS CURRENTLY PROGRAMMED, CONVERGENCE IS CHECKED BY EXAMINING C THE SIZES OF ALL OF THE COMPUTED ERROR ESTIMATES ON ALL OF THE C INTERVALS SPECIFIED BY THE USER. IDEALLY CONVERGENCE SHOULD C BE CHECKED ONLY ON THOSE EIGENVALUES OF INTEREST AND C ONCE THE EIGENVALUES ON SUB-INTERVALS OF THESE INTERVALS HAVE C CONVERGED, ANY SUBSEQUENT EIGENVALUE COMPUTATIONS SHOULD BE C MADE ONLY ON THE UNCONVERGED PORTIONS. OBVIOUSLY, IT WOULD BE C DIFFICULT TO INCORPORATE CODE TO DO THE ABOVE WITHOUT KNOWING C A PRIORI PRECISELY WHAT THE USER IS TRYING TO COMPUTE. C THEREFORE, WE DID NOT ATTEMPT TO DO THIS. IF ONE WISHES TO C MAKE SUCH A MODIFICATION THEN ONE MUST ALSO MODIFY THE PROGRAM C SO THAT IT CREATES AN OVERALL LIST OF THE CONVERGED 'GOOD' C T-EIGENVALUES AS THEY ARE COMPUTED, SINCE CONVERGED 'GOOD' C T-EIGENVALUES WHICH WERE COMPUTED AT A PARTICULAR VALUE OF MEV C WOULD NO LONGER BE RECOMPUTED AT LARGER VALUES OF MEV. C C IF ONLY A FEW EIGENVALUES ARE TO BE COMPUTED THEN SUCH CHANGES C WOULD NOT MAKE MUCH DIFFERENCE IN THE RUNNING TIME. C C C-----ARRAYS REQUIRED BY THE EIGENVALUE PROGRAMS------------------------ C C C ALL 4 CASES C C ALPHA(J) = REAL*8 ARRAY. ITS DIMENSION MUST BE AT LEAST KMAX, C THE LENGTH OF THE LARGEST T-MATRIX ALLOWED. THIS C ARRAY CONTAINS THE DIAGONAL ENTRIES OF THE T-MATRICES. C C BETA(J) = REAL*8 ARRAY. ITS DIMENSION MUST BE AT LEAST KMAX+1. C THIS ARRAY CONTAINS THE SUBDIAGONAL ENTRIES OF THE C T-MATRICES. C C THE ALPHA AND BETA VECTORS ARE NOT ALTERED C DURING THE CALCULATIONS. C C V1(J),V2(J),VS(J) = REAL*8 ARRAYS IN REAL SYMMETRIC CASES. C V1 AND V2 ARE COMPLEX*16 IN HERMITIAN CASE. C IN CASES (1) AND (2) VS MUST BE OF C DIMENSION AT LEAST KMAX. IN CASES (3) AND C (4) VS MUST BE OF DIMENSION AT LEAST C MAX(N,KMAX). IN REAL SYMMETRIC CASES C V1 MUST BE OF DIMENSION AT LEAST C MAX(KMAX+1,N) AND V2 MUST BE OF DIMENSION C MAX(KMAX,N). IN HERMITIAN CASES V1 C MUST BE OF DIMENSION MAX(N,(KMAX+1)/2) C AND V2 OF DIMENSION AT LEAST MAX(N,KMAX/2). C IN ALL CASES HOWEVER, THE ABOVE DIMENSIONS C FOR V2 ARE VALID ONLY IF NO MORE C THAN KMAX/2 EIGENVALUES OF THE GIVEN C T-MATRICES ARE TO BE COMPUTED IN ANY GIVEN C SUBINTERVAL. V2 IS USED IN THE SUBROUTINE C BISEC TO HOLD THE UPPER AND LOWER C ENDPOINTS OF THE SUBINTERVALS GENERATED C DURING THE BISECTIONS. THEREFORE, ITS C REAL*8 DIMENSION MUST ALWAYS BE AT LEAST C 2*Q WHERE Q IS THE MAXIMUM NUMBER OF C EIGENVALUES OF THE SPECIFIED T-MATRIX IN ANY C ONE OF THE SPECIFIED INTERVALS. C NOTE THAT IN THE HERMITIAN CASE, V1 AND V2 C ARE DEFINED AS COMPLEX*16 IN THE MAIN PROGRAM C AND IN THE LANCZS SUBROUTINE BUT ARE C REDEFINED AS REAL*8 IN OTHER SUBROUTINES. C C LB(J),UB(J) = REAL*8 ARRAYS. EACH MUST BE OF DIMENSION AT LEAST C NINT, THE NUMBER OF SUBINTERVALS TO BE CONSIDERED. C LB CONTAINS THE LEFT-END POINTS OF THE INTERVALS C ON WHICH EIGENVALUES ARE TO BE COMPUTED. UB C CONTAINS THE RIGHT-END POINTS. C C EXPLAN(J) = REAL*4 ARRAY. ITS DIMENSION IS 20. THIS ARRAY IS C USED TO ALLOW EXPLANATORY COMMENTS IN THE INPUT FILES. C C G(J) = REAL*4 ARRAY. ITS DIMENSION MUST BE >= MAX(2*KMAX,N) C IT IS USED FOR HOLDING THE RANDOM VECTORS GENERATED, C HOLDING THE COMPUTED ERROR ESTIMATES AND THE COMPUTED C MINIMAL GAPS FOR THE GOOD T-EIGENVALUES. C C MP(J) = INTEGER*4 ARRAY. ITS DIMENSION MUST BE AT LEAST KMAX, C THE MAXIMUM SIZE OF THE T-MATRICES ALLOWED. IT CONTAINS C THE T-MULTIPLICITIES OF THE COMPUTED EIGENVALUES. NOTE C THAT 'SPURIOUS' T-EIGENVALUES ARE DENOTED BY A C T-MULTIPLICITY OF 0. T-EIGENVALUES THAT THE SUBROUTINE C PRTEST HAS IDENTIFIED AS 'GOOD' BUT HIDDEN ARE IDENTIFIED C BY A T-MULTIPLICITY OF -10. C C NMEV(J) = INTEGER*4 ARRAY. ITS DIMENSION MUST BE AT LEAST THE C NUMBER OF T-MATRICES ALLOWED. IT CONTAINS THE ORDERS C OF THE T-MATRICES TO BE CONSIDERED. C C C FOR CASE (3) ONLY: C GR(J),GC(J) = REAL*8 ARRAYS. USED ONLY IN THE HERMITIAN CASE. C GR AND GC MUST EACH BE OF DIMENSION AT LEAST N. C BOTH ARE USED IN THE RANDOM VECTOR GENERATION. C GR IS ALSO USED TO STORE MINIMAL GAPS BETWEEN C 'GOOD' T-EIGENVALUES. C C FOR CASES (3) AND (4) FOR THE PERMUTATION: C C IPR(J), IPT(J) = INTEGER*4 ARRAYS. EACH OF DIMENSION AT LEAST N. C USED TO STORE THE REORDERING OF THE GIVEN MATRIX C OR MATRICES. C C C OTHER ARRAYS C C THE USER MUST SPECIFY IN THE SUBROUTINE USPEC (A OR B) WHATEVER C ARRAYS ARE REQUIRED TO DEFINE THE MATRIX OR MATRICES BEING USED. C ALSO IN CASES (3) AND (4) ONLY, WHEN WORKING WITH INVERSES C OF SPARSE SYMMETRIC MATRICES, IN THE OPTIONAL, PREPROCESSING C PROGRAMS PERMUT, LFACT, LORDER, AND LTEST IT IS NECESSARY TO C SPECIFY ADDITIONAL ARRAYS JUST FOR THESE COMPUTATIONS. THE USER C IS REFERRED TO THOSE PROGRAMS FOR DETAILS. C C---------------------------------------------------------------------- C C-----SUBROUTINES INCLUDED---------------------------------------------- C C C LANCZS = COMPUTES THE ALPHA/BETA HISTORY. IN CASES (1) AND (2) C REAL SYMMETRIC AND HERMITIAN MATRICES, USES SUBROUTINES C FINPRO, GENRAN AND CMATV. IN CASE (3), INVERSES OF C REAL SYMMETRIC MATRICES, USES SUBROUTINES FINPRO, C GENRAN AND BSOLV. IN CASE (4), GENERALIZED EIGENVALUE C PROBLEM, USES SUBROUTINES FINPRO, GENRAN, AMATV AND C LSOLV. C C BISEC = COMPUTES EIGENVALUES OF THE SPECIFIED T-MATRIX C USING STURM SEQUENCING, ON SEQUENCE OF INTERVALS C SPECIFIED BY THE USER. EACH SUBINTERVAL IS TREATED C AS OPEN ON THE LEFT AND CLOSED ON THE RIGHT. C EIGENVALUES ARE COMPUTED WITH SIMULTANEOUS DETERMINATION C OF THE T-MULTIPLICITIES AND OF SPURIOUS T-EIGENVALUES. C C INVERR = USES INVERSE ITERATION ON T-MATRICES TO COMPUTE ERROR C ESTIMATES ON COMPUTED GOOD T-EIGENVALUES. (USES GENRAN) C C LUMP = 'COMBINES' EIGENVALUES OF T-MATRIX USING THE RELATIVE C TOLERANCE RELTOL. C C ISOEV = CALCULATES GAPS BETWEEN DISTINCT EIGENVALUES OF T-MATRIX C AND THEN USES THESE GAPS TO LABEL THOSE 'GOOD' C T-EIGENVALUES FOR WHICH ERROR ESTIMATES ARE NOT COMPUTED. C C TNORM = COMPUTES THE SCALE TKMAX USED IN DETERMINING THE C TOLERANCES FOR THE SPURIOUS, T-MULTIPLICITY AND PRTESTS. C IT ALSO CHECKS FOR LOCAL ORTHOGONALITY OF THE LANCZOS C VECTORS BY TESTING THE RELATIVE SIZE OF THE BETAS USING C THE RELATIVE TOLERANCE BTOL. C C PRTEST = LOOKS FOR GOOD T-EIGENVALUES THAT HAVE BEEN MISLABELLED C BY THE SPURIOUS TEST BECAUSE THEY HAD 'TOO SMALL' A C PROJECTION ON THE STARTING LANCZOS VECTOR. C (LESS THAN SINGLE PRECISION) C TESTS INDICATE THAT SUCH EIGENVALUES ARE RARE. C PRTEST SHOULD BE CALLED ONLY AFTER CONVERGENCE C HAS BEEN ESTABLISHED. C C INVERM = USED TO COMPUTE ERROR ESTIMATES FOR ANY T-EIGENVALUES C WHICH PRTEST INDICATES MAY HAVE BEEN MISLABELLED. C SUCH T-EIGENVALUES ARE RELABELLED ONLY IF THEIR ERROR C ESTIMATES ARE SUFFICIENTLY SMALL. PRIMARY USE OF C INVERM IS IN THE CORRESPONDING EIGENVECTOR COMPUTATIONS. C C CASES (3) AND (4) ONLY, FACTORED INVERSES: C C FOR OPTIONAL, PRELIMINARY PROCESSING: C PERMUT (PROGRAM CALLS SPARSPAK PACKAGE) : C USES THE NONZERO STRUCTURE OF A GIVEN MATRIX A. C CAN BE USED TO OBTAIN A REORDERING OF A THAT PRESERVES C SPARSITY UNDER FACTORIZATION. PERMUT CALLS C THE SPARSPAK PROGRAMS, (A. GEORGE, J. LIU, E.NG, C U. WATERLOO). PERMUT ALSO TAKES THE USER-SPECIFIED MATRIX, C APPLIES THE SCALE S0 AND THE SHIFT TO IT, AND THEN WRITES C OUT THE CORRESPONDING SPARSE MATRIX DATA FILE FOR THE C RESULTING MATRIX C = S0*A + SHIFT*I. SEE THE PERMUT FORTRAN C CODE FOR DETAILS. C C LORDER (STAND-ALONE PROGRAM): C GIVEN A MATRIX C IN SPARSE FORMAT AND A PERMUTATION P, C COMPUTES THE REORDERED MATRIX B = P*C*P' AND WRITES IT C TO FILE 9 IN SPARSE FORMAT. SEE THE LORDER FORTRAN CODE C FOR DETAILS. C C LFACT (STAND-ALONE PROGRAM) : C GIVEN A POSITIVE DEFINITE MATRIX B IN SPARSE FORMAT C COMPUTES THE SPARSE CHOLESKY FACTOR L OF B AND WRITES IT C TO FILE 7 IN SPARSE FORMAT. THUS, B = L*L'. C SEE THE LFACT FORTRAN CODE FOR DETAILS. C C LTEST (CALLS 3 USER-SUPPLIED PROGRAMS CMATV, CMATS, AND BSOLV): C GIVEN THE FACTORIZATION OF A MATRIX B, LTEST COMPUTES C THE SOLUTION OF THE EQUATION B*U = B*V1 FOR A SPECIFIC RANDOMLY- C GENERATED V1, WITH AND WITHOUT ITERATIVE REFINEMENT, TO C OBTAIN A ROUGH CHECK ON THE NUMERICAL CONDITION OF THE MATRIX B. C THIS PROGRAM USES 3 SUBROUTINES CMATV, CMATS, AND BLSOLV. C SEE THE LTEST FORTRAN PROGRAM FOR DETAILS. C C C-----OTHER PROGRAMS PROVIDED------------------------------------------- C C LECOMPAC (STAND ALONE PROGRAM): C TRANSLATES A REAL SYMMETRIC MATRIX PROVIDED IN THE C FORMAT I, J, A(I,J) INTO THE SPARSE MATRIX C FORMAT USED IN THE SAMPLE USPEC, CMATV, BSOLV AND C LSOLV SUBROUTINES PROVIDED. IT ASSUMES THAT THE C MATRIX ENTRIES ARE GIVEN EITHER COLUMN BY COLUMN OR C ROW BY ROW. THE DATA SET CREATED IS WRITTEN TO C FILE 8. C C C-----COMMENTS ON THE STORAGE REQUIRED FOR EIGENVALUE PROGRAMS---------- C C C CASE (1), REAL SYMMETRIC MATRICES: C C THE ARRAYS IN THE REAL SYMMETRIC EIGENVALUE PROGRAM REQUIRE C APPROXIMATELY THE EQUIVALENT OF ONE REAL*8 ARRAY OF DIMENSION C C 3.5*KMAX + 2*MAX(KMAX,N) + .5* MAX(2*KMAX,N) C C PLUS WHATEVER IS NEEDED TO GENERATE A*X FOR THE GIVEN MATRIX A. C THE ARRAYS ALPHA, BETA, VS AND MP CONSUME 3.5*KMAX*8 BYTES. C THE ARRAYS V1 AND V2 CONSUME 2*MAXIMUM(KMAX,N)*8 BYTES, WITH THE C QUALIFICATION STATED ABOVE WHERE V2 IS DEFINED. THE G-ARRAY C CONSUMES .5*MAX(2*KMAX,N)*8 BYTES. C C CASE (2), HERMITIAN MATRICES: C C THE ARRAYS IN THE HERMITIAN EIGENVALUE PROGRAMS REQUIRE C THE EQUIVALENT OF ONE REAL*8 ARRAY OF DIMENSION C C 3.5*KMAX + 4*MAX(KMAX/2,N) + .5*MAX(2*KMAX,N) + 2*N C C PLUS WHATEVER IS NEEDED TO GENERATE A*X FOR THE GIVEN MATRIX A. C THE ARRAYS ALPHA, BETA, VS, AND MP CONSUME 3.5*KMAX*8 BYTES. C THE ARRAYS V1 AND V2 CONSUME 4*MAXIMUM(KMAX/2,N)*8 BYTES, WITH THE C QUALIFICATION STATED ABOVE WHERE V2 IS DEFINED. THE G-ARRAY C CONSUMES .5*MAX(2*KMAX,N)*8 BYTES. GR REQUIRES C AND GC REQUIRE 2*N*8BYTES. C C C CASE (3), INVERSES OF REAL SYMMETRIC MATRICES: C C THE ARRAYS IN THE EIGENVALUE PROGRAMS DESIGNED FOR C CASE (3), INVERSES OF REAL SYMMETRIC MATRICES USING C REORDERING AND FACTORIZATION, REQUIRE C THE EQUIVALENT OF ONE REAL*8 ARRAY OF DIMENSION C C 3*KMAX + 3*MAX(KMAX,N) + .5*MAX(2*KMAX,N) C C PLUS WHATEVER IS NEEDED TO GENERATE B(INVERSE)*X FOR THE C SCALED, SHIFTED AND PERMUTED VERSION OF A WHICH WE DENOTE C BY B. THE ARRAYS ALPHA, BETA, MP, AND MP2 CONSUME 3*KMAX*8 C BYTES. THE ARRAYS V1, V2, AND VS CONSUME 3*MAX(KMAX,N)*8 BYTES, C WITH THE QUALIFICATION STATED ABOVE WHERE V2 IS DEFINED. C THE G ARRAY CONSUMES .5*MAX(2*KMAX,N)*8 BYTES. THESE NUMBERS C DO NOT INCLUDE THE STORAGE REQUIRED BY THE PREPROCESSING PROGRAMS C PERMUT, LORDER, LFACT, AND LTEST. C C C A SYMMETRIC, SPARSE MATRIX OF ORDER N WITH NZS NONZERO ELEMENTS C BELOW THE MAIN DIAGONAL WOULD REQUIRE THE EQUIVALENT OF ONE C REAL*8 ARRAY OF DIMENSION 1.5*(NZS + N) IF THE POINTERS USED C ARE INTEGER*4. C C SOME OF THE ARRAY STORAGE IS NOT ESSENTIAL AND COULD BE C ELIMINATED IF STORAGE IS A PROBLEM. C THE FOLLOWING COMMENTS APPLY DIRECTLY ONLY TO CASE (1), C THE PROGRAMS FOR REAL SYMMETRIC MATRICES, HOWEVER, SIMILAR C STATEMENTS COULD BE MADE ABOUT THE OTHER CASES. C C CASE (1), REAL SYMMETRIC PROGRAMS: C THE G ARRAY COULD BE REMOVED IF THE USER IS WILLING TO C C (1) REGENERATE THE RANDOM STARTING VECTOR IN INVERR C FOR EACH ERROR ESTIMATE C (2) WRITE OUT THE ERROR ESTIMATES AND VARIOUS GAPS AS C THEY ARE GENERATED RATHER THAN STORING THEM IN G FOR C LATER PRINTOUT C (3) CHECK CONVERGENCE WITHIN INVERR C C CLEARLY THE INDEX VECTOR MP COULD BE AN INTEGER*2 ARRAY AS COULD C THE POINTERS USED TO DEFINE THE USER'S MATRIX. C C THE USER SHOULD NOTE THAT WITH AN EIGENVALUE SUBROUTINE THAT C USES BISECTION (LIKE BISEC) IF MORE THAN 25% OF THE C EIGENVALUES ARE TO BE COMPUTED, THEN IT MAY BE MORE C ECONOMICAL TO USE THE EISPACK SUBROUTINE IMTQL1. C (SEE MATRIX EIGENSYTEM ROUTINES-EISPACK GUIDE (2ND EDITION) C B.T. SMITH ET AL, SPRINGER-VERLAG, NEW YORK, 1976, P213.). C HOWEVER, IF THE SUBROUTINE IMTQL1 IS TO BE USED IN PLACE C OF BISEC, THEN NONTRIVIAL CHANGES IN THE LANCZOS CODE MUST BE C MADE. FOR DETAILS OF ONE SUCH IMPLEMENTATION SEE C IBM RESEARCH REPORT 8298, COMPUTING C EIGENVALUES OF LARGE SYMMETRIC MATRICES - AN IMPLEMENTATION OF A C LANCZOS ALGORITHM WITH NO REORTHOGONALIZATION. PART II. COMPUTER C PROGRAMS., DECEMBER 1980, WHICH CONTAINS A GENERAL C LANCZOS CODE WHICH INCLUDES AN IMTQL1 OPTION OR C PREFERABLY CONTACT THE AUTHORS. C C THE BISEC SUBROUTINE WHICH IS INCLUDED IS A MODIFIED FORM OF C THE BISECT SUBROUTINE IN EISPACK. BISEC ASSUMES THAT THE C VECTOR V2 IS LONG ENOUGH TO HOLD BOTH THE UPPER AND THE C LOWER BOUNDS ON THE BISECTION INTERVALS USED TO COMPUTE C THE EIGENVALUES OF THE T-MATRICES. THEREFORE, IF THE C LENGTH OF V2 IS ONLY KMAX, BISEC CAN COMPUTE ONLY AT MOST C KMAX/2 EIGENVALUES OF THE GIVEN T-MATRIX IN ANY GIVEN C SUBINTERVAL. C C AS PROGRAMMED BISEC USES THE ARRAYS ALPHA,BETA,V1,V2,VS AND MP. C HOWEVER, V1 IS USED ONLY TO STORE BETA(J)**2 SO THAT THEY DO NOT C HAVE TO BE REGENERATED ON EACH STURM. IF THE USER IS WILLING TO C COMPUTE THE BETA(J)**2 AS NEEDED, THEN V1 COULD BE ELIMINATED C FROM BISEC. BISEC STORAGE THEN BECOMES A REAL*8 ARRAY OF DIMENSION C 4.25*KMAX IF WE ALSO REDUCE MP TO INTEGER*2. FURTHERMORE, C IF ONE KNEW THAT ONLY Q*MEV EIGENVALUES OF T(1,MEV) WERE TO BE C COMPUTED AT EACH STAGE FOR SOME Q<.5 THEN FURTHER REDUCTIONS IN C STORAGE COULD BE MADE IN BISEC. C C AS PROGRAMMED INVERR USES ALPHA, BETA,V1,V2,VS,G AND MP. C VS CONTAINS THE COMPUTED EIGENVALUES OF T(1,MEV). MP GIVES C THEIR T-MULTIPLICITIES AND FLAGS WHICH EIGENVALUES ARE TO HAVE C ERROR ESTIMATES COMPUTED. V2 IS USED FOR THE SOLUTION C VECTOR IN THE INVERSE ITERATION AND V1 FOR THE FACTORIZATION. C G CONTAINS THE RANDOMLY-GENERATED STARTING VECTOR FOR THE C INVERSE ITERATION. THE BASIC STORAGE FOR INVERR IS THEREFORE C A REAL*8 ARRAY OF DIMENSION 4*KMAX PLUS THE STORAGE NEEDED FOR C THE COMPUTED T-EIGENVALUES AND THEIR T-MULTIPLICITIES. C C VS COULD BE USED TO STORE ONLY THOSE COMPUTED EIGENVALUES OF C T(1,MEV) THAT ARE OF INTEREST. IN THAT CASE THE DIMENSIONS OF VS C AND OF MP NEED NOT BE ANY LONGER THAN THE NUMBER OF SUCH C EIGENVALUES. AS PROGRAMMED, ALL THE COMPUTED DISTINCT EIGENVALUES C OF T(1,MEV) ARE STORED IN VS. THEREFORE TO TAKE ADVANTAGE OF C SUCH A REDUCTION IN STORAGE THE USER WOULD HAVE TO MODIFY THAT C PART OF THE PROGRAM AND ALSO COMMENT OUT THE CALL TO THE C SUBROUTINE PRTEST. C C C----------------------------------------------------------------------- C C COMMENTS FOR EIGENVECTOR COMPUTATIONS C C----------------------------------------------------------------------- C C C THE EIGENVALUES WHOSE EIGENVECTORS ARE TO BE COMPUTED MUST C HAVE BEEN COMPUTED USING THE CORRESPONDING LANCZOS EIGENVALUE C PROGRAMS BECAUSE THE EIGENVECTOR PROGRAMS WILL USE THE SAME C FAMILY OF LANCZOS TRIDIAGONAL MATRICES THAT WAS USED IN THE C CORRESPONDING EIGENVALUE COMPUTATIONS. C C THESE PROGRAMS ASSUME THAT THE EIGENVALUES SUPPLIED TO IT C HAVE BEEN COMPUTED ACCURATELY, AS MEASURED BY THE C ERROR ESTIMATES COMPUTED IN THE CORRESPONDING LANCZOS C EIGENVALUE COMPUTATIONS, ALTHOUGH THESE ESTIMATES ARE C TYPICALLY CONSERVATIVE. IN CASES (1), (2) AND (4), THE C EIGENVALUES OF INTEREST ARE STORED IN THE ARRAY GOODEV(J), C J=1,NGOOD. IN CASE (3) THE PROGRAM WORKS WITH THE C EIGENVALUES OF B(INVERSE) WHICH ARE STORED IN THE ARRAY C GOODBI(J), J=1,NGOOD. THE CORRESPONDING EIGENVALUES C OF A ARE STORED IN GOODA(J), J=1,NGOOD. C C FOR EACH GOODEV(J), THE SUBROUTINE STURMI COMPUTES THE C SMALLEST SIZE LANCZOS TRIDIAGONAL MATRIX, T(1,M1(J)), FOR C WHICH GOODEV(J) IS AN EIGENVALUE TO WITHIN A SPECIFIED C TOLERANCE. IT ALSO ATTEMPTS TO COMPUTE THE SIZE, M2(J), C BY WHICH THE GIVEN EIGENVALUE BECOMES A DOUBLE EIGENVALUE C TO WITHIN THE GIVEN TOLERANCE. THESE VALUES ARE USED C TO DETERMINE 1ST GUESSES AT SIZES FOR THE T-EIGENVECTORS C THAT WILL BE USED IN THE RITZ VECTOR COMPUTATIONS. C SUBROUTINE INVERM SUCCESSIVELY COMPUTES CORRESPONDING C EIGENVECTORS OF ENLARGED T-MATRICES UNTIL A SUITABLE C SIZE T-MATRIX IS DETERMINED FOR EACH J. UP TO 10 SUCH C EIGENVECTOR COMPUTATIONS ARE ALLOWED FOR EACH EIGENVALUE. C C AFTER APPROPRIATE T-EIGENVECTORS HAVE BEEN COMPUTED, C RITZ VECTOR CORRESPONDING TO THESE T-EIGENVECTORS ARE THEN C COMPUTED AND TAKEN AS APPROXIMATE EIGENVECTORS OF A FOR THE C GIVEN EIGENVALUES, GOODEV(J), J = 1, ..., NGOOD. C C THIS IMPLEMENTATION FIRST COMPUTES ALL OF THE RELEVANT C EIGENVECTORS OF THE SYMMETRIC TRIDIAGONAL MATRICES C IN THE VECTOR, TVEC. C C THEN, AS EACH OF THE LANCZOS VECTORS IS REGENERATED, ALL C OF THE RITZ VECTORS CORRESPONDING TO THESE C T-EIGENVECTORS ARE UPDATED USING THE CURRENTLY-GENERATED C LANCZOS VECTOR. LANCZOS VECTORS ARE GENERATED (NOTE C THAT THEY ARE NOT BEING KEPT), UNTIL ENOUGH HAVE C BEEN GENERATED TO MAP THE LONGEST T-EIGENVECTOR INTO ITS C CORRESPONDING RITZ VECTOR. THE ARRAY RITVEC CONTAINS THE C SUCCESSIVE RITZ VECTORS WHICH ARE THE APPROXIMATE C EIGENVECTORS OF A. C C C-----PARAMETER CONTROLS FOR EIGENVECTOR PROGRAMS----------------------- C C C IN CASES (3) AND (4) WHERE A SPARSE FACTORIZATION OF A C SPECIFIED MATRIX IS USED, THE USER SPECIFIES USING THE FLAG C JPERM WHETHER OR NOT THE FACTORIZATION SUPPLIED CORRESPONDS C TO THE ORIGINAL MATRIX OR TO A PERMUTATION OF THE ORIGINAL C MATRIX. C C JPERM = (0,1) MEANS C 0 NO PERMUTATION C 1 MATRIX HAS BEEN PERMUTED. NOTE THAT IN C CASE (4) THE PROGRAM WILL ASSUME THAT THE C DATA SUPPLIED FOR THE A-MATRIX CORRESPONDS TO THE C CORRESPONDING PERMUTATION OF THE ORIGINAL A-MATRIX. C IN BOTH CASES THE LANCZS CODES WILL WORK WITH THE C PERMUTED MATRICES AND THE PERMUTATION WILL BE C UNDONE ONLY IN THE EIGENVECTOR PROGRAM AFTER C THE RITZ VECTORS FOR THE PERMUTED PROBLEM HAVE C BEEN COMPUTED. C C OTHER PARAMETER CONTROLS ARE INTRODUCED TO ALLOW SEGMENTATION C OF THE EIGENVECTOR COMPUTATIONS AND TO ALLOW VARIOUS COMBINATIONS C OF READ/WRITES. C C THE FLAG MBOUND ALLOWS THE USER TO DETERMINE A FIRST GUESS ON THE C STORAGE THAT WILL BE REQUIRED BY THE T-EIGENVECTORS FOR THE C EIGENVALUES WHOSE EIGENVECTORS ARE TO BE COMPUTED. C THIS CAN BE USED TO ESTIMATE THE REQUIRED SIZE OF THE TVEC ARRAY. C C MBOUND = (0,1) MEANS C C (0) PROGRAM COMPUTES FIRST GUESSES AT THE SIZES C OF THE T-MATRICES REQUIRED BY EACH OF THE C EIGENVALUES SUPPLIED AND THEN CONTINUES WITH C THE CORRESPONDING T-EIGENVECTOR COMPUTATIONS. C C (1) PROGRAM COMPUTES FIRST GUESSES AT THE SIZES C OF THE T-MATRICES REQUIRED BY EACH OF THE C EIGENVALUES SUPPLIED, STORES THESE IN FILE 10 C AND THEN TERMINATES. THE USER CAN USE THESE C SIZES TO ESTIMATE THE SIZE TVEC ARRAY NEEDED C FOR THE DESIRED T-EIGENVECTOR COMPUTATIONS. C C THE FLAGS NTVCON, TVSTOP, LVCONT, AND ERCONT CONTROL THE STOPPING C CRITERIA FOR INTERMEDIATE POINTS IN THE LANCZOS PROCEDURE. THEY C TERMINATE THE PROCEDURE IF VARIOUS QUANTITIES COULD NOT BE C COMPUTED AS DESIRED. C C NTVCON = (0,1) MEANS C C (0) IF THE ESTIMATED STORAGE FOR THE T-EIGENVECTORS C EXCEEDS THE USER-SPECIFIED DIMENSION OF THE C TVEC ARRAY PROGRAM DOES NOT CONTINUE WITH THE C T-EIGENVECTOR COMPUTATIONS. TERMINATION OCCURS. C C (1) CONTINUE WITH THE T-EIGENVECTOR COMPUTATIONS C EVEN IF THE ESTIMATED STORAGE FOR TVEC EXCEEDS C THE USER-SPECIFIED DIMENSION OF THE TVEC ARRAY. C IN THIS SITUATION THE PROGRAM COMPUTES AS MANY C T-EIGENVECTORS AS IT HAS ROOM FOR, IN THE SAME C ORDER IN WHICH THE EIGENVALUES ARE PROVIDED. C C SVTVEC = (0,1) MEANS C C (0) DO NOT STORE THE COMPUTED T-EIGENVECTORS ON C FILE 11 UNLESS ALSO HAVE THE FLAG TVSTOP = 1, C IN WHICH CASE THE T-EIGENVECTORS ARE ALWAYS C WRITTEN TO FILE 11. C C (1) STORE THE COMPUTED T-EIGENVECTORS ON FILE 11. C C TVSTOP = (0,1) MEANS C C (0) ATTEMPT TO CONTINUE ON TO THE COMPUTATION C OF THE RITZVECTORS AFTER COMPLETING THE C COMPUTATION OF THE T-EIGENVECTORS. C C (1) TERMINATE AFTER COMPUTING THE C T-EIGENVECTORS AND STORING THEM ON FILE 11. C C LVCONT = (0,1) MEANS C C (0) IF SOME OF THE T-EIGENVECTORS THAT WERE C REQUESTED WERE NOT COMPUTED, EXIT C FROM THE PROGRAM WITHOUT COMPUTING THE C CORRESPONDING RITZ VECTORS. C C (1) CONTINUE ON TO THE RITZ VECTOR COMPUTATIONS C EVEN IF NOT ALL OF THE T-EIGENVECTORS C REQUESTED WERE COMPUTED. C C ERCONT = (0,1) MEANS C C (0) PROCEDURE WILL NOT COMPUTE A RITZ VECTOR FOR C ANY EIGENVALUE FOR WHICH NO T-EIGENVECTOR WHICH C SATISFIES THE ERROR ESTIMATE TEST (ERTOL) HAS C BEEN IDENTIFIED. C C (1) A RITZ VECTOR WILL BE COMPUTED FOR EVERY C EIGENVALUE FOR WHICH A T-EIGENVECTOR HAS BEEN C COMPUTED REGARDLESS OF WHETHER OR NOT THAT C T-EIGENVECTOR SATISFIED THE ERROR ESTIMATE TEST. C C C-----INPUT/OUTPUT FILES FOR THE EIGENVECTOR COMPUTATIONS--------------- C C C ANY INPUT DATA OTHER THAN THE T-MATRIX HISTORY FILE AND THE C PREVIOUSLY COMPUTED EIGENVALUES AND CORRESPONDING ERROR C ESTIMATES SHOULD BE STORED ON FILE 5 IN FREE FORMAT. C SEE SAMPLE INPUT/OUTPUT FOR TYPICAL INPUT FILE. 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. ADDITIONAL PRINTOUT IS GENERATED WHEN C THE FLAG IWRITE = 1. C C C DESCRIPTION OF OTHER I/O FILES C C FILE (K) CONTAINS: C C (2) INPUT FILE: C PREVIOUSLY-GENERATED T-MATRICES (ALPHA/BETA ARRAYS) C AND THE FINAL TWO LANCZOS VECTORS USED ON THAT C COMPUTATION. THIS PROGRAM ALLOWS ENLARGEMENT C OF ANY T-MATRICES PROVIDED ON FILE 2. NOTE THAT C IN CASE (4), THREE 'LANCZOS' VECTORS ARE ON FILE 2. C C (3) INPUT FILE: C THE GOOD EIGENVALUES OF THE T-MATRIX T(1,MEV) C FOR WHICH RITZ VECTORS ARE REQUESTED. C FILE 3 ALSO CONTAINS THE T-MULTIPLICITIES OF THESE C EIGENVALUES AND THEIR COMPUTED GAPS IN THE C T-MATRICES AND IN THE USER-SUPPLIED MATRIX. IN C CASE (3) FILE 3 CONTAINS THE EIGENVALUES OF THE C B-INVERSE MATRIX AND 3 SETS OF CORRESPONDING GAPS. C THIS FILE IS CREATED IN THE LANCZOS EIGENVALUE C COMPUTATIONS. C C (4) INPUT FILE: C ERROR ESTIMATES FOR THE ISOLATED GOOD T-EIGENVALUES C ON FILE 3. THIS FILE IS CREATED DURING THE LANCZOS C EIGENVALUE COMPUTATIONS. C C (7) INPUT FILE: C IN CASE (3) ((4)), C CONTAINS SPARSE MATRIX REPRESENTATION OF FACTORIZATION C OF MATRIX (B-MATRIX) USED BY LANCZS SUBROUTINE. C C (8) INPUT FILE: C IN CASES (1),(2) AND (4), USPEC SUBROUTINE ASSUMES C THAT USER-SUPPLIED A-MATRIX IS ON FILE 8. IN CASE (3) C A-MATRIX CAN BE STORED ON FILE 8, BUT IT IS NOT C USED BY THE LANCZOS PROGRAMS. C C (9) OUTPUT FILE: C ERROR ESTIMATES FOR THE COMPUTED RITZ VECTORS CONSIDERED C AS EIGENVECTORS OF THE MATRIX USED BY THE LANCZS C SUBROUTINE. THESE ESTIMATES ARE OF THE FORM C AERROR = || A*RITVEC - EVAL*RITVEC || C WHERE A DENOTES THE MATRIX USED BY LANCZS, EVAL DENOTES C THE EIGENVALUE BEING CONSIDERED AND RITVEC DENOTES C THE COMPUTED RITZ VECTOR. C C (10) OUTPUT FILE: C GUESSES AT APPROPRIATE SIZE T-MATRICES FOR THE C T-EIGENVECTORS FOR EACH SUPPLIED EIGENVALUE, GOODEV(J). C C (11) OUTPUT FILE: C COMPUTED T-EIGENVECTORS CORRESPONDING TO EIGENVALUES C IN THE GOODEV ARRAY. NOTE THAT IT IS POSSIBLE IN C CERTAIN SITUATIONS THAT FOR SOME EIGENVALUES IN THE C GOODEV ARRAY A T-EIGENVECTOR WILL NOT BE COMPUTED. C C (12) OUTPUT FILE: C CONTAINS COMPUTED RITZ VECTORS CORRESPONDING TO C THE T-EIGENVECTORS ON FILE 11. NOTE THAT IN C SOME SITUATIONS THAT FOR SOME EIGENVALUES IN C THE GOODEV ARRAY FOR WHICH T-EIGENVECTORS HAVE C BEEN COMPUTED NO RITZ VECTOR WILL HAVE BEEN C COMPUTED. C C (13) OUTPUT FILE: C ADDITIONAL INFORMATION ABOUT THE BOUNDS AND ERROR C ESTIMATES OBTAINED. C C C-----SEEDS FOR EIGENVECTOR PROGRAMS------------------------------------ C C SEEDS FOR RANDOM NUMBER GENERATOR GENRAN C (1) SVSEED = INTEGER*4 SCALAR USED IN THE SUBROUTINE C GENRAN TO GENERATE THE STARTING VECTOR FOR C THE REGENERATION OF THE LANCZOS VECTORS. C C (2) RHSEED = INTEGER*4 SCALAR USED IN THE SUBROUTINE C GENRAN TO GENERATE A RANDOM VECTOR FOR C USE IN SUBROUTINE INVERM. C C USER SHOULD NOTE THAT SVSEED MUST BE THE SAME SEED THAT C WAS USED TO GENERATE THE T-MATRICES THAT WERE USED TO C COMPUTE THE EIGENVALUES WHOSE EIGENVECTORS ARE TO BE COMPUTED. C SVSEED IS READ IN FROM FILE 3. C C C-----USER-SPECIFIED PARAMETERS FOR THE EIGENVECTOR PROGRAMS------------ C C C NGOOD = NUMBER OF EIGENVALUES READ INTO THE GOODEV ARRAY C READ FROM FILE 3. C C N = SIZE OF THE USER-SUPPLIED MATRIX. C C MEV = SIZE OF THE T-MATRIX THAT WAS USED TO COMPUTE C THE EIGENVALUES WHOSE EIGENVECTORS ARE REQUESTED. C MEV IS READ IN FROM FILE 3. C C KMAX = SIZE OF THE T-MATRIX PROVIDED ON FILE 2. C C MDIMTV = MAXIMUM CUMULATIVE SIZE OF THE TVEC ARRAY ALLOWED C FOR ALL OF THE T-EIGENVECTORS REQUIRED. MDIMTV C MUST NOT EXCEED THE USER-SPECIFIED DIMENSION OF C THE TVEC ARRAY. PROGRAM CAN BE RUN WITH THE FLAG C MBOUND = 1 TO DETERMINE AN EDUCATED GUESS ON AN C APPROPRIATE DIMENSION FOR THE TVEC ARRAY. C C MDIMRV = MAXIMUM CUMULATIVE SIZE OF THE RITVEC ARRAY ALLOWED C FOR ALL OF THE RITZ VECTORS TO BE COMPUTED. MDIMRV C MUST NOT EXCEED THE USER-SPECIFIED DIMENSION OF C THE RITVEC ARRAY. MUST BE SELECTED SO THAT C THERE IS ENOUGH ROOM FOR A RITZ VECTOR FOR EVERY C GOODEV(J) READ INTO PROGRAM. (>= NGOOD*N) C C C-----ARRAYS REQUIRED BY THE EIGENVECTOR PROGRAMS--------------------- C C C ALL 4 CASES C C ALPHA(J) = REAL*8 ARRAY. ITS DIMENSION MUST BE AT LEAST C KMAXN, THE LARGEST SIZE T-MATRIX CONSIDERED BY C THE PROGRAM. NOTE THAT KMAXN IS THE LARGER OF C THE SIZE OF THE ALPHA, BETA HISTORY PROVIDED C ON FILE 2 (IF ANY ) AND THE SIZE WHICH THE PROGRAM C SPECIFIES INTERNALLY, THIS LATTER IS ALWAYS C < = 11*MEV / 8 + 12, WHERE MEV IS THE SIZE C T-MATRIX THAT WAS USED IN THE CORRESPONDING EIGENVALUE C COMPUTATIONS. ALPHA CONTAINS THE DIAGONAL ENTRIES C OF THE LANCZOS T-MATRICES. ALPHA IS NOT DESTROYED C IN THE COMPUTATIONS. C C BETA(J) = REAL*8 ARRAY. ITS DIMENSION MUST BE AT LEAST 1 C MORE THAN THAT OF ALPHA. DIMENSION COMMENTS ABOVE C ABOUT ALPHA APPLY ALSO TO THE BETA ARRAY. BETA C CONTAINS THE SUBDIAGONAL ENTRIES OF THE T-MATRICES. C BETA IS NOT DESTROYED IN THE COMPUTATIONS. C C RITVEC(J) = REAL*8 ARRAY IN REAL SYMMETRIC AND INVERSE OF C REAL SYMMETRIC CASES. COMPLEX*16 IN CASE (2), C HERMITIAN MATRICES. IN EACH CASE ITS DIMENSION >= C NGOOD*N WHERE N IS THE ORDER OF THE USER-SUPPLIED C MATRIX AND NGOOD IS THE NUMBER OF EIGENVALUES WHOSE C EIGENVECTORS ARE TO BE COMPUTED. IT CONTAINS THE C COMPUTED APPROXIMATE EIGENVECTORS OF A. THESE C COMPUTED RITZ VECTORS ARE STORED ON FILE 12. C C TVEC(J) = REAL*8 ARRAY. ITS DIMENSION MUST BE AT LEAST C MTOL = |MA(1)| + |MA(2)| + ... + |MA(NGOOD)| C WHERE NGOOD IS THE NUMBER OF EIGENVALUES BEING C CONSIDERED AND |MA(J)| IS THE SIZE OF THE C T-MATRIX BEING USED FOR THE RITZ VECTOR C COMPUTATION FOR GOODEV(J). THESE SIZES C ARE COMPUTED BY THE PROGRAM. AN ESTIMATE OF C MTOL CAN BE OBTAINED BY SETTING MBOUND = 1, C RUNNING THE PROGRAM, AND MULTIPLYING THE RESULTING C TOTAL OF THE T-SIZES SPECIFIED BY 5/4. THE ARRAY C TVEC CONTAINS THE COMPUTED T-EIGENVECTORS. IF THE C FLAG SVTVEC = 1 OR THE FLAG TVSTOP = 1, THEN C THESE VECTORS ARE SAVED ON FILE 11. C C V1(J) = REAL*8 ARRAY IN REAL SYMMETRIC AND INVERSE OF C REAL SYMMETRIC CASES. COMPLEX*16 IN CASE (2), C HERMITIAN MATRICES. IN THE REAL CASES ITS C DIMENSION MUST BE THE MAXIMUM OF KMAX AND N. C IN THE HERMITIAN CASE ITS DIMENSION MUST BE C THE MAXIMUM OF KMAX/2 AND N WHERE KMAX IS THE C LARGEST SIZE T-MATRIX THAT IS TO BE CONSIDERED C IN THE T-EIGENVECTOR COMPUTATIONS. V1 IS USED C IN THE SUBROUTINE INVERM AND IN THE REGENERATION C OF THE LANCZOS VECTORS. C C V2(J) = REAL*8 ARRAY IN THE REAL SYMMETRIC AND INVERSE C OF REAL SYMMETRIC CASES. COMPLEX*16 IN CASE (2), C HERMITIAN MATRICES. IN CASES (1),(3) AND (4), ITS C DIMENSION MUST BE > = MAX(KMAX,N); IN CASE (2) C > = MAX(KMAX/2,N). IT IS USED IN THE REGENERATION C OF THE LANCZOS VECTORS AND IN SUBROUTINE INVERM. C C GOODEV(J), = REAL*8 ARRAYS EACH OF DIMENSION AT LEAST NGOOD. C EVNEW(J) CONTAIN THE EIGENVALUES FOR WHICH EIGENVECTORS C ARE REQUESTED. EIGENVALUES IN GOODEV ARE READ C IN FROM FILE 3. IN CASE (3) GOODEV IS REPLACED C BY GOODA AND GOODBI ARRAYS, SEE BELOW. C C AMINGP(J), = REAL*4 ARRAYS OF DIMENSION AT LEAST NGOOD. C TMINGP(J) CONTAIN, RESPECTIVELY, THE MINIMAL GAPS FOR C CORRESPONDING EIGENVALUES IN GOODEV ARRAY IN C A-MATRIX AND IN T-MATRIX. C C TERR(J), ERR(J), = REAL*4 ARRAYS (EXCEPT TLAST WHICH IS C ERRDGP(J), TLAST(J) REAL*8). EACH MUST BE OF DIMENSION C RNORM(J), TBETA(J) AT LEAST NGOOD. USED TO STORE QUANTITIES C GENERATED DURING THE COMPUTATIONS FOR C LATER PRINTOUT. C C G(J) = REAL*4 ARRAY WHOSE DIMENSION MUST BE AT LEAST C MAX(KMAX,N). USED IN SUBROUTINE GENRAN TO HOLD C RANDOM NUMBERS NEEDED FOR THE LANCZOS VECTOR C REGENERATION AND FOR THE INVERSE ITERATION C COMPUTATIONS IN THE SUBROUTINE INVERM. C C MP(J) = INTEGER*4 ARRAY WHOSE DIMENSION IS AT LEAST NGOOD. C INITIALLY CONTAINS THE MULTIPLICITY OF THE EIGENVALUE C GOODEV(J) AS AN EIGENVALUE OF THE T-MATRIX T(1,MEV). C USED TO FLAG EIGENVALUES FOR WHICH NO T-EIGENVECTOR C OR NO RITZ VECTOR IS TO BE COMPUTED. C C MA(J) = INTEGER*4 ARRAYS EACH OF WHOSE DIMENSIONS C IS AT LEAST NGOOD. USED IN DETERMINING C AN APPROPRIATE T-MATRIX FOR EACH EIGENVALUE C IN GOODEV ARRAY. C C MINT(J),MFIN(J) = INTEGER*4 ARRAYS WHOSE DIMENSIONS MUST BE AT C LEAST NGOOD. USED TO POINT TO THE BEGINNINGS C AND THE ENDS OF THE COMPUTED EIGENVECTOR C OF THE T-MATRIX, T(1,|MA(J)|). C C IDELTA(J) = INTEGER*4 ARRAY WHOSE DIMENSION MUST BE AT C LEAST NGOOD. CONTAINS INCREMENTS USED IN LOOPS C ON APPROPRIATE SIZE T-MATRIX FOR THE T-EIGENVECTOR C COMPUTATIONS. C C C CASE (2) ONLY, HERMITIAN MATRICES: C C GR(J),GC(J) = REAL*8 ARRAYS USED ONLY IN CASE (2), C HERMITIAN MATRICES. EACH MUST BE AT C LEAST MAX(N,KMAX). USED TO HOLD C STARTING VECTORS FOR LANCZS C COMPUTATIONS AND FOR INVERM SUBROUTINES. C C CASES (3) AND (4) ONLY, FACTORED INVERSES OF REAL SYMMETRIC C MATRICES AND GENERALIZED EIGENVALUE PROBLEMS: C C VS(J) = REAL*8 ARRAY WHOSE DIMENSION MUST BE AT LEAST N. C USED IN REGENERATION OF THE LANCZOS VECTORS. C C IPR(J), IPT(J) = INTEGER*4 ARRAYS. EACH MUST BE OF DIMENSION C AT LEAST N, THE ORDER OF A. USED TO STORE C THE REORDERING OF THE GIVEN MATRIX. C C CASE (3) ONLY, INVERSES OF REAL SYMMETRIC MATRICES: C C GOODA(J), GOODBI(J) = REAL*8 ARRAYS. EACH MUST BE OF DIMENSION C AT LEAST NGOOD, THE NUMBER OF EIGENVALUES C BEING CONSIDERED. GOODA CONTAINS THE C EIGENVALUES OF A AND GOODBI CONTAINS THE C EIGENVALUES OF B(INVERSE). THE PROGRAM C WORKS DIRECTLY WITH THE GOODBI ARRAY. C C C-----SUBROUTINES INCLUDED FOR THE EIGENVECTOR COMPUTATIONS------------- C C C STURMI = FOR EACH GIVEN EIGENVALUE GOODEV(J) DETERMINES C THE SMALLEST SIZE T-MATRIX FOR WHICH GOODEV(J) IS C AN EIGENVALUE (TO WITHIN A GIVEN TOLERANCE) AND IF C POSSIBLE THE SMALLEST SIZE T-MATRIX FOR WHICH C IT IS A DOUBLE EIGENVALUE (TO WITHIN THE SAME C TOLERANCE). THE SIZE T-MATRIX USED IN THE C EIGENVECTOR COMPUTATIONS IS THEN DETERMINED BY C STARTING WITH AN INITIAL GUESS BASED UPON THE C INFORMATION FROM STURMI, AND LOOPING ON THE SIZE C OF THE T-EIGENVECTOR COMPUTATIONS. C C LBISEC = RECOMPUTES THE VALUE OF THE GIVEN EIGENVALUE AT THE C T-SIZE SPECIFIED FOR THE T-EIGENVECTOR COMPUTATION. C LBISEC IS A SIMPLIFICATION OF THE BISEC SUBROUTINE C USED IN THE LANCZOS EIGENVALUE COMPUTATIONS. C C INVERM = FOR THE T-SIZES CONSIDERED BY THE PROGRAM COMPUTES C THE CORRESPONDING EIGENVECTORS OF THESE T-MATRICES C CORRESPONDING TO THE USER-SUPPLIED EIGENVALUES IN C THE GOODEV ARRAY. C C THE LANCZS, TNORM , AND CINPRD (CASE (2) ONLY) SUBROUTINES C ARE USED HERE AS WELL AS IN THE EIGENVALUE COMPUTATIONS. C C IN CASES (3) AND (4) ONLY AND THEN ONLY IF THE ORIGINAL MATRIX C (MATRICES) HAS (HAVE) BEEN PERMUTED: C C LPERM = (USED IN CASE (3) AND (4) ONLY). GIVEN A B-MATRIX AND C A PERMUTATION P DEFINED IN THE VECTORS IPR AND IPT, C AND A VECTOR X COMPUTE EITHER (P-TRANSPOSE)*X OR PX. C C-----------------------------------------------------------------------