SUBROUTINE ZGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(1) COMPLEX*16 A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZGECO FACTORS A COMPLEX*16 MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZGECO BY ZGESL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZGECO BY ZGESL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZGECO BY ZGEDI. C TO COMPUTE INVERSE(A) , FOLLOW ZGECO BY ZGEDI. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZGEFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(1),INFO COMPLEX*16 A(LDA,1) C C ZGEFA FACTORS A COMPLEX*16 MATRIX BY GAUSSIAN ELIMINATION. C C ZGEFA IS USUALLY CALLED BY ZGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR ZGECO) = (1 + 9/N)*(TIME FOR ZGEFA) . C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT ZGESL OR ZGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN ZGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL,IZAMAX C FORTRAN DABS C C INTERNAL VARIABLES C SUBROUTINE ZGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB COMPLEX*16 A(LDA,1),B(1) C C ZGESL SOLVES THE COMPLEX*16 SYSTEM C A * X = B OR CTRANS(A) * X = B C USING THE FACTORS COMPUTED BY ZGECO OR ZGEFA. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE OUTPUT FROM ZGECO OR ZGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM ZGECO OR ZGEFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE CTRANS(A)*X = B WHERE C CTRANS(A) IS THE CONJUGATE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF ZGECO HAS SET RCOND .GT. 0.0 C OR ZGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL ZGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(1),JOB COMPLEX*16 A(LDA,1),DET(2),WORK(1) C C ZGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY ZGECO OR ZGEFA. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE OUTPUT FROM ZGECO OR ZGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM ZGECO OR ZGEFA. C C WORK COMPLEX*16(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE UNCHANGED. C C DET COMPLEX*16(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. CABS1(DET(1)) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND IF ZGECO HAS SET RCOND .GT. 0.0 OR ZGEFA HAS SET C INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL,ZSWAP C FORTRAN DABS,DCMPLX,MOD C C INTERNAL VARIABLES C SUBROUTINE ZGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) INTEGER LDA,N,ML,MU,IPVT(1) COMPLEX*16 ABD(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZGBCO FACTORS A COMPLEX*16 BAND MATRIX BY GAUSSIAN C ELIMINATION AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZGBFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZGBCO BY ZGBSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZGBCO BY ZGBSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZGBCO BY ZGBDI. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS C ML+1 THROUGH 2*ML+MU+1 OF ABD . C SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. 2*ML + MU + 1 . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C 0 .LE. ML .LT. N . C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. MU .LT. N . C MORE EFFICIENT IF ML .LE. MU . C C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C BAND STORAGE C C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT C WILL SET UP THE INPUT. C C ML = (BAND WIDTH BELOW THE DIAGONAL) C MU = (BAND WIDTH ABOVE THE DIAGONAL) C M = ML + MU + 1 C DO 20 J = 1, N C I1 = MAX0(1, J-MU) C I2 = MIN0(N, J+ML) C DO 10 I = I1, I2 C K = I - J + M C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD . C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR C ELEMENTS GENERATED DURING THE TRIANGULARIZATION. C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 . C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED. C C EXAMPLE.. IF THE ORIGINAL MATRIX IS C C 11 12 13 0 0 0 C 21 22 23 24 0 0 C 0 32 33 34 35 0 C 0 0 43 44 45 46 C 0 0 0 54 55 56 C 0 0 0 0 65 66 C C THEN N = 6, ML = 1, MU = 2, LDA .GE. 5 AND ABD SHOULD CONTAIN C C * * * + + + , * = NOT USED C * * 13 24 35 46 , + = USED FOR PIVOTING C * 12 23 34 45 56 C 11 22 33 44 55 66 C 21 32 43 54 65 * C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZGBFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,MAX0,MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZGBFA(ABD,LDA,N,ML,MU,IPVT,INFO) INTEGER LDA,N,ML,MU,IPVT(1),INFO COMPLEX*16 ABD(LDA,1) C C ZGBFA FACTORS A COMPLEX*16 BAND MATRIX BY ELIMINATION. C C ZGBFA IS USUALLY CALLED BY ZGBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS C ML+1 THROUGH 2*ML+MU+1 OF ABD . C SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. 2*ML + MU + 1 . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C 0 .LE. ML .LT. N . C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. MU .LT. N . C MORE EFFICIENT IF ML .LE. MU . C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT ZGBSL WILL DIVIDE BY ZERO IF C CALLED. USE RCOND IN ZGBCO FOR A RELIABLE C INDICATION OF SINGULARITY. C C BAND STORAGE C C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT C WILL SET UP THE INPUT. C C ML = (BAND WIDTH BELOW THE DIAGONAL) C MU = (BAND WIDTH ABOVE THE DIAGONAL) C M = ML + MU + 1 C DO 20 J = 1, N C I1 = MAX0(1, J-MU) C I2 = MIN0(N, J+ML) C DO 10 I = I1, I2 C K = I - J + M C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD . C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR C ELEMENTS GENERATED DURING THE TRIANGULARIZATION. C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 . C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL,IZAMAX C FORTRAN DABS,MAX0,MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB) INTEGER LDA,N,ML,MU,IPVT(1),JOB COMPLEX*16 ABD(LDA,1),B(1) C C ZGBSL SOLVES THE COMPLEX*16 BAND SYSTEM C A * X = B OR CTRANS(A) * X = B C USING THE FACTORS COMPUTED BY ZGBCO OR ZGBFA. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE OUTPUT FROM ZGBCO OR ZGBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM ZGBCO OR ZGBFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE CTRANS(A)*X = B , WHERE C CTRANS(A) IS THE CONJUGATE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF ZGBCO HAS SET RCOND .GT. 0.0 C OR ZGBFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL ZGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN DCONJG,MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZGBDI(ABD,LDA,N,ML,MU,IPVT,DET) INTEGER LDA,N,ML,MU,IPVT(1) COMPLEX*16 ABD(LDA,1),DET(2) C C ZGBDI COMPUTES THE DETERMINANT OF A BAND MATRIX C USING THE FACTORS COMPUTED BY ZGBCO OR ZGBFA. C IF THE INVERSE IS NEEDED, USE ZGBSL N TIMES. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE OUTPUT FROM ZGBCO OR ZGBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM ZGBCO OR ZGBFA. C C ON RETURN C C DET COMPLEX*16(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. CABS1(DET(1)) .LT. 10.0 C OR DET(1) = 0.0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C FORTRAN DABS,DCMPLX C C INTERNAL VARIABLES C SUBROUTINE ZPOCO(A,LDA,N,RCOND,Z,INFO) INTEGER LDA,N,INFO COMPLEX*16 A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZPOCO FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZPOFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZPOCO BY ZPOSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZPOCO BY ZPOSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZPOCO BY ZPODI. C TO COMPUTE INVERSE(A) , FOLLOW ZPOCO BY ZPODI. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE HERMITIAN MATRIX TO BE FACTORED. ONLY THE C DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX R SO THAT A = C CTRANS(R)*R WHERE CTRANS(R) IS THE CONJUGATE C TRANSPOSE. THE STRICT LOWER TRIANGLE IS UNALTERED. C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. IF INFO .NE. 0 , RCOND IS UNCHANGED. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C IF INFO .NE. 0 , Z IS UNCHANGED. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZPOFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZPOFA(A,LDA,N,INFO) INTEGER LDA,N,INFO COMPLEX*16 A(LDA,1) C C ZPOFA FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX. C C ZPOFA IS USUALLY CALLED BY ZPOCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR ZPOCO) = (1 + 18/N)*(TIME FOR ZPOFA) . C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE HERMITIAN MATRIX TO BE FACTORED. ONLY THE C DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX R SO THAT A = C CTRANS(R)*R WHERE CTRANS(R) IS THE CONJUGATE C TRANSPOSE. THE STRICT LOWER TRIANGLE IS UNALTERED. C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZDOTC C FORTRAN DCMPLX,DCONJG,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZPOSL(A,LDA,N,B) INTEGER LDA,N COMPLEX*16 A(LDA,1),B(1) C C ZPOSL SOLVES THE COMPLEX*16 HERMITIAN POSITIVE DEFINITE SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZPOCO OR ZPOFA. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE OUTPUT FROM ZPOCO OR ZPOFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES C SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE C ARGUMENTS. IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED C CORRECTLY AND INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZPOCO(A,LDA,N,RCOND,Z,INFO) C IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZPOSL(A,LDA,N,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C C INTERNAL VARIABLES C SUBROUTINE ZPODI(A,LDA,N,DET,JOB) INTEGER LDA,N,JOB COMPLEX*16 A(LDA,1) DOUBLE PRECISION DET(2) C C ZPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN C COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX (SEE BELOW) C USING THE FACTORS COMPUTED BY ZPOCO, ZPOFA OR ZQRDC. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE OUTPUT A FROM ZPOCO OR ZPOFA C OR THE OUTPUT X FROM ZQRDC. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C A IF ZPOCO OR ZPOFA WAS USED TO FACTOR A THEN C ZPODI PRODUCES THE UPPER HALF OF INVERSE(A) . C IF ZQRDC WAS USED TO DECOMPOSE X THEN C ZPODI PRODUCES THE UPPER HALF OF INVERSE(CTRANS(X)*X) C WHERE CTRANS(X) IS THE CONJUGATE TRANSPOSE. C ELEMENTS OF A BELOW THE DIAGONAL ARE UNCHANGED. C IF THE UNITS DIGIT OF JOB IS ZERO, A IS UNCHANGED. C C DET DOUBLE PRECISION(2) C DETERMINANT OF A OR OF CTRANS(X)*X IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DET(1) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND IF ZPOCO OR ZPOFA HAS SET INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL C FORTRAN DCONJG,MOD C C INTERNAL VARIABLES C SUBROUTINE ZPPCO(AP,N,RCOND,Z,INFO) INTEGER N,INFO COMPLEX*16 AP(1),Z(1) DOUBLE PRECISION RCOND C C ZPPCO FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C STORED IN PACKED FORM C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZPPFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZPPCO BY ZPPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZPPCO BY ZPPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZPPCO BY ZPPDI. C TO COMPUTE INVERSE(A) , FOLLOW ZPPCO BY ZPPDI. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A HERMITIAN MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C AP AN UPPER TRIANGULAR MATRIX R , STORED IN PACKED C FORM, SO THAT A = CTRANS(R)*R . C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. IF INFO .NE. 0 , RCOND IS UNCHANGED. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS SINGULAR TO WORKING PRECISION, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C IF INFO .NE. 0 , Z IS UNCHANGED. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A HERMITIAN MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZPPFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZPPFA(AP,N,INFO) INTEGER N,INFO COMPLEX*16 AP(1) C C ZPPFA FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C STORED IN PACKED FORM. C C ZPPFA IS USUALLY CALLED BY ZPPCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR ZPPCO) = (1 + 18/N)*(TIME FOR ZPPFA) . C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A HERMITIAN MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C AP AN UPPER TRIANGULAR MATRIX R , STORED IN PACKED C FORM, SO THAT A = CTRANS(R)*R . C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K IF THE LEADING MINOR OF ORDER K IS NOT C POSITIVE DEFINITE. C C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A HERMITIAN MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZDOTC C FORTRAN DCMPLX,DCONJG,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZPPSL(AP,N,B) INTEGER N COMPLEX*16 AP(1),B(1) C C ZPPSL SOLVES THE COMPLEX*16 HERMITIAN POSITIVE DEFINITE SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZPPCO OR ZPPFA. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE OUTPUT FROM ZPPCO OR ZPPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES C SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE C ARGUMENTS. IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED C CORRECTLY AND INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZPPCO(AP,N,RCOND,Z,INFO) C IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZPPSL(AP,N,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C C INTERNAL VARIABLES C SUBROUTINE ZPPDI(AP,N,DET,JOB) INTEGER N,JOB COMPLEX*16 AP(1) DOUBLE PRECISION DET(2) C C ZPPDI COMPUTES THE DETERMINANT AND INVERSE C OF A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C USING THE FACTORS COMPUTED BY ZPPCO OR ZPPFA . C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE OUTPUT FROM ZPPCO OR ZPPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C AP THE UPPER TRIANGULAR HALF OF THE INVERSE . C THE STRICT LOWER TRIANGLE IS UNALTERED. C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DET(1) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND IF ZPOCO OR ZPOFA HAS SET INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL C FORTRAN DCONJG,MOD C C INTERNAL VARIABLES C SUBROUTINE ZPBCO(ABD,LDA,N,M,RCOND,Z,INFO) INTEGER LDA,N,M,INFO COMPLEX*16 ABD(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZPBCO FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C STORED IN BAND FORM AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZPBFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZPBCO BY ZPBSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZPBCO BY ZPBSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZPBCO BY ZPBDI. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE MATRIX TO BE FACTORED. THE COLUMNS OF THE UPPER C TRIANGLE ARE STORED IN THE COLUMNS OF ABD AND THE C DIAGONALS OF THE UPPER TRIANGLE ARE STORED IN THE C ROWS OF ABD . SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. M + 1 . C C N INTEGER C THE ORDER OF THE MATRIX A . C C M INTEGER C THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. M .LT. N . C C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX R , STORED IN BAND C FORM, SO THAT A = CTRANS(R)*R . C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. IF INFO .NE. 0 , RCOND IS UNCHANGED. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS SINGULAR TO WORKING PRECISION, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C IF INFO .NE. 0 , Z IS UNCHANGED. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. C C BAND STORAGE C C IF A IS A HERMITIAN POSITIVE DEFINITE BAND MATRIX, C THE FOLLOWING PROGRAM SEGMENT WILL SET UP THE INPUT. C C M = (BAND WIDTH ABOVE DIAGONAL) C DO 20 J = 1, N C I1 = MAX0(1, J-M) C DO 10 I = I1, J C K = I-J+M+1 C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C THIS USES M + 1 ROWS OF A , EXCEPT FOR THE M BY M C UPPER LEFT TRIANGLE, WHICH IS IGNORED. C C EXAMPLE.. IF THE ORIGINAL MATRIX IS C C 11 12 13 0 0 0 C 12 22 23 24 0 0 C 13 23 33 34 35 0 C 0 24 34 44 45 46 C 0 0 35 45 55 56 C 0 0 0 46 56 66 C C THEN N = 6 , M = 2 AND ABD SHOULD CONTAIN C C * * 13 24 35 46 C * 12 23 34 45 56 C 11 22 33 44 55 66 C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZPBFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,MAX0,MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZPBFA(ABD,LDA,N,M,INFO) INTEGER LDA,N,M,INFO COMPLEX*16 ABD(LDA,1) C C ZPBFA FACTORS A COMPLEX*16 HERMITIAN POSITIVE DEFINITE MATRIX C STORED IN BAND FORM. C C ZPBFA IS USUALLY CALLED BY ZPBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE MATRIX TO BE FACTORED. THE COLUMNS OF THE UPPER C TRIANGLE ARE STORED IN THE COLUMNS OF ABD AND THE C DIAGONALS OF THE UPPER TRIANGLE ARE STORED IN THE C ROWS OF ABD . SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. M + 1 . C C N INTEGER C THE ORDER OF THE MATRIX A . C C M INTEGER C THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. M .LT. N . C C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX R , STORED IN BAND C FORM, SO THAT A = CTRANS(R)*R . C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K IF THE LEADING MINOR OF ORDER K IS NOT C POSITIVE DEFINITE. C C BAND STORAGE C C IF A IS A HERMITIAN POSITIVE DEFINITE BAND MATRIX, C THE FOLLOWING PROGRAM SEGMENT WILL SET UP THE INPUT. C C M = (BAND WIDTH ABOVE DIAGONAL) C DO 20 J = 1, N C I1 = MAX0(1, J-M) C DO 10 I = I1, J C K = I-J+M+1 C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZDOTC C FORTRAN DCMPLX,DCONJG,MAX0,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZPBSL(ABD,LDA,N,M,B) INTEGER LDA,N,M COMPLEX*16 ABD(LDA,1),B(1) C C ZPBSL SOLVES THE COMPLEX*16 HERMITIAN POSITIVE DEFINITE BAND C SYSTEM A*X = B C USING THE FACTORS COMPUTED BY ZPBCO OR ZPBFA. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE OUTPUT FROM ZPBCO OR ZPBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE MATRIX A . C C M INTEGER C THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES C SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE C ARGUMENTS. IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED C CORRECTLY AND INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZPBCO(ABD,LDA,N,RCOND,Z,INFO) C IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZPBSL(ABD,LDA,N,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZPBDI(ABD,LDA,N,M,DET) INTEGER LDA,N,M COMPLEX*16 ABD(LDA,1) DOUBLE PRECISION DET(2) C C ZPBDI COMPUTES THE DETERMINANT C OF A COMPLEX*16 HERMITIAN POSITIVE DEFINITE BAND MATRIX C USING THE FACTORS COMPUTED BY ZPBCO OR ZPBFA. C IF THE INVERSE IS NEEDED, USE ZPBSL N TIMES. C C ON ENTRY C C ABD COMPLEX*16(LDA, N) C THE OUTPUT FROM ZPBCO OR ZPBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE MATRIX A . C C M INTEGER C THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C ON RETURN C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX IN THE FORM C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DET(1) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C C INTERNAL VARIABLES C SUBROUTINE ZSICO(A,LDA,N,KPVT,RCOND,Z) INTEGER LDA,N,KPVT(1) COMPLEX*16 A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZSICO FACTORS A COMPLEX*16 SYMMETRIC MATRIX BY ELIMINATION WITH C SYMMETRIC PIVOTING AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZSIFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZSICO BY ZSISL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZSICO BY ZSISL. C TO COMPUTE INVERSE(A) , FOLLOW ZSICO BY ZSIDI. C TO COMPUTE DETERMINANT(A) , FOLLOW ZSICO BY ZSIDI. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE SYMMETRIC MATRIX TO BE FACTORED. C ONLY THE DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = U*D*TRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , TRANS(U) IS THE C TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZSIFA C BLAS ZAXPY,ZDOTU,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,IABS C C INTERNAL VARIABLES C SUBROUTINE ZSIFA(A,LDA,N,KPVT,INFO) INTEGER LDA,N,KPVT(1),INFO COMPLEX*16 A(LDA,1) C C ZSIFA FACTORS A COMPLEX*16 SYMMETRIC MATRIX BY ELIMINATION C WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW ZSIFA BY ZSISL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZSIFA BY ZSISL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZSIFA BY ZSIDI. C TO COMPUTE INVERSE(A) , FOLLOW ZSIFA BY ZSIDI. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE SYMMETRIC MATRIX TO BE FACTORED. C ONLY THE DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = U*D*TRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , TRANS(U) IS THE C TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF THE K-TH PIVOT BLOCK IS SINGULAR. THIS IS C NOT AN ERROR CONDITION FOR THIS SUBROUTINE, C BUT IT DOES INDICATE THAT ZSISL OR ZSIDI MAY C DIVIDE BY ZERO IF CALLED. C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSWAP,IZAMAX C FORTRAN DABS,DMAX1,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZSISL(A,LDA,N,KPVT,B) INTEGER LDA,N,KPVT(1) COMPLEX*16 A(LDA,1),B(1) C C ZSISL SOLVES THE COMPLEX*16 SYMMETRIC SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZSIFA. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE OUTPUT FROM ZSIFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZSIFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF ZSICO HAS SET RCOND .EQ. 0.0 C OR ZSIFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZSIFA(A,LDA,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZSISL(A,LDA,N,KPVT,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTU C FORTRAN IABS C C INTERNAL VARIABLES. C SUBROUTINE ZSIDI(A,LDA,N,KPVT,DET,WORK,JOB) INTEGER LDA,N,JOB COMPLEX*16 A(LDA,1),DET(2),WORK(1) INTEGER KPVT(1) C C ZSIDI COMPUTES THE DETERMINANT AND INVERSE C OF A COMPLEX*16 SYMMETRIC MATRIX USING THE FACTORS FROM ZSIFA. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE OUTPUT FROM ZSIFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A. C C N INTEGER C THE ORDER OF THE MATRIX A. C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZSIFA. C C WORK COMPLEX*16(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C JOB HAS THE DECIMAL EXPANSION AB WHERE C IF B .NE. 0, THE INVERSE IS COMPUTED, C IF A .NE. 0, THE DETERMINANT IS COMPUTED, C C FOR EXAMPLE, JOB = 11 GIVES BOTH. C C ON RETURN C C VARIABLES NOT REQUESTED BY JOB ARE NOT USED. C C A CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF C THE ORIGINAL MATRIX. THE STRICT LOWER TRIANGLE C IS NEVER REFERENCED. C C DET COMPLEX*16(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) = 0.0. C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF THE INVERSE IS REQUESTED C AND ZSICO HAS SET RCOND .EQ. 0.0 C OR ZSIFA HAS SET INFO .NE. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZCOPY,ZDOTU,ZSWAP C FORTRAN DABS,DCMPLX,IABS,MOD C C INTERNAL VARIABLES. C SUBROUTINE ZSPCO(AP,N,KPVT,RCOND,Z) INTEGER N,KPVT(1) COMPLEX*16 AP(1),Z(1) DOUBLE PRECISION RCOND C C ZSPCO FACTORS A COMPLEX*16 SYMMETRIC MATRIX STORED IN PACKED C FORM BY ELIMINATION WITH SYMMETRIC PIVOTING AND ESTIMATES C THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZSPFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZSPCO BY ZSPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZSPCO BY ZSPSL. C TO COMPUTE INVERSE(A) , FOLLOW ZSPCO BY ZSPDI. C TO COMPUTE DETERMINANT(A) , FOLLOW ZSPCO BY ZSPDI. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A SYMMETRIC MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C AP A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT STORED IN PACKED FORM. C THE FACTORIZATION CAN BE WRITTEN A = U*D*TRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , TRANS(U) IS THE C TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A SYMMETRIC MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZSPFA C BLAS ZAXPY,ZDOTU,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,IABS C C INTERNAL VARIABLES C SUBROUTINE ZSPFA(AP,N,KPVT,INFO) INTEGER N,KPVT(1),INFO COMPLEX*16 AP(1) C C ZSPFA FACTORS A COMPLEX*16 SYMMETRIC MATRIX STORED IN C PACKED FORM BY ELIMINATION WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW ZSPFA BY ZSPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZSPFA BY ZSPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZSPFA BY ZSPDI. C TO COMPUTE INVERSE(A) , FOLLOW ZSPFA BY ZSPDI. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A SYMMETRIC MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C AP A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT STORED IN PACKED FORM. C THE FACTORIZATION CAN BE WRITTEN A = U*D*TRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , TRANS(U) IS THE C TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF THE K-TH PIVOT BLOCK IS SINGULAR. THIS IS C NOT AN ERROR CONDITION FOR THIS SUBROUTINE, C BUT IT DOES INDICATE THAT ZSPSL OR ZSPDI MAY C DIVIDE BY ZERO IF CALLED. C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A SYMMETRIC MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSWAP,IZAMAX C FORTRAN DABS,DMAX1,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZSPSL(AP,N,KPVT,B) INTEGER N,KPVT(1) COMPLEX*16 AP(1),B(1) C C ZSISL SOLVES THE COMPLEX*16 SYMMETRIC SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZSPFA. C C ON ENTRY C C AP COMPLEX*16(N*(N+1)/2) C THE OUTPUT FROM ZSPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZSPFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF ZSPCO HAS SET RCOND .EQ. 0.0 C OR ZSPFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZSPFA(AP,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZSPSL(AP,N,KPVT,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTU C FORTRAN IABS C C INTERNAL VARIABLES. C SUBROUTINE ZSPDI(AP,N,KPVT,DET,WORK,JOB) INTEGER N,JOB COMPLEX*16 AP(1),WORK(1),DET(2) INTEGER KPVT(1) C C ZSPDI COMPUTES THE DETERMINANT AND INVERSE C OF A COMPLEX*16 SYMMETRIC MATRIX USING THE FACTORS FROM ZSPFA, C WHERE THE MATRIX IS STORED IN PACKED FORM. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE OUTPUT FROM ZSPFA. C C N INTEGER C THE ORDER OF THE MATRIX A. C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZSPFA. C C WORK COMPLEX*16(N) C WORK VECTOR. CONTENTS IGNORED. C C JOB INTEGER C JOB HAS THE DECIMAL EXPANSION AB WHERE C IF B .NE. 0, THE INVERSE IS COMPUTED, C IF A .NE. 0, THE DETERMINANT IS COMPUTED, C C FOR EXAMPLE, JOB = 11 GIVES BOTH. C C ON RETURN C C VARIABLES NOT REQUESTED BY JOB ARE NOT USED. C C AP CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF C THE ORIGINAL MATRIX, STORED IN PACKED FORM. C THE COLUMNS OF THE UPPER TRIANGLE ARE STORED C SEQUENTIALLY IN A ONE-DIMENSIONAL ARRAY. C C DET COMPLEX*16(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) = 0.0. C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INVERSE IS REQUESTED C AND ZSPCO HAS SET RCOND .EQ. 0.0 C OR ZSPFA HAS SET INFO .NE. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZCOPY,ZDOTU,ZSWAP C FORTRAN DABS,DCMPLX,IABS,MOD C C INTERNAL VARIABLES. C SUBROUTINE ZHICO(A,LDA,N,KPVT,RCOND,Z) INTEGER LDA,N,KPVT(1) COMPLEX*16 A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C ZHICO FACTORS A COMPLEX*16 HERMITIAN MATRIX BY ELIMINATION WITH C SYMMETRIC PIVOTING AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZHIFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZHICO BY ZHISL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZHICO BY ZHISL. C TO COMPUTE INVERSE(A) , FOLLOW ZHICO BY ZHIDI. C TO COMPUTE DETERMINANT(A) , FOLLOW ZHICO BY ZHIDI. C TO COMPUTE INERTIA(A), FOLLOW ZHICO BY ZHIDI. C C ON ENTRY C C A COMPLEX*16(LDA, N) C THE HERMITIAN MATRIX TO BE FACTORED. C ONLY THE DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C OUTPUT C C A A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = U*D*CTRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , CTRANS(U) IS THE C CONJUGATE TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZHIFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,IABS C C INTERNAL VARIABLES C SUBROUTINE ZHIFA(A,LDA,N,KPVT,INFO) INTEGER LDA,N,KPVT(1),INFO COMPLEX*16 A(LDA,1) C C ZHIFA FACTORS A COMPLEX*16 HERMITIAN MATRIX BY ELIMINATION C WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW ZHIFA BY ZHISL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZHIFA BY ZHISL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZHIFA BY ZHIDI. C TO COMPUTE INERTIA(A) , FOLLOW ZHIFA BY ZHIDI. C TO COMPUTE INVERSE(A) , FOLLOW ZHIFA BY ZHIDI. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE HERMITIAN MATRIX TO BE FACTORED. C ONLY THE DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = U*D*CTRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , CTRANS(U) IS THE C CONJUGATE TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF THE K-TH PIVOT BLOCK IS SINGULAR. THIS IS C NOT AN ERROR CONDITION FOR THIS SUBROUTINE, C BUT IT DOES INDICATE THAT ZHISL OR ZHIDI MAY C DIVIDE BY ZERO IF CALLED. C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSWAP,IZAMAX C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZHISL(A,LDA,N,KPVT,B) INTEGER LDA,N,KPVT(1) COMPLEX*16 A(LDA,1),B(1) C C ZHISL SOLVES THE COMPLEX*16 HERMITIAN SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZHIFA. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE OUTPUT FROM ZHIFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZHIFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF ZHICO HAS SET RCOND .EQ. 0.0 C OR ZHIFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZHIFA(A,LDA,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZHISL(A,LDA,N,KPVT,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN DCONJG,IABS C C INTERNAL VARIABLES. C SUBROUTINE ZHIDI(A,LDA,N,KPVT,DET,INERT,WORK,JOB) INTEGER LDA,N,JOB COMPLEX*16 A(LDA,1),WORK(1) DOUBLE PRECISION DET(2) INTEGER KPVT(1),INERT(3) C C ZHIDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE C OF A COMPLEX*16 HERMITIAN MATRIX USING THE FACTORS FROM ZHIFA. C C ON ENTRY C C A COMPLEX*16(LDA,N) C THE OUTPUT FROM ZHIFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A. C C N INTEGER C THE ORDER OF THE MATRIX A. C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZHIFA. C C WORK COMPLEX*16(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C JOB HAS THE DECIMAL EXPANSION ABC WHERE C IF C .NE. 0, THE INVERSE IS COMPUTED, C IF B .NE. 0, THE DETERMINANT IS COMPUTED, C IF A .NE. 0, THE INERTIA IS COMPUTED. C C FOR EXAMPLE, JOB = 111 GIVES ALL THREE. C C ON RETURN C C VARIABLES NOT REQUESTED BY JOB ARE NOT USED. C C A CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF C THE ORIGINAL MATRIX. THE STRICT LOWER TRIANGLE C IS NEVER REFERENCED. C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) = 0.0. C C INERT INTEGER(3) C THE INERTIA OF THE ORIGINAL MATRIX. C INERT(1) = NUMBER OF POSITIVE EIGENVALUES. C INERT(2) = NUMBER OF NEGATIVE EIGENVALUES. C INERT(3) = NUMBER OF ZERO EIGENVALUES. C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF THE INVERSE IS REQUESTED C AND ZHICO HAS SET RCOND .EQ. 0.0 C OR ZHIFA HAS SET INFO .NE. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZCOPY,ZDOTC,ZSWAP C FORTRAN DABS,CDABS,DCMPLX,DCONJG,IABS,MOD C C INTERNAL VARIABLES. C SUBROUTINE ZHPCO(AP,N,KPVT,RCOND,Z) INTEGER N,KPVT(1) COMPLEX*16 AP(1),Z(1) DOUBLE PRECISION RCOND C C ZHPCO FACTORS A COMPLEX*16 HERMITIAN MATRIX STORED IN PACKED C FORM BY ELIMINATION WITH SYMMETRIC PIVOTING AND ESTIMATES C THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, ZHPFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW ZHPCO BY ZHPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZHPCO BY ZHPSL. C TO COMPUTE INVERSE(A) , FOLLOW ZHPCO BY ZHPDI. C TO COMPUTE DETERMINANT(A) , FOLLOW ZHPCO BY ZHPDI. C TO COMPUTE INERTIA(A), FOLLOW ZHPCO BY ZHPDI. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A HERMITIAN MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C OUTPUT C C AP A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT STORED IN PACKED FORM. C THE FACTORIZATION CAN BE WRITTEN A = U*D*CTRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , CTRANS(U) IS THE C CONJUGATE TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A HERMITIAN MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK ZHPFA C BLAS ZAXPY,ZDOTC,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,IABS C C INTERNAL VARIABLES C SUBROUTINE ZHPFA(AP,N,KPVT,INFO) INTEGER N,KPVT(1),INFO COMPLEX*16 AP(1) C C ZHPFA FACTORS A COMPLEX*16 HERMITIAN MATRIX STORED IN C PACKED FORM BY ELIMINATION WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW ZHPFA BY ZHPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW ZHPFA BY ZHPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW ZHPFA BY ZHPDI. C TO COMPUTE INERTIA(A) , FOLLOW ZHPFA BY ZHPDI. C TO COMPUTE INVERSE(A) , FOLLOW ZHPFA BY ZHPDI. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE PACKED FORM OF A HERMITIAN MATRIX A . THE C COLUMNS OF THE UPPER TRIANGLE ARE STORED SEQUENTIALLY C IN A ONE-DIMENSIONAL ARRAY OF LENGTH N*(N+1)/2 . C SEE COMMENTS BELOW FOR DETAILS. C C N INTEGER C THE ORDER OF THE MATRIX A . C C OUTPUT C C AP A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT STORED IN PACKED FORM. C THE FACTORIZATION CAN BE WRITTEN A = U*D*CTRANS(U) C WHERE U IS A PRODUCT OF PERMUTATION AND UNIT C UPPER TRIANGULAR MATRICES , CTRANS(U) IS THE C CONJUGATE TRANSPOSE OF U , AND D IS BLOCK DIAGONAL C WITH 1 BY 1 AND 2 BY 2 BLOCKS. C C KPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF THE K-TH PIVOT BLOCK IS SINGULAR. THIS IS C NOT AN ERROR CONDITION FOR THIS SUBROUTINE, C BUT IT DOES INDICATE THAT ZHPSL OR ZHPDI MAY C DIVIDE BY ZERO IF CALLED. C C PACKED STORAGE C C THE FOLLOWING PROGRAM SEGMENT WILL PACK THE UPPER C TRIANGLE OF A HERMITIAN MATRIX. C C K = 0 C DO 20 J = 1, N C DO 10 I = 1, J C K = K + 1 C AP(K) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSWAP,IZAMAX C FORTRAN DABS,DMAX1,DCMPLX,DCONJG,DSQRT C C INTERNAL VARIABLES C SUBROUTINE ZHPSL(AP,N,KPVT,B) INTEGER N,KPVT(1) COMPLEX*16 AP(1),B(1) C C ZHISL SOLVES THE COMPLEX*16 HERMITIAN SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY ZHPFA. C C ON ENTRY C C AP COMPLEX*16(N*(N+1)/2) C THE OUTPUT FROM ZHPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZHPFA. C C B COMPLEX*16(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO MAY OCCUR IF ZHPCO HAS SET RCOND .EQ. 0.0 C OR ZHPFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL ZHPFA(AP,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL ZHPSL(AP,N,KPVT,C(1,J)) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN DCONJG,IABS C C INTERNAL VARIABLES. C SUBROUTINE ZHPDI(AP,N,KPVT,DET,INERT,WORK,JOB) INTEGER N,JOB COMPLEX*16 AP(1),WORK(1) DOUBLE PRECISION DET(2) INTEGER KPVT(1),INERT(3) C C ZHPDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE C OF A COMPLEX*16 HERMITIAN MATRIX USING THE FACTORS FROM ZHPFA, C WHERE THE MATRIX IS STORED IN PACKED FORM. C C ON ENTRY C C AP COMPLEX*16 (N*(N+1)/2) C THE OUTPUT FROM ZHPFA. C C N INTEGER C THE ORDER OF THE MATRIX A. C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM ZHPFA. C C WORK COMPLEX*16(N) C WORK VECTOR. CONTENTS IGNORED. C C JOB INTEGER C JOB HAS THE DECIMAL EXPANSION ABC WHERE C IF C .NE. 0, THE INVERSE IS COMPUTED, C IF B .NE. 0, THE DETERMINANT IS COMPUTED, C IF A .NE. 0, THE INERTIA IS COMPUTED. C C FOR EXAMPLE, JOB = 111 GIVES ALL THREE. C C ON RETURN C C VARIABLES NOT REQUESTED BY JOB ARE NOT USED. C C AP CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF C THE ORIGINAL MATRIX, STORED IN PACKED FORM. C THE COLUMNS OF THE UPPER TRIANGLE ARE STORED C SEQUENTIALLY IN A ONE-DIMENSIONAL ARRAY. C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) = 0.0. C C INERT INTEGER(3) C THE INERTIA OF THE ORIGINAL MATRIX. C INERT(1) = NUMBER OF POSITIVE EIGENVALUES. C INERT(2) = NUMBER OF NEGATIVE EIGENVALUES. C INERT(3) = NUMBER OF ZERO EIGENVALUES. C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INVERSE IS REQUESTED C AND ZHPCO HAS SET RCOND .EQ. 0.0 C OR ZHPFA HAS SET INFO .NE. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C JAMES BUNCH, UNIV. CALIF. SAN DIEGO, ARGONNE NAT. LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZCOPY,ZDOTC,ZSWAP C FORTRAN DABS,CDABS,DCMPLX,DCONJG,IABS,MOD C C INTERNAL VARIABLES. C SUBROUTINE ZTRCO(T,LDT,N,RCOND,Z,JOB) INTEGER LDT,N,JOB COMPLEX*16 T(LDT,1),Z(1) DOUBLE PRECISION RCOND C C ZTRCO ESTIMATES THE CONDITION OF A COMPLEX*16 TRIANGULAR MATRIX. C C ON ENTRY C C T COMPLEX*16(LDT,N) C T CONTAINS THE TRIANGULAR MATRIX. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C C N INTEGER C N IS THE ORDER OF THE SYSTEM. C C JOB INTEGER C = 0 T IS LOWER TRIANGULAR. C = NONZERO T IS UPPER TRIANGULAR. C C ON RETURN C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF T . C FOR THE SYSTEM T*X = B , RELATIVE PERTURBATIONS C IN T AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN T MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX*16(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF T IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDSCAL,DZASUM C FORTRAN DABS,DMAX1,DCMPLX,DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZTRSL(T,LDT,N,B,JOB,INFO) INTEGER LDT,N,JOB,INFO COMPLEX*16 T(LDT,1),B(1) C C C ZTRSL SOLVES SYSTEMS OF THE FORM C C T * X = B C OR C CTRANS(T) * X = B C C WHERE T IS A TRIANGULAR MATRIX OF ORDER N. HERE CTRANS(T) C DENOTES THE CONJUGATE TRANSPOSE OF THE MATRIX T. C C ON ENTRY C C T COMPLEX*16(LDT,N) C T CONTAINS THE MATRIX OF THE SYSTEM. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C C N INTEGER C N IS THE ORDER OF THE SYSTEM. C C B COMPLEX*16(N). C B CONTAINS THE RIGHT HAND SIDE OF THE SYSTEM. C C JOB INTEGER C JOB SPECIFIES WHAT KIND OF SYSTEM IS TO BE SOLVED. C IF JOB IS C C 00 SOLVE T*X=B, T LOWER TRIANGULAR, C 01 SOLVE T*X=B, T UPPER TRIANGULAR, C 10 SOLVE CTRANS(T)*X=B, T LOWER TRIANGULAR, C 11 SOLVE CTRANS(T)*X=B, T UPPER TRIANGULAR. C C ON RETURN C C B B CONTAINS THE SOLUTION, IF INFO .EQ. 0. C OTHERWISE B IS UNALTERED. C C INFO INTEGER C INFO CONTAINS ZERO IF THE SYSTEM IS NONSINGULAR. C OTHERWISE INFO CONTAINS THE INDEX OF C THE FIRST ZERO DIAGONAL ELEMENT OF T. C C LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZDOTC C FORTRAN DABS,DCONJG,MOD C C INTERNAL VARIABLES C SUBROUTINE ZTRDI(T,LDT,N,DET,JOB,INFO) INTEGER LDT,N,JOB,INFO COMPLEX*16 T(LDT,1),DET(2) C C ZTRDI COMPUTES THE DETERMINANT AND INVERSE OF A COMPLEX*16 C TRIANGULAR MATRIX. C C ON ENTRY C C T COMPLEX*16(LDT,N) C T CONTAINS THE TRIANGULAR MATRIX. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C C N INTEGER C N IS THE ORDER OF THE SYSTEM. C C JOB INTEGER C = 010 NO DET, INVERSE OF LOWER TRIANGULAR. C = 011 NO DET, INVERSE OF UPPER TRIANGULAR. C = 100 DET, NO INVERSE. C = 110 DET, INVERSE OF LOWER TRIANGULAR. C = 111 DET, INVERSE OF UPPER TRIANGULAR. C C ON RETURN C C T INVERSE OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE UNCHANGED. C C DET COMPLEX*16(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. CABS1(DET(1)) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C INFO INTEGER C INFO CONTAINS ZERO IF THE SYSTEM IS NONSINGULAR C AND THE INVERSE IS REQUESTED. C OTHERWISE INFO CONTAINS THE INDEX OF C A ZERO DIAGONAL ELEMENT OF T. C C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS ZAXPY,ZSCAL C FORTRAN DABS,DCMPLX,MOD C C INTERNAL VARIABLES C SUBROUTINE ZGTSL(N,C,D,E,B,INFO) INTEGER N,INFO COMPLEX*16 C(1),D(1),E(1),B(1) C C ZGTSL GIVEN A GENERAL TRIDIAGONAL MATRIX AND A RIGHT HAND C SIDE WILL FIND THE SOLUTION. C C ON ENTRY C C N INTEGER C IS THE ORDER OF THE TRIDIAGONAL MATRIX. C C C COMPLEX*16(N) C IS THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX. C C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL. C ON OUTPUT C IS DESTROYED. C C D COMPLEX*16(N) C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX. C ON OUTPUT D IS DESTROYED. C C E COMPLEX*16(N) C IS THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX. C E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPERDIAGONAL. C ON OUTPUT E IS DESTROYED. C C B COMPLEX*16(N) C IS THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B IS THE SOLUTION VECTOR. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF THE K-TH ELEMENT OF THE DIAGONAL BECOMES C EXACTLY ZERO. THE SUBROUTINE RETURNS WHEN C THIS IS DETECTED. C C LINPACK. THIS VERSION DATED 08/14/78 . C JACK DONGARRA, ARGONNE NATIONAL LABORATORY. C C NO EXTERNALS C FORTRAN DABS C C INTERNAL VARIABLES C SUBROUTINE ZPTSL(N,D,E,B) INTEGER N COMPLEX*16 D(1),E(1),B(1) C C ZPTSL GIVEN A POSITIVE DEFINITE TRIDIAGONAL MATRIX AND A RIGHT C HAND SIDE WILL FIND THE SOLUTION. C C ON ENTRY C C N INTEGER C IS THE ORDER OF THE TRIDIAGONAL MATRIX. C C D COMPLEX*16(N) C IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX. C ON OUTPUT D IS DESTROYED. C C E COMPLEX*16(N) C IS THE OFFDIAGONAL OF THE TRIDIAGONAL MATRIX. C E(1) THROUGH E(N-1) SHOULD CONTAIN THE C OFFDIAGONAL. C C B COMPLEX*16(N) C IS THE RIGHT HAND SIDE VECTOR. C C ON RETURN C C B CONTAINS THE SOULTION. C C LINPACK. THIS VERSION DATED 08/14/78 . C JACK DONGARRA, ARGONNE NATIONAL LABORATORY. C C NO EXTERNALS C FORTRAN DCONJG,MOD C C INTERNAL VARIABLES C SUBROUTINE ZQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) INTEGER LDX,N,P,JOB INTEGER JPVT(1) COMPLEX*16 X(LDX,1),QRAUX(1),WORK(1) C C ZQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE C PERFORMED AT THE USERS OPTION. C C ON ENTRY C C X COMPLEX*16(LDX,P), WHERE LDX .GE. N. C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE C COMPUTED. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C JPVT INTEGER(P). C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE C VALUE OF JPVT(K). C C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C COLUMN. C C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN. C C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN. C C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST C REDUCED NORM. JPVT IS NOT REFERENCED IF C JOB .EQ. 0. C C WORK COMPLEX*16(P). C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF C JOB .EQ. 0. C C JOB INTEGER. C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING. C IF JOB .EQ. 0, NO PIVOTING IS DONE. C IF JOB .NE. 0, PIVOTING IS DONE. C C ON RETURN C C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER C TRIANGULAR MATRIX R OF THE QR FACTORIZATION. C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM C WHICH THE UNITARY PART OF THE DECOMPOSITION C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT C OF THE ORIGINAL MATRIX X BUT THAT OF X C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT. C C QRAUX COMPLEX*16(P). C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER C THE UNITARY PART OF THE DECOMPOSITION. C C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS ZAXPY,ZDOTC,ZSCAL,ZSWAP,DZNRM2 C FORTRAN DABS,DMAX1,CDABS,DCMPLX,CDSQRT,MIN0 C C INTERNAL VARIABLES C SUBROUTINE ZQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) INTEGER LDX,N,K,JOB,INFO COMPLEX*16 X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1) C C ZQRSL APPLIES THE OUTPUT OF ZQRDC TO COMPUTE COORDINATE C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS. C FOR K .LE. MIN(N,P), LET XK BE THE MATRIX C C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C C FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL C N X P MATRIX X THAT WAS INPUT TO ZQRDC (IF NO PIVOTING WAS C DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR C ORIGINAL ORDER). ZQRDC PRODUCES A FACTORED UNITARY MATRIX Q C AND AN UPPER TRIANGULAR MATRIX R SUCH THAT C C XK = Q * (R) C (0) C C THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS C X AND QRAUX. C C ON ENTRY C C X COMPLEX*16(LDX,P). C X CONTAINS THE OUTPUT OF ZQRDC. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX XK. IT MUST C HAVE THE SAME VALUE AS N IN ZQRDC. C C K INTEGER. C K IS THE NUMBER OF COLUMNS OF THE MATRIX XK. K C MUST NNOT BE GREATER THAN MIN(N,P), WHERE P IS THE C SAME AS IN THE CALLING SEQUENCE TO ZQRDC. C C QRAUX COMPLEX*16(P). C QRAUX CONTAINS THE AUXILIARY OUTPUT FROM ZQRDC. C C Y COMPLEX*16(N) C Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED C BY ZQRSL. C C JOB INTEGER. C JOB SPECIFIES WHAT IS TO BE COMPUTED. JOB HAS C THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING C MEANING. C C IF A.NE.0, COMPUTE QY. C IF B,C,D, OR E .NE. 0, COMPUTE QTY. C IF C.NE.0, COMPUTE B. C IF D.NE.0, COMPUTE RSD. C IF E.NE.0, COMPUTE XB. C C NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB C AUTOMATICALLY TRIGGERS THE COMPUTATION OF QTY, FOR C WHICH AN ARRAY MUST BE PROVIDED IN THE CALLING C SEQUENCE. C C ON RETURN C C QY COMPLEX*16(N). C QY CONNTAINS Q*Y, IF ITS COMPUTATION HAS BEEN C REQUESTED. C C QTY COMPLEX*16(N). C QTY CONTAINS CTRANS(Q)*Y, IF ITS COMPUTATION HAS C BEEN REQUESTED. HERE CTRANS(Q) IS THE CONJUGATE C TRANSPOSE OF THE MATRIX Q. C C B COMPLEX*16(K) C B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM C C MINIMIZE NORM2(Y - XK*B), C C IF ITS COMPUTATION HAS BEEN REQUESTED. (NOTE THAT C IF PIVOTING WAS REQUESTED IN ZQRDC, THE J-TH C COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J) C OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO ZQRDC.) C C RSD COMPLEX*16(N). C RSD CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. RSD IS C ALSO THE ORTHOGONAL PROJECTION OF Y ONTO THE C ORTHOGONAL COMPLEMENT OF THE COLUMN SPACE OF XK. C C XB COMPLEX*16(N). C XB CONTAINS THE LEAST SQUARES APPROXIMATION XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. XB IS ALSO C THE ORTHOGONAL PROJECTION OF Y ONTO THE COLUMN SPACE C OF X. C C INFO INTEGER. C INFO IS ZERO UNLESS THE COMPUTATION OF B HAS C BEEN REQUESTED AND R IS EXACTLY SINGULAR. IN C THIS CASE, INFO IS THE INDEX OF THE FIRST ZERO C DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED. C C THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED C IF THEIR COMPUTATION IS NOT REQUESTED AND IN THIS CASE C CAN BE REPLACED BY DUMMY VARIABLES IN THE CALLING PROGRAM. C TO SAVE STORAGE, THE USER MAY IN SOME CASES USE THE SAME C ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE. A C FREQUENTLY OCCURING EXAMPLE IS WHEN ONE WISHES TO COMPUTE C ANY OF B, RSD, OR XB AND DOES NOT NEED Y OR QTY. IN THIS C CASE ONE MAY IDENTIFY Y, QTY, AND ONE OF B, RSD, OR XB, WHILE C PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE THAT IS TO BE C COMPUTED. THUS THE CALLING SEQUENCE C C CALL ZQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) C C WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD C OVERWRITING Y. MORE GENERALLY, EACH ITEM IN THE FOLLOWING C LIST CONTAINS GROUPS OF PERMISSIBLE IDENTIFICATIONS FOR C A SINGLE CALLINNG SEQUENCE. C C 1. (Y,QTY,B) (RSD) (XB) (QY) C C 2. (Y,QTY,RSD) (B) (XB) (QY) C C 3. (Y,QTY,XB) (B) (RSD) (QY) C C 4. (Y,QY) (QTY,B) (RSD) (XB) C C 5. (Y,QY) (QTY,RSD) (B) (XB) C C 6. (Y,QY) (QTY,XB) (B) (RSD) C C IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO C THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZQRSL USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS ZAXPY,ZCOPY,ZDOTC C FORTRAN DABS,MIN0,MOD C C INTERNAL VARIABLES C SUBROUTINE ZCHDC(A,LDA,P,WORK,JPVT,JOB,INFO) INTEGER LDA,P,JPVT(1),JOB,INFO COMPLEX*16 A(LDA,1),WORK(1) C C ZCHDC COMPUTES THE CHOLESKY DECOMPOSITION OF A POSITIVE DEFINITE C MATRIX. A PIVOTING OPTION ALLOWS THE USER TO ESTIMATE THE C CONDITION OF A POSITIVE DEFINITE MATRIX OR DETERMINE THE RANK C OF A POSITIVE SEMIDEFINITE MATRIX. C C ON ENTRY C C A COMPLEX*16(LDA,P). C A CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO C BE COMPUTED. ONLT THE UPPER HALF OF A NEED BE STORED. C THE LOWER PART OF THE ARRAY A IS NOT REFERENCED. C C LDA INTEGER. C LDA IS THE LEADING DIMENSION OF THE ARRAY A. C C P INTEGER. C P IS THE ORDER OF THE MATRIX. C C WORK COMPLEX*16. C WORK IS A WORK ARRAY. C C JPVT INTEGER(P). C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION C OF THE PIVOT ELEMENTS, IF PIVOTING HAS BEEN REQUESTED. C EACH DIAGONAL ELEMENT A(K,K) C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE C VALUE OF JPVT(K). C C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C ELEMENT. C C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE ELEMENT. C C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL ELEMENT. C C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL ELEMENTS C ARE MOVED BY SYMMETRIC ROW AND COLUMN INTERCHANGES TO C THE BEGINNING OF THE ARRAY A AND FINAL C ELEMENTS TO THE END. BOTH INITIAL AND FINAL ELEMENTS C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY C FREE ELEMENTS ARE MOVED. AT THE K-TH STAGE OF THE C REDUCTION, IF A(K,K) IS OCCUPIED BY A FREE ELEMENT C IT IS INTERCHANGED WITH THE LARGEST FREE ELEMENT C A(L,L) WITH L .GE. K. JPVT IS NOT REFERENCED IF C JOB .EQ. 0. C C JOB INTEGER. C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING. C IF JOB .EQ. 0, NO PIVOTING IS DONE. C IF JOB .NE. 0, PIVOTING IS DONE. C C ON RETURN C C A A CONTAINS IN ITS UPPER HALF THE CHOLESKY FACTOR C OF THE MATRIX A AS IT HAS BEEN PERMUTED BY PIVOTING. C C JPVT JPVT(J) CONTAINS THE INDEX OF THE DIAGONAL ELEMENT C OF A THAT WAS MOVED INTO THE J-TH POSITION, C PROVIDED PIVOTING WAS REQUESTED. C C INFO CONTAINS THE INDEX OF THE LAST POSITIVE DIAGONAL C ELEMENT OF THE CHOLESKY FACTOR. C C FOR POSITIVE DEFINITE MATRICES INFO = P IS THE NORMAL RETURN. C FOR PIVOTING WITH POSITIVE SEMIDEFINITE MATRICES INFO WILL C IN GENERAL BE LESS THAN P. HOWEVER, INFO MAY BE GREATER THAN C THE RANK OF A, SINCE ROUNDING ERROR CAN CAUSE AN OTHERWISE ZERO C ELEMENT TO BE POSITIVE. INDEFINITE SYSTEMS WILL ALWAYS CAUSE C INFO TO BE LESS THAN P. C C LINPACK. THIS VERSION DATED 03/19/79 . C J.J. DONGARRA AND G.W. STEWART, ARGONNE NATIONAL LABORATORY AND C UNIVERSITY OF MARYLAND. C C C BLAS ZAXPY,ZSWAP C FORTRAN DSQRT,DCONJG C C INTERNAL VARIABLES C SUBROUTINE ZCHDD(R,LDR,P,X,Z,LDZ,NZ,Y,RHO,C,S,INFO) INTEGER LDR,P,LDZ,NZ,INFO COMPLEX*16 R(LDR,1),X(1),Z(LDZ,1),Y(1),S(1) DOUBLE PRECISION RHO(1),C(1) C C ZCHDD DOWNDATES AN AUGMENTED CHOLESKY DECOMPOSITION OR THE C TRIANGULAR FACTOR OF AN AUGMENTED QR DECOMPOSITION. C SPECIFICALLY, GIVEN AN UPPER TRIANGULAR MATRIX R OF ORDER P, A C ROW VECTOR X, A COLUMN VECTOR Z, AND A SCALAR Y, ZCHDD C DETERMINEDS A UNITARY MATRIX U AND A SCALAR ZETA SUCH THAT C C (R Z ) (RR ZZ) C U * ( ) = ( ) , C (0 ZETA) ( X Y) C C WHERE RR IS UPPER TRIANGULAR. IF R AND Z HAVE BEEN OBTAINED C FROM THE FACTORIZATION OF A LEAST SQUARES PROBLEM, THEN C RR AND ZZ ARE THE FACTORS CORRESPONDING TO THE PROBLEM C WITH THE OBSERVATION (X,Y) REMOVED. IN THIS CASE, IF RHO C IS THE NORM OF THE RESIDUAL VECTOR, THEN THE NORM OF C THE RESIDUAL VECTOR OF THE DOWNDATED PROBLEM IS C DSQRT(RHO**2 - ZETA**2). ZCHDD WILL SIMULTANEOUSLY DOWNDATE C SEVERAL TRIPLETS (Z,Y,RHO) ALONG WITH R. C FOR A LESS TERSE DESCRIPTION OF WHAT ZCHDD DOES AND HOW C IT MAY BE APPLIED, SEE THE LINPACK GUIDE. C C THE MATRIX U IS DETERMINED AS THE PRODUCT U(1)*...*U(P) C WHERE U(I) IS A ROTATION IN THE (P+1,I)-PLANE OF THE C FORM C C ( C(I) -DCONJG(S(I)) ) C ( ) . C ( S(I) C(I) ) C C THE ROTATIONS ARE CHOSEN SO THAT C(I) IS DOUBLE PRECISION. C C THE USER IS WARNED THAT A GIVEN DOWNDATING PROBLEM MAY C BE IMPOSSIBLE TO ACCOMPLISH OR MAY PRODUCE C INACCURATE RESULTS. FOR EXAMPLE, THIS CAN HAPPEN C IF X IS NEAR A VECTOR WHOSE REMOVAL WILL REDUCE THE C RANK OF R. BEWARE. C C ON ENTRY C C R COMPLEX*16(LDR,P), WHERE LDR .GE. P. C R CONTAINS THE UPPER TRIANGULAR MATRIX C THAT IS TO BE DOWNDATED. THE PART OF R C BELOW THE DIAGONAL IS NOT REFERENCED. C C LDR INTEGER. C LDR IS THE LEADING DIMENSION FO THE ARRAY R. C C P INTEGER. C P IS THE ORDER OF THE MATRIX R. C C X COMPLEX*16(P). C X CONTAINS THE ROW VECTOR THAT IS TO C BE REMOVED FROM R. X IS NOT ALTERED BY ZCHDD. C C Z COMPLEX*16(LDZ,NZ), WHERE LDZ .GE. P. C Z IS AN ARRAY OF NZ P-VECTORS WHICH C ARE TO BE DOWNDATED ALONG WITH R. C C LDZ INTEGER. C LDZ IS THE LEADING DIMENSION OF THE ARRAY Z. C C NZ INTEGER. C NZ IS THE NUMBER OF VECTORS TO BE DOWNDATED C NZ MAY BE ZERO, IN WHICH CASE Z, Y, AND RHO C ARE NOT REFERENCED. C C Y COMPLEX*16(NZ). C Y CONTAINS THE SCALARS FOR THE DOWNDATING C OF THE VECTORS Z. Y IS NOT ALTERED BY ZCHDD. C C RHO DOUBLE PRECISION(NZ). C RHO CONTAINS THE NORMS OF THE RESIDUAL C VECTORS THAT ARE TO BE DOWNDATED. C C ON RETURN C C R C Z CONTAIN THE DOWNDATED QUANTITIES. C RHO C C C DOUBLE PRECISION(P). C C CONTAINS THE COSINES OF THE TRANSFORMING C ROTATIONS. C C S COMPLEX*16(P). C S CONTAINS THE SINES OF THE TRANSFORMING C ROTATIONS. C C INFO INTEGER. C INFO IS SET AS FOLLOWS. C C INFO = 0 IF THE ENTIRE DOWNDATING C WAS SUCCESSFUL. C C INFO =-1 IF R COULD NOT BE DOWNDATED. C IN THIS CASE, ALL QUANTITIES C ARE LEFT UNALTERED. C C INFO = 1 IF SOME RHO COULD NOT BE C DOWNDATED. THE OFFENDING RHOS ARE C SET TO -1. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZCHDD USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C FORTRAN CDABS,DCONJG C BLAS ZDOTC, DZNRM2 C SUBROUTINE ZCHEX(R,LDR,P,K,L,Z,LDZ,NZ,C,S,JOB) INTEGER LDR,P,K,L,LDZ,NZ,JOB COMPLEX*16 R(LDR,1),Z(LDZ,1),S(1) DOUBLE PRECISION C(1) C C ZCHEX UPDATES THE CHOLESKY FACTORIZATION C C A = CTRANS(R)*R C C OF A POSITIVE DEFINITE MATRIX A OF ORDER P UNDER DIAGONAL C PERMUTATIONS OF THE FORM C C TRANS(E)*A*E C C WHERE E IS A PERMUTATION MATRIX. SPECIFICALLY, GIVEN C AN UPPER TRIANGULAR MATRIX R AND A PERMUTATION MATRIX C E (WHICH IS SPECIFIED BY K, L, AND JOB), ZCHEX DETERMINES C A UNITARY MATRIX U SUCH THAT C C U*R*E = RR, C C WHERE RR IS UPPER TRIANGULAR. AT THE USERS OPTION, THE C TRANSFORMATION U WILL BE MULTIPLIED INTO THE ARRAY Z. C IF A = CTRANS(X)*X, SO THAT R IS THE TRIANGULAR PART OF THE C QR FACTORIZATION OF X, THEN RR IS THE TRIANGULAR PART OF THE C QR FACTORIZATION OF X*E, I.E. X WITH ITS COLUMNS PERMUTED. C FOR A LESS TERSE DESCRIPTION OF WHAT ZCHEX DOES AND HOW C IT MAY BE APPLIED, SEE THE LINPACK GUIDE. C C THE MATRIX Q IS DETERMINED AS THE PRODUCT U(L-K)*...*U(1) C OF PLANE ROTATIONS OF THE FORM C C ( C(I) S(I) ) C ( ) , C ( -DCONJG(S(I)) C(I) ) C C WHERE C(I) IS DOUBLE PRECISION, THE ROWS THESE ROTATIONS OPERATE C ON ARE DESCRIBED BELOW. C C THERE ARE TWO TYPES OF PERMUTATIONS, WHICH ARE DETERMINED C BY THE VALUE OF JOB. C C 1. RIGHT CIRCULAR SHIFT (JOB = 1). C C THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER. C C 1,...,K-1,L,K,K+1,...,L-1,L+1,...,P. C C U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I) C ACTS IN THE (L-I,L-I+1)-PLANE. C C 2. LEFT CIRCULAR SHIFT (JOB = 2). C THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER C C 1,...,K-1,K+1,K+2,...,L,K,L+1,...,P. C C U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I) C ACTS IN THE (K+I-1,K+I)-PLANE. C C ON ENTRY C C R COMPLEX*16(LDR,P), WHERE LDR.GE.P. C R CONTAINS THE UPPER TRIANGULAR FACTOR C THAT IS TO BE UPDATED. ELEMENTS OF R C BELOW THE DIAGONAL ARE NOT REFERENCED. C C LDR INTEGER. C LDR IS THE LEADING DIMENSION OF THE ARRAY R. C C P INTEGER. C P IS THE ORDER OF THE MATRIX R. C C K INTEGER. C K IS THE FIRST COLUMN TO BE PERMUTED. C C L INTEGER. C L IS THE LAST COLUMN TO BE PERMUTED. C L MUST BE STRICTLY GREATER THAN K. C C Z COMPLEX*16(LDZ,NZ), WHERE LDZ.GE.P. C Z IS AN ARRAY OF NZ P-VECTORS INTO WHICH THE C TRANSFORMATION U IS MULTIPLIED. Z IS C NOT REFERENCED IF NZ = 0. C C LDZ INTEGER. C LDZ IS THE LEADING DIMENSION OF THE ARRAY Z. C C NZ INTEGER. C NZ IS THE NUMBER OF COLUMNS OF THE MATRIX Z. C C JOB INTEGER. C JOB DETERMINES THE TYPE OF PERMUTATION. C JOB = 1 RIGHT CIRCULAR SHIFT. C JOB = 2 LEFT CIRCULAR SHIFT. C C ON RETURN C C R CONTAINS THE UPDATED FACTOR. C C Z CONTAINS THE UPDATED MATRIX Z. C C C DOUBLE PRECISION(P). C C CONTAINS THE COSINES OF THE TRANSFORMING ROTATIONS. C C S COMPLEX*16(P). C S CONTAINS THE SINES OF THE TRANSFORMING ROTATIONS. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZCHEX USES THE FOLLOWING FUNCTIONS AND SUBROUTINES. C C EXTENDED BLAS ZROTG C FORTRAN MIN0 C SUBROUTINE ZCHUD(R,LDR,P,X,Z,LDZ,NZ,Y,RHO,C,S) INTEGER LDR,P,LDZ,NZ DOUBLE PRECISION RHO(1),C(1) COMPLEX*16 R(LDR,1),X(1),Z(LDZ,1),Y(1),S(1) C C ZCHUD UPDATES AN AUGMENTED CHOLESKY DECOMPOSITION OF THE C TRIANGULAR PART OF AN AUGMENTED QR DECOMPOSITION. SPECIFICALLY, C GIVEN AN UPPER TRIANGULAR MATRIX R OF ORDER P, A ROW VECTOR C X, A COLUMN VECTOR Z, AND A SCALAR Y, ZCHUD DETERMINES A C UNTIARY MATRIX U AND A SCALAR ZETA SUCH THAT C C C (R Z) (RR ZZ ) C U * ( ) = ( ) , C (X Y) ( 0 ZETA) C C WHERE RR IS UPPER TRIANGULAR. IF R AND Z HAVE BEEN C OBTAINED FROM THE FACTORIZATION OF A LEAST SQUARES C PROBLEM, THEN RR AND ZZ ARE THE FACTORS CORRESPONDING TO C THE PROBLEM WITH THE OBSERVATION (X,Y) APPENDED. IN THIS C CASE, IF RHO IS THE NORM OF THE RESIDUAL VECTOR, THEN THE C NORM OF THE RESIDUAL VECTOR OF THE UPDATED PROBLEM IS C DSQRT(RHO**2 + ZETA**2). ZCHUD WILL SIMULTANEOUSLY UPDATE C SEVERAL TRIPLETS (Z,Y,RHO). C FOR A LESS TERSE DESCRIPTION OF WHAT ZCHUD DOES AND HOW C IT MAY BE APPLIED, SEE THE LINPACK GUIDE. C C THE MATRIX U IS DETERMINED AS THE PRODUCT U(P)*...*U(1), C WHERE U(I) IS A ROTATION IN THE (I,P+1) PLANE OF THE C FORM C C ( C(I) S(I) ) C ( ) . C ( -DCONJG(S(I)) C(I) ) C C THE ROTATIONS ARE CHOSEN SO THAT C(I) IS DOUBLE PRECISION. C C ON ENTRY C C R COMPLEX*16(LDR,P), WHERE LDR .GE. P. C R CONTAINS THE UPPER TRIANGULAR MATRIX C THAT IS TO BE UPDATED. THE PART OF R C BELOW THE DIAGONAL IS NOT REFERENCED. C C LDR INTEGER. C LDR IS THE LEADING DIMENSION OF THE ARRAY R. C C P INTEGER. C P IS THE ORDER OF THE MATRIX R. C C X COMPLEX*16(P). C X CONTAINS THE ROW TO BE ADDED TO R. X IS C NOT ALTERED BY ZCHUD. C C Z COMPLEX*16(LDZ,NZ), WHERE LDZ .GE. P. C Z IS AN ARRAY CONTAINING NZ P-VECTORS TO C BE UPDATED WITH R. C C LDZ INTEGER. C LDZ IS THE LEADING DIMENSION OF THE ARRAY Z. C C NZ INTEGER. C NZ IS THE NUMBER OF VECTORS TO BE UPDATED C NZ MAY BE ZERO, IN WHICH CASE Z, Y, AND RHO C ARE NOT REFERENCED. C C Y COMPLEX*16(NZ). C Y CONTAINS THE SCALARS FOR UPDATING THE VECTORS C Z. Y IS NOT ALTERED BY ZCHUD. C C RHO DOUBLE PRECISION(NZ). C RHO CONTAINS THE NORMS OF THE RESIDUAL C VECTORS THAT ARE TO BE UPDATED. IF RHO(J) C IS NEGATIVE, IT IS LEFT UNALTERED. C C ON RETURN C C RC C RHO CONTAIN THE UPDATED QUANTITIES. C Z C C C DOUBLE PRECISION(P). C C CONTAINS THE COSINES OF THE TRANSFORMING C ROTATIONS. C C S COMPLEX*16(P). C S CONTAINS THE SINES OF THE TRANSFORMING C ROTATIONS. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZCHUD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES. C C EXTENDED BLAS ZROTG C FORTRAN DCONJG,DSQRT C SUBROUTINE ZSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO) INTEGER LDX,N,P,LDU,LDV,JOB,INFO COMPLEX*16 X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1) C C C ZSVDC IS A SUBROUTINE TO REDUCE A COMPLEX*16 NXP MATRIX X BY C UNITARY TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS, C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS. C C ON ENTRY C C X COMPLEX*16(LDX,P), WHERE LDX.GE.N. C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE C DECOMPOSITION IS TO BE COMPUTED. X IS C DESTROYED BY ZSVDC. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF ROWS OF THE MATRIX X. C C LDU INTEGER. C LDU IS THE LEADING DIMENSION OF THE ARRAY U C (SEE BELOW). C C LDV INTEGER. C LDV IS THE LEADING DIMENSION OF THE ARRAY V C (SEE BELOW). C C WORK COMPLEX*16(N). C WORK IS A SCRATCH ARRAY. C C JOB INTEGER. C JOB CONTROLS THE COMPUTATION OF THE SINGULAR C VECTORS. IT HAS THE DECIMAL EXPANSION AB C WITH THE FOLLOWING MEANING C C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR C VECTORS. C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS C IN U. C A.GE.2 RETURNS THE FIRST MIN(N,P) C LEFT SINGULAR VECTORS IN U. C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR C VECTORS. C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS C IN V. C C ON RETURN C C S COMPLEX*16(MM), WHERE MM=MIN(N+1,P). C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE C SINGULAR VALUES OF X ARRANGED IN DESCENDING C ORDER OF MAGNITUDE. C C E COMPLEX*16(P). C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE C DISCUSSION OF INFO FOR EXCEPTIONS. C C U COMPLEX*16(LDU,K), WHERE LDU.GE.N. IF JOBA.EQ.1 C THEN K.EQ.N, IF JOBA.GE.2 THEN C C K.EQ.MIN(N,P). C U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P C OR IF JOBA.GT.2, THEN U MAY BE IDENTIFIED WITH X C IN THE SUBROUTINE CALL. C C V COMPLEX*16(LDV,P), WHERE LDV.GE.P. C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. C V IS NOT REFERENCED IF JOBB.EQ.0. IF P.LE.N, C THEN V MAY BE IDENTIFIED WHTH X IN THE C SUBROUTINE CALL. C C INFO INTEGER. C THE SINGULAR VALUES (AND THEIR CORRESPONDING C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M) C ARE CORRECT (HERE M=MIN(N,P)). THUS IF C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX C B = CTRANS(U)*X*V IS THE BIDIAGONAL MATRIX C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE C ELEMENTS OF E ON ITS SUPER-DIAGONAL (CTRANS(U) C IS THE CONJUGATE-TRANSPOSE OF U). THUS THE C SINGULAR VALUES OF X AND B ARE THE SAME. C C LINPACK. THIS VERSION DATED 03/19/79 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C ZSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C EXTERNAL ZDROT C BLAS ZAXPY,ZDOTC,ZSCAL,ZSWAP,DZNRM2,DROTG C FORTRAN DABS,DMAX1,CDABS,DCMPLX C FORTRAN DCONJG,MAX0,MIN0,MOD,DSQRT C C INTERNAL VARIABLES C .