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