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