SUBROUTINE SPPCO(AP,N,RCOND,Z,INFO) INTEGER N,INFO REAL AP(1),Z(1) REAL RCOND C C SPPCO FACTORS A REAL 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 SPPCO BY SPPSL. C TO COMPUTE INVERSE(A)*C , FOLLOW SPPCO BY SPPSL. C TO COMPUTE DETERMINANT(A) , FOLLOW SPPCO BY SPPDI. C TO COMPUTE INVERSE(A) , FOLLOW SPPCO BY SPPDI. C C ON ENTRY C C AP REAL (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 REAL 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 REAL(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 ABS,AMAX1,REAL,SIGN C C INTERNAL VARIABLES C REAL SDOT,EK,T,WK,WKM REAL 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) + ABS(AP(IJ)) IJ = IJ + 1 10 CONTINUE 20 CONTINUE 30 CONTINUE ANORM = 0.0E0 DO 40 J = 1, N ANORM = AMAX1(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.0E0 DO 50 J = 1, N Z(J) = 0.0E0 50 CONTINUE KK = 0 DO 110 K = 1, N KK = KK + K IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K)) IF (ABS(EK-Z(K)) .LE. AP(KK)) GO TO 60 S = AP(KK)/ABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 60 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = ABS(WK) SM = ABS(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 + ABS(Z(J)+WKM*AP(KJ)) Z(J) = Z(J) + WK*AP(KJ) S = S + ABS(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.0E0/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 (ABS(Z(K)) .LE. AP(KK)) GO TO 120 S = AP(KK)/ABS(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.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0E0 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 (ABS(Z(K)) .LE. AP(KK)) GO TO 140 S = AP(KK)/ABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 140 CONTINUE Z(K) = Z(K)/AP(KK) 150 CONTINUE S = 1.0E0/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 (ABS(Z(K)) .LE. AP(KK)) GO TO 160 S = AP(KK)/ABS(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.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0 180 CONTINUE RETURN END