1 SUBROUTINE CSICO(A,LDA,N,KPVT,RCOND,Z) 2C***BEGIN PROLOGUE CSICO 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 SYMMETRIC 10C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) 11C***PURPOSE Factors a COMPLEX SYMMETRIC matrix by elimination with 12C symmetric pivoting and estimates the condition of the 13C matrix. 14C***DESCRIPTION 15C 16C CSICO factors a complex symmetric matrix by elimination with 17C symmetric pivoting and estimates the condition of the matrix. 18C 19C If RCOND is not needed, CSIFA is slightly faster. 20C To solve A*X = B , follow CSICO by CSISL. 21C To compute INVERSE(A)*C , follow CSICO by CSISL. 22C To compute INVERSE(A) , follow CSICO by CSIDI. 23C To compute DETERMINANT(A) , follow CSICO by CSIDI. 24C 25C On Entry 26C 27C A COMPLEX(LDA, N) 28C the symmetric matrix to be factored. 29C Only the diagonal and upper triangle are used. 30C 31C LDA INTEGER 32C the leading dimension of the array A . 33C 34C N INTEGER 35C the order of the matrix A . 36C 37C On Return 38C 39C A a block diagonal matrix and the multipliers which 40C were used to obtain it. 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 LINPACK. This version dated 08/14/78 . 69C Cleve Moler, University of New Mexico, Argonne National Lab. 70C 71C Subroutines and Functions 72C 73C LINPACK CSIFA 74C BLAS CAXPY,CDOTU,CSSCAL,SCASUM 75C Fortran ABS,AIMAG,AMAX1,CMPLX,IABS,REAL 76C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., 77C *LINPACK USERS GUIDE*, SIAM, 1979. 78C***ROUTINES CALLED CAXPY,CDOTU,CSIFA,CSSCAL,SCASUM 79C***END PROLOGUE CSICO 80 INTEGER LDA,N,KPVT(*) 81 COMPLEX A(LDA,*),Z(*) 82 REAL RCOND 83C 84 COMPLEX AK,AKM1,BK,BKM1,CDOTU,DENOM,EK,T 85 REAL ANORM,S,SCASUM,YNORM 86 INTEGER I,INFO,J,JM1,K,KP,KPS,KS 87 COMPLEX ZDUM,ZDUM2,CSIGN1 88 REAL CABS1 89 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) 90 CSIGN1(ZDUM,ZDUM2) = CABS1(ZDUM)*(ZDUM2/CABS1(ZDUM2)) 91C 92C FIND NORM OF A USING ONLY UPPER HALF 93C 94C***FIRST EXECUTABLE STATEMENT CSICO 95 DO 30 J = 1, N 96 Z(J) = CMPLX(SCASUM(J,A(1,J),1),0.0E0) 97 JM1 = J - 1 98 IF (JM1 .LT. 1) GO TO 20 99 DO 10 I = 1, JM1 100 Z(I) = CMPLX(REAL(Z(I))+CABS1(A(I,J)),0.0E0) 101 10 CONTINUE 102 20 CONTINUE 103 30 CONTINUE 104 ANORM = 0.0E0 105 DO 40 J = 1, N 106 ANORM = AMAX1(ANORM,REAL(Z(J))) 107 40 CONTINUE 108C 109C FACTOR 110C 111 CALL CSIFA(A,LDA,N,KPVT,INFO) 112C 113C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . 114C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E . 115C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL 116C GROWTH IN THE ELEMENTS OF W WHERE U*D*W = E . 117C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 118C 119C SOLVE U*D*W = E 120C 121 EK = (1.0E0,0.0E0) 122 DO 50 J = 1, N 123 Z(J) = (0.0E0,0.0E0) 124 50 CONTINUE 125 K = N 126 60 IF (K .EQ. 0) GO TO 120 127 KS = 1 128 IF (KPVT(K) .LT. 0) KS = 2 129 KP = IABS(KPVT(K)) 130 KPS = K + 1 - KS 131 IF (KP .EQ. KPS) GO TO 70 132 T = Z(KPS) 133 Z(KPS) = Z(KP) 134 Z(KP) = T 135 70 CONTINUE 136 IF (CABS1(Z(K)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K)) 137 Z(K) = Z(K) + EK 138 CALL CAXPY(K-KS,Z(K),A(1,K),1,Z(1),1) 139 IF (KS .EQ. 1) GO TO 80 140 IF (CABS1(Z(K-1)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K-1)) 141 Z(K-1) = Z(K-1) + EK 142 CALL CAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1) 143 80 CONTINUE 144 IF (KS .EQ. 2) GO TO 100 145 IF (CABS1(Z(K)) .LE. CABS1(A(K,K))) GO TO 90 146 S = CABS1(A(K,K))/CABS1(Z(K)) 147 CALL CSSCAL(N,S,Z,1) 148 EK = CMPLX(S,0.0E0)*EK 149 90 CONTINUE 150 IF (CABS1(A(K,K)) .NE. 0.0E0) Z(K) = Z(K)/A(K,K) 151 IF (CABS1(A(K,K)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0) 152 GO TO 110 153 100 CONTINUE 154 AK = A(K,K)/A(K-1,K) 155 AKM1 = A(K-1,K-1)/A(K-1,K) 156 BK = Z(K)/A(K-1,K) 157 BKM1 = Z(K-1)/A(K-1,K) 158 DENOM = AK*AKM1 - 1.0E0 159 Z(K) = (AKM1*BK - BKM1)/DENOM 160 Z(K-1) = (AK*BKM1 - BK)/DENOM 161 110 CONTINUE 162 K = K - KS 163 GO TO 60 164 120 CONTINUE 165 S = 1.0E0/SCASUM(N,Z,1) 166 CALL CSSCAL(N,S,Z,1) 167C 168C SOLVE TRANS(U)*Y = W 169C 170 K = 1 171 130 IF (K .GT. N) GO TO 160 172 KS = 1 173 IF (KPVT(K) .LT. 0) KS = 2 174 IF (K .EQ. 1) GO TO 150 175 Z(K) = Z(K) + CDOTU(K-1,A(1,K),1,Z(1),1) 176 IF (KS .EQ. 2) 177 1 Z(K+1) = Z(K+1) + CDOTU(K-1,A(1,K+1),1,Z(1),1) 178 KP = IABS(KPVT(K)) 179 IF (KP .EQ. K) GO TO 140 180 T = Z(K) 181 Z(K) = Z(KP) 182 Z(KP) = T 183 140 CONTINUE 184 150 CONTINUE 185 K = K + KS 186 GO TO 130 187 160 CONTINUE 188 S = 1.0E0/SCASUM(N,Z,1) 189 CALL CSSCAL(N,S,Z,1) 190C 191 YNORM = 1.0E0 192C 193C SOLVE U*D*V = Y 194C 195 K = N 196 170 IF (K .EQ. 0) GO TO 230 197 KS = 1 198 IF (KPVT(K) .LT. 0) KS = 2 199 IF (K .EQ. KS) GO TO 190 200 KP = IABS(KPVT(K)) 201 KPS = K + 1 - KS 202 IF (KP .EQ. KPS) GO TO 180 203 T = Z(KPS) 204 Z(KPS) = Z(KP) 205 Z(KP) = T 206 180 CONTINUE 207 CALL CAXPY(K-KS,Z(K),A(1,K),1,Z(1),1) 208 IF (KS .EQ. 2) CALL CAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1) 209 190 CONTINUE 210 IF (KS .EQ. 2) GO TO 210 211 IF (CABS1(Z(K)) .LE. CABS1(A(K,K))) GO TO 200 212 S = CABS1(A(K,K))/CABS1(Z(K)) 213 CALL CSSCAL(N,S,Z,1) 214 YNORM = S*YNORM 215 200 CONTINUE 216 IF (CABS1(A(K,K)) .NE. 0.0E0) Z(K) = Z(K)/A(K,K) 217 IF (CABS1(A(K,K)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0) 218 GO TO 220 219 210 CONTINUE 220 AK = A(K,K)/A(K-1,K) 221 AKM1 = A(K-1,K-1)/A(K-1,K) 222 BK = Z(K)/A(K-1,K) 223 BKM1 = Z(K-1)/A(K-1,K) 224 DENOM = AK*AKM1 - 1.0E0 225 Z(K) = (AKM1*BK - BKM1)/DENOM 226 Z(K-1) = (AK*BKM1 - BK)/DENOM 227 220 CONTINUE 228 K = K - KS 229 GO TO 170 230 230 CONTINUE 231 S = 1.0E0/SCASUM(N,Z,1) 232 CALL CSSCAL(N,S,Z,1) 233 YNORM = S*YNORM 234C 235C SOLVE TRANS(U)*Z = V 236C 237 K = 1 238 240 IF (K .GT. N) GO TO 270 239 KS = 1 240 IF (KPVT(K) .LT. 0) KS = 2 241 IF (K .EQ. 1) GO TO 260 242 Z(K) = Z(K) + CDOTU(K-1,A(1,K),1,Z(1),1) 243 IF (KS .EQ. 2) 244 1 Z(K+1) = Z(K+1) + CDOTU(K-1,A(1,K+1),1,Z(1),1) 245 KP = IABS(KPVT(K)) 246 IF (KP .EQ. K) GO TO 250 247 T = Z(K) 248 Z(K) = Z(KP) 249 Z(KP) = T 250 250 CONTINUE 251 260 CONTINUE 252 K = K + KS 253 GO TO 240 254 270 CONTINUE 255C MAKE ZNORM = 1.0 256 S = 1.0E0/SCASUM(N,Z,1) 257 CALL CSSCAL(N,S,Z,1) 258 YNORM = S*YNORM 259C 260 IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM 261 IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0 262 RETURN 263 END 264