.QA EISPACK .ts +5,+8,+10,+2,+2,+2,+2 .p EISPACK is a systematized collection of subroutines which compute the eigenvalues and/or eigenvectors of three special classes of matrix problems. The three problem classes are : .list 0 .le;- real symmetric band matrices, .le;- generalized real symmetric matrices, and .le;- generalized real matrices. .els 1 .p The singular value decomposition of an arbitrary matrix can also be computed, which, in turn, enables the solution of certain linear least squares problems. .QB Subroutines .QC parameters .nf Temporary variable Types: i = integer f = floating s = scaler v = vector m = matrix Array parameters: nm,mm,mb - defines dimension of variable A(nm,mb) N = order of matrix nm <= N mb <= N TRANSFORM: MATZ = .true. - accumulate transfromation matrix Z(nm,N) = transformation matrix ERRORS: IERR = 0 no error .QC BAKVEC .nf .f.P forms the eigenvectors of a certain real non-symmetric tridiagonal matrix from the eigenvectors of that symmetric tridiagonal matrix determined by FIGI (F280). .QC BALANC .nf .f.P balances a real general matrix and isolates eigenvalues whenever possible. Sums of the magnitudes of elements in corresponding rows and columns are made nearly equal by exact similarity transformations, and eigenvalues are isolated by permutation similarity transformations. Balancing reduces the 1-norm of the original matrix whenever sums of the magnitudes of elements in some row and corresponding column are markedly different, while at the same time leaving the eigenvalues unchanged. Reducing the norm in this way can improve the accuracy of the computed eigenvalues and/or eigenvectors. .QC BALBAK .nf .f.P Forms the eigenvectors of a real general matrix from the eigenvectors of that matrix transformed by BALANC (F269). .QC BANDR .nf .br;BANDR(nm,N,mb,A,D,E,E2,MATZ,Z) Input: nm,N,mb,A,MATZ Output: D,E,E2,Z D(N) = Diagonal elements E(N) = subdiagonal elements E2(N) = squares of subdiagonal elem .nf .f.P reduces a real symmetric band matrix to a symmetric tridiagonal matrix using and optionally accumulating orthogonal similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC BANDV .nf BANDV(nm,N,mb,A,E21,M,W,Z,IERR,NV,RV,RV6) Input: nm,N,MBW,A,E21,M,W,NV mb = Number of cols in A used for band. E21 = 0. eigenvalues in ascending order = 2. descending order (eigen problem) = 1. coef matrix is symm. (lin eq.) = -1. coef matrix is not. W(M) = eigenvalues = define coef matrix with A Output: IERR,Z,RV,RV6 Z(nm,m) = eigenvectors RV(N*(2*MB-1) = upper triangular matrix RV6(N) = approx eigenvec = solution vectors .f.P determines those eigenvectors of a real symmetric band matrix corresponding to a set of ordered approximate eigenvalues, using inverse iteration. It can also be used to solve a system of linear equations with a symmetric or non-symmetric band coefficient matrix, or a succession of related systems. .QC BISECT .nf BISECT(N,EPS1,D,E,Ex,LB,UB,MM,M,W,IND,IERR,RV4,RV5) Input: N,EPS1,D,E,E2,LB,UB,MM Output: E2,W,IND,IERR .f.P determines those eigenvalues of a symmetric tridiagonal matrix in a specified interval using Sturm sequencing. .QC BQR .nf BQR(nm,N,mb,A,T,R,IERR,NV,RV) Input: nm,N,mb,A,T,R,NV T = shift for diagonal elements of A R = 0 for first call to BQR Output: R,IERR,RV NV = 2*mb**2 + 4*MB -3 RV(NV) = aux storage + half lengths + vectors .f.P determines some of the eigenvalues of a real symmetric band matrix using the QR algorithm with shifts of origin. It determines one eigenvalue for each call, deflating the matrix before returning. The subroutine could be used to determine all of the eigenvalues; however, more efficient routines exist in this case (see section 2C). .QC CBABK2 .nf .f.P forms the eigenvectors of a complex general matrix from the eigenvectors of that matrix transformed by CBAL (F271). .QC CBAL .nf .f.P balances a complex general matrix and isolates eigenvalues whenever possible. Sums of the magnitudes of elements in corresponding rows and columns are made nearly equal by exact similarity transformations, and eigenvalues are isolated by permutation similarity transformations. Balancing reduces the 1-norm of the original matrix whenever sums of the magnitudes of elements in some row and corresponding column are markedly different, while at the same time leaving the eigenvalues unchanged. Reducing the norm in this way can improve the accuracy of the computed eigenvalues and/or eigenvectors. .QC CG .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a complex general matrix. .QC CH .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a complex Hermitian matrix. .QC CINVIT .nf .f.P computes the eigenvectors of a complex upper Hessenberg matrix corresponding to specified eigenvalues using inverse iteration. .QC COMBAK .nf .f.P forms the eigenvectors of a complex general matrix from the eigenvectors of that upper Hessenberg matrix determined by COMHES (F282). .QC COMHES .nf .f.P reduces a complex general matrix to complex upper Hessenberg form using stabilized elementary similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC COMLR .nf .f.P computes the eigenvalues of a complex upper Hessenberg matrix using the modified LR method. .QC COMLR2 .nf .f.P computes the eigenvalues and eigenvectors of a complex upper Hessenberg matrix. COMLR2 uses the modified LR method to compute the eigenvalues and accumulates the LR transformations to compute the eigenvectors. The eigenvectors of a complex general matrix can also be computed directly by COMLR2, if COMHES (F282) has been used to reduce this matrix to Hessenberg form. .QC COMQR .nf .f.P computes the eigenvalues of a complex upper Hessenberg matrix using the QR method. .QC COMQR2 .nf .f.P computes the eigenvalues and eigenvectors of a complex upper Hessenberg matrix. COMQR2 uses the QR method to compute the eigenvalues and accumulates the QR transformations to compute the eigenvectors. The eigenvectors of a complex general matrix can also be computed directly by COMQR2, if CORTH (F244) has been used to reduce this matrix to Hessenberg form. .QC CORTB .nf .f.P forms the eigenvectors of a complex general matrix from the eigenvectors of that upper Hessenberg matrix determined by CORTH (F244). .QC CORTH .nf .f.P reduces a complex general matrix to complex upper Hessenberg form using unitary similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC ELMBAK .nf .f.P forms the eigenvectors of a real general matrix from the eigenvectors of that upper Hessenberg matrix determined by ELMHES (F273). .QC ELMHES .nf .f.P reduces a real general matrix to upper Hessenberg form using stabilized elementary similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC ELTRAN .nf .f.P accumulates the stabilized elementary similarity transformations used in the reduction of a real general matrix to upper Hessenberg form by ELMHES (F273). .QC FIGI .nf .f.P transforms a certain real non-symmetric tridiagonal matrix to a symmetric tridiagonal matrix. The property of the non-symmetric matrix required for use of this subroutine is that the products of pairs of corresponding off-diagonal elements be all non-negative. The transformed matrix is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC FIGI2 .nf .f.P transforms a certain real non-symmetric tridiagonal matrix to a symmetric tridiagonal matrix using and accumulating diagonal similarity transformations. The property of the non-symmetric matrix required for use of this subroutine is that the products of pairs of corresponding off-diagonal elements be all non- negative, and zero only when both factors are zero. The symmetrized matrix and the transformation matrix are used by subroutine IMTQL2 (F292) to find the eigenvalues and eigenvectors of the original matrix. .QC HQR .nf .f.P computes the eigenvalues of a real upper Hessenberg matrix using the QR method. .QC HQR2 .nf .f.P computes the eigenvalues and eigenvectors of a real upper Hessenberg matrix using the QR method. The eigenvectors of a real general matrix can also be computed if ELMHES (F273) and ELTRAN (F220) or ORTHES (F275) and ORTRAN (F221) have been used to reduce this general matrix to Hessenberg form and to accumulate the transformations. .QC HTRIBK .nf .f.P forms the eigenvectors of a complex Hermitian matrix from the eigenvectors of that real symmetric tridiagonal matrix determined by HTRIDI (F284). .QC HTRIB3 .nf .f.P forms the eigenvectors of a complex Hermitian matrix from the eigenvectors of that real symmetric tridiagonal matrix determined by HTRID3 (F242). .QC HTRIDI .nf .f.P reduces a complex Hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC HTRID3 .nf .f.P reduces a complex Hermitian matrix, stored as a single square array, to a real symmetric tridiagonal matrix using unitary similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC IMTQLV .nf .f.P is a variant of IMTQL1 (F291) that determines the eigenvalues of a symmetric tridiagonal matrix using the implicit QL method, and associates with them their corresponding submatrix indices. This feature and the preservation of the input matrix make it possible to later determine the eigenvectors of the matrix. .QC IMTQL1 .nf .f.P determines the eigenvalues of a symmetric tridiagonal matrix using the implicit QL method. .QC IMTQL2 .nf .f.P determines the eigenvalues and eigenvectors of a symmetric tridiagonal matrix. IMTQL2 uses the implicit QL method to compute the eigenvalues and accumulates the QL transformations to compute the eigenvectors. The eigenvectors of a full symmetric matrix can also be computed directly by IMTQL2, if TRED2 (F278) has been used to reduce this matrix to tridiagonal form. .QC INVIT .nf .f.P computes the eigenvectors of a real upper Hessenberg matrix corresponding to specified eigenvalues using inverse iteration. .QC MINFIT .nf MINFIT(NM,M,N,A,W,IP,B,IERR,RV1) Input: nm,M,M,A,IP,B M = number of rows of A,B <= nm B(NM,IP) = constant matrix IP = 0 (B is not used) Output: A,W,B A(nm,N) = orthogonal matrix W(N) = singular values B(NM,IP) = U B .f.P computes the singular values and complete orthogonal decomposition of the real rectangular coefficient matrix A of the linear system AX = B. A is decomposed into T U * DIAG(S) * V T T T with U U = V V = I. U B is formed rather than U itself, enabling a determination of the solution X of minimal norm (see section 2C). The diagonal elements of S are the singular values of A, equal to the non-negative T square roots of the eigenvalues of A A. .QC ORTBAK .nf .f.P forms the eigenvectors of a real general matrix from the eigenvectors of that upper Hessenberg matrix determined by ORTHES (F275). .QC ORTHES .nf .f.P reduces a real general matrix to upper Hessenberg form using orthogonal similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC ORTRAN .nf .f.P accumulates the orthogonal similarity transformations used in the reduction of a real general matrix to upper Hessenberg form by ORTHES (F275). .QC QZHES .nf QZHES(nm,M,A,B,MATZ,A) Input: nm,N,A,B,MATZ A,B = input real matrices Output: A,B,Z A = upper Heisenberg matrix B = upper triangular matrix .f.P reduces the generalized real eigenproblem AX = (LAMBDA)BX to an equivalent problem with A in upper Hessenberg form and B in upper triangular form, using and optionally accumulating orthogonal transformations. These reduced forms are used by other subroutines to find the generalized eigenvalues and/or eigenvectors of the original matrix system. See section 2C for the specific routines. .QC QZIT .nf QZIT(nm,N,A,B,EPS1,MATZ,Z,IERR) Input: nm,N,A,B,EPS1 A = upper heisenberg matrix B = upper triangular mat. EPS1 = tolerance for negligible elements Output: A,B,Z,IERR A = quasi-upper triangular B = new upper tri. .f.P reduces the generalized real eigenproblem AX = (LAMBDA)BX, where A is in upper Hessenberg form and B is in upper triangular form, to an equivalent problem with A now in quasi-upper triangular form and B still in upper triangular form, using and optionally accumulating orthogonal transformations. These reduced forms are used by other subroutines to find the generalized eigenvalues and/or eigenvectors of this matrix system or of an original real general matrix system from which this system was itself reduced. See section 2C for the specific routines. .QC QZVAL .nf QZVAL(nm,N,A,B,ALFR,ALFI,BETA,MATZ,Z) Input: nm,N,A,B,MATZ,Z A(nm,N) = quasi upper tri. matrix B(nm,N) = upper triang matrix. Output: ALFR,ALFI,BETA,Z ALFR(N) = real of diagonal elements of A ALFI(N) = imaginary of diagonal elements of A BETA(N) = real non neg elemen of B... .f.P reduces the generalized real eigenproblem AX = (LAMBDA)BX, where A is in quasi- upper triangular form and B is in upper triangular form, to an equivalent problem with A further reduced and B still in upper triangular form, using and optionally accumulating orthogonal transformations. The reduced A is upper triangular except for possible 2x2 diagonal blocks corresponding to pairs of complex eigenvalues. The subroutine returns quantities whose ratios give the generalized eigenvalues. The reduced forms may then be used to find the generalized eigenvectors of this matrix system or of an original real general matrix system from which this system was itself reduced. See section 2C. .QC QZVEC .nf QZVEC(nm,N,A,B,ALFR,ALFI,BETA,Z) Input: nm,N,A,B,ALFR,ALFI,BETA,Z Output: B,Z Z(nm,N) = Real and imaginary eigen values .f.P determines the eigenvectors corresponding to a set of specified eigenvalues for the generalized real eigenproblem AX = (LAMBDA)BX, where A is in quasi-upper triangular form and B is in upper triangular form, using back substitution. Quasi-triangular, in this case, means triangular except for possible 2x2 diagonal blocks corresponding to pairs of complex eigenvalues. The eigenvectors of a real general matrix system can also be computed if QZHES (F238), QZIT (F239), and QZVAL (F240) have been used to reduce the matrices of the system to quasi- triangular and triangular form and to accumulate the transformations. .QC RATQR .nf .f.P determines the algebraically smallest or largest eigenvalues of a symmetric tridiagonal matrix using the rational QR method with Newton corrections. .QC REBAKB .nf .f.P forms the eigenvectors of a generalized symmetric eigensystem from the eigenvectors of that derived symmetric matrix determined by REDUC2 (F230). .br;SEE - REBAK .QC REBAK .nf REBAK(nm,N,B,DL,M,Z) Input: nm,N,B,DL,M,Z M = number of cols of Z to back transform Z(nm,N) = eigenvect to transform Ouptut: Z Z(nm,M) = transformed .f.P forms the eigenvectors of a generalized symmetric eigensystem from the eigenvectors of that derived symmetric matrix determined by REDUC (F224) or REDUC2 (F230). .QC REDUC .nf REDUC(nm,N,A,B,DL,IERR) Input: nm,N,A,B,DL DL(N) = diag. elem. of Cholesky factor if N<0 Output: A,B,DL,IERR DL = diag. elem. of Cholescky factor L of B .f.P reduces the generalized symmetric eigenproblem AX = (LAMBDA)BX, where B is positive definite, to the standard symmetric eigenproblem using the Cholesky factorization of B. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix system. See section 2C for the specific routines. .QC REDUC2 .nf .f.P reduces either of the generalized symmetric eigenproblems ABX = (LAMBDA)X or BAX = (LAMBDA)X, where B is positive definite, to a standard symmetric eigenproblem using the Cholesky factorization of B. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix system. See section 2C for the specific routines. .br;SEE - REDUC .QC RG .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a real general matrix. .QC RGG .nf RGG(nm,N,A,B,ALFR,ALFI,BETA,MATZ,Z,IERR) Input: nm,N,A,B,MATZ Output: A,B,ALFR,ALFI,BETA,Z,IERR ALFR(N) = real numerator of eigenvalues ALFI(N) = imaginary "" BETA(N) = denominators of eigenvalues .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) for the real general generalized eigenproblem A*X = (LAMBDA)*B*X. .QC RSB .nf RSB(nm,N,mb,A,W,MATZ,Z,fv1,fv2,IERR) Input: nm,N,mb,A,MATZ A = lower triangle of symm matrix Output: A,W,Z,fv1,fv2,IERR A(nm,mb) W(N) = eigenvalues in ascending order .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a real symmetric band matrix. .QC RS .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a real symmetric matrix. .QC RSGAB .nf RSGAB(nm,N,A,B,W,MATZ,Z,fv1,fv2,IERR) Input: nm,N,A,B,MATZ A = upper triangle of symm matrix B(nm,N) = pos. def. upper triangle. Output: W,Z,fv1,fv2,IERR W(N) = eigenvalues in ascending order .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) for the real symmetric generalized eigenproblem A*B*X = (LAMBDA)*X. .QC RSGBA .nf RSGBA(nm,N,A,B,W,MATZ,Z,fv1,fv2,IERR) Input: nm,N,A,B,MATZ A = upper triangle of symm matrix B(nm,N) = pos. def. upper triangle. Output: W,Z,fv1,fv2,IERR W(N) = eigenvalues in ascending order .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) for the real symmetric generalized eigenproblem B*A*X = (LAMBDA)*X. .QC RSG .nf RSG(nm,N,A,B,W,MATZ,Z,fv1,fv2,IERR) Input: nm,N,A,B,MATZ A = upper triangle of symm matrix B(nm,N) = pos. def. upper triangle. Output: W,Z,fv1,fv2,IERR W(N) = eigenvalues in ascending order .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) for the real symmetric generalized eigenproblem A*X = (LAMBDA)*B*X. .QC RSP .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a real symmetric packed matrix. .QC RST .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a real symmetric tridiagonal matrix. .QC RT .nf .f.P calls the recommended sequence of subroutines from the eigensystem subroutine package EISPACK to determine the eigenvalues and eigenvectors (if desired) of a certain real tridiagonal matrix. The property of the matrix required for use of this subroutine is that the products of pairs of corresponding off-diagonal elements be all non-negative, and further if eigenvectors are desired, no product be zero unless both factors are zero. .QC SVD .nf SVD(nm,M,N,A,W,MATU,U,MATV,V,IERR,rv1) Input: nm,M,N,A,MATU,MATV Ouptut: W,U,V,IERR,rv1 U(nm,N) = orthogonal columns if MATU=.true. = temp storage if .false. V(nm,N) = V matrix in decomp. if MATV=.true. = dummy if .false. .f.P computes the singular values and complete orthogonal decomposition of a real rectangular matrix A. A is decomposed into T U * DIAG(S) * V T T with U U = V V = I. The diagonal elements of S are the singular values of A, equal to the non-negative square T roots of the eigenvalues of A A. .QC TINVIT .nf TINVIT(nm,N,D,E,E2,M,W,IND,Z,IERR,rv1,rv2,rv3,rv4,rv6) Input: nm,N,D,E,E2,M,W D(N) = diagonal elements E(N) = subdiagonal elem. E2(N) = sq. of elemets of E corr to negligible M = numnber of eigenvalues W(M) = ordered eigenvalues Output: Z,IERR,rvn... Z(nm,M) = eigenvectors .f.P determines those eigenvectors of a symmetric tridiagonal matrix corresponding to a set of ordered approximate eigenvalues, using inverse iteration. .QC TQLRAT .nf TQLRAT(N,D,E2,IERR) Input: N,D,E2 D(N) = diagonal elem E2(N) = sq of subdiagonal Output: D,E2,IERR D = ascending eigenvalues .f.P determines the eigenvalues of a symmetric tridiagonal matrix using a rational variant of the QL method. .QC TQL1 .nf .f.P determines the eigenvalues of a symmetric tridiagonal matrix using the QL method. .QC TQL2 .nf TQL2(nm,N,D,E,Z,IERR) Input: nm,N,D,E,Z D(N) = diagonal elements E(N) = subdiagonal Z = indentity matrix or = transform matrix from TRED2 Output: D = ascending eigenvalues Z(nm,N) = eigenvectors of this trid mat. = eigenv. of full symm matrix .f.P determines the eigenvalues and eigenvectors of a symmetric tridiagonal matrix. TQL2 uses the QL method to compute the eigenvalues and accumulates the QL transformations to compute the eigenvectors. The eigenvectors of a full symmetric matrix can also be computed directly by TQL2, if TRED2 (F278) has been used to reduce this matrix to tridiagonal form. .QC TRBAK1 .nf TRBAK1(nm,N,A,E,M,Z) Input: nm,N,A,E,M Z(nm,M) = eigenvectors to transform Output: Z Z(nm,M) = back transformed eigenv .f.P forms the eigenvectors of a real symmetric matrix from the eigenvectors of that symmetric tridiagonal matrix determined by TRED1 (F277). .QC TRBAK3 .nf .f.P forms the eigenvectors of a real symmetric matrix from the eigenvectors of that symmetric tridiagonal matrix determined by TRED3 (F228). .QC TRED1 .nf TRED1(nm,N,A,D,E,E2) Input: nm,N,A A(nm,N) = lower tri. symm. matrix Output: A,D,E,E2 A = info about trans in lower tri. D(N) = diagonal elem E(N) = subdiagonal E2(N) = sq of subdiag. .f.P reduces a real symmetric matrix to a symmetric tridiagonal matrix using orthogonal similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC TRED2 .nf TRED2(nm,N,A,D,E,E2) Input: nm,N,A A(nm,N) = lower tri. symm. matrix Output: A,D,E,E2 A = info about trans in lower tri. D(N) = diagonal elem E(N) = subdiagonal E2(N) = sq of subdiag. .f.P reduces a real symmetric matrix to a symmetric tridiagonal matrix using and accumulating orthogonal similarity transformations. This reduced form and the transformation matrix are used by subroutine TQL2 (F290) or IMTQL2 (F292) to find the eigenvalues and eigenvectors of the original matrix. .QC TRED3 .nf .f.P reduces a real symmetric matrix, stored as a one-dimensional array, to a symmetric tridiagonal matrix using orthogonal similarity transformations. This reduced form is used by other subroutines to find the eigenvalues and/or eigenvectors of the original matrix. See section 2C for the specific routines. .QC TRIDIB .nf .f.P determines those eigenvalues of a symmetric tridiagonal matrix between specified boundary indices using Sturm sequencing. .QC TSTURM .nf .f.P determines those eigenvalues of a symmetric tridiagonal matrix in a specified interval and their corresponding eigenvectors, using Sturm sequencing and inverse iteration. ******************************* .nf .QB Comments .BR This is a listing of comments from each subroutine .nf .QC CDIV__ SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI) REAL AR,AI,BR,BI,CR,CI C C COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) C .QC CSROOT__ SUBROUTINE CSROOT(XR,XI,YR,YI) REAL XR,XI,YR,YI C C (YR,YI) = COMPLEX SQRT(XR,XI) C BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI) C REAL FUNCTION EPSLON (X) REAL X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C REAL FUNCTION PYTHAG(A,B) REAL A,B C C FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C .QC BAKVEC__ SUBROUTINE BAKVEC(NM,N,T,E,M,Z,IERR) C INTEGER I,J,M,N,NM,IERR REAL 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 .QC BALANC__ SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL A(NM,N),SCALE(N) REAL 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 .QC BALBAK__ SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z) C INTEGER I,J,K,M,N,II,NM,IGH,LOW REAL SCALE(N),Z(NM,M) REAL 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 .QC BANDR__ 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 REAL A(NM,MB),D(N),E(N),E2(N),Z(NM,N) REAL 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 .QC BANDV__ 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 REAL A(NM,MBW),W(M),Z(NM,M),RV(NV),RV6(N) REAL 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.0E0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR C 2.0E0 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.0E0 IF THE COEFFICIENT C MATRIX IS SYMMETRIC AND TO -1.0E0 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 SQRT(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 .QC BISECT__ 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 REAL D(N),E(N),E2(N),W(MM),RV4(N),RV5(N) REAL 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 .QC BQR__ 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 REAL A(NM,MB),RV(NV) REAL 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 SQRT(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 .QC CBABK2__ SUBROUTINE CBABK2(NM,N,LOW,IGH,SCALE,M,ZR,ZI) C INTEGER I,J,K,M,N,II,NM,IGH,LOW REAL SCALE(N),ZR(NM,M),ZI(NM,M) REAL 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 .QC CBAL__ SUBROUTINE CBAL(NM,N,AR,AI,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL AR(NM,N),AI(NM,N),SCALE(N) REAL 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 .QC CG__ SUBROUTINE CG(NM,N,AR,AI,WR,WI,MATZ,ZR,ZI,FV1,FV2,FV3,IERR) C INTEGER N,NM,IS1,IS2,IERR,MATZ REAL 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 .QC CH__ SUBROUTINE CH(NM,N,AR,AI,W,MATZ,ZR,ZI,FV1,FV2,FM1,IERR) C INTEGER I,J,N,NM,IERR,MATZ REAL 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 .QC CINVIT__ 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 REAL 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) REAL 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 SQRT(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 .QC COMBAK__ SUBROUTINE COMBAK(NM,LOW,IGH,AR,AI,INT,M,ZR,ZI) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 REAL AR(NM,IGH),AI(NM,IGH),ZR(NM,M),ZI(NM,M) REAL 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 .QC COMHES__ SUBROUTINE COMHES(NM,N,LOW,IGH,AR,AI,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 REAL AR(NM,N),AI(NM,N) REAL 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 .QC COMLR__ 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 REAL HR(NM,N),HI(NM,N),WR(N),WI(N) REAL 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 .QC COMLR2__ 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 REAL HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N) REAL 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 .QC COMQR__ 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 REAL HR(NM,N),HI(NM,N),WR(N),WI(N) REAL 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 SQRT(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 .QC COMQR2__ 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 REAL HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N), X ORTR(IGH),ORTI(IGH) REAL 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 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.0E0 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 SQRT(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 .QC CORTB__ 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 REAL AR(NM,IGH),AI(NM,IGH),ORTR(IGH),ORTI(IGH), X ZR(NM,M),ZI(NM,M) REAL 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 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 .QC CORTH__ SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI) C INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW REAL AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH) REAL 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 SQRT(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 .QC ELMBAK__ SUBROUTINE ELMBAK(NM,LOW,IGH,A,INT,M,Z) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 REAL A(NM,IGH),Z(NM,M) REAL 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 .QC ELMHES__ SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 REAL A(NM,N) REAL 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 .QC ELTRAN__ SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 REAL 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 .QC FIGI__ SUBROUTINE FIGI(NM,N,T,D,E,E2,IERR) C INTEGER I,N,NM,IERR REAL 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 .QC FIGI2__ SUBROUTINE FIGI2(NM,N,T,D,E,Z,IERR) C INTEGER I,J,N,NM,IERR REAL T(NM,3),D(N),E(N),Z(NM,N) REAL 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 .QC HQR__ 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 REAL H(NM,N),WR(N),WI(N) REAL 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 .QC HQR2__ 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 REAL H(NM,N),WR(N),WI(N),Z(NM,N) REAL 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 .QC HTRIB3__ SUBROUTINE HTRIB3(NM,N,A,TAU,M,ZR,ZI) C INTEGER I,J,K,L,M,N,NM REAL A(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) REAL 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 .QC HTRIBK__ SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI) C INTEGER I,J,K,L,M,N,NM REAL AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M) REAL 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 .QC HTRID3__ SUBROUTINE HTRID3(NM,N,A,D,E,E2,TAU) C INTEGER I,J,K,L,N,II,NM,JM1,JP1 REAL A(NM,N),D(N),E(N),E2(N),TAU(2,N) REAL 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 SQRT(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 .QC HTRIDI__ SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU) C INTEGER I,J,K,L,N,II,NM,JP1 REAL AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N) REAL 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 SQRT(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 .QC IMTQL1__ SUBROUTINE IMTQL1(N,D,E,IERR) C INTEGER I,J,L,M,N,II,MML,IERR REAL D(N),E(N) REAL 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 SQRT(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 .QC IMTQL2__ SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,NM,MML,IERR REAL D(N),E(N),Z(NM,N) REAL 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 SQRT(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 .QC IMTQLV__ SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1) C INTEGER I,J,K,L,M,N,II,MML,TAG,IERR REAL D(N),E(N),E2(N),W(N),RV1(N) REAL 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 SQRT(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 .QC INVIT__ 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 REAL A(NM,N),WR(N),WI(N),Z(NM,MM),RM1(N,N), X RV1(N),RV2(N) REAL 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 SQRT(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 .QC MINFIT__ 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 REAL A(NM,N),W(N),B(NM,IP),RV1(N) REAL 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 SQRT(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 .QC ORTBAK__ SUBROUTINE ORTBAK(NM,LOW,IGH,A,ORT,M,Z) C INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1 REAL A(NM,IGH),ORT(IGH),Z(NM,M) REAL 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 .QC OTRHES__ SUBROUTINE ORTHES(NM,N,LOW,IGH,A,ORT) C INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW REAL A(NM,N),ORT(IGH) REAL 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 .QC ORTRAN__ SUBROUTINE ORTRAN(NM,N,LOW,IGH,A,ORT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 REAL A(NM,IGH),ORT(IGH),Z(NM,N) REAL 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 .QC QZHES__ SUBROUTINE QZHES(NM,N,A,B,MATZ,Z) C INTEGER I,J,K,L,N,LB,L1,NM,NK1,NM1,NM2 REAL A(NM,N),B(NM,N),Z(NM,N) REAL 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 .QC QZIT__ 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 REAL A(NM,N),B(NM,N),Z(NM,N) REAL 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 .QC QZVAL__ SUBROUTINE QZVAL(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z) C INTEGER I,J,N,EN,NA,NM,NN,ISW REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) REAL 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 .QC QZVEC__ 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 REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) REAL 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 .QC RATQR__ 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 REAL D(N),E(N),E2(N),W(N),BD(N) REAL 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.0E0 IF THE SMALLEST EIGENVALUES HAVE BEEN C FOUND, AND TO 2.0E0 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 .QC REBAK__ SUBROUTINE REBAK(NM,N,B,DL,M,Z) C INTEGER I,J,K,M,N,I1,II,NM REAL B(NM,N),DL(N),Z(NM,M) REAL 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 .QC REBAKB__ SUBROUTINE REBAKB(NM,N,B,DL,M,Z) C INTEGER I,J,K,M,N,I1,II,NM REAL B(NM,N),DL(N),Z(NM,M) REAL 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 .QC REDUC__ SUBROUTINE REDUC(NM,N,A,B,DL,IERR) C INTEGER I,J,K,N,I1,J1,NM,NN,IERR REAL A(NM,N),B(NM,N),DL(N) REAL 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 .QC REDUC2__ SUBROUTINE REDUC2(NM,N,A,B,DL,IERR) C INTEGER I,J,K,N,I1,J1,NM,NN,IERR REAL A(NM,N),B(NM,N),DL(N) REAL 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 .QC RG__ SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR) C INTEGER N,NM,IS1,IS2,IERR,MATZ REAL 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 .QC RGG__ SUBROUTINE RGG(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC RS__ SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC RSB__ SUBROUTINE RSB(NM,N,MB,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,MB,NM,IERR,MATZ REAL 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 .QC RSG__ SUBROUTINE RSG(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC RSGAB__ SUBROUTINE RSGAB(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC RSGBA__ SUBROUTINE RSGBA(NM,N,A,B,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC RSM__ SUBROUTINE RSM(NM,N,A,W,M,Z,FWORK,IWORK,IERR) C INTEGER N,NM,M,IWORK(N),IERR REAL 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 .QC RSP__ SUBROUTINE RSP(NM,N,NV,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER I,J,N,NM,NV,IERR,MATZ REAL 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 .QC RST__ SUBROUTINE RST(NM,N,W,E,MATZ,Z,IERR) C INTEGER I,J,N,NM,IERR,MATZ REAL 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 .QC RT__ SUBROUTINE RT(NM,N,A,W,MATZ,Z,FV1,IERR) C INTEGER N,NM,IERR,MATZ REAL 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 .QC SVD__ 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 REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N) REAL 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 SQRT(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 .QC TINVIT__ 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 REAL D(N),E(N),E2(N),W(M),Z(NM,M), X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) REAL 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.0E0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0E0 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 SQRT(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 .QC TQL1__ SUBROUTINE TQL1(N,D,E,IERR) C INTEGER I,J,L,M,N,II,L1,L2,MML,IERR REAL D(N),E(N) REAL 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 SQRT(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 .QC TQL2__ SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR REAL D(N),E(N),Z(NM,N) REAL 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 SQRT(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 .QC TQLRAT__ SUBROUTINE TQLRAT(N,D,E2,IERR) C INTEGER I,J,L,M,N,II,L1,MML,IERR REAL D(N),E2(N) REAL 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 SQRT(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 .QC TRBAK1__ SUBROUTINE TRBAK1(NM,N,A,E,M,Z) C INTEGER I,J,K,L,M,N,NM REAL A(NM,N),E(N),Z(NM,M) REAL 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 .QC TRBAK3__ SUBROUTINE TRBAK3(NM,N,NV,A,M,Z) C INTEGER I,J,K,L,M,N,IK,IZ,NM,NV REAL A(NV),Z(NM,M) REAL 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 .QC TRED1__ SUBROUTINE TRED1(NM,N,A,D,E,E2) C INTEGER I,J,K,L,N,II,NM,JP1 REAL A(NM,N),D(N),E(N),E2(N) REAL 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 .QC TRED2__ SUBROUTINE TRED2(NM,N,A,D,E,Z) C INTEGER I,J,K,L,N,II,NM,JP1 REAL A(NM,N),D(N),E(N),Z(NM,N) REAL 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 .QC TRED3__ SUBROUTINE TRED3(N,NV,A,D,E,E2) C INTEGER I,J,K,L,N,II,IZ,JK,NV,JM1 REAL A(NV),D(N),E(N),E2(N) REAL 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 .QC TRIDIB__ 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 REAL D(N),E(N),E2(N),W(M),RV4(N),RV5(N) REAL 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 .QC TSTURM__ 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 REAL D(N),E(N),E2(N),W(MM),Z(NM,MM), X RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N) REAL 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 SQRT(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 .f .QB Library__Name(s) .br sys$userlib:EISPACK - Single precision sys$userlib:EISPACKD - Double precision .QB Documentation .i5;Documentation is available in the VAX computer room. .P Additional information may be obtained from : B.S. Garbow et al. " Matrix Eigensystem Routines - EISPACK Guide Extension ", Springer-Verlag, 1977. (Paperback, available at the Bookstore) .QB How__To__Use .list 0 .le; - Compile your program: .i5; _$ FOR MYPROG .le; - Link your object module with the appropriate library that can be found in the directory sys_$userlib: .i5;( e.g. _$#LINK#MYPROG,sys_$userlib:EISPACK/LIB ) .le; - Run the executable image: .i5;_$ RUN MYPROG .els 0 .QB Authors .c text B.S. Garbow, J.M. Boyle J.J. Dongarra, C.B. Moler Argonne National Laboratory July 1975 .end center