1*DECK DSPFA 2 SUBROUTINE DSPFA (AP, N, KPVT, INFO) 3C***BEGIN PROLOGUE DSPFA 4C***PURPOSE Factor a real symmetric matrix stored in packed form by 5C elimination with symmetric pivoting. 6C***LIBRARY SLATEC (LINPACK) 7C***CATEGORY D2B1A 8C***TYPE DOUBLE PRECISION (SSPFA-S, DSPFA-D, CHPFA-C, CSPFA-C) 9C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED, 10C SYMMETRIC 11C***AUTHOR Bunch, J., (UCSD) 12C***DESCRIPTION 13C 14C DSPFA factors a double precision symmetric matrix stored in 15C packed form by elimination with symmetric pivoting. 16C 17C To solve A*X = B , follow DSPFA by DSPSL. 18C To compute INVERSE(A)*C , follow DSPFA by DSPSL. 19C To compute DETERMINANT(A) , follow DSPFA by DSPDI. 20C To compute INERTIA(A) , follow DSPFA by DSPDI. 21C To compute INVERSE(A) , follow DSPFA by DSPDI. 22C 23C On Entry 24C 25C AP DOUBLE PRECISION (N*(N+1)/2) 26C the packed form of a symmetric matrix A . The 27C columns of the upper triangle are stored sequentially 28C in a one-dimensional array of length N*(N+1)/2 . 29C See comments below for details. 30C 31C N INTEGER 32C the order of the matrix A . 33C 34C Output 35C 36C AP a block diagonal matrix and the multipliers which 37C were used to obtain it stored in packed form. 38C The factorization can be written A = U*D*TRANS(U) 39C where U is a product of permutation and unit 40C upper triangular matrices, TRANS(U) is the 41C transpose of U , and D is block diagonal 42C with 1 by 1 and 2 by 2 blocks. 43C 44C KPVT INTEGER(N) 45C an integer vector of pivot indices. 46C 47C INFO INTEGER 48C = 0 normal value. 49C = K if the K-th pivot block is singular. This is 50C not an error condition for this subroutine, 51C but it does indicate that DSPSL or DSPDI may 52C divide by zero if called. 53C 54C Packed Storage 55C 56C The following program segment will pack the upper 57C triangle of a symmetric matrix. 58C 59C K = 0 60C DO 20 J = 1, N 61C DO 10 I = 1, J 62C K = K + 1 63C AP(K) = A(I,J) 64C 10 CONTINUE 65C 20 CONTINUE 66C 67C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. 68C Stewart, LINPACK Users' Guide, SIAM, 1979. 69C***ROUTINES CALLED DAXPY, DSWAP, IDAMAX 70C***REVISION HISTORY (YYMMDD) 71C 780814 DATE WRITTEN 72C 890531 Changed all specific intrinsics to generic. (WRB) 73C 890831 Modified array declarations. (WRB) 74C 891107 Modified routine equivalence list. (WRB) 75C 891107 REVISION DATE from Version 3.2 76C 891214 Prologue converted to Version 4.0 format. (BAB) 77C 900326 Removed duplicate information from DESCRIPTION section. 78C (WRB) 79C 920501 Reformatted the REFERENCES section. (WRB) 80C***END PROLOGUE DSPFA 81 INTEGER N,KPVT(*),INFO 82 DOUBLE PRECISION AP(*) 83C 84 DOUBLE PRECISION AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T 85 DOUBLE PRECISION ABSAKK,ALPHA,COLMAX,ROWMAX 86 INTEGER IDAMAX,IJ,IK,IKM1,IM,IMAX,IMAXP1,IMIM,IMJ,IMK 87 INTEGER J,JJ,JK,JKM1,JMAX,JMIM,K,KK,KM1,KM1K,KM1KM1,KM2,KSTEP 88 LOGICAL SWAP 89C***FIRST EXECUTABLE STATEMENT DSPFA 90C 91C INITIALIZE 92C 93C ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE. 94C 95 ALPHA = (1.0D0 + SQRT(17.0D0))/8.0D0 96C 97 INFO = 0 98C 99C MAIN LOOP ON K, WHICH GOES FROM N TO 1. 100C 101 K = N 102 IK = (N*(N - 1))/2 103 10 CONTINUE 104C 105C LEAVE THE LOOP IF K=0 OR K=1. 106C 107 IF (K .EQ. 0) GO TO 200 108 IF (K .GT. 1) GO TO 20 109 KPVT(1) = 1 110 IF (AP(1) .EQ. 0.0D0) INFO = 1 111 GO TO 200 112 20 CONTINUE 113C 114C THIS SECTION OF CODE DETERMINES THE KIND OF 115C ELIMINATION TO BE PERFORMED. WHEN IT IS COMPLETED, 116C KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND 117C SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS 118C REQUIRED. 119C 120 KM1 = K - 1 121 KK = IK + K 122 ABSAKK = ABS(AP(KK)) 123C 124C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN 125C COLUMN K. 126C 127 IMAX = IDAMAX(K-1,AP(IK+1),1) 128 IMK = IK + IMAX 129 COLMAX = ABS(AP(IMK)) 130 IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30 131 KSTEP = 1 132 SWAP = .FALSE. 133 GO TO 90 134 30 CONTINUE 135C 136C DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN 137C ROW IMAX. 138C 139 ROWMAX = 0.0D0 140 IMAXP1 = IMAX + 1 141 IM = IMAX*(IMAX - 1)/2 142 IMJ = IM + 2*IMAX 143 DO 40 J = IMAXP1, K 144 ROWMAX = MAX(ROWMAX,ABS(AP(IMJ))) 145 IMJ = IMJ + J 146 40 CONTINUE 147 IF (IMAX .EQ. 1) GO TO 50 148 JMAX = IDAMAX(IMAX-1,AP(IM+1),1) 149 JMIM = JMAX + IM 150 ROWMAX = MAX(ROWMAX,ABS(AP(JMIM))) 151 50 CONTINUE 152 IMIM = IMAX + IM 153 IF (ABS(AP(IMIM)) .LT. ALPHA*ROWMAX) GO TO 60 154 KSTEP = 1 155 SWAP = .TRUE. 156 GO TO 80 157 60 CONTINUE 158 IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70 159 KSTEP = 1 160 SWAP = .FALSE. 161 GO TO 80 162 70 CONTINUE 163 KSTEP = 2 164 SWAP = IMAX .NE. KM1 165 80 CONTINUE 166 90 CONTINUE 167 IF (MAX(ABSAKK,COLMAX) .NE. 0.0D0) GO TO 100 168C 169C COLUMN K IS ZERO. SET INFO AND ITERATE THE LOOP. 170C 171 KPVT(K) = K 172 INFO = K 173 GO TO 190 174 100 CONTINUE 175 IF (KSTEP .EQ. 2) GO TO 140 176C 177C 1 X 1 PIVOT BLOCK. 178C 179 IF (.NOT.SWAP) GO TO 120 180C 181C PERFORM AN INTERCHANGE. 182C 183 CALL DSWAP(IMAX,AP(IM+1),1,AP(IK+1),1) 184 IMJ = IK + IMAX 185 DO 110 JJ = IMAX, K 186 J = K + IMAX - JJ 187 JK = IK + J 188 T = AP(JK) 189 AP(JK) = AP(IMJ) 190 AP(IMJ) = T 191 IMJ = IMJ - (J - 1) 192 110 CONTINUE 193 120 CONTINUE 194C 195C PERFORM THE ELIMINATION. 196C 197 IJ = IK - (K - 1) 198 DO 130 JJ = 1, KM1 199 J = K - JJ 200 JK = IK + J 201 MULK = -AP(JK)/AP(KK) 202 T = MULK 203 CALL DAXPY(J,T,AP(IK+1),1,AP(IJ+1),1) 204 AP(JK) = MULK 205 IJ = IJ - (J - 1) 206 130 CONTINUE 207C 208C SET THE PIVOT ARRAY. 209C 210 KPVT(K) = K 211 IF (SWAP) KPVT(K) = IMAX 212 GO TO 190 213 140 CONTINUE 214C 215C 2 X 2 PIVOT BLOCK. 216C 217 KM1K = IK + K - 1 218 IKM1 = IK - (K - 1) 219 IF (.NOT.SWAP) GO TO 160 220C 221C PERFORM AN INTERCHANGE. 222C 223 CALL DSWAP(IMAX,AP(IM+1),1,AP(IKM1+1),1) 224 IMJ = IKM1 + IMAX 225 DO 150 JJ = IMAX, KM1 226 J = KM1 + IMAX - JJ 227 JKM1 = IKM1 + J 228 T = AP(JKM1) 229 AP(JKM1) = AP(IMJ) 230 AP(IMJ) = T 231 IMJ = IMJ - (J - 1) 232 150 CONTINUE 233 T = AP(KM1K) 234 AP(KM1K) = AP(IMK) 235 AP(IMK) = T 236 160 CONTINUE 237C 238C PERFORM THE ELIMINATION. 239C 240 KM2 = K - 2 241 IF (KM2 .EQ. 0) GO TO 180 242 AK = AP(KK)/AP(KM1K) 243 KM1KM1 = IKM1 + K - 1 244 AKM1 = AP(KM1KM1)/AP(KM1K) 245 DENOM = 1.0D0 - AK*AKM1 246 IJ = IK - (K - 1) - (K - 2) 247 DO 170 JJ = 1, KM2 248 J = KM1 - JJ 249 JK = IK + J 250 BK = AP(JK)/AP(KM1K) 251 JKM1 = IKM1 + J 252 BKM1 = AP(JKM1)/AP(KM1K) 253 MULK = (AKM1*BK - BKM1)/DENOM 254 MULKM1 = (AK*BKM1 - BK)/DENOM 255 T = MULK 256 CALL DAXPY(J,T,AP(IK+1),1,AP(IJ+1),1) 257 T = MULKM1 258 CALL DAXPY(J,T,AP(IKM1+1),1,AP(IJ+1),1) 259 AP(JK) = MULK 260 AP(JKM1) = MULKM1 261 IJ = IJ - (J - 1) 262 170 CONTINUE 263 180 CONTINUE 264C 265C SET THE PIVOT ARRAY. 266C 267 KPVT(K) = 1 - K 268 IF (SWAP) KPVT(K) = -IMAX 269 KPVT(K-1) = KPVT(K) 270 190 CONTINUE 271 IK = IK - (K - 1) 272 IF (KSTEP .EQ. 2) IK = IK - (K - 2) 273 K = K - KSTEP 274 GO TO 10 275 200 CONTINUE 276 RETURN 277 END 278