SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(1) DOUBLE PRECISION A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. C TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. C TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. C TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. C C ON ENTRY C C A DOUBLE PRECISION(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 DOUBLE PRECISION(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 SGEFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,SASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C C C COMPUTE 1-NORM OF A C ANORM = 0.0D0 DO 10 J = 1, N ANORM = DMAX1(ANORM,SASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL SGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID C OVERFLOW. C C SOLVE TRANS(U)*W = E C EK = 1.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 S = DABS(A(K,K))/DABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) IF (A(K,K) .EQ. 0.0D0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + DABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + DABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM DO 70 J = KP1, N Z(J) = Z(J) + T*A(K,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + SDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 S = 1.0D0/DABS(Z(K)) CALL SSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE L*V = Y C DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) CALL SAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 S = 1.0D0/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 S = DABS(A(K,K))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 T = -Z(K) CALL SAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(1),INFO DOUBLE PRECISION A(LDA,1) C C SGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. C C SGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR SGEFA) . C C ON ENTRY C C A DOUBLE PRECISION(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 DGESL OR DGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN DGECO 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 SAXPY,SSCAL,ISAMAX C C INTERNAL VARIABLES C DOUBLE PRECISION T INTEGER ISAMAX,J,K,KP1,L,NM1 C C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = ISAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K) .EQ. 0.0D0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0D0/A(K,K) CALL SSCAL(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0D0) INFO = N RETURN END SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),B(1) C C DGESL SOLVES THE DOUBLE PRECISION SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGECO OR SGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR SGEFA. 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 DGECO OR SGEFA. C C B DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE 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 DGECO HAS SET RCOND .GT. 0.0 C OR SGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGESL(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 SAXPY,SDOT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T INTEGER K,KB,L,NM1 C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE SGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),DET(2),WORK(1) C C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY DGECO OR SGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR SGEFA. 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 DGECO OR SGEFA. C C WORK DOUBLE PRECISION(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 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. DABS(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 DGECO HAS SET RCOND .GT. 0.0 OR SGEFA 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 SAXPY,SSCAL,SSWAP C FORTRAN DABS,MOD C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION TEN INTEGER I,J,K,KB,KP1,L,NM1 C C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(U) C IF (MOD(JOB,10) .EQ. 0) GO TO 150 DO 100 K = 1, N A(K,K) = 1.0D0/A(K,K) T = -A(K,K) CALL SSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0D0 CALL SAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(U)*INVERSE(L) C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 140 DO 130 KB = 1, NM1 K = N - KB KP1 = K + 1 DO 110 I = KP1, N WORK(I) = A(I,K) A(I,K) = 0.0D0 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL SAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL SSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END SUBROUTINE SGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) INTEGER LDA,N,ML,MU,IPVT(1) DOUBLE PRECISION ABD(LDA,1),Z(1) DOUBLE PRECISION RCOND C C DGBCO FACTORS A DOUBLE PRECISION BAND MATRIX BY GAUSSIAN C ELIMINATION AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SGBFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DGBCO BY DGBSL. C TO COMPUTE INVERSE(A)*C , FOLLOW DGBCO BY DGBSL. C TO COMPUTE DETERMINANT(A) , FOLLOW DGBCO BY DGBDI. C C ON ENTRY C C ABD DOUBLE PRECISION(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 DOUBLE PRECISION(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 SGBFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,MAX0,MIN0,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,SASUM,SM,YNORM INTEGER IS,INFO,J,JU,K,KB,KP1,L,LA,LM,LZ,M,MM C C C COMPUTE 1-NORM OF A C ANORM = 0.0D0 L = ML + 1 IS = L + MU DO 10 J = 1, N ANORM = DMAX1(ANORM,SASUM(L,ABD(IS,J),1)) IF (IS .GT. ML + 1) IS = IS - 1 IF (J .LE. MU) L = L + 1 IF (J .GE. N - ML) L = L - 1 10 CONTINUE C C FACTOR C CALL SGBFA(ABD,LDA,N,ML,MU,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID C OVERFLOW. C C SOLVE TRANS(U)*W = E C EK = 1.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE M = ML + MU + 1 JU = 0 DO 100 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. DABS(ABD(M,K))) GO TO 30 S = DABS(ABD(M,K))/DABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) IF (ABD(M,K) .EQ. 0.0D0) GO TO 40 WK = WK/ABD(M,K) WKM = WKM/ABD(M,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE KP1 = K + 1 JU = MIN0(MAX0(JU,MU+IPVT(K)),N) MM = M IF (KP1 .GT. JU) GO TO 90 DO 60 J = KP1, JU MM = MM - 1 SM = SM + DABS(Z(J)+WKM*ABD(MM,J)) Z(J) = Z(J) + WK*ABD(MM,J) S = S + DABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM MM = M DO 70 J = KP1, JU MM = MM - 1 Z(J) = Z(J) + T*ABD(MM,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB LM = MIN0(ML,N-K) IF (K .LT. N) Z(K) = Z(K) + SDOT(LM,ABD(M+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 S = 1.0D0/DABS(Z(K)) CALL SSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE L*V = Y C DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T LM = MIN0(ML,N-K) IF (K .LT. N) CALL SAXPY(LM,T,ABD(M+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 S = 1.0D0/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = W C DO 160 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. DABS(ABD(M,K))) GO TO 150 S = DABS(ABD(M,K))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (ABD(M,K) .NE. 0.0D0) Z(K) = Z(K)/ABD(M,K) IF (ABD(M,K) .EQ. 0.0D0) Z(K) = 1.0D0 LM = MIN0(K,M) - 1 LA = M - LM LZ = K - LM T = -Z(K) CALL SAXPY(LM,T,ABD(LA,K),1,Z(LZ),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END SUBROUTINE SGBFA(ABD,LDA,N,ML,MU,IPVT,INFO) INTEGER LDA,N,ML,MU,IPVT(1),INFO DOUBLE PRECISION ABD(LDA,1) C C SGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION. C C SGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD DOUBLE PRECISION(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 DGBSL WILL DIVIDE BY ZERO IF C CALLED. USE RCOND IN DGBCO 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 SAXPY,SSCAL,ISAMAX C FORTRAN MAX0,MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION T INTEGER I,ISAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 C C M = ML + MU + 1 INFO = 0 C C ZERO INITIAL FILL-IN COLUMNS C J0 = MU + 2 J1 = MIN0(N,M) - 1 IF (J1 .LT. J0) GO TO 30 DO 20 JZ = J0, J1 I0 = M + 1 - JZ DO 10 I = I0, ML ABD(I,JZ) = 0.0D0 10 CONTINUE 20 CONTINUE 30 CONTINUE JZ = J1 JU = 0 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 130 DO 120 K = 1, NM1 KP1 = K + 1 C C ZERO NEXT FILL-IN COLUMN C JZ = JZ + 1 IF (JZ .GT. N) GO TO 50 IF (ML .LT. 1) GO TO 50 DO 40 I = 1, ML ABD(I,JZ) = 0.0D0 40 CONTINUE 50 CONTINUE C C FIND L = PIVOT INDEX C LM = MIN0(ML,N-K) L = ISAMAX(LM+1,ABD(M,K),1) + M - 1 IPVT(K) = L + K - M C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (ABD(L,K) .EQ. 0.0D0) GO TO 100 C C INTERCHANGE IF NECESSARY C IF (L .EQ. M) GO TO 60 T = ABD(L,K) ABD(L,K) = ABD(M,K) ABD(M,K) = T 60 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0D0/ABD(M,K) CALL SSCAL(LM,T,ABD(M+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C JU = MIN0(MAX0(JU,MU+IPVT(K)),N) MM = M IF (JU .LT. KP1) GO TO 90 DO 80 J = KP1, JU L = L - 1 MM = MM - 1 T = ABD(L,J) IF (L .EQ. MM) GO TO 70 ABD(L,J) = ABD(MM,J) ABD(MM,J) = T 70 CONTINUE CALL SAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1) 80 CONTINUE 90 CONTINUE GO TO 110 100 CONTINUE INFO = K 110 CONTINUE 120 CONTINUE 130 CONTINUE IPVT(N) = N IF (ABD(M,N) .EQ. 0.0D0) INFO = N RETURN END SUBROUTINE SGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB) INTEGER LDA,N,ML,MU,IPVT(1),JOB DOUBLE PRECISION ABD(LDA,1),B(1) C C DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGBCO OR SGBFA. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGBCO OR SGBFA. 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 DGBCO OR SGBFA. C C B DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B , WHERE C TRANS(A) IS THE 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 DGBCO HAS SET RCOND .GT. 0.0 C OR SGBFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGBSL(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 SAXPY,SDOT C FORTRAN MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T INTEGER K,KB,L,LA,LB,LM,M,NM1 C M = MU + ML + 1 NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (ML .EQ. 0) GO TO 30 IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 LM = MIN0(ML,N-K) L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(LM,T,ABD(M+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/ABD(M,K) LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = -B(K) CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = SDOT(LM,ABD(LA,K),1,B(LB),1) B(K) = (B(K) - T)/ABD(M,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (ML .EQ. 0) GO TO 90 IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB LM = MIN0(ML,N-K) B(K) = B(K) + SDOT(LM,ABD(M+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE SGBDI(ABD,LDA,N,ML,MU,IPVT,DET) INTEGER LDA,N,ML,MU,IPVT(1) DOUBLE PRECISION ABD(LDA,1),DET(2) C C DGBDI COMPUTES THE DETERMINANT OF A BAND MATRIX C USING THE FACTORS COMPUTED BY DGBCO OR SGBFA. C IF THE INVERSE IS NEEDED, USE DGBSL N TIMES. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGBCO OR SGBFA. 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 DGBCO OR SGBFA. C C ON RETURN 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 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 C C INTERNAL VARIABLES C DOUBLE PRECISION TEN INTEGER I,M C C M = ML + MU + 1 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = ABD(M,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE SPOCO(A,LDA,N,RCOND,Z,INFO) INTEGER LDA,N,INFO DOUBLE PRECISION A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C DPOCO FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SPOFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DPOCO BY DPOSL. C TO COMPUTE INVERSE(A)*C , FOLLOW DPOCO BY DPOSL. C TO COMPUTE DETERMINANT(A) , FOLLOW DPOCO BY DPODI. C TO COMPUTE INVERSE(A) , FOLLOW DPOCO BY DPODI. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE SYMMETRIC 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 = TRANS(R)*R C WHERE TRANS(R) IS THE TRANSPOSE. C 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 DOUBLE PRECISION(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 SPOFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,DREAL,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,SASUM,SM,YNORM INTEGER I,J,JM1,K,KB,KP1 C C C FIND NORM OF A USING ONLY UPPER HALF C DO 30 J = 1, N Z(J) = SASUM(J,A(1,J),1) JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 I = 1, JM1 Z(I) = Z(I) + DABS(A(I,J)) 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0D0 DO 40 J = 1, N ANORM = DMAX1(ANORM,Z(J)) 40 CONTINUE C C FACTOR C CALL SPOFA(A,LDA,N,INFO) IF (INFO .NE. 0) GO TO 180 C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE TRANS(R)*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE TRANS(R)*W = E C EK = 1.0D0 DO 50 J = 1, N Z(J) = 0.0D0 50 CONTINUE DO 110 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. A(K,K)) GO TO 60 S = A(K,K)/DABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 60 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) WK = WK/A(K,K) WKM = WKM/A(K,K) KP1 = K + 1 IF (KP1 .GT. N) GO TO 100 DO 70 J = KP1, N SM = SM + DABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + DABS(Z(J)) 70 CONTINUE IF (S .GE. SM) GO TO 90 T = WKM - WK WK = WKM DO 80 J = KP1, N Z(J) = Z(J) + T*A(K,J) 80 CONTINUE 90 CONTINUE 100 CONTINUE Z(K) = WK 110 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE R*Y = W C DO 130 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. A(K,K)) GO TO 120 S = A(K,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) 120 CONTINUE Z(K) = Z(K)/A(K,K) T = -Z(K) CALL SAXPY(K-1,T,A(1,K),1,Z(1),1) 130 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE TRANS(R)*V = Y C DO 150 K = 1, N Z(K) = Z(K) - SDOT(K-1,A(1,K),1,Z(1),1) IF (DABS(Z(K)) .LE. A(K,K)) GO TO 140 S = A(K,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 140 CONTINUE Z(K) = Z(K)/A(K,K) 150 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE R*Z = V C DO 170 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. A(K,K)) GO TO 160 S = A(K,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 160 CONTINUE Z(K) = Z(K)/A(K,K) T = -Z(K) CALL SAXPY(K-1,T,A(1,K),1,Z(1),1) 170 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 180 CONTINUE RETURN END SUBROUTINE SPOFA(A,LDA,N,INFO) INTEGER LDA,N,INFO DOUBLE PRECISION A(LDA,1) C C SPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX. C C SPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR SPOFA) . C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE SYMMETRIC 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 = TRANS(R)*R C WHERE TRANS(R) IS THE TRANSPOSE. C 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 SDOT C FORTRAN DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T DOUBLE PRECISION S INTEGER J,JM1,K C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1, N INFO = J S = 0.0D0 JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 K = 1, JM1 T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) A(K,J) = T S = S + T*T 10 CONTINUE 20 CONTINUE S = A(J,J) - S C ......EXIT IF (S .LE. 0.0D0) GO TO 40 A(J,J) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE SPOSL(A,LDA,N,B) INTEGER LDA,N DOUBLE PRECISION A(LDA,1),B(1) C C DPOSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C SYSTEM A * X = B C USING THE FACTORS COMPUTED BY DPOCO OR SPOFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DPOCO OR SPOFA. 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 DOUBLE PRECISION(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 DPOCO(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 DPOSL(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 SAXPY,SDOT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T INTEGER K,KB C C SOLVE TRANS(R)*Y = B C DO 10 K = 1, N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 10 CONTINUE C C SOLVE R*X = Y C DO 20 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1,T,A(1,K),1,B(1),1) 20 CONTINUE RETURN END SUBROUTINE SPODI(A,LDA,N,DET,JOB) INTEGER LDA,N,JOB DOUBLE PRECISION A(LDA,1) DOUBLE PRECISION DET(2) C C DPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN C DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX (SEE BELOW) C USING THE FACTORS COMPUTED BY DPOCO, SPOFA OR DQRDC. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT A FROM DPOCO OR SPOFA C OR THE OUTPUT X FROM DQRDC. 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 DPOCO OR SPOFA WAS USED TO FACTOR A THEN C DPODI PRODUCES THE UPPER HALF OF INVERSE(A) . C IF DQRDC WAS USED TO DECOMPOSE X THEN C DPODI PRODUCES THE UPPER HALF OF INVERSE(TRANS(X)*X) C WHERE TRANS(X) IS THE 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 TRANS(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 DPOCO OR SPOFA 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 SAXPY,SSCAL C FORTRAN MOD C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION S INTEGER I,J,JM1,K,KP1 C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 S = 10.0D0 DO 50 I = 1, N DET(1) = A(I,I)**2*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DET(1) .GE. 1.0D0) GO TO 20 DET(1) = S*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DET(1) .LT. S) GO TO 40 DET(1) = DET(1)/S DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(R) C IF (MOD(JOB,10) .EQ. 0) GO TO 140 DO 100 K = 1, N A(K,K) = 1.0D0/A(K,K) T = -A(K,K) CALL SSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0D0 CALL SAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(R) * TRANS(INVERSE(R)) C DO 130 J = 1, N JM1 = J - 1 IF (JM1 .LT. 1) GO TO 120 DO 110 K = 1, JM1 T = A(K,J) CALL SAXPY(K,T,A(1,J),1,A(1,K),1) 110 CONTINUE 120 CONTINUE T = A(J,J) CALL SSCAL(J,T,A(1,J),1) 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE SPPCO(AP,N,RCOND,Z,INFO) INTEGER N,INFO DOUBLE PRECISION AP(1),Z(1) DOUBLE PRECISION RCOND C C DPPCO FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX STORED IN PACKED FORM C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SPPFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DPPCO BY DPPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW DPPCO BY DPPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW DPPCO BY DPPDI. C TO COMPUTE INVERSE(A) , FOLLOW DPPCO BY DPPDI. C C ON ENTRY C C AP DOUBLE PRECISION (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 AN UPPER TRIANGULAR MATRIX R , STORED IN PACKED C FORM, SO THAT A = TRANS(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 DOUBLE PRECISION(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 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 SPPFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,DREAL,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,SASUM,SM,YNORM INTEGER I,IJ,J,JM1,J1,K,KB,KJ,KK,KP1 C C C FIND NORM OF A C J1 = 1 DO 30 J = 1, N Z(J) = SASUM(J,AP(J1),1) IJ = J1 J1 = J1 + J JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 I = 1, JM1 Z(I) = Z(I) + DABS(AP(IJ)) IJ = IJ + 1 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0D0 DO 40 J = 1, N ANORM = DMAX1(ANORM,Z(J)) 40 CONTINUE C C FACTOR C CALL SPPFA(AP,N,INFO) IF (INFO .NE. 0) GO TO 180 C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE TRANS(R)*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE TRANS(R)*W = E C EK = 1.0D0 DO 50 J = 1, N Z(J) = 0.0D0 50 CONTINUE KK = 0 DO 110 K = 1, N KK = KK + K IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. AP(KK)) GO TO 60 S = AP(KK)/DABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 60 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) WK = WK/AP(KK) WKM = WKM/AP(KK) KP1 = K + 1 KJ = KK + K IF (KP1 .GT. N) GO TO 100 DO 70 J = KP1, N SM = SM + DABS(Z(J)+WKM*AP(KJ)) Z(J) = Z(J) + WK*AP(KJ) S = S + DABS(Z(J)) KJ = KJ + J 70 CONTINUE IF (S .GE. SM) GO TO 90 T = WKM - WK WK = WKM KJ = KK + K DO 80 J = KP1, N Z(J) = Z(J) + T*AP(KJ) KJ = KJ + J 80 CONTINUE 90 CONTINUE 100 CONTINUE Z(K) = WK 110 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE R*Y = W C DO 130 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. AP(KK)) GO TO 120 S = AP(KK)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) 120 CONTINUE Z(K) = Z(K)/AP(KK) KK = KK - K T = -Z(K) CALL SAXPY(K-1,T,AP(KK+1),1,Z(1),1) 130 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE TRANS(R)*V = Y C DO 150 K = 1, N Z(K) = Z(K) - SDOT(K-1,AP(KK+1),1,Z(1),1) KK = KK + K IF (DABS(Z(K)) .LE. AP(KK)) GO TO 140 S = AP(KK)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 140 CONTINUE Z(K) = Z(K)/AP(KK) 150 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE R*Z = V C DO 170 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. AP(KK)) GO TO 160 S = AP(KK)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 160 CONTINUE Z(K) = Z(K)/AP(KK) KK = KK - K T = -Z(K) CALL SAXPY(K-1,T,AP(KK+1),1,Z(1),1) 170 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 180 CONTINUE RETURN END SUBROUTINE SPPFA(AP,N,INFO) INTEGER N,INFO DOUBLE PRECISION AP(1) C C SPPFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX STORED IN PACKED FORM. C C SPPFA IS USUALLY CALLED BY DPPCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DPPCO) = (1 + 18/N)*(TIME FOR SPPFA) . C C ON ENTRY C C AP DOUBLE PRECISION (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 AN UPPER TRIANGULAR MATRIX R , STORED IN PACKED C FORM, SO THAT A = TRANS(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 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 BLAS SDOT C FORTRAN DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T DOUBLE PRECISION S INTEGER J,JJ,JM1,K,KJ,KK C BEGIN BLOCK WITH ...EXITS TO 40 C C JJ = 0 DO 30 J = 1, N INFO = J S = 0.0D0 JM1 = J - 1 KJ = JJ KK = 0 IF (JM1 .LT. 1) GO TO 20 DO 10 K = 1, JM1 KJ = KJ + 1 T = AP(KJ) - SDOT(K-1,AP(KK+1),1,AP(JJ+1),1) KK = KK + K T = T/AP(KK) AP(KJ) = T S = S + T*T 10 CONTINUE 20 CONTINUE JJ = JJ + J S = AP(JJ) - S C ......EXIT IF (S .LE. 0.0D0) GO TO 40 AP(JJ) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE SPPSL(AP,N,B) INTEGER N DOUBLE PRECISION AP(1),B(1) C C DPPSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C SYSTEM A * X = B C USING THE FACTORS COMPUTED BY DPPCO OR SPPFA. C C ON ENTRY C C AP DOUBLE PRECISION (N*(N+1)/2) C THE OUTPUT FROM DPPCO OR SPPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C B DOUBLE PRECISION(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 DPPCO(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 DPPSL(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 SAXPY,SDOT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T INTEGER K,KB,KK C KK = 0 DO 10 K = 1, N T = SDOT(K-1,AP(KK+1),1,B(1),1) KK = KK + K B(K) = (B(K) - T)/AP(KK) 10 CONTINUE DO 20 KB = 1, N K = N + 1 - KB B(K) = B(K)/AP(KK) KK = KK - K T = -B(K) CALL SAXPY(K-1,T,AP(KK+1),1,B(1),1) 20 CONTINUE RETURN END SUBROUTINE SPPDI(AP,N,DET,JOB) INTEGER N,JOB DOUBLE PRECISION AP(1) DOUBLE PRECISION DET(2) C C DPPDI COMPUTES THE DETERMINANT AND INVERSE C OF A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX C USING THE FACTORS COMPUTED BY DPPCO OR SPPFA . C C ON ENTRY C C AP DOUBLE PRECISION (N*(N+1)/2) C THE OUTPUT FROM DPPCO OR SPPFA. 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 DPOCO OR SPOFA 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 SAXPY,SSCAL C FORTRAN MOD C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION S INTEGER I,II,J,JJ,JM1,J1,K,KJ,KK,KP1,K1 C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 S = 10.0D0 II = 0 DO 50 I = 1, N II = II + I DET(1) = AP(II)**2*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DET(1) .GE. 1.0D0) GO TO 20 DET(1) = S*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DET(1) .LT. S) GO TO 40 DET(1) = DET(1)/S DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(R) C IF (MOD(JOB,10) .EQ. 0) GO TO 140 KK = 0 DO 100 K = 1, N K1 = KK + 1 KK = KK + K AP(KK) = 1.0D0/AP(KK) T = -AP(KK) CALL SSCAL(K-1,T,AP(K1),1) KP1 = K + 1 J1 = KK + 1 KJ = KK + K IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = AP(KJ) AP(KJ) = 0.0D0 CALL SAXPY(K,T,AP(K1),1,AP(J1),1) J1 = J1 + J KJ = KJ + J 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(R) * TRANS(INVERSE(R)) C JJ = 0 DO 130 J = 1, N J1 = JJ + 1 JJ = JJ + J JM1 = J - 1 K1 = 1 KJ = J1 IF (JM1 .LT. 1) GO TO 120 DO 110 K = 1, JM1 T = AP(KJ) CALL SAXPY(K,T,AP(J1),1,AP(K1),1) K1 = K1 + K KJ = KJ + 1 110 CONTINUE 120 CONTINUE T = AP(JJ) CALL SSCAL(J,T,AP(J1),1) 130 CONTINUE 140 CONTINUE RETURN END SUBROUTINE SPBCO(ABD,LDA,N,M,RCOND,Z,INFO) INTEGER LDA,N,M,INFO DOUBLE PRECISION ABD(LDA,1),Z(1) DOUBLE PRECISION RCOND C C DPBCO FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX STORED IN BAND FORM AND ESTIMATES THE CONDITION OF THE C MATRIX. C C IF RCOND IS NOT NEEDED, SPBFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DPBCO BY DPBSL. C TO COMPUTE INVERSE(A)*C , FOLLOW DPBCO BY DPBSL. C TO COMPUTE DETERMINANT(A) , FOLLOW DPBCO BY DPBDI. C C ON ENTRY C C ABD DOUBLE PRECISION(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 = TRANS(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 DOUBLE PRECISION(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 SYMMETRIC 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 SPBFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,MAX0,MIN0,DREAL,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,SASUM,SM,YNORM INTEGER I,J,J2,K,KB,KP1,L,LA,LB,LM,MU C C C FIND NORM OF A C DO 30 J = 1, N L = MIN0(J,M+1) MU = MAX0(M+2-J,1) Z(J) = SASUM(L,ABD(MU,J),1) K = J - L IF (M .LT. MU) GO TO 20 DO 10 I = MU, M K = K + 1 Z(K) = Z(K) + DABS(ABD(I,J)) 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0D0 DO 40 J = 1, N ANORM = DMAX1(ANORM,Z(J)) 40 CONTINUE C C FACTOR C CALL SPBFA(ABD,LDA,N,M,INFO) IF (INFO .NE. 0) GO TO 180 C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE TRANS(R)*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE TRANS(R)*W = E C EK = 1.0D0 DO 50 J = 1, N Z(J) = 0.0D0 50 CONTINUE DO 110 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. ABD(M+1,K)) GO TO 60 S = ABD(M+1,K)/DABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 60 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) WK = WK/ABD(M+1,K) WKM = WKM/ABD(M+1,K) KP1 = K + 1 J2 = MIN0(K+M,N) I = M + 1 IF (KP1 .GT. J2) GO TO 100 DO 70 J = KP1, J2 I = I - 1 SM = SM + DABS(Z(J)+WKM*ABD(I,J)) Z(J) = Z(J) + WK*ABD(I,J) S = S + DABS(Z(J)) 70 CONTINUE IF (S .GE. SM) GO TO 90 T = WKM - WK WK = WKM I = M + 1 DO 80 J = KP1, J2 I = I - 1 Z(J) = Z(J) + T*ABD(I,J) 80 CONTINUE 90 CONTINUE 100 CONTINUE Z(K) = WK 110 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE R*Y = W C DO 130 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. ABD(M+1,K)) GO TO 120 S = ABD(M+1,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) 120 CONTINUE Z(K) = Z(K)/ABD(M+1,K) LM = MIN0(K-1,M) LA = M + 1 - LM LB = K - LM T = -Z(K) CALL SAXPY(LM,T,ABD(LA,K),1,Z(LB),1) 130 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE TRANS(R)*V = Y C DO 150 K = 1, N LM = MIN0(K-1,M) LA = M + 1 - LM LB = K - LM Z(K) = Z(K) - SDOT(LM,ABD(LA,K),1,Z(LB),1) IF (DABS(Z(K)) .LE. ABD(M+1,K)) GO TO 140 S = ABD(M+1,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 140 CONTINUE Z(K) = Z(K)/ABD(M+1,K) 150 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE R*Z = W C DO 170 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. ABD(M+1,K)) GO TO 160 S = ABD(M+1,K)/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 160 CONTINUE Z(K) = Z(K)/ABD(M+1,K) LM = MIN0(K-1,M) LA = M + 1 - LM LB = K - LM T = -Z(K) CALL SAXPY(LM,T,ABD(LA,K),1,Z(LB),1) 170 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 180 CONTINUE RETURN END SUBROUTINE SPBFA(ABD,LDA,N,M,INFO) INTEGER LDA,N,M,INFO DOUBLE PRECISION ABD(LDA,1) C C SPBFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX STORED IN BAND FORM. C C SPBFA IS USUALLY CALLED BY DPBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD DOUBLE PRECISION(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 = TRANS(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 SYMMETRIC 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 SDOT C FORTRAN MAX0,DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T DOUBLE PRECISION S INTEGER IK,J,JK,K,MU C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1, N INFO = J S = 0.0D0 IK = M + 1 JK = MAX0(J-M,1) MU = MAX0(M+2-J,1) IF (M .LT. MU) GO TO 20 DO 10 K = MU, M T = ABD(K,J) - SDOT(K-MU,ABD(IK,JK),1,ABD(MU,J),1) T = T/ABD(M+1,JK) ABD(K,J) = T S = S + T*T IK = IK - 1 JK = JK + 1 10 CONTINUE 20 CONTINUE S = ABD(M+1,J) - S C ......EXIT IF (S .LE. 0.0D0) GO TO 40 ABD(M+1,J) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE SPBSL(ABD,LDA,N,M,B) INTEGER LDA,N,M DOUBLE PRECISION ABD(LDA,1),B(1) C C DPBSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C BAND SYSTEM A*X = B C USING THE FACTORS COMPUTED BY DPBCO OR SPBFA. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DPBCO OR SPBFA. 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 DOUBLE PRECISION(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 DPBCO(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 DPBSL(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 SAXPY,SDOT C FORTRAN MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T INTEGER K,KB,LA,LB,LM C C SOLVE TRANS(R)*Y = B C DO 10 K = 1, N LM = MIN0(K-1,M) LA = M + 1 - LM LB = K - LM T = SDOT(LM,ABD(LA,K),1,B(LB),1) B(K) = (B(K) - T)/ABD(M+1,K) 10 CONTINUE C C SOLVE R*X = Y C DO 20 KB = 1, N K = N + 1 - KB LM = MIN0(K-1,M) LA = M + 1 - LM LB = K - LM B(K) = B(K)/ABD(M+1,K) T = -B(K) CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1) 20 CONTINUE RETURN END SUBROUTINE SPBDI(ABD,LDA,N,M,DET) INTEGER LDA,N,M DOUBLE PRECISION ABD(LDA,1) DOUBLE PRECISION DET(2) C C DPBDI COMPUTES THE DETERMINANT C OF A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE BAND MATRIX C USING THE FACTORS COMPUTED BY DPBCO OR SPBFA. C IF THE INVERSE IS NEEDED, USE DPBSL N TIMES. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DPBCO OR SPBFA. 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 DOUBLE PRECISION S INTEGER I C C COMPUTE DETERMINANT C DET(1) = 1.0D0 DET(2) = 0.0D0 S = 10.0D0 DO 50 I = 1, N DET(1) = ABD(M+1,I)**2*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DET(1) .GE. 1.0D0) GO TO 20 DET(1) = S*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DET(1) .LT. S) GO TO 40 DET(1) = DET(1)/S DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE SSICO(A,LDA,N,KPVT,RCOND,Z) INTEGER LDA,N,KPVT(1) DOUBLE PRECISION A(LDA,1),Z(1) DOUBLE PRECISION RCOND C C DSICO FACTORS A DOUBLE PRECISION SYMMETRIC MATRIX BY ELIMINATION C WITH SYMMETRIC PIVOTING AND ESTIMATES THE CONDITION OF THE C MATRIX. C C IF RCOND IS NOT NEEDED, SSIFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DSICO BY DSISL. C TO COMPUTE INVERSE(A)*C , FOLLOW DSICO BY DSISL. C TO COMPUTE INVERSE(A) , FOLLOW DSICO BY DSIDI. C TO COMPUTE DETERMINANT(A) , FOLLOW DSICO BY DSIDI. C TO COMPUTE INERTIA(A), FOLLOW DSICO BY DSIDI. C C ON ENTRY C C A DOUBLE PRECISION(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 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*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 DOUBLE PRECISION(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 SSIFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,IABS,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T DOUBLE PRECISION ANORM,S,SASUM,YNORM INTEGER I,INFO,J,JM1,K,KP,KPS,KS C C C FIND NORM OF A USING ONLY UPPER HALF C DO 30 J = 1, N Z(J) = SASUM(J,A(1,J),1) JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 I = 1, JM1 Z(I) = Z(I) + DABS(A(I,J)) 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0D0 DO 40 J = 1, N ANORM = DMAX1(ANORM,Z(J)) 40 CONTINUE C C FACTOR C CALL SSIFA(A,LDA,N,KPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE U*D*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE U*D*W = E C EK = 1.0D0 DO 50 J = 1, N Z(J) = 0.0D0 50 CONTINUE K = N 60 IF (K .EQ. 0) GO TO 120 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 KP = IABS(KPVT(K)) KPS = K + 1 - KS IF (KP .EQ. KPS) GO TO 70 T = Z(KPS) Z(KPS) = Z(KP) Z(KP) = T 70 CONTINUE IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,Z(K)) Z(K) = Z(K) + EK CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1) IF (KS .EQ. 1) GO TO 80 IF (Z(K-1) .NE. 0.0D0) EK = DSIGN(EK,Z(K-1)) Z(K-1) = Z(K-1) + EK CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1) 80 CONTINUE IF (KS .EQ. 2) GO TO 100 IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 90 S = DABS(A(K,K))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 90 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 GO TO 110 100 CONTINUE AK = A(K,K)/A(K-1,K) AKM1 = A(K-1,K-1)/A(K-1,K) BK = Z(K)/A(K-1,K) BKM1 = Z(K-1)/A(K-1,K) DENOM = AK*AKM1 - 1.0D0 Z(K) = (AKM1*BK - BKM1)/DENOM Z(K-1) = (AK*BKM1 - BK)/DENOM 110 CONTINUE K = K - KS GO TO 60 120 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(U)*Y = W C K = 1 130 IF (K .GT. N) GO TO 160 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. 1) GO TO 150 Z(K) = Z(K) + SDOT(K-1,A(1,K),1,Z(1),1) IF (KS .EQ. 2) * Z(K+1) = Z(K+1) + SDOT(K-1,A(1,K+1),1,Z(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 140 T = Z(K) Z(K) = Z(KP) Z(KP) = T 140 CONTINUE 150 CONTINUE K = K + KS GO TO 130 160 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE U*D*V = Y C K = N 170 IF (K .EQ. 0) GO TO 230 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. KS) GO TO 190 KP = IABS(KPVT(K)) KPS = K + 1 - KS IF (KP .EQ. KPS) GO TO 180 T = Z(KPS) Z(KPS) = Z(KP) Z(KP) = T 180 CONTINUE CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1) IF (KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1) 190 CONTINUE IF (KS .EQ. 2) GO TO 210 IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 200 S = DABS(A(K,K))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 200 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 GO TO 220 210 CONTINUE AK = A(K,K)/A(K-1,K) AKM1 = A(K-1,K-1)/A(K-1,K) BK = Z(K)/A(K-1,K) BKM1 = Z(K-1)/A(K-1,K) DENOM = AK*AKM1 - 1.0D0 Z(K) = (AKM1*BK - BKM1)/DENOM Z(K-1) = (AK*BKM1 - BK)/DENOM 220 CONTINUE K = K - KS GO TO 170 230 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE TRANS(U)*Z = V C K = 1 240 IF (K .GT. N) GO TO 270 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. 1) GO TO 260 Z(K) = Z(K) + SDOT(K-1,A(1,K),1,Z(1),1) IF (KS .EQ. 2) * Z(K+1) = Z(K+1) + SDOT(K-1,A(1,K+1),1,Z(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 250 T = Z(K) Z(K) = Z(KP) Z(KP) = T 250 CONTINUE 260 CONTINUE K = K + KS GO TO 240 270 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END SUBROUTINE SSIFA(A,LDA,N,KPVT,INFO) INTEGER LDA,N,KPVT(1),INFO DOUBLE PRECISION A(LDA,1) C C SSIFA FACTORS A DOUBLE PRECISION SYMMETRIC MATRIX BY ELIMINATION C WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW SSIFA BY DSISL. C TO COMPUTE INVERSE(A)*C , FOLLOW SSIFA BY DSISL. C TO COMPUTE DETERMINANT(A) , FOLLOW SSIFA BY DSIDI. C TO COMPUTE INERTIA(A) , FOLLOW SSIFA BY DSIDI. C TO COMPUTE INVERSE(A) , FOLLOW SSIFA BY DSIDI. C C ON ENTRY C C A DOUBLE PRECISION(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 DSISL OR DSIDI 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 SAXPY,SSWAP,ISAMAX C FORTRAN DABS,DMAX1,DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T DOUBLE PRECISION ABSAKK,ALPHA,COLMAX,ROWMAX INTEGER IMAX,IMAXP1,J,JJ,JMAX,K,KM1,KM2,KSTEP,ISAMAX LOGICAL SWAP C C C INITIALIZE C C ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE. ALPHA = (1.0D0 + DSQRT(17.0D0))/8.0D0 C INFO = 0 C C MAIN LOOP ON K, WHICH GOES FROM N TO 1. C K = N 10 CONTINUE C C LEAVE THE LOOP IF K=0 OR K=1. C C ...EXIT IF (K .EQ. 0) GO TO 200 IF (K .GT. 1) GO TO 20 KPVT(1) = 1 IF (A(1,1) .EQ. 0.0D0) INFO = 1 C ......EXIT GO TO 200 20 CONTINUE C C THIS SECTION OF CODE DETERMINES THE KIND OF C ELIMINATION TO BE PERFORMED. WHEN IT IS COMPLETED, C KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND C SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS C REQUIRED. C KM1 = K - 1 ABSAKK = DABS(A(K,K)) C C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN C COLUMN K. C IMAX = ISAMAX(K-1,A(1,K),1) COLMAX = DABS(A(IMAX,K)) IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30 KSTEP = 1 SWAP = .FALSE. GO TO 90 30 CONTINUE C C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN C ROW IMAX. C ROWMAX = 0.0D0 IMAXP1 = IMAX + 1 DO 40 J = IMAXP1, K ROWMAX = DMAX1(ROWMAX,DABS(A(IMAX,J))) 40 CONTINUE IF (IMAX .EQ. 1) GO TO 50 JMAX = ISAMAX(IMAX-1,A(1,IMAX),1) ROWMAX = DMAX1(ROWMAX,DABS(A(JMAX,IMAX))) 50 CONTINUE IF (DABS(A(IMAX,IMAX)) .LT. ALPHA*ROWMAX) GO TO 60 KSTEP = 1 SWAP = .TRUE. GO TO 80 60 CONTINUE IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70 KSTEP = 1 SWAP = .FALSE. GO TO 80 70 CONTINUE KSTEP = 2 SWAP = IMAX .NE. KM1 80 CONTINUE 90 CONTINUE IF (DMAX1(ABSAKK,COLMAX) .NE. 0.0D0) GO TO 100 C C COLUMN K IS ZERO. SET INFO AND ITERATE THE LOOP. C KPVT(K) = K INFO = K GO TO 190 100 CONTINUE IF (KSTEP .EQ. 2) GO TO 140 C C 1 X 1 PIVOT BLOCK. C IF (.NOT.SWAP) GO TO 120 C C PERFORM AN INTERCHANGE. C CALL SSWAP(IMAX,A(1,IMAX),1,A(1,K),1) DO 110 JJ = IMAX, K J = K + IMAX - JJ T = A(J,K) A(J,K) = A(IMAX,J) A(IMAX,J) = T 110 CONTINUE 120 CONTINUE C C PERFORM THE ELIMINATION. C DO 130 JJ = 1, KM1 J = K - JJ MULK = -A(J,K)/A(K,K) T = MULK CALL SAXPY(J,T,A(1,K),1,A(1,J),1) A(J,K) = MULK 130 CONTINUE C C SET THE PIVOT ARRAY. C KPVT(K) = K IF (SWAP) KPVT(K) = IMAX GO TO 190 140 CONTINUE C C 2 X 2 PIVOT BLOCK. C IF (.NOT.SWAP) GO TO 160 C C PERFORM AN INTERCHANGE. C CALL SSWAP(IMAX,A(1,IMAX),1,A(1,K-1),1) DO 150 JJ = IMAX, KM1 J = KM1 + IMAX - JJ T = A(J,K-1) A(J,K-1) = A(IMAX,J) A(IMAX,J) = T 150 CONTINUE T = A(K-1,K) A(K-1,K) = A(IMAX,K) A(IMAX,K) = T 160 CONTINUE C C PERFORM THE ELIMINATION. C KM2 = K - 2 IF (KM2 .EQ. 0) GO TO 180 AK = A(K,K)/A(K-1,K) AKM1 = A(K-1,K-1)/A(K-1,K) DENOM = 1.0D0 - AK*AKM1 DO 170 JJ = 1, KM2 J = KM1 - JJ BK = A(J,K)/A(K-1,K) BKM1 = A(J,K-1)/A(K-1,K) MULK = (AKM1*BK - BKM1)/DENOM MULKM1 = (AK*BKM1 - BK)/DENOM T = MULK CALL SAXPY(J,T,A(1,K),1,A(1,J),1) T = MULKM1 CALL SAXPY(J,T,A(1,K-1),1,A(1,J),1) A(J,K) = MULK A(J,K-1) = MULKM1 170 CONTINUE 180 CONTINUE C C SET THE PIVOT ARRAY. C KPVT(K) = 1 - K IF (SWAP) KPVT(K) = -IMAX KPVT(K-1) = KPVT(K) 190 CONTINUE K = K - KSTEP GO TO 10 200 CONTINUE RETURN END SUBROUTINE SSISL(A,LDA,N,KPVT,B) INTEGER LDA,N,KPVT(1) DOUBLE PRECISION A(LDA,1),B(1) C C DSISL SOLVES THE DOUBLE PRECISION SYMMETRIC SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY SSIFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA,N) C THE OUTPUT FROM SSIFA. 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 SSIFA. C C B DOUBLE PRECISION(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 DSICO HAS SET RCOND .EQ. 0.0 C OR SSIFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL SSIFA(A,LDA,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL DSISL(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 SAXPY,SDOT C FORTRAN IABS C C INTERNAL VARIABLES. C DOUBLE PRECISION AK,AKM1,BK,BKM1,SDOT,DENOM,TEMP INTEGER K,KP C C LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND C D INVERSE TO B. C K = N 10 IF (K .EQ. 0) GO TO 80 IF (KPVT(K) .LT. 0) GO TO 40 C C 1 X 1 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 30 KP = KPVT(K) IF (KP .EQ. K) GO TO 20 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 20 CONTINUE C C APPLY THE TRANSFORMATION. C CALL SAXPY(K-1,B(K),A(1,K),1,B(1),1) 30 CONTINUE C C APPLY D INVERSE. C B(K) = B(K)/A(K,K) K = K - 1 GO TO 70 40 CONTINUE C C 2 X 2 PIVOT BLOCK. C IF (K .EQ. 2) GO TO 60 KP = IABS(KPVT(K)) IF (KP .EQ. K - 1) GO TO 50 C C INTERCHANGE. C TEMP = B(K-1) B(K-1) = B(KP) B(KP) = TEMP 50 CONTINUE C C APPLY THE TRANSFORMATION. C CALL SAXPY(K-2,B(K),A(1,K),1,B(1),1) CALL SAXPY(K-2,B(K-1),A(1,K-1),1,B(1),1) 60 CONTINUE C C APPLY D INVERSE. C AK = A(K,K)/A(K-1,K) AKM1 = A(K-1,K-1)/A(K-1,K) BK = B(K)/A(K-1,K) BKM1 = B(K-1)/A(K-1,K) DENOM = AK*AKM1 - 1.0D0 B(K) = (AKM1*BK - BKM1)/DENOM B(K-1) = (AK*BKM1 - BK)/DENOM K = K - 2 70 CONTINUE GO TO 10 80 CONTINUE C C LOOP FORWARD APPLYING THE TRANSFORMATIONS. C K = 1 90 IF (K .GT. N) GO TO 160 IF (KPVT(K) .LT. 0) GO TO 120 C C 1 X 1 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 110 C C APPLY THE TRANSFORMATION. C B(K) = B(K) + SDOT(K-1,A(1,K),1,B(1),1) KP = KPVT(K) IF (KP .EQ. K) GO TO 100 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 100 CONTINUE 110 CONTINUE K = K + 1 GO TO 150 120 CONTINUE C C 2 X 2 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 140 C C APPLY THE TRANSFORMATION. C B(K) = B(K) + SDOT(K-1,A(1,K),1,B(1),1) B(K+1) = B(K+1) + SDOT(K-1,A(1,K+1),1,B(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 130 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 130 CONTINUE 140 CONTINUE K = K + 2 150 CONTINUE GO TO 90 160 CONTINUE RETURN END SUBROUTINE SSIDI(A,LDA,N,KPVT,DET,INERT,WORK,JOB) INTEGER LDA,N,JOB DOUBLE PRECISION A(LDA,1),WORK(1) DOUBLE PRECISION DET(2) INTEGER KPVT(1),INERT(3) C C DSIDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE C OF A DOUBLE PRECISION SYMMETRIC MATRIX USING THE FACTORS FROM C SSIFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA,N) C THE OUTPUT FROM SSIFA. 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 SSIFA. C C WORK DOUBLE PRECISION(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 DSICO HAS SET RCOND .EQ. 0.0 C OR SSIFA 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 SAXPY,SCOPY,SDOT,SSWAP C FORTRAN DABS,IABS,MOD C C INTERNAL VARIABLES. C DOUBLE PRECISION AKKP1,SDOT,TEMP DOUBLE PRECISION TEN,D,T,AK,AKP1 INTEGER J,JB,K,KM1,KS,KSTEP LOGICAL NOINV,NODET,NOERT C NOINV = MOD(JOB,10) .EQ. 0 NODET = MOD(JOB,100)/10 .EQ. 0 NOERT = MOD(JOB,1000)/100 .EQ. 0 C IF (NODET .AND. NOERT) GO TO 140 IF (NOERT) GO TO 10 INERT(1) = 0 INERT(2) = 0 INERT(3) = 0 10 CONTINUE IF (NODET) GO TO 20 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 20 CONTINUE T = 0.0D0 DO 130 K = 1, N D = A(K,K) C C CHECK IF 1 BY 1 C IF (KPVT(K) .GT. 0) GO TO 50 C C 2 BY 2 BLOCK C USE DET (D S) = (D/T * C - T) * T , T = DABS(S) C (S C) C TO AVOID UNDERFLOW/OVERFLOW TROUBLES. C TAKE TWO PASSES THROUGH SCALING. USE T FOR FLAG. C IF (T .NE. 0.0D0) GO TO 30 T = DABS(A(K,K+1)) D = (D/T)*A(K+1,K+1) - T GO TO 40 30 CONTINUE D = T T = 0.0D0 40 CONTINUE 50 CONTINUE C IF (NOERT) GO TO 60 IF (D .GT. 0.0D0) INERT(1) = INERT(1) + 1 IF (D .LT. 0.0D0) INERT(2) = INERT(2) + 1 IF (D .EQ. 0.0D0) INERT(3) = INERT(3) + 1 60 CONTINUE C IF (NODET) GO TO 120 DET(1) = D*DET(1) IF (DET(1) .EQ. 0.0D0) GO TO 110 70 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 80 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 70 80 CONTINUE 90 IF (DABS(DET(1)) .LT. TEN) GO TO 100 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 90 100 CONTINUE 110 CONTINUE 120 CONTINUE 130 CONTINUE 140 CONTINUE C C COMPUTE INVERSE(A) C IF (NOINV) GO TO 270 K = 1 150 IF (K .GT. N) GO TO 260 KM1 = K - 1 IF (KPVT(K) .LT. 0) GO TO 180 C C 1 BY 1 C A(K,K) = 1.0D0/A(K,K) IF (KM1 .LT. 1) GO TO 170 CALL SCOPY(KM1,A(1,K),1,WORK,1) DO 160 J = 1, KM1 A(J,K) = SDOT(J,A(1,J),1,WORK,1) CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1) 160 CONTINUE A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1) 170 CONTINUE KSTEP = 1 GO TO 220 180 CONTINUE C C 2 BY 2 C T = DABS(A(K,K+1)) AK = A(K,K)/T AKP1 = A(K+1,K+1)/T AKKP1 = A(K,K+1)/T D = T*(AK*AKP1 - 1.0D0) A(K,K) = AKP1/D A(K+1,K+1) = AK/D A(K,K+1) = -AKKP1/D IF (KM1 .LT. 1) GO TO 210 CALL SCOPY(KM1,A(1,K+1),1,WORK,1) DO 190 J = 1, KM1 A(J,K+1) = SDOT(J,A(1,J),1,WORK,1) CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K+1),1) 190 CONTINUE A(K+1,K+1) = A(K+1,K+1) + SDOT(KM1,WORK,1,A(1,K+1),1) A(K,K+1) = A(K,K+1) + SDOT(KM1,A(1,K),1,A(1,K+1),1) CALL SCOPY(KM1,A(1,K),1,WORK,1) DO 200 J = 1, KM1 A(J,K) = SDOT(J,A(1,J),1,WORK,1) CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1) 200 CONTINUE A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1) 210 CONTINUE KSTEP = 2 220 CONTINUE C C SWAP C KS = IABS(KPVT(K)) IF (KS .EQ. K) GO TO 250 CALL SSWAP(KS,A(1,KS),1,A(1,K),1) DO 230 JB = KS, K J = K + KS - JB TEMP = A(J,K) A(J,K) = A(KS,J) A(KS,J) = TEMP 230 CONTINUE IF (KSTEP .EQ. 1) GO TO 240 TEMP = A(KS,K+1) A(KS,K+1) = A(K,K+1) A(K,K+1) = TEMP 240 CONTINUE 250 CONTINUE K = K + KSTEP GO TO 150 260 CONTINUE 270 CONTINUE RETURN END SUBROUTINE SSPCO(AP,N,KPVT,RCOND,Z) INTEGER N,KPVT(1) DOUBLE PRECISION AP(1),Z(1) DOUBLE PRECISION RCOND C C DSPCO FACTORS A DOUBLE PRECISION SYMMETRIC MATRIX STORED IN C PACKED FORM BY ELIMINATION WITH SYMMETRIC PIVOTING AND ESTIMATES C THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SSPFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DSPCO BY DSPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW DSPCO BY DSPSL. C TO COMPUTE INVERSE(A) , FOLLOW DSPCO BY DSPDI. C TO COMPUTE DETERMINANT(A) , FOLLOW DSPCO BY DSPDI. C TO COMPUTE INERTIA(A), FOLLOW DSPCO BY DSPDI. C C ON ENTRY C C AP DOUBLE PRECISION (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 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*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 DOUBLE PRECISION(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 SSPFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN DABS,DMAX1,IABS,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T DOUBLE PRECISION ANORM,S,SASUM,YNORM INTEGER I,IJ,IK,IKM1,IKP1,INFO,J,JM1,J1 INTEGER K,KK,KM1K,KM1KM1,KP,KPS,KS C C C FIND NORM OF A USING ONLY UPPER HALF C J1 = 1 DO 30 J = 1, N Z(J) = SASUM(J,AP(J1),1) IJ = J1 J1 = J1 + J JM1 = J - 1 IF (JM1 .LT. 1) GO TO 20 DO 10 I = 1, JM1 Z(I) = Z(I) + DABS(AP(IJ)) IJ = IJ + 1 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0D0 DO 40 J = 1, N ANORM = DMAX1(ANORM,Z(J)) 40 CONTINUE C C FACTOR C CALL SSPFA(AP,N,KPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE U*D*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE U*D*W = E C EK = 1.0D0 DO 50 J = 1, N Z(J) = 0.0D0 50 CONTINUE K = N IK = (N*(N - 1))/2 60 IF (K .EQ. 0) GO TO 120 KK = IK + K IKM1 = IK - (K - 1) KS = 1 IF (KPVT(K) .LT. 0) KS = 2 KP = IABS(KPVT(K)) KPS = K + 1 - KS IF (KP .EQ. KPS) GO TO 70 T = Z(KPS) Z(KPS) = Z(KP) Z(KP) = T 70 CONTINUE IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,Z(K)) Z(K) = Z(K) + EK CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1) IF (KS .EQ. 1) GO TO 80 IF (Z(K-1) .NE. 0.0D0) EK = DSIGN(EK,Z(K-1)) Z(K-1) = Z(K-1) + EK CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1) 80 CONTINUE IF (KS .EQ. 2) GO TO 100 IF (DABS(Z(K)) .LE. DABS(AP(KK))) GO TO 90 S = DABS(AP(KK))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 90 CONTINUE IF (AP(KK) .NE. 0.0D0) Z(K) = Z(K)/AP(KK) IF (AP(KK) .EQ. 0.0D0) Z(K) = 1.0D0 GO TO 110 100 CONTINUE KM1K = IK + K - 1 KM1KM1 = IKM1 + K - 1 AK = AP(KK)/AP(KM1K) AKM1 = AP(KM1KM1)/AP(KM1K) BK = Z(K)/AP(KM1K) BKM1 = Z(K-1)/AP(KM1K) DENOM = AK*AKM1 - 1.0D0 Z(K) = (AKM1*BK - BKM1)/DENOM Z(K-1) = (AK*BKM1 - BK)/DENOM 110 CONTINUE K = K - KS IK = IK - K IF (KS .EQ. 2) IK = IK - (K + 1) GO TO 60 120 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(U)*Y = W C K = 1 IK = 0 130 IF (K .GT. N) GO TO 160 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. 1) GO TO 150 Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1) IKP1 = IK + K IF (KS .EQ. 2) * Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 140 T = Z(K) Z(K) = Z(KP) Z(KP) = T 140 CONTINUE 150 CONTINUE IK = IK + K IF (KS .EQ. 2) IK = IK + (K + 1) K = K + KS GO TO 130 160 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE U*D*V = Y C K = N IK = N*(N - 1)/2 170 IF (K .EQ. 0) GO TO 230 KK = IK + K IKM1 = IK - (K - 1) KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. KS) GO TO 190 KP = IABS(KPVT(K)) KPS = K + 1 - KS IF (KP .EQ. KPS) GO TO 180 T = Z(KPS) Z(KPS) = Z(KP) Z(KP) = T 180 CONTINUE CALL SAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1) IF (KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1) 190 CONTINUE IF (KS .EQ. 2) GO TO 210 IF (DABS(Z(K)) .LE. DABS(AP(KK))) GO TO 200 S = DABS(AP(KK))/DABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 200 CONTINUE IF (AP(KK) .NE. 0.0D0) Z(K) = Z(K)/AP(KK) IF (AP(KK) .EQ. 0.0D0) Z(K) = 1.0D0 GO TO 220 210 CONTINUE KM1K = IK + K - 1 KM1KM1 = IKM1 + K - 1 AK = AP(KK)/AP(KM1K) AKM1 = AP(KM1KM1)/AP(KM1K) BK = Z(K)/AP(KM1K) BKM1 = Z(K-1)/AP(KM1K) DENOM = AK*AKM1 - 1.0D0 Z(K) = (AKM1*BK - BKM1)/DENOM Z(K-1) = (AK*BKM1 - BK)/DENOM 220 CONTINUE K = K - KS IK = IK - K IF (KS .EQ. 2) IK = IK - (K + 1) GO TO 170 230 CONTINUE S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE TRANS(U)*Z = V C K = 1 IK = 0 240 IF (K .GT. N) GO TO 270 KS = 1 IF (KPVT(K) .LT. 0) KS = 2 IF (K .EQ. 1) GO TO 260 Z(K) = Z(K) + SDOT(K-1,AP(IK+1),1,Z(1),1) IKP1 = IK + K IF (KS .EQ. 2) * Z(K+1) = Z(K+1) + SDOT(K-1,AP(IKP1+1),1,Z(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 250 T = Z(K) Z(K) = Z(KP) Z(KP) = T 250 CONTINUE 260 CONTINUE IK = IK + K IF (KS .EQ. 2) IK = IK + (K + 1) K = K + KS GO TO 240 270 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END SUBROUTINE SSPFA(AP,N,KPVT,INFO) INTEGER N,KPVT(1),INFO DOUBLE PRECISION AP(1) C C SSPFA FACTORS A DOUBLE PRECISION SYMMETRIC MATRIX STORED IN C PACKED FORM BY ELIMINATION WITH SYMMETRIC PIVOTING. C C TO SOLVE A*X = B , FOLLOW SSPFA BY DSPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW SSPFA BY DSPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW SSPFA BY DSPDI. C TO COMPUTE INERTIA(A) , FOLLOW SSPFA BY DSPDI. C TO COMPUTE INVERSE(A) , FOLLOW SSPFA BY DSPDI. C C ON ENTRY C C AP DOUBLE PRECISION (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 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*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 DSPSL OR DSPDI 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 SAXPY,SSWAP,ISAMAX C FORTRAN DABS,DMAX1,DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T DOUBLE PRECISION ABSAKK,ALPHA,COLMAX,ROWMAX INTEGER ISAMAX,IJ,IJJ,IK,IKM1,IM,IMAX,IMAXP1,IMIM,IMJ,IMK INTEGER J,JJ,JK,JKM1,JMAX,JMIM,K,KK,KM1,KM1K,KM1KM1,KM2,KSTEP LOGICAL SWAP C C C INITIALIZE C C ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE. ALPHA = (1.0D0 + DSQRT(17.0D0))/8.0D0 C INFO = 0 C C MAIN LOOP ON K, WHICH GOES FROM N TO 1. C K = N IK = (N*(N - 1))/2 10 CONTINUE C C LEAVE THE LOOP IF K=0 OR K=1. C C ...EXIT IF (K .EQ. 0) GO TO 200 IF (K .GT. 1) GO TO 20 KPVT(1) = 1 IF (AP(1) .EQ. 0.0D0) INFO = 1 C ......EXIT GO TO 200 20 CONTINUE C C THIS SECTION OF CODE DETERMINES THE KIND OF C ELIMINATION TO BE PERFORMED. WHEN IT IS COMPLETED, C KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND C SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS C REQUIRED. C KM1 = K - 1 KK = IK + K ABSAKK = DABS(AP(KK)) C C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN C COLUMN K. C IMAX = ISAMAX(K-1,AP(IK+1),1) IMK = IK + IMAX COLMAX = DABS(AP(IMK)) IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30 KSTEP = 1 SWAP = .FALSE. GO TO 90 30 CONTINUE C C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN C ROW IMAX. C ROWMAX = 0.0D0 IMAXP1 = IMAX + 1 IM = IMAX*(IMAX - 1)/2 IMJ = IM + 2*IMAX DO 40 J = IMAXP1, K ROWMAX = DMAX1(ROWMAX,DABS(AP(IMJ))) IMJ = IMJ + J 40 CONTINUE IF (IMAX .EQ. 1) GO TO 50 JMAX = ISAMAX(IMAX-1,AP(IM+1),1) JMIM = JMAX + IM ROWMAX = DMAX1(ROWMAX,DABS(AP(JMIM))) 50 CONTINUE IMIM = IMAX + IM IF (DABS(AP(IMIM)) .LT. ALPHA*ROWMAX) GO TO 60 KSTEP = 1 SWAP = .TRUE. GO TO 80 60 CONTINUE IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70 KSTEP = 1 SWAP = .FALSE. GO TO 80 70 CONTINUE KSTEP = 2 SWAP = IMAX .NE. KM1 80 CONTINUE 90 CONTINUE IF (DMAX1(ABSAKK,COLMAX) .NE. 0.0D0) GO TO 100 C C COLUMN K IS ZERO. SET INFO AND ITERATE THE LOOP. C KPVT(K) = K INFO = K GO TO 190 100 CONTINUE IF (KSTEP .EQ. 2) GO TO 140 C C 1 X 1 PIVOT BLOCK. C IF (.NOT.SWAP) GO TO 120 C C PERFORM AN INTERCHANGE. C CALL SSWAP(IMAX,AP(IM+1),1,AP(IK+1),1) IMJ = IK + IMAX DO 110 JJ = IMAX, K J = K + IMAX - JJ JK = IK + J T = AP(JK) AP(JK) = AP(IMJ) AP(IMJ) = T IMJ = IMJ - (J - 1) 110 CONTINUE 120 CONTINUE C C PERFORM THE ELIMINATION. C IJ = IK - (K - 1) DO 130 JJ = 1, KM1 J = K - JJ JK = IK + J MULK = -AP(JK)/AP(KK) T = MULK CALL SAXPY(J,T,AP(IK+1),1,AP(IJ+1),1) IJJ = IJ + J AP(JK) = MULK IJ = IJ - (J - 1) 130 CONTINUE C C SET THE PIVOT ARRAY. C KPVT(K) = K IF (SWAP) KPVT(K) = IMAX GO TO 190 140 CONTINUE C C 2 X 2 PIVOT BLOCK. C KM1K = IK + K - 1 IKM1 = IK - (K - 1) IF (.NOT.SWAP) GO TO 160 C C PERFORM AN INTERCHANGE. C CALL SSWAP(IMAX,AP(IM+1),1,AP(IKM1+1),1) IMJ = IKM1 + IMAX DO 150 JJ = IMAX, KM1 J = KM1 + IMAX - JJ JKM1 = IKM1 + J T = AP(JKM1) AP(JKM1) = AP(IMJ) AP(IMJ) = T IMJ = IMJ - (J - 1) 150 CONTINUE T = AP(KM1K) AP(KM1K) = AP(IMK) AP(IMK) = T 160 CONTINUE C C PERFORM THE ELIMINATION. C KM2 = K - 2 IF (KM2 .EQ. 0) GO TO 180 AK = AP(KK)/AP(KM1K) KM1KM1 = IKM1 + K - 1 AKM1 = AP(KM1KM1)/AP(KM1K) DENOM = 1.0D0 - AK*AKM1 IJ = IK - (K - 1) - (K - 2) DO 170 JJ = 1, KM2 J = KM1 - JJ JK = IK + J BK = AP(JK)/AP(KM1K) JKM1 = IKM1 + J BKM1 = AP(JKM1)/AP(KM1K) MULK = (AKM1*BK - BKM1)/DENOM MULKM1 = (AK*BKM1 - BK)/DENOM T = MULK CALL SAXPY(J,T,AP(IK+1),1,AP(IJ+1),1) T = MULKM1 CALL SAXPY(J,T,AP(IKM1+1),1,AP(IJ+1),1) AP(JK) = MULK AP(JKM1) = MULKM1 IJJ = IJ + J IJ = IJ - (J - 1) 170 CONTINUE 180 CONTINUE C C SET THE PIVOT ARRAY. C KPVT(K) = 1 - K IF (SWAP) KPVT(K) = -IMAX KPVT(K-1) = KPVT(K) 190 CONTINUE IK = IK - (K - 1) IF (KSTEP .EQ. 2) IK = IK - (K - 2) K = K - KSTEP GO TO 10 200 CONTINUE RETURN END SUBROUTINE SSPSL(AP,N,KPVT,B) INTEGER N,KPVT(1) DOUBLE PRECISION AP(1),B(1) C C DSISL SOLVES THE DOUBLE PRECISION SYMMETRIC SYSTEM C A * X = B C USING THE FACTORS COMPUTED BY SSPFA. C C ON ENTRY C C AP DOUBLE PRECISION(N*(N+1)/2) C THE OUTPUT FROM SSPFA. C C N INTEGER C THE ORDER OF THE MATRIX A . C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM SSPFA. C C B DOUBLE PRECISION(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 DSPCO HAS SET RCOND .EQ. 0.0 C OR SSPFA HAS SET INFO .NE. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL SSPFA(AP,N,KPVT,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL DSPSL(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 SAXPY,SDOT C FORTRAN IABS C C INTERNAL VARIABLES. C DOUBLE PRECISION AK,AKM1,BK,BKM1,SDOT,DENOM,TEMP INTEGER IK,IKM1,IKP1,K,KK,KM1K,KM1KM1,KP C C LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND C D INVERSE TO B. C K = N IK = (N*(N - 1))/2 10 IF (K .EQ. 0) GO TO 80 KK = IK + K IF (KPVT(K) .LT. 0) GO TO 40 C C 1 X 1 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 30 KP = KPVT(K) IF (KP .EQ. K) GO TO 20 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 20 CONTINUE C C APPLY THE TRANSFORMATION. C CALL SAXPY(K-1,B(K),AP(IK+1),1,B(1),1) 30 CONTINUE C C APPLY D INVERSE. C B(K) = B(K)/AP(KK) K = K - 1 IK = IK - K GO TO 70 40 CONTINUE C C 2 X 2 PIVOT BLOCK. C IKM1 = IK - (K - 1) IF (K .EQ. 2) GO TO 60 KP = IABS(KPVT(K)) IF (KP .EQ. K - 1) GO TO 50 C C INTERCHANGE. C TEMP = B(K-1) B(K-1) = B(KP) B(KP) = TEMP 50 CONTINUE C C APPLY THE TRANSFORMATION. C CALL SAXPY(K-2,B(K),AP(IK+1),1,B(1),1) CALL SAXPY(K-2,B(K-1),AP(IKM1+1),1,B(1),1) 60 CONTINUE C C APPLY D INVERSE. C KM1K = IK + K - 1 KK = IK + K AK = AP(KK)/AP(KM1K) KM1KM1 = IKM1 + K - 1 AKM1 = AP(KM1KM1)/AP(KM1K) BK = B(K)/AP(KM1K) BKM1 = B(K-1)/AP(KM1K) DENOM = AK*AKM1 - 1.0D0 B(K) = (AKM1*BK - BKM1)/DENOM B(K-1) = (AK*BKM1 - BK)/DENOM K = K - 2 IK = IK - (K + 1) - K 70 CONTINUE GO TO 10 80 CONTINUE C C LOOP FORWARD APPLYING THE TRANSFORMATIONS. C K = 1 IK = 0 90 IF (K .GT. N) GO TO 160 IF (KPVT(K) .LT. 0) GO TO 120 C C 1 X 1 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 110 C C APPLY THE TRANSFORMATION. C B(K) = B(K) + SDOT(K-1,AP(IK+1),1,B(1),1) KP = KPVT(K) IF (KP .EQ. K) GO TO 100 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 100 CONTINUE 110 CONTINUE IK = IK + K K = K + 1 GO TO 150 120 CONTINUE C C 2 X 2 PIVOT BLOCK. C IF (K .EQ. 1) GO TO 140 C C APPLY THE TRANSFORMATION. C B(K) = B(K) + SDOT(K-1,AP(IK+1),1,B(1),1) IKP1 = IK + K B(K+1) = B(K+1) + SDOT(K-1,AP(IKP1+1),1,B(1),1) KP = IABS(KPVT(K)) IF (KP .EQ. K) GO TO 130 C C INTERCHANGE. C TEMP = B(K) B(K) = B(KP) B(KP) = TEMP 130 CONTINUE 140 CONTINUE IK = IK + K + K + 1 K = K + 2 150 CONTINUE GO TO 90 160 CONTINUE RETURN END SUBROUTINE SSPDI(AP,N,KPVT,DET,INERT,WORK,JOB) INTEGER N,JOB DOUBLE PRECISION AP(1),WORK(1) DOUBLE PRECISION DET(2) INTEGER KPVT(1),INERT(3) C C DSPDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE C OF A DOUBLE PRECISION SYMMETRIC MATRIX USING THE FACTORS FROM C SSPFA, WHERE THE MATRIX IS STORED IN PACKED FORM. C C ON ENTRY C C AP DOUBLE PRECISION (N*(N+1)/2) C THE OUTPUT FROM SSPFA. C C N INTEGER C THE ORDER OF THE MATRIX A. C C KPVT INTEGER(N) C THE PIVOT VECTOR FROM SSPFA. C C WORK DOUBLE PRECISION(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 DSPCO HAS SET RCOND .EQ. 0.0 C OR SSPFA 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 SAXPY,SCOPY,SDOT,SSWAP C FORTRAN DABS,IABS,MOD C C INTERNAL VARIABLES. C DOUBLE PRECISION AKKP1,SDOT,TEMP DOUBLE PRECISION TEN,D,T,AK,AKP1 INTEGER IJ,IK,IKP1,IKS,J,JB,JK,JKP1 INTEGER K,KK,KKP1,KM1,KS,KSJ,KSKP1,KSTEP LOGICAL NOINV,NODET,NOERT C NOINV = MOD(JOB,10) .EQ. 0 NODET = MOD(JOB,100)/10 .EQ. 0 NOERT = MOD(JOB,1000)/100 .EQ. 0 C IF (NODET .AND. NOERT) GO TO 140 IF (NOERT) GO TO 10 INERT(1) = 0 INERT(2) = 0 INERT(3) = 0 10 CONTINUE IF (NODET) GO TO 20 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 20 CONTINUE T = 0.0D0 IK = 0 DO 130 K = 1, N KK = IK + K D = AP(KK) C C CHECK IF 1 BY 1 C IF (KPVT(K) .GT. 0) GO TO 50 C C 2 BY 2 BLOCK C USE DET (D S) = (D/T * C - T) * T , T = DABS(S) C (S C) C TO AVOID UNDERFLOW/OVERFLOW TROUBLES. C TAKE TWO PASSES THROUGH SCALING. USE T FOR FLAG. C IF (T .NE. 0.0D0) GO TO 30 IKP1 = IK + K KKP1 = IKP1 + K T = DABS(AP(KKP1)) D = (D/T)*AP(KKP1+1) - T GO TO 40 30 CONTINUE D = T T = 0.0D0 40 CONTINUE 50 CONTINUE C IF (NOERT) GO TO 60 IF (D .GT. 0.0D0) INERT(1) = INERT(1) + 1 IF (D .LT. 0.0D0) INERT(2) = INERT(2) + 1 IF (D .EQ. 0.0D0) INERT(3) = INERT(3) + 1 60 CONTINUE C IF (NODET) GO TO 120 DET(1) = D*DET(1) IF (DET(1) .EQ. 0.0D0) GO TO 110 70 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 80 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 70 80 CONTINUE 90 IF (DABS(DET(1)) .LT. TEN) GO TO 100 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 90 100 CONTINUE 110 CONTINUE 120 CONTINUE IK = IK + K 130 CONTINUE 140 CONTINUE C C COMPUTE INVERSE(A) C IF (NOINV) GO TO 270 K = 1 IK = 0 150 IF (K .GT. N) GO TO 260 KM1 = K - 1 KK = IK + K IKP1 = IK + K KKP1 = IKP1 + K IF (KPVT(K) .LT. 0) GO TO 180 C C 1 BY 1 C AP(KK) = 1.0D0/AP(KK) IF (KM1 .LT. 1) GO TO 170 CALL SCOPY(KM1,AP(IK+1),1,WORK,1) IJ = 0 DO 160 J = 1, KM1 JK = IK + J AP(JK) = SDOT(J,AP(IJ+1),1,WORK,1) CALL SAXPY(J-1,WORK(J),AP(IJ+1),1,AP(IK+1),1) IJ = IJ + J 160 CONTINUE AP(KK) = AP(KK) + SDOT(KM1,WORK,1,AP(IK+1),1) 170 CONTINUE KSTEP = 1 GO TO 220 180 CONTINUE C C 2 BY 2 C T = DABS(AP(KKP1)) AK = AP(KK)/T AKP1 = AP(KKP1+1)/T AKKP1 = AP(KKP1)/T D = T*(AK*AKP1 - 1.0D0) AP(KK) = AKP1/D AP(KKP1+1) = AK/D AP(KKP1) = -AKKP1/D IF (KM1 .LT. 1) GO TO 210 CALL SCOPY(KM1,AP(IKP1+1),1,WORK,1) IJ = 0 DO 190 J = 1, KM1 JKP1 = IKP1 + J AP(JKP1) = SDOT(J,AP(IJ+1),1,WORK,1) CALL SAXPY(J-1,WORK(J),AP(IJ+1),1,AP(IKP1+1),1) IJ = IJ + J 190 CONTINUE AP(KKP1+1) = AP(KKP1+1) * + SDOT(KM1,WORK,1,AP(IKP1+1),1) AP(KKP1) = AP(KKP1) * + SDOT(KM1,AP(IK+1),1,AP(IKP1+1),1) CALL SCOPY(KM1,AP(IK+1),1,WORK,1) IJ = 0 DO 200 J = 1, KM1 JK = IK + J AP(JK) = SDOT(J,AP(IJ+1),1,WORK,1) CALL SAXPY(J-1,WORK(J),AP(IJ+1),1,AP(IK+1),1) IJ = IJ + J 200 CONTINUE AP(KK) = AP(KK) + SDOT(KM1,WORK,1,AP(IK+1),1) 210 CONTINUE KSTEP = 2 220 CONTINUE C C SWAP C KS = IABS(KPVT(K)) IF (KS .EQ. K) GO TO 250 IKS = (KS*(KS - 1))/2 CALL SSWAP(KS,AP(IKS+1),1,AP(IK+1),1) KSJ = IK + KS DO 230 JB = KS, K J = K + KS - JB JK = IK + J TEMP = AP(JK) AP(JK) = AP(KSJ) AP(KSJ) = TEMP KSJ = KSJ - (J - 1) 230 CONTINUE IF (KSTEP .EQ. 1) GO TO 240 KSKP1 = IKP1 + KS TEMP = AP(KSKP1) AP(KSKP1) = AP(KKP1) AP(KKP1) = TEMP 240 CONTINUE 250 CONTINUE IK = IK + K IF (KSTEP .EQ. 2) IK = IK + K + 1 K = K + KSTEP GO TO 150 260 CONTINUE 270 CONTINUE RETURN END SUBROUTINE STRCO(T,LDT,N,RCOND,Z,JOB) INTEGER LDT,N,JOB DOUBLE PRECISION T(LDT,1),Z(1) DOUBLE PRECISION RCOND C C DTRCO ESTIMATES THE CONDITION OF A DOUBLE PRECISION TRIANGULAR C MATRIX. C C ON ENTRY C C T DOUBLE PRECISION(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