SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI) DOUBLE PRECISION AR,AI,BR,BI,CR,CI C C COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) C SUBROUTINE CSROOT(XR,XI,YR,YI) DOUBLE PRECISION XR,XI,YR,YI C C (YR,YI) = COMPLEX DSQRT(XR,XI) C BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI) C DOUBLE PRECISION FUNCTION EPSLON (X) DOUBLE PRECISION X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C DOUBLE PRECISION FUNCTION PYTHAG(A,B) DOUBLE PRECISION A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C SUBROUTINE BAKVEC(NM,N,T,E,M,Z,IERR) C INTEGER I,J,M,N,NM,IERR DOUBLE PRECISION T(NM,3),E(N),Z(NM,M) C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A NONSYMMETRIC C TRIDIAGONAL MATRIX BY BACK TRANSFORMING THOSE OF THE C CORRESPONDING SYMMETRIC MATRIX DETERMINED BY FIGI. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C T CONTAINS THE NONSYMMETRIC MATRIX. ITS SUBDIAGONAL IS C STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, C ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN, C AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF C THE THIRD COLUMN. T(1,1) AND T(N,3) ARE ARBITRARY. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C T IS UNALTERED. C C E IS DESTROYED. C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 2*N+I IF E(I) IS ZERO WITH T(I,1) OR T(I-1,3) NON-ZERO. C IN THIS CASE, THE SYMMETRIC MATRIX IS NOT SIMILAR C TO THE ORIGINAL MATRIX, AND THE EIGENVECTORS C CANNOT BE FOUND BY THIS PROGRAM. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC DOUBLE PRECISION A(NM,N),SCALE(N) DOUBLE PRECISION C,F,G,R,S,B2,RADIX LOGICAL NOCONV C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES C EIGENVALUES WHENEVER POSSIBLE. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE INPUT MATRIX TO BE BALANCED. C C ON OUTPUT C C A CONTAINS THE BALANCED MATRIX. C C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) C IS EQUAL TO ZERO IF C (1) I IS GREATER THAN J AND C (2) J=1,...,LOW-1 OR I=IGH+1,...,N. C C SCALE CONTAINS INFORMATION DETERMINING THE C PERMUTATIONS AND SCALING FACTORS USED. C C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN C SCALE(J) = P(J), FOR J = 1,...,LOW-1 C = D(J,J), J = LOW,...,IGH C = P(J) J = IGH+1,...,N. C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, C THEN 1 TO LOW-1. C C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. C C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS C K,L HAVE BEEN REVERSED.) C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z) C INTEGER I,J,K,M,N,II,NM,IGH,LOW DOUBLE PRECISION SCALE(N),Z(NM,M) DOUBLE PRECISION S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C BALANCED MATRIX DETERMINED BY BALANC. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY BALANC. C C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS C AND SCALING FACTORS USED BY BALANC. C C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BANDR(NM,N,MB,A,D,E,E2,MATZ,Z) C INTEGER J,K,L,N,R,I1,I2,J1,J2,KR,MB,MR,M1,NM,N2,R1,UGL,MAXL,MAXR DOUBLE PRECISION A(NM,MB),D(N),E(N),E2(N),Z(NM,N) DOUBLE PRECISION G,U,B1,B2,C2,F1,F2,S2,DMIN,DMINRT LOGICAL MATZ C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, C NUM. MATH. 12, 231-241(1968) BY SCHWARZ. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX C TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY C ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE C NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL C DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE C LOWER TRIANGLE OF THE MATRIX. C C A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT C MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL C IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, C ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY C ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN. C CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. C C MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS C TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE. C C ON OUTPUT C C A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH C CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN C THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z C IS NOT REFERENCED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BANDV(NM,N,MBW,A,E21,M,W,Z,IERR,NV,RV,RV6) C INTEGER I,J,K,M,N,R,II,IJ,JJ,KJ,MB,M1,NM,NV,IJ1,ITS,KJ1,MBW,M21, X IERR,MAXJ,MAXK,GROUP DOUBLE PRECISION A(NM,MBW),W(M),Z(NM,M),RV(NV),RV6(N) DOUBLE PRECISION U,V,UK,XU,X0,X1,E21,EPS2,EPS3,EPS4,NORM,ORDER, X EPSLON,PYTHAG C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC C BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE C ITERATION. THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS C OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND C COEFFICIENT MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE C BAND MATRIX. IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) C BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT C DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO C SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE C MATRIX. IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS C OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT C SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT C DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS C CASE, MBW=2*MB-1. C C A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT C MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL C IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, C ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY C ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. C IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR C EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS C N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH C ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF C COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2 C POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY, C AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB C POSITIONS OF THE LAST COLUMN. C CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. C C E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS C 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR C 2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER. C IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR C EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT C MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT. C C M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF C SYSTEMS OF LINEAR EQUATIONS. C C W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. C IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR C EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY C MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M. C C Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF C THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS. C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. C C ON OUTPUT C C A AND W ARE UNALTERED. C C Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. IF THE C SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, C Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M). C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH C EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH C SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR. C C RV AND RV6 ARE TEMPORARY STORAGE ARRAYS. NOTE THAT RV IS C OF DIMENSION AT LEAST N*(2*MB-1). IF THE SUBROUTINE C IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE C DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON C RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5) C INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM DOUBLE PRECISION D(N),E(N),E2(N),W(MM),RV4(N),RV5(N) DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON INTEGER IND(MM) C C THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE C IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL C SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL, C USING BISECTION. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED C EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, C IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, C NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE C PRECISION AND THE 1-NORM OF THE SUBMATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY. C C LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. C IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. C C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF C EIGENVALUES IN THE INTERVAL. WARNING. IF MORE THAN C MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, C AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND. C C ON OUTPUT C C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS C (LAST) DEFAULT VALUE. C C D AND E ARE UNALTERED. C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS ALSO SET TO ZERO. C C M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). C C W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER. C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 3*N+1 IF M EXCEEDS MM. C C RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. C C THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM C APPEARS IN BISECT IN-LINE. C C NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN C BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE BQR(NM,N,MB,A,T,R,IERR,NV,RV) C INTEGER I,J,K,L,M,N,II,IK,JK,JM,KJ,KK,KM,LL,MB,MK,MN,MZ, X M1,M2,M3,M4,NI,NM,NV,ITS,KJ1,M21,M31,IERR,IMULT DOUBLE PRECISION A(NM,MB),RV(NV) DOUBLE PRECISION F,G,Q,R,S,T,TST1,TST2,SCALE,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR, C NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY) C MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE C QR ALGORITHM WITH SHIFTS OF ORIGIN. CONSECUTIVE CALLS C CAN BE MADE TO FIND FURTHER EIGENVALUES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE C NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL C DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE C LOWER TRIANGLE OF THE MATRIX. C C A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT C MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL C IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, C ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY C ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN. C CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. C ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS C CALL SHOULD BE PASSED. C C T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL C OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED C IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST C TO T. ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE C PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE C IS SOUGHT. C C R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS C OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL. C IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF C THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE. C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. C C ON OUTPUT C C A CONTAINS THE TRANSFORMED BAND MATRIX. THE MATRIX A+TI C DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE C INPUT A+TI TO WITHIN ROUNDING ERRORS. ITS LAST ROW AND C COLUMN ARE NULL (IF IERR IS ZERO). C C T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO). C C R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE C LAST COLUMN OF THE INPUT MATRIX A. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C N IF THE EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST C (2*MB**2+4*MB-3). THE FIRST (3*MB-2) LOCATIONS CORRESPOND C TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND C TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS C CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U. C C NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT C MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CBABK2(NM,N,LOW,IGH,SCALE,M,ZR,ZI) C INTEGER I,J,K,M,N,II,NM,IGH,LOW DOUBLE PRECISION SCALE(N),ZR(NM,M),ZI(NM,M) DOUBLE PRECISION S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE C CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C BALANCED MATRIX DETERMINED BY CBAL. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY CBAL. C C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS C AND SCALING FACTORS USED BY CBAL. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS TO BE C BACK TRANSFORMED IN THEIR FIRST M COLUMNS. C C ON OUTPUT C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CBAL(NM,N,AR,AI,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC DOUBLE PRECISION AR(NM,N),AI(NM,N),SCALE(N) DOUBLE PRECISION C,F,G,R,S,B2,RADIX LOGICAL NOCONV C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE C CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES C EIGENVALUES WHENEVER POSSIBLE. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED. C C ON OUTPUT C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE BALANCED MATRIX. C C LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J) C ARE EQUAL TO ZERO IF C (1) I IS GREATER THAN J AND C (2) J=1,...,LOW-1 OR I=IGH+1,...,N. C C SCALE CONTAINS INFORMATION DETERMINING THE C PERMUTATIONS AND SCALING FACTORS USED. C C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN C SCALE(J) = P(J), FOR J = 1,...,LOW-1 C = D(J,J) J = LOW,...,IGH C = P(J) J = IGH+1,...,N. C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, C THEN 1 TO LOW-1. C C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. C C THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN C CBAL IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS C K,L HAVE BEEN REVERSED.) C C ARITHMETIC IS REAL THROUGHOUT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CG(NM,N,AR,AI,WR,WI,MATZ,ZR,ZI,FV1,FV2,FV3,IERR) C INTEGER N,NM,IS1,IS2,IERR,MATZ DOUBLE PRECISION AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N), X FV1(N),FV2(N),FV3(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A COMPLEX GENERAL MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A=(AR,AI). C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX GENERAL MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR COMQR C AND COMQR2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1, FV2, AND FV3 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CH(NM,N,AR,AI,W,MATZ,ZR,ZI,FV1,FV2,FM1,IERR) C INTEGER I,J,N,NM,IERR,MATZ DOUBLE PRECISION AR(NM,N),AI(NM,N),W(N),ZR(NM,N),ZI(NM,N), X FV1(N),FV2(N),FM1(2,N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A COMPLEX HERMITIAN MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A=(AR,AI). C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX HERMITIAN MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1, FV2, AND FM1 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CINVIT(NM,N,AR,AI,WR,WI,SELECT,MM,M,ZR,ZI, X IERR,RM1,RM2,RV1,RV2) C INTEGER I,J,K,M,N,S,II,MM,MP,NM,UK,IP1,ITS,KM1,IERR DOUBLE PRECISION AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,MM), X ZI(NM,MM),RM1(N,N),RM2(N,N),RV1(N),RV2(N) DOUBLE PRECISION X,Y,EPS3,NORM,NORMV,EPSLON,GROWTO,ILAMBD,PYTHAG, X RLAMBD,UKROOT LOGICAL SELECT(N) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT C BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER C HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, C USING INVERSE ITERATION. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE HESSENBERG MATRIX. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, C OF THE EIGENVALUES OF THE MATRIX. THE EIGENVALUES MUST BE C STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE COMLR, C WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. C C SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE C EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS C SPECIFIED BY SETTING SELECT(J) TO .TRUE.. C C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF C EIGENVECTORS TO BE FOUND. C C ON OUTPUT C C AR, AI, WI, AND SELECT ARE UNALTERED. C C WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED C SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. C C M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, C OF THE EIGENVECTORS. THE EIGENVECTORS ARE NORMALIZED C SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. C ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -(2*N+1) IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED, C -K IF THE ITERATION CORRESPONDING TO THE K-TH C VALUE FAILS, C -(N+K) IF BOTH ERROR SITUATIONS OCCUR. C C RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. C C THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE. C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMBAK(NM,LOW,IGH,AR,AI,INT,M,ZR,ZI) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 DOUBLE PRECISION AR(NM,IGH),AI(NM,IGH),ZR(NM,M),ZI(NM,M) DOUBLE PRECISION XR,XI INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C UPPER HESSENBERG MATRIX DETERMINED BY COMHES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE C REDUCTION BY COMHES IN THEIR LOWER TRIANGLES C BELOW THE SUBDIAGONAL. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION BY COMHES. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS TO BE C BACK TRANSFORMED IN THEIR FIRST M COLUMNS. C C ON OUTPUT C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMHES(NM,N,LOW,IGH,AR,AI,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 DOUBLE PRECISION AR(NM,N),AI(NM,N) DOUBLE PRECISION XR,XI,YR,YI INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY C STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. C C ON OUTPUT C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE HESSENBERG MATRIX. THE C MULTIPLIERS WHICH WERE USED IN THE REDUCTION C ARE STORED IN THE REMAINING TRIANGLES UNDER THE C HESSENBERG MATRIX. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C CALLS CDIV FOR COMPLEX DIVISION. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMLR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR) C INTEGER I,J,L,M,N,EN,LL,MM,NM,IGH,IM1,ITN,ITS,LOW,MP1,ENM1,IERR DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N) DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,TST1,TST2 C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, C NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX C UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE C MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, C IF PERFORMED. C C ON OUTPUT C C THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN C DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE C CALLING COMLR IF SUBSEQUENT CALCULATION OF C EIGENVECTORS IS TO BE PERFORMED. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS CSROOT FOR COMPLEX SQUARE ROOT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMLR2(NM,N,LOW,IGH,INT,HR,HI,WR,WI,ZR,ZI,IERR) C INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NM,NN,IGH,IM1,IP1, X ITN,ITS,LOW,MP1,ENM1,IEND,IERR DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N) DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2 INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR2, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A COMPLEX UPPER HESSENBERG MATRIX BY THE MODIFIED LR C METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX C CAN ALSO BE FOUND IF COMHES HAS BEEN USED TO REDUCE C THIS GENERAL MATRIX TO HESSENBERG FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS INTERCHANGED C IN THE REDUCTION BY COMHES, IF PERFORMED. ONLY ELEMENTS C LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS OF THE HESSEN- C BERG MATRIX ARE DESIRED, SET INT(J)=J FOR THESE ELEMENTS. C C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE C MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, C IF PERFORMED. IF THE EIGENVECTORS OF THE HESSENBERG C MATRIX ARE DESIRED, THESE ELEMENTS MUST BE SET TO ZERO. C C ON OUTPUT C C THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN C DESTROYED, BUT THE LOCATION HR(1,1) CONTAINS THE NORM C OF THE TRIANGULARIZED MATRIX. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS C ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF C THE EIGENVECTORS HAS BEEN FOUND. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS CSROOT FOR COMPLEX SQUARE ROOT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR) C INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N) DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2, X PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE C ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN C AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX C UPPER HESSENBERG MATRIX BY THE QR METHOD. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN C INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN C THE REDUCTION BY CORTH, IF PERFORMED. C C ON OUTPUT C C THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN C DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE C CALLING COMQR IF SUBSEQUENT CALCULATION OF C EIGENVECTORS IS TO BE PERFORMED. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS CSROOT FOR COMPLEX SQUARE ROOT. C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE COMQR2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR) C INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1, X ITN,ITS,LOW,LP1,ENM1,IEND,IERR DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N), X ORTR(IGH),ORTI(IGH) DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2, X PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE C ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS C AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR C METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX C CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE C THIS GENERAL MATRIX TO HESSENBERG FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- C FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS C OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND C ORTI(J) TO 0.0D0 FOR THESE ELEMENTS. C C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER C INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE C REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF C THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE C ARBITRARY. C C ON OUTPUT C C ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI C HAVE BEEN DESTROYED. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS C ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF C THE EIGENVECTORS HAS BEEN FOUND. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS CSROOT FOR COMPLEX SQUARE ROOT. C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CORTB(NM,LOW,IGH,AR,AI,ORTR,ORTI,M,ZR,ZI) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 DOUBLE PRECISION AR(NM,IGH),AI(NM,IGH),ORTR(IGH),ORTI(IGH), X ZR(NM,M),ZI(NM,M) DOUBLE PRECISION H,GI,GR C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968) C BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C UPPER HESSENBERG MATRIX DETERMINED BY CORTH. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY C TRANSFORMATIONS USED IN THE REDUCTION BY CORTH C IN THEIR STRICT LOWER TRIANGLES. C C ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE C TRANSFORMATIONS USED IN THE REDUCTION BY CORTH. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED. C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVECTORS TO BE C BACK TRANSFORMED IN THEIR FIRST M COLUMNS. C C ON OUTPUT C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C ORTR AND ORTI HAVE BEEN ALTERED. C C NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI) C INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW DOUBLE PRECISION AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH) DOUBLE PRECISION F,G,H,FI,FR,SCALE,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968) C BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY C UNITARY SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. C C ON OUTPUT C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE HESSENBERG MATRIX. INFORMATION C ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION C IS STORED IN THE REMAINING TRIANGLES UNDER THE C HESSENBERG MATRIX. C C ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE C TRANSFORMATIONS. ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ELMBAK(NM,LOW,IGH,A,INT,M,Z) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 DOUBLE PRECISION A(NM,IGH),Z(NM,M) DOUBLE PRECISION X INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMBAK, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C UPPER HESSENBERG MATRIX DETERMINED BY ELMHES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. C C A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE C REDUCTION BY ELMHES IN ITS LOWER TRIANGLE C BELOW THE SUBDIAGONAL. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION BY ELMHES. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 DOUBLE PRECISION A(NM,N) DOUBLE PRECISION X,Y INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY C STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C A CONTAINS THE INPUT MATRIX. C C ON OUTPUT C C A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS C WHICH WERE USED IN THE REDUCTION ARE STORED IN THE C REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 DOUBLE PRECISION A(NM,IGH),Z(NM,N) INTEGER INT(IGH) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY C SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A C REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY ELMHES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE C REDUCTION BY ELMHES IN ITS LOWER TRIANGLE C BELOW THE SUBDIAGONAL. C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION BY ELMHES. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY ELMHES. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C C SUBROUTINE FIGI(NM,N,T,D,E,E2,IERR) C INTEGER I,N,NM,IERR DOUBLE PRECISION T(NM,3),D(N),E(N),E2(N) C C GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS C OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL C NON-NEGATIVE, THIS SUBROUTINE REDUCES IT TO A SYMMETRIC C TRIDIAGONAL MATRIX WITH THE SAME EIGENVALUES. IF, FURTHER, C A ZERO PRODUCT ONLY OCCURS WHEN BOTH FACTORS ARE ZERO, C THE REDUCED MATRIX IS SIMILAR TO THE ORIGINAL MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C T CONTAINS THE INPUT MATRIX. ITS SUBDIAGONAL IS C STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, C ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN, C AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF C THE THIRD COLUMN. T(1,1) AND T(N,3) ARE ARBITRARY. C C ON OUTPUT C C T IS UNALTERED. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS NOT SET. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C N+I IF T(I,1)*T(I-1,3) IS NEGATIVE, C -(3*N+I) IF T(I,1)*T(I-1,3) IS ZERO WITH ONE FACTOR C NON-ZERO. IN THIS CASE, THE EIGENVECTORS OF C THE SYMMETRIC MATRIX ARE NOT SIMPLY RELATED C TO THOSE OF T AND SHOULD NOT BE SOUGHT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE FIGI2(NM,N,T,D,E,Z,IERR) C INTEGER I,J,N,NM,IERR DOUBLE PRECISION T(NM,3),D(N),E(N),Z(NM,N) DOUBLE PRECISION H C C GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS C OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL C NON-NEGATIVE, AND ZERO ONLY WHEN BOTH FACTORS ARE ZERO, THIS C SUBROUTINE REDUCES IT TO A SYMMETRIC TRIDIAGONAL MATRIX C USING AND ACCUMULATING DIAGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C T CONTAINS THE INPUT MATRIX. ITS SUBDIAGONAL IS C STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, C ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN, C AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF C THE THIRD COLUMN. T(1,1) AND T(N,3) ARE ARBITRARY. C C ON OUTPUT C C T IS UNALTERED. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS NOT SET. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN C THE REDUCTION. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C N+I IF T(I,1)*T(I-1,3) IS NEGATIVE, C 2*N+I IF T(I,1)*T(I-1,3) IS ZERO WITH C ONE FACTOR NON-ZERO. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) C INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR DOUBLE PRECISION H(NM,N),WR(N),WI(N) DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2 LOGICAL NOTLAS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL C UPPER HESSENBERG MATRIX BY THE QR METHOD. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. C C ON OUTPUT C C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR) C INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN, X IGH,ITN,ITS,LOW,MP2,ENM2,IERR DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N) DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2 LOGICAL NOTLAS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE C EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND C IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE C BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM C AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C H CONTAINS THE UPPER HESSENBERG MATRIX. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN C AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE C REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS C OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE C IDENTITY MATRIX. C C ON OUTPUT C C H HAS BEEN DESTROYED. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. C IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z C CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX C WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH C COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS C EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN C ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C CALLS CDIV FOR COMPLEX DIVISION. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HTRIB3(NM,N,A,TAU,M,ZR,ZI) C INTEGER I,J,K,L,M,N,NM DOUBLE PRECISION A(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) DOUBLE PRECISION H,S,SI C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRBAK3, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY HTRID3. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS C USED IN THE REDUCTION BY HTRID3. C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI) C INTEGER I,J,K,L,M,N,NM DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) DOUBLE PRECISION H,S,SI C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY HTRIDI. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- C FORMATIONS USED IN THE REDUCTION BY HTRIDI IN THEIR C FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR. C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS C IN THEIR FIRST M COLUMNS. C C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HTRID3(NM,N,A,D,E,E2,TAU) C INTEGER I,J,K,L,N,II,NM,JM1,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N),TAU(2,N) DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRED3, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX, STORED AS C A SINGLE SQUARE ARRAY, TO A REAL SYMMETRIC TRIDIAGONAL MATRIX C USING UNITARY SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE LOWER TRIANGLE OF THE COMPLEX HERMITIAN INPUT C MATRIX. THE REAL PARTS OF THE MATRIX ELEMENTS ARE STORED C IN THE FULL LOWER TRIANGLE OF A, AND THE IMAGINARY PARTS C ARE STORED IN THE TRANSPOSED POSITIONS OF THE STRICT UPPER C TRIANGLE OF A. NO STORAGE IS REQUIRED FOR THE ZERO C IMAGINARY PARTS OF THE DIAGONAL ELEMENTS. C C ON OUTPUT C C A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS C USED IN THE REDUCTION. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N) DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF C THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) C BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX C TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING C UNITARY SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. C ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT C C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- C FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER C TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE C DIAGONAL OF AR ARE UNALTERED. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE IMTQL1(N,D,E,IERR) C INTEGER I,J,L,M,N,II,MML,IERR DOUBLE PRECISION D(N),E(N) DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E HAS BEEN DESTROYED. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1. C C E HAS BEEN DESTROYED. C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1) C INTEGER I,J,K,L,M,N,II,MML,TAG,IERR DOUBLE PRECISION D(N),E(N),E2(N),W(N),RV1(N) DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG INTEGER IND(N) C C THIS SUBROUTINE IS A VARIANT OF IMTQL1 WHICH IS A TRANSLATION OF C ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND C WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL C MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM C THEIR CORRESPONDING SUBMATRIX INDICES. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY. C C ON OUTPUT C C D AND E ARE UNALTERED. C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS ALSO SET TO ZERO. C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE C CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES C BELONGING TO THE FIRST SUBMATRIX FROM THE TOP, C 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C RV1 IS A TEMPORARY STORAGE ARRAY. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM1,RV1,RV2) C INTEGER I,J,K,L,M,N,S,II,IP,MM,MP,NM,NS,N1,UK,IP1,ITS,KM1,IERR DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,MM),RM1(N,N), X RV1(N),RV2(N) DOUBLE PRECISION T,W,X,Y,EPS3,NORM,NORMV,EPSLON,GROWTO,ILAMBD, X PYTHAG,RLAMBD,UKROOT LOGICAL SELECT(N) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE INVIT C BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL UPPER C HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, C USING INVERSE ITERATION. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE HESSENBERG MATRIX. C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, C OF THE EIGENVALUES OF THE MATRIX. THE EIGENVALUES MUST BE C STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE HQR, C WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. C C SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE C EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS C SPECIFIED BY SETTING SELECT(J) TO .TRUE.. C C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF C COLUMNS REQUIRED TO STORE THE EIGENVECTORS TO BE FOUND. C NOTE THAT TWO COLUMNS ARE REQUIRED TO STORE THE C EIGENVECTOR CORRESPONDING TO A COMPLEX EIGENVALUE. C C ON OUTPUT C C A AND WI ARE UNALTERED. C C WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED C SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. C C SELECT MAY HAVE BEEN ALTERED. IF THE ELEMENTS CORRESPONDING C TO A PAIR OF CONJUGATE COMPLEX EIGENVALUES WERE EACH C INITIALLY SET TO .TRUE., THE PROGRAM RESETS THE SECOND OF C THE TWO ELEMENTS TO .FALSE.. C C M IS THE NUMBER OF COLUMNS ACTUALLY USED TO STORE C THE EIGENVECTORS. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. C IF THE NEXT SELECTED EIGENVALUE IS REAL, THE NEXT COLUMN C OF Z CONTAINS ITS EIGENVECTOR. IF THE EIGENVALUE IS C COMPLEX, THE NEXT TWO COLUMNS OF Z CONTAIN THE REAL AND C IMAGINARY PARTS OF ITS EIGENVECTOR. THE EIGENVECTORS ARE C NORMALIZED SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. C ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -(2*N+1) IF MORE THAN MM COLUMNS OF Z ARE NECESSARY C TO STORE THE EIGENVECTORS CORRESPONDING TO C THE SPECIFIED EIGENVALUES. C -K IF THE ITERATION CORRESPONDING TO THE K-TH C VALUE FAILS, C -(N+K) IF BOTH ERROR SITUATIONS OCCUR. C C RM1, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. NOTE THAT RM1 C IS SQUARE OF DIMENSION N BY N AND, AUGMENTED BY TWO COLUMNS C OF Z, IS THE TRANSPOSE OF THE CORRESPONDING ALGOL B ARRAY. C C THE ALGOL PROCEDURE GUESSVEC APPEARS IN INVIT IN LINE. C C CALLS CDIV FOR COMPLEX DIVISION. C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE MINFIT(NM,M,N,A,W,IP,B,IERR,RV1) C INTEGER I,J,K,L,M,N,II,IP,I1,KK,K1,LL,L1,M1,NM,ITS,IERR DOUBLE PRECISION A(NM,N),W(N),B(NM,IP),RV1(N) DOUBLE PRECISION C,F,G,H,S,X,Y,Z,TST1,TST2,SCALE,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT, C NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). C C THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR C T C SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV OF A REAL C T C M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U. HOUSEHOLDER C BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. NOTE THAT NM MUST BE AT LEAST C AS LARGE AS THE MAXIMUM OF M AND N. C C M IS THE NUMBER OF ROWS OF A AND B. C C N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V. C C A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM. C C IP IS THE NUMBER OF COLUMNS OF B. IP CAN BE ZERO. C C B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM C IF IP IS NOT ZERO. OTHERWISE B IS NOT REFERENCED. C C ON OUTPUT C C A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE C DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS. IF AN C ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO C INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. C C W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE C DIAGONAL ELEMENTS OF S). THEY ARE UNORDERED. IF AN C ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT C FOR INDICES IERR+1,IERR+2,...,N. C C T C B HAS BEEN OVERWRITTEN BY U B. IF AN ERROR EXIT IS MADE, C T C THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT C SINGULAR VALUES SHOULD BE CORRECT. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C K IF THE K-TH SINGULAR VALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C RV1 IS A TEMPORARY STORAGE ARRAY. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ORTBAK(NM,LOW,IGH,A,ORT,M,Z) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 DOUBLE PRECISION A(NM,IGH),ORT(IGH),Z(NM,M) DOUBLE PRECISION G C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTBAK, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C UPPER HESSENBERG MATRIX DETERMINED BY ORTHES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES C IN ITS STRICT LOWER TRIANGLE. C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. C C ORT HAS BEEN ALTERED. C C NOTE THAT ORTBAK PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ORTHES(NM,N,LOW,IGH,A,ORT) C INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW DOUBLE PRECISION A(NM,N),ORT(IGH) DOUBLE PRECISION F,G,H,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C A CONTAINS THE INPUT MATRIX. C C ON OUTPUT C C A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT C THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION C IS STORED IN THE REMAINING TRIANGLE UNDER THE C HESSENBERG MATRIX. C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE ORTRAN(NM,N,LOW,IGH,A,ORT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 DOUBLE PRECISION A(NM,IGH),ORT(IGH),Z(NM,N) DOUBLE PRECISION G C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY C TRANSFORMATIONS USED IN THE REDUCTION OF A REAL GENERAL C MATRIX TO UPPER HESSENBERG FORM BY ORTHES. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N. C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES C IN ITS STRICT LOWER TRIANGLE. C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY ORTHES. C C ORT HAS BEEN ALTERED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C C SUBROUTINE QZHES(NM,N,A,B,MATZ,Z) C INTEGER I,J,K,L,N,LB,L1,NM,NK1,NM1,NM2 DOUBLE PRECISION A(NM,N),B(NM,N),Z(NM,N) DOUBLE PRECISION R,S,T,U1,U2,V1,V2,RHO LOGICAL MATZ C C THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND C REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER C TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. C IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES. C C A CONTAINS A REAL GENERAL MATRIX. C C B CONTAINS A REAL GENERAL MATRIX. C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE. C C ON OUTPUT C C A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS C BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. C C B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS C BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF C MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C C SUBROUTINE QZIT(NM,N,A,B,EPS1,MATZ,Z,IERR) C INTEGER I,J,K,L,N,EN,K1,K2,LD,LL,L1,NA,NM,ISH,ITN,ITS,KM1,LM1, X ENM2,IERR,LOR1,ENORN DOUBLE PRECISION A(NM,N),B(NM,N),Z(NM,N) DOUBLE PRECISION R,S,T,A1,A2,A3,EP,SH,U1,U2,U3,V1,V2,V3,ANI,A11, X A12,A21,A22,A33,A34,A43,A44,BNI,B11,B12,B22,B33,B34, X B44,EPSA,EPSB,EPS1,ANORM,BNORM,EPSLON LOGICAL MATZ,NOTLAS C C THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, C AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM C IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. C IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING C ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM C OF THE OTHER MATRIX. IT IS USUALLY PRECEDED BY QZHES AND C FOLLOWED BY QZVAL AND, POSSIBLY, QZVEC. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES. C C A CONTAINS A REAL UPPER HESSENBERG MATRIX. C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. C C EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. C EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN C ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF C ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS C POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE C IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A C POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, C BUT LESS ACCURATE RESULTS. C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE. C C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION C BY QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. C C ON OUTPUT C C A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM. THE ELEMENTS C BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO C CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. C C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS C HAVE BEEN ALTERED. THE LOCATION B(N,1) IS USED TO STORE C EPS1 TIMES THE NORM OF B FOR LATER USE BY QZVAL AND QZVEC. C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS C (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE.. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE QZVAL(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z) C INTEGER I,J,N,EN,NA,NM,NN,ISW DOUBLE PRECISION A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) DOUBLE PRECISION C,D,E,R,S,T,AN,A1,A2,BN,CQ,CZ,DI,DR,EI,TI,TR,U1, X U2,V1,V2,A1I,A11,A12,A2I,A21,A22,B11,B12,B22,SQI,SQR, X SSI,SSR,SZI,SZR,A11I,A11R,A12I,A12R,A22I,A22R,EPSB LOGICAL MATZ C C THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM C IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. C IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY C REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX C EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE C GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES C AND QZIT AND MAY BE FOLLOWED BY QZVEC. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES. C C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) C COMPUTED AND SAVED IN QZIT. C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE. C C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES C AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. C C ON OUTPUT C C A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX C IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO C PAIRS OF COMPLEX EIGENVALUES. C C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS C HAVE BEEN ALTERED. B(N,1) IS UNALTERED. C C ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE C DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE C OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM C BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR C IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE. C C BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, C NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED C EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS C (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE QZVEC(NM,N,A,B,ALFR,ALFI,BETA,Z) C INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN,ISW,ENM2 DOUBLE PRECISION A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) DOUBLE PRECISION D,Q,R,S,T,W,X,Y,DI,DR,RA,RR,SA,TI,TR,T1,T2,W1,X1, X ZZ,Z1,ALFM,ALMI,ALMR,BETM,EPSB C C THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM IN C QUASI-TRIANGULAR FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO C A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR C FORM. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND C TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. C IT IS USUALLY PRECEDED BY QZHES, QZIT, AND QZVAL. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES. C C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) C COMPUTED AND SAVED IN QZIT. C C ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE C RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED C EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVAL. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTIONS BY QZHES, QZIT, AND QZVAL, IF PERFORMED. C IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE C DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. C C ON OUTPUT C C A IS UNALTERED. ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION C ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS. C C B HAS BEEN DESTROYED. C C ALFR, ALFI, AND BETA ARE UNALTERED. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. C IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND C THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. C IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX. C IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF C A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS C OF Z CONTAIN ITS EIGENVECTOR. C IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF C A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS C OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR. C EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS C OF ITS LARGEST COMPONENT IS 1.0 . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RATQR(N,EPS1,D,E,E2,M,W,IND,BD,TYPE,IDEF,IERR) C INTEGER I,J,K,M,N,II,JJ,K1,IDEF,IERR,JDEF DOUBLE PRECISION D(N),E(N),E2(N),W(N),BD(N) DOUBLE PRECISION F,P,Q,R,S,EP,QP,ERR,TOT,EPS1,DELTA,EPSLON INTEGER IND(N) LOGICAL TYPE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RATQR, C NUM. MATH. 11, 264-272(1968) BY REINSCH AND BAUER. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 257-265(1971). C C THIS SUBROUTINE FINDS THE ALGEBRAICALLY SMALLEST OR LARGEST C EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE C RATIONAL QR METHOD WITH NEWTON CORRECTIONS. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C EPS1 IS A THEORETICAL ABSOLUTE ERROR TOLERANCE FOR THE C COMPUTED EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, C OR INDEED SMALLER THAN ITS DEFAULT VALUE, IT IS RESET C AT EACH ITERATION TO THE RESPECTIVE DEFAULT VALUE, C NAMELY, THE PRODUCT OF THE RELATIVE MACHINE PRECISION C AND THE MAGNITUDE OF THE CURRENT EIGENVALUE ITERATE. C THE THEORETICAL ABSOLUTE ERROR IN THE K-TH EIGENVALUE C IS USUALLY NOT GREATER THAN K TIMES EPS1. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY. C C M IS THE NUMBER OF EIGENVALUES TO BE FOUND. C C IDEF SHOULD BE SET TO 1 IF THE INPUT MATRIX IS KNOWN TO BE C POSITIVE DEFINITE, TO -1 IF THE INPUT MATRIX IS KNOWN TO C BE NEGATIVE DEFINITE, AND TO 0 OTHERWISE. C C TYPE SHOULD BE SET TO .TRUE. IF THE SMALLEST EIGENVALUES C ARE TO BE FOUND, AND TO .FALSE. IF THE LARGEST EIGENVALUES C ARE TO BE FOUND. C C ON OUTPUT C C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS C (LAST) DEFAULT VALUE. C C D AND E ARE UNALTERED (UNLESS W OVERWRITES D). C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS SET TO 0.0D0 IF THE SMALLEST EIGENVALUES HAVE BEEN C FOUND, AND TO 2.0D0 IF THE LARGEST EIGENVALUES HAVE BEEN C FOUND. E2 IS OTHERWISE UNALTERED (UNLESS OVERWRITTEN BY BD). C C W CONTAINS THE M ALGEBRAICALLY SMALLEST EIGENVALUES IN C ASCENDING ORDER, OR THE M LARGEST EIGENVALUES IN C DESCENDING ORDER. IF AN ERROR EXIT IS MADE BECAUSE OF C AN INCORRECT SPECIFICATION OF IDEF, NO EIGENVALUES C ARE FOUND. IF THE NEWTON ITERATES FOR A PARTICULAR C EIGENVALUE ARE NOT MONOTONE, THE BEST ESTIMATE OBTAINED C IS RETURNED AND IERR IS SET. W MAY COINCIDE WITH D. C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. C C BD CONTAINS REFINED BOUNDS FOR THE THEORETICAL ERRORS OF THE C CORRESPONDING EIGENVALUES IN W. THESE BOUNDS ARE USUALLY C WITHIN THE TOLERANCE SPECIFIED BY EPS1. BD MAY COINCIDE C WITH E2. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 6*N+1 IF IDEF IS SET TO 1 AND TYPE TO .TRUE. C WHEN THE MATRIX IS NOT POSITIVE DEFINITE, OR C IF IDEF IS SET TO -1 AND TYPE TO .FALSE. C WHEN THE MATRIX IS NOT NEGATIVE DEFINITE, C 5*N+K IF SUCCESSIVE ITERATES TO THE K-TH EIGENVALUE C ARE NOT MONOTONE INCREASING, WHERE K REFERS C TO THE LAST SUCH OCCURRENCE. C C NOTE THAT SUBROUTINE TRIDIB IS GENERALLY FASTER AND MORE C ACCURATE THAN RATQR IF THE EIGENVALUES ARE CLUSTERED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE REBAK(NM,N,B,DL,M,Z) C INTEGER I,J,K,M,N,I1,II,NM DOUBLE PRECISION B(NM,N),DL(N),Z(NM,M) DOUBLE PRECISION X C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REBAKA, C NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A GENERALIZED C SYMMETRIC EIGENSYSTEM BY BACK TRANSFORMING THOSE OF THE C DERIVED SYMMETRIC MATRIX DETERMINED BY REDUC. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX SYSTEM. C C B CONTAINS INFORMATION ABOUT THE SIMILARITY TRANSFORMATION C (CHOLESKY DECOMPOSITION) USED IN THE REDUCTION BY REDUC C IN ITS STRICT LOWER TRIANGLE. C C DL CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATION. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE REBAKB(NM,N,B,DL,M,Z) C INTEGER I,J,K,M,N,I1,II,NM DOUBLE PRECISION B(NM,N),DL(N),Z(NM,M) DOUBLE PRECISION X C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REBAKB, C NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A GENERALIZED C SYMMETRIC EIGENSYSTEM BY BACK TRANSFORMING THOSE OF THE C DERIVED SYMMETRIC MATRIX DETERMINED BY REDUC2. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX SYSTEM. C C B CONTAINS INFORMATION ABOUT THE SIMILARITY TRANSFORMATION C (CHOLESKY DECOMPOSITION) USED IN THE REDUCTION BY REDUC2 C IN ITS STRICT LOWER TRIANGLE. C C DL CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATION. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE REDUC(NM,N,A,B,DL,IERR) C INTEGER I,J,K,N,I1,J1,NM,NN,IERR DOUBLE PRECISION A(NM,N),B(NM,N),DL(N) DOUBLE PRECISION X,Y C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REDUC1, C NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971). C C THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEM C AX=(LAMBDA)BX, WHERE B IS POSITIVE DEFINITE, TO THE STANDARD C SYMMETRIC EIGENPROBLEM USING THE CHOLESKY FACTORIZATION OF B. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. IF THE CHOLESKY C FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED C WITH A MINUS SIGN. C C A AND B CONTAIN THE REAL SYMMETRIC INPUT MATRICES. ONLY THE C FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED. IF C N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS, C INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L. C C DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L. C C ON OUTPUT C C A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE C OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE C STANDARD FORM. THE STRICT UPPER TRIANGLE OF A IS UNALTERED. C C B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER C TRIANGLE OF ITS CHOLESKY FACTOR L. THE FULL UPPER C TRIANGLE OF B IS UNALTERED. C C DL CONTAINS THE DIAGONAL ELEMENTS OF L. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 7*N+1 IF B IS NOT POSITIVE DEFINITE. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE REDUC2(NM,N,A,B,DL,IERR) C INTEGER I,J,K,N,I1,J1,NM,NN,IERR DOUBLE PRECISION A(NM,N),B(NM,N),DL(N) DOUBLE PRECISION X,Y C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REDUC2, C NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971). C C THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEMS C ABX=(LAMBDA)X OR BAY=(LAMBDA)Y, WHERE B IS POSITIVE DEFINITE, C TO THE STANDARD SYMMETRIC EIGENPROBLEM USING THE CHOLESKY C FACTORIZATION OF B. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. IF THE CHOLESKY C FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED C WITH A MINUS SIGN. C C A AND B CONTAIN THE REAL SYMMETRIC INPUT MATRICES. ONLY THE C FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED. IF C N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS, C INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L. C C DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L. C C ON OUTPUT C C A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE C OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE C STANDARD FORM. THE STRICT UPPER TRIANGLE OF A IS UNALTERED. C C B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER C TRIANGLE OF ITS CHOLESKY FACTOR L. THE FULL UPPER C TRIANGLE OF B IS UNALTERED. C C DL CONTAINS THE DIAGONAL ELEMENTS OF L. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 7*N+1 IF B IS NOT POSITIVE DEFINITE. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR) C INTEGER N,NM,IS1,IS2,IERR,MATZ DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,N),FV1(N) INTEGER IV1(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL GENERAL MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C A CONTAINS THE REAL GENERAL MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. COMPLEX CONJUGATE C PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE C EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS C IF MATZ IS NOT ZERO. IF THE J-TH EIGENVALUE IS REAL, THE C J-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF THE J-TH C EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE C J-TH AND (J+1)-TH COLUMNS OF Z CONTAIN THE REAL AND C IMAGINARY PARTS OF ITS EIGENVECTOR. THE CONJUGATE OF THIS C VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR C AND HQR2. THE NORMAL COMPLETION CODE IS ZERO. C C IV1 AND FV1 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RGG(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) LOGICAL TF C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C FOR THE REAL GENERAL GENERALIZED EIGENPROBLEM AX = (LAMBDA)BX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS A REAL GENERAL MATRIX. C C B CONTAINS A REAL GENERAL MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE NUMERATORS OF THE EIGENVALUES. C C BETA CONTAINS THE DENOMINATORS OF THE EIGENVALUES, C WHICH ARE THUS GIVEN BY THE RATIOS (ALFR+I*ALFI)/BETA. C COMPLEX CONJUGATE PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY C WITH THE EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST. C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS C IF MATZ IS NOT ZERO. IF THE J-TH EIGENVALUE IS REAL, THE C J-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF THE J-TH C EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE C J-TH AND (J+1)-TH COLUMNS OF Z CONTAIN THE REAL AND C IMAGINARY PARTS OF ITS EIGENVECTOR. THE CONJUGATE OF THIS C VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR QZIT. C THE NORMAL COMPLETION CODE IS ZERO. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C A CONTAINS THE REAL SYMMETRIC MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSB(NM,N,MB,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,MB,NM,IERR,MATZ DOUBLE PRECISION A(NM,MB),W(N),Z(NM,N),FV1(N),FV2(N) LOGICAL TF C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC BAND MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C MB IS THE HALF BAND WIDTH OF THE MATRIX, DEFINED AS THE C NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL C DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE C LOWER TRIANGLE OF THE MATRIX. C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C BAND MATRIX. ITS LOWEST SUBDIAGONAL IS STORED IN THE C LAST N+1-MB POSITIONS OF THE FIRST COLUMN, ITS NEXT C SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND C FINALLY ITS PRINCIPAL DIAGONAL IN THE N POSITIONS C OF THE LAST COLUMN. CONTENTS OF STORAGES NOT PART C OF THE MATRIX ARE ARBITRARY. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSG(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM AX = (LAMBDA)BX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS A REAL SYMMETRIC MATRIX. C C B CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSGAB(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM ABX = (LAMBDA)X. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS A REAL SYMMETRIC MATRIX. C C B CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSGBA(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),B(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C FOR THE REAL SYMMETRIC GENERALIZED EIGENPROBLEM BAX = (LAMBDA)X. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS A REAL SYMMETRIC MATRIX. C C B CONTAINS A POSITIVE DEFINITE REAL SYMMETRIC MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSM(NM,N,A,W,M,Z,FWORK,IWORK,IERR) C INTEGER N,NM,M,IWORK(N),IERR DOUBLE PRECISION A(NM,N),W(N),Z(NM,M),FWORK(1) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND ALL OF THE EIGENVALUES AND SOME OF THE EIGENVECTORS C OF A REAL SYMMETRIC MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C A CONTAINS THE REAL SYMMETRIC MATRIX. C C M THE EIGENVECTORS CORRESPONDING TO THE FIRST M EIGENVALUES C ARE TO BE COMPUTED. C IF M = 0 THEN NO EIGENVECTORS ARE COMPUTED. C IF M = N THEN ALL OF THE EIGENVECTORS ARE COMPUTED. C C ON OUTPUT C C W CONTAINS ALL N EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE ORTHONORMAL EIGENVECTORS ASSOCIATED WITH C THE FIRST M EIGENVALUES. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT, C IMTQLV AND TINVIT. THE NORMAL COMPLETION CODE IS ZERO. C C FWORK IS A TEMPORARY STORAGE ARRAY OF DIMENSION 8*N. C C IWORK IS AN INTEGER TEMPORARY STORAGE ARRAY OF DIMENSION N. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RSP(NM,N,NV,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER I,J,N,NM,NV,IERR,MATZ DOUBLE PRECISION A(NV),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC PACKED MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C NV IS AN INTEGER VARIABLE SET EQUAL TO THE C DIMENSION OF THE ARRAY A AS SPECIFIED FOR C A IN THE CALLING PROGRAM. NV MUST NOT BE C LESS THAN N*(N+1)/2. C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C PACKED MATRIX STORED ROW-WISE. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RST(NM,N,W,E,MATZ,Z,IERR) C INTEGER I,J,N,NM,IERR,MATZ DOUBLE PRECISION W(N),E(N),Z(NM,N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC TRIDIAGONAL MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C W CONTAINS THE DIAGONAL ELEMENTS OF THE REAL C SYMMETRIC TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE MATRIX IN C ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR IMTQL1 C AND IMTQL2. THE NORMAL COMPLETION CODE IS ZERO. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE RT(NM,N,A,W,MATZ,Z,FV1,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,3),W(N),Z(NM,N),FV1(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A SPECIAL REAL TRIDIAGONAL MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C A CONTAINS THE SPECIAL REAL TRIDIAGONAL MATRIX IN ITS C FIRST THREE COLUMNS. THE SUBDIAGONAL ELEMENTS ARE STORED C IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, THE C DIAGONAL ELEMENTS IN THE SECOND COLUMN, AND THE SUPERDIAGONAL C ELEMENTS IN THE FIRST N-1 POSITIONS OF THE THIRD COLUMN. C ELEMENTS A(1,1) AND A(N,3) ARE ARBITRARY. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR IMTQL1 C AND IMTQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 IS A TEMPORARY STORAGE ARRAY. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1) C INTEGER I,J,K,L,M,N,II,I1,KK,K1,LL,L1,MN,NM,ITS,IERR DOUBLE PRECISION A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N) DOUBLE PRECISION C,F,G,H,S,X,Y,Z,TST1,TST2,SCALE,PYTHAG LOGICAL MATU,MATV C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD, C NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). C C THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION C T C A=USV OF A REAL M BY N RECTANGULAR MATRIX. HOUSEHOLDER C BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. NOTE THAT NM MUST BE AT LEAST C AS LARGE AS THE MAXIMUM OF M AND N. C C M IS THE NUMBER OF ROWS OF A (AND U). C C N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V. C C A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED. C C MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE C DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. C C MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE C DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. C C ON OUTPUT C C A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V). C C W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE C DIAGONAL ELEMENTS OF S). THEY ARE UNORDERED. IF AN C ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT C FOR INDICES IERR+1,IERR+2,...,N. C C U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE C DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE. OTHERWISE C U IS USED AS A TEMPORARY ARRAY. U MAY COINCIDE WITH A. C IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING C TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. C C V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF C MATV HAS BEEN SET TO .TRUE. OTHERWISE V IS NOT REFERENCED. C V MAY ALSO COINCIDE WITH A IF U IS NOT NEEDED. IF AN ERROR C EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF C CORRECT SINGULAR VALUES SHOULD BE CORRECT. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C K IF THE K-TH SINGULAR VALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C RV1 IS A TEMPORARY STORAGE ARRAY. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TINVIT(NM,N,D,E,E2,M,W,IND,Z, X IERR,RV1,RV2,RV3,RV4,RV6) C INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M), X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) DOUBLE PRECISION U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,EPSLON, X PYTHAG INTEGER IND(M) C C THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH- C NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL C SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, C USING INVERSE ITERATION. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E, C WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. C E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN C THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM C OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN C 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0D0 C IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT, C TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES, C THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE. C C M IS THE NUMBER OF SPECIFIED EIGENVALUES. C C W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC. C C ON OUTPUT C C ALL INPUT ARRAYS ARE UNALTERED. C C Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH C EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. C C RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TQL1(N,D,E,IERR) C INTEGER I,J,L,M,N,II,L1,L2,MML,IERR DOUBLE PRECISION D(N),E(N) DOUBLE PRECISION C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E HAS BEEN DESTROYED. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1. C C E HAS BEEN DESTROYED. C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TQLRAT(N,D,E2,IERR) C INTEGER I,J,L,M,N,II,L1,MML,IERR DOUBLE PRECISION D(N),E2(N) DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E2 HAS BEEN DESTROYED. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TRBAK1(NM,N,A,E,M,Z) C INTEGER I,J,K,L,M,N,NM DOUBLE PRECISION A(NM,N),E(N),Z(NM,M) DOUBLE PRECISION S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK1, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED1. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION BY TRED1 C IN ITS STRICT LOWER TRIANGLE. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C NOTE THAT TRBAK1 PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TRBAK3(NM,N,NV,A,M,Z) C INTEGER I,J,K,L,M,N,IK,IZ,NM,NV DOUBLE PRECISION A(NV),Z(NM,M) DOUBLE PRECISION H,S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS C USED IN THE REDUCTION BY TRED3 IN ITS FIRST C N*(N+1)/2 POSITIONS. C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TRED1(NM,N,A,D,E,E2) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N) DOUBLE PRECISION F,G,H,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX C TO A SYMMETRIC TRIDIAGONAL MATRIX USING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TRED2(NM,N,A,D,E,Z) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N) DOUBLE PRECISION F,G,H,HH,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX C PRODUCED IN THE REDUCTION. C C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TRED3(N,NV,A,D,E,E2) C INTEGER I,J,K,L,N,II,IZ,JK,NV,JM1 DOUBLE PRECISION A(NV),D(N),E(N),E2(N) DOUBLE PRECISION F,G,H,HH,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS C A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX C USING ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL C ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS. C C ON OUTPUT C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL C TRANSFORMATIONS USED IN THE REDUCTION. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C C SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5) C INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM DOUBLE PRECISION D(N),E(N),E2(N),W(M),RV4(N),RV5(N) DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON INTEGER IND(M) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BISECT, C NUM. MATH. 9, 386-393(1967) BY BARTH, MARTIN, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL C SYMMETRIC MATRIX BETWEEN SPECIFIED BOUNDARY INDICES, C USING BISECTION. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED C EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, C IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, C NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE C PRECISION AND THE 1-NORM OF THE SUBMATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY. C C M11 SPECIFIES THE LOWER BOUNDARY INDEX FOR THE DESIRED C EIGENVALUES. C C M SPECIFIES THE NUMBER OF EIGENVALUES DESIRED. THE UPPER C BOUNDARY INDEX M22 IS THEN OBTAINED AS M22=M11+M-1. C C ON OUTPUT C C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS C (LAST) DEFAULT VALUE. C C D AND E ARE UNALTERED. C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS ALSO SET TO ZERO. C C LB AND UB DEFINE AN INTERVAL CONTAINING EXACTLY THE DESIRED C EIGENVALUES. C C W CONTAINS, IN ITS FIRST M POSITIONS, THE EIGENVALUES C BETWEEN INDICES M11 AND M22 IN ASCENDING ORDER. C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 3*N+1 IF MULTIPLE EIGENVALUES AT INDEX M11 MAKE C UNIQUE SELECTION IMPOSSIBLE, C 3*N+2 IF MULTIPLE EIGENVALUES AT INDEX M22 MAKE C UNIQUE SELECTION IMPOSSIBLE. C C RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. C C NOTE THAT SUBROUTINE TQL1, IMTQL1, OR TQLRAT IS GENERALLY FASTER C THAN TRIDIB, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C SUBROUTINE TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z, X IERR,RV1,RV2,RV3,RV4,RV5,RV6) C INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS, X IERR,GROUP,ISTURM DOUBLE PRECISION D(N),E(N),E2(N),W(MM),Z(NM,MM), X RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N) DOUBLE PRECISION U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4, X NORM,TST1,TST2,EPSLON,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRISTURM C BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL C SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL AND THEIR C ASSOCIATED EIGENVECTORS, USING BISECTION AND INVERSE ITERATION. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED C EIGENVALUES. IT SHOULD BE CHOSEN COMMENSURATE WITH C RELATIVE PERTURBATIONS IN THE MATRIX ELEMENTS OF THE C ORDER OF THE RELATIVE MACHINE PRECISION. IF THE C INPUT EPS1 IS NON-POSITIVE, IT IS RESET FOR EACH C SUBMATRIX TO A DEFAULT VALUE, NAMELY, MINUS THE C PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE C 1-NORM OF THE SUBMATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY. C C LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. C IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. C C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF C EIGENVALUES IN THE INTERVAL. WARNING. IF MORE THAN C MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, C AN ERROR RETURN IS MADE WITH NO VALUES OR VECTORS FOUND. C C ON OUTPUT C C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS C (LAST) DEFAULT VALUE. C C D AND E ARE UNALTERED. C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS ALSO SET TO ZERO. C C M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). C C W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER IF THE MATRIX C DOES NOT SPLIT. IF THE MATRIX SPLITS, THE EIGENVALUES ARE C IN ASCENDING ORDER FOR EACH SUBMATRIX. IF A VECTOR ERROR C EXIT IS MADE, W CONTAINS THOSE VALUES ALREADY FOUND. C C Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. C IF AN ERROR EXIT IS MADE, Z CONTAINS THOSE VECTORS C ALREADY FOUND. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 3*N+1 IF M EXCEEDS MM. C 4*N+R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH C EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. C C RV1, RV2, RV3, RV4, RV5, AND RV6 ARE TEMPORARY STORAGE ARRAYS. C C THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM C APPEARS IN TSTURM IN-LINE. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C .