1 SUBROUTINE CHPCO(AP,N,KPVT,RCOND,Z) 2C***BEGIN PROLOGUE CHPCO 3C***DATE WRITTEN 780814 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C***CATEGORY NO. D2D1A 8C***KEYWORDS COMPLEX,CONDITION,FACTOR,HERMITIAN,LINEAR ALGEBRA,LINPACK, 9C MATRIX,PACKED 10C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) 11C***PURPOSE Factors a COMPLEX HERMITIAN matrix (packed form) by 12C elimination with symmetric pivoting and estimates the 13C condition of the matrix. 14C***DESCRIPTION 15C 16C CHPCO factors a complex Hermitian matrix stored in packed 17C form by elimination with symmetric pivoting and estimates 18C the condition of the matrix. 19C 20C if RCOND is not needed, CHPFA is slightly faster. 21C To solve A*X = B , follow CHPCO by CHPSL. 22C To compute INVERSE(A)*C , follow CHPCO by CHPSL. 23C To compute INVERSE(A) , follow CHPCO by CHPDI. 24C To compute DETERMINANT(A) , follow CHPCO by CHPDI. 25C To compute INERTIA(A), follow CHPCO by CHPDI. 26C 27C On Entry 28C 29C AP COMPLEX (N*(N+1)/2) 30C the packed form of a Hermitian matrix A . The 31C columns of the upper triangle are stored sequentially 32C in a one-dimensional array of length N*(N+1)/2 . 33C See comments below for details. 34C 35C N INTEGER 36C the order of the matrix A . 37C 38C Output 39C 40C AP a block diagonal matrix and the multipliers which 41C were used to obtain it stored in packed form. 42C The factorization can be written A = U*D*CTRANS(U) 43C where U is a product of permutation and unit 44C upper triangular matrices , CTRANS(U) is the 45C conjugate transpose of U , and D is block diagonal 46C with 1 by 1 and 2 by 2 blocks. 47C 48C KVPT INTEGER(N) 49C an integer vector of pivot indices. 50C 51C RCOND REAL 52C an estimate of the reciprocal condition of A . 53C For the system A*X = B , relative perturbations 54C in A and B of size EPSILON may cause 55C relative perturbations in X of size EPSILON/RCOND . 56C If RCOND is so small that the logical expression 57C 1.0 + RCOND .EQ. 1.0 58C is true, then A may be singular to working 59C precision. In particular, RCOND is zero if 60C exact singularity is detected or the estimate 61C underflows. 62C 63C Z COMPLEX(N) 64C a work vector whose contents are usually unimportant. 65C If A is close to a singular matrix, then Z is 66C an approximate null vector in the sense that 67C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . 68C 69C Packed Storage 70C 71C The following program segment will pack the upper 72C triangle of a Hermitian matrix. 73C 74C K = 0 75C DO 20 J = 1, N 76C DO 10 I = 1, J 77C K = K + 1 78C AP(K) = A(I,J) 79C 10 CONTINUE 80C 20 CONTINUE 81C 82C LINPACK. This version dated 08/14/78 . 83C Cleve Moler, University of New Mexico, Argonne National Lab. 84C 85C Subroutines and Functions 86C 87C LINPACK CHPFA 88C BLAS CAXPY,CDOTC,CSSCAL,SCASUM 89C Fortran ABS,AIMAG,AMAX1,CMPLX,CONJG,IABS,REAL 90C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., 91C *LINPACK USERS GUIDE*, SIAM, 1979. 92C***ROUTINES CALLED CAXPY,CDOTC,CHPFA,CSSCAL,SCASUM 93C***END PROLOGUE CHPCO 94 INTEGER N,KPVT(*) 95 COMPLEX AP(*),Z(*) 96 REAL RCOND 97C 98 COMPLEX AK,AKM1,BK,BKM1,CDOTC,DENOM,EK,T 99 REAL ANORM,S,SCASUM,YNORM 100 INTEGER I,IJ,IK,IKM1,IKP1,INFO,J,JM1,J1 101 INTEGER K,KK,KM1K,KM1KM1,KP,KPS,KS 102 COMPLEX ZDUM,ZDUM2,CSIGN1 103 REAL CABS1 104 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) 105 CSIGN1(ZDUM,ZDUM2) = CABS1(ZDUM)*(ZDUM2/CABS1(ZDUM2)) 106C 107C FIND NORM OF A USING ONLY UPPER HALF 108C 109C***FIRST EXECUTABLE STATEMENT CHPCO 110 J1 = 1 111 DO 30 J = 1, N 112 Z(J) = CMPLX(SCASUM(J,AP(J1),1),0.0E0) 113 IJ = J1 114 J1 = J1 + J 115 JM1 = J - 1 116 IF (JM1 .LT. 1) GO TO 20 117 DO 10 I = 1, JM1 118 Z(I) = CMPLX(REAL(Z(I))+CABS1(AP(IJ)),0.0E0) 119 IJ = IJ + 1 120 10 CONTINUE 121 20 CONTINUE 122 30 CONTINUE 123 ANORM = 0.0E0 124 DO 40 J = 1, N 125 ANORM = AMAX1(ANORM,REAL(Z(J))) 126 40 CONTINUE 127C 128C FACTOR 129C 130 CALL CHPFA(AP,N,KPVT,INFO) 131C 132C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . 133C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . 134C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL 135C GROWTH IN THE ELEMENTS OF W WHERE U*D*W = E . 136C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 137C 138C SOLVE U*D*W = E 139C 140 EK = (1.0E0,0.0E0) 141 DO 50 J = 1, N 142 Z(J) = (0.0E0,0.0E0) 143 50 CONTINUE 144 K = N 145 IK = (N*(N - 1))/2 146 60 IF (K .EQ. 0) GO TO 120 147 KK = IK + K 148 IKM1 = IK - (K - 1) 149 KS = 1 150 IF (KPVT(K) .LT. 0) KS = 2 151 KP = IABS(KPVT(K)) 152 KPS = K + 1 - KS 153 IF (KP .EQ. KPS) GO TO 70 154 T = Z(KPS) 155 Z(KPS) = Z(KP) 156 Z(KP) = T 157 70 CONTINUE 158 IF (CABS1(Z(K)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K)) 159 Z(K) = Z(K) + EK 160 CALL CAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1) 161 IF (KS .EQ. 1) GO TO 80 162 IF (CABS1(Z(K-1)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K-1)) 163 Z(K-1) = Z(K-1) + EK 164 CALL CAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1) 165 80 CONTINUE 166 IF (KS .EQ. 2) GO TO 100 167 IF (CABS1(Z(K)) .LE. CABS1(AP(KK))) GO TO 90 168 S = CABS1(AP(KK))/CABS1(Z(K)) 169 CALL CSSCAL(N,S,Z,1) 170 EK = CMPLX(S,0.0E0)*EK 171 90 CONTINUE 172 IF (CABS1(AP(KK)) .NE. 0.0E0) Z(K) = Z(K)/AP(KK) 173 IF (CABS1(AP(KK)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0) 174 GO TO 110 175 100 CONTINUE 176 KM1K = IK + K - 1 177 KM1KM1 = IKM1 + K - 1 178 AK = AP(KK)/CONJG(AP(KM1K)) 179 AKM1 = AP(KM1KM1)/AP(KM1K) 180 BK = Z(K)/CONJG(AP(KM1K)) 181 BKM1 = Z(K-1)/AP(KM1K) 182 DENOM = AK*AKM1 - 1.0E0 183 Z(K) = (AKM1*BK - BKM1)/DENOM 184 Z(K-1) = (AK*BKM1 - BK)/DENOM 185 110 CONTINUE 186 K = K - KS 187 IK = IK - K 188 IF (KS .EQ. 2) IK = IK - (K + 1) 189 GO TO 60 190 120 CONTINUE 191 S = 1.0E0/SCASUM(N,Z,1) 192 CALL CSSCAL(N,S,Z,1) 193C 194C SOLVE CTRANS(U)*Y = W 195C 196 K = 1 197 IK = 0 198 130 IF (K .GT. N) GO TO 160 199 KS = 1 200 IF (KPVT(K) .LT. 0) KS = 2 201 IF (K .EQ. 1) GO TO 150 202 Z(K) = Z(K) + CDOTC(K-1,AP(IK+1),1,Z(1),1) 203 IKP1 = IK + K 204 IF (KS .EQ. 2) 205 1 Z(K+1) = Z(K+1) + CDOTC(K-1,AP(IKP1+1),1,Z(1),1) 206 KP = IABS(KPVT(K)) 207 IF (KP .EQ. K) GO TO 140 208 T = Z(K) 209 Z(K) = Z(KP) 210 Z(KP) = T 211 140 CONTINUE 212 150 CONTINUE 213 IK = IK + K 214 IF (KS .EQ. 2) IK = IK + (K + 1) 215 K = K + KS 216 GO TO 130 217 160 CONTINUE 218 S = 1.0E0/SCASUM(N,Z,1) 219 CALL CSSCAL(N,S,Z,1) 220C 221 YNORM = 1.0E0 222C 223C SOLVE U*D*V = Y 224C 225 K = N 226 IK = N*(N - 1)/2 227 170 IF (K .EQ. 0) GO TO 230 228 KK = IK + K 229 IKM1 = IK - (K - 1) 230 KS = 1 231 IF (KPVT(K) .LT. 0) KS = 2 232 IF (K .EQ. KS) GO TO 190 233 KP = IABS(KPVT(K)) 234 KPS = K + 1 - KS 235 IF (KP .EQ. KPS) GO TO 180 236 T = Z(KPS) 237 Z(KPS) = Z(KP) 238 Z(KP) = T 239 180 CONTINUE 240 CALL CAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1) 241 IF (KS .EQ. 2) CALL CAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1) 242 190 CONTINUE 243 IF (KS .EQ. 2) GO TO 210 244 IF (CABS1(Z(K)) .LE. CABS1(AP(KK))) GO TO 200 245 S = CABS1(AP(KK))/CABS1(Z(K)) 246 CALL CSSCAL(N,S,Z,1) 247 YNORM = S*YNORM 248 200 CONTINUE 249 IF (CABS1(AP(KK)) .NE. 0.0E0) Z(K) = Z(K)/AP(KK) 250 IF (CABS1(AP(KK)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0) 251 GO TO 220 252 210 CONTINUE 253 KM1K = IK + K - 1 254 KM1KM1 = IKM1 + K - 1 255 AK = AP(KK)/CONJG(AP(KM1K)) 256 AKM1 = AP(KM1KM1)/AP(KM1K) 257 BK = Z(K)/CONJG(AP(KM1K)) 258 BKM1 = Z(K-1)/AP(KM1K) 259 DENOM = AK*AKM1 - 1.0E0 260 Z(K) = (AKM1*BK - BKM1)/DENOM 261 Z(K-1) = (AK*BKM1 - BK)/DENOM 262 220 CONTINUE 263 K = K - KS 264 IK = IK - K 265 IF (KS .EQ. 2) IK = IK - (K + 1) 266 GO TO 170 267 230 CONTINUE 268 S = 1.0E0/SCASUM(N,Z,1) 269 CALL CSSCAL(N,S,Z,1) 270 YNORM = S*YNORM 271C 272C SOLVE CTRANS(U)*Z = V 273C 274 K = 1 275 IK = 0 276 240 IF (K .GT. N) GO TO 270 277 KS = 1 278 IF (KPVT(K) .LT. 0) KS = 2 279 IF (K .EQ. 1) GO TO 260 280 Z(K) = Z(K) + CDOTC(K-1,AP(IK+1),1,Z(1),1) 281 IKP1 = IK + K 282 IF (KS .EQ. 2) 283 1 Z(K+1) = Z(K+1) + CDOTC(K-1,AP(IKP1+1),1,Z(1),1) 284 KP = IABS(KPVT(K)) 285 IF (KP .EQ. K) GO TO 250 286 T = Z(K) 287 Z(K) = Z(KP) 288 Z(KP) = T 289 250 CONTINUE 290 260 CONTINUE 291 IK = IK + K 292 IF (KS .EQ. 2) IK = IK + (K + 1) 293 K = K + KS 294 GO TO 240 295 270 CONTINUE 296C MAKE ZNORM = 1.0 297 S = 1.0E0/SCASUM(N,Z,1) 298 CALL CSSCAL(N,S,Z,1) 299 YNORM = S*YNORM 300C 301 IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM 302 IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0 303 RETURN 304 END 305