![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: CHARACTER UPLO 10: INTEGER INFO, N 11: * .. 12: * .. Array Arguments .. 13: INTEGER IPIV( * ) 14: DOUBLE PRECISION AP( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DSPTRF computes the factorization of a real symmetric matrix A stored 21: * in packed format using the Bunch-Kaufman diagonal pivoting method: 22: * 23: * A = U*D*U**T or A = L*D*L**T 24: * 25: * where U (or L) is a product of permutation and unit upper (lower) 26: * triangular matrices, and D is symmetric and block diagonal with 27: * 1-by-1 and 2-by-2 diagonal blocks. 28: * 29: * Arguments 30: * ========= 31: * 32: * UPLO (input) CHARACTER*1 33: * = 'U': Upper triangle of A is stored; 34: * = 'L': Lower triangle of A is stored. 35: * 36: * N (input) INTEGER 37: * The order of the matrix A. N >= 0. 38: * 39: * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 40: * On entry, the upper or lower triangle of the symmetric matrix 41: * A, packed columnwise in a linear array. The j-th column of A 42: * is stored in the array AP as follows: 43: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 44: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 45: * 46: * On exit, the block diagonal matrix D and the multipliers used 47: * to obtain the factor U or L, stored as a packed triangular 48: * matrix overwriting A (see below for further details). 49: * 50: * IPIV (output) INTEGER array, dimension (N) 51: * Details of the interchanges and the block structure of D. 52: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 53: * interchanged and D(k,k) is a 1-by-1 diagonal block. 54: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 55: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 56: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 57: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 58: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 59: * 60: * INFO (output) INTEGER 61: * = 0: successful exit 62: * < 0: if INFO = -i, the i-th argument had an illegal value 63: * > 0: if INFO = i, D(i,i) is exactly zero. The factorization 64: * has been completed, but the block diagonal matrix D is 65: * exactly singular, and division by zero will occur if it 66: * is used to solve a system of equations. 67: * 68: * Further Details 69: * =============== 70: * 71: * 5-96 - Based on modifications by J. Lewis, Boeing Computer Services 72: * Company 73: * 74: * If UPLO = 'U', then A = U*D*U', where 75: * U = P(n)*U(n)* ... *P(k)U(k)* ..., 76: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to 77: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 78: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 79: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 80: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 81: * 82: * ( I v 0 ) k-s 83: * U(k) = ( 0 I 0 ) s 84: * ( 0 0 I ) n-k 85: * k-s s n-k 86: * 87: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 88: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 89: * and A(k,k), and v overwrites A(1:k-2,k-1:k). 90: * 91: * If UPLO = 'L', then A = L*D*L', where 92: * L = P(1)*L(1)* ... *P(k)*L(k)* ..., 93: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 94: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 95: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 96: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 97: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 98: * 99: * ( I 0 0 ) k-1 100: * L(k) = ( 0 I 0 ) s 101: * ( 0 v I ) n-k-s+1 102: * k-1 s n-k-s+1 103: * 104: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 105: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 106: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 107: * 108: * ===================================================================== 109: * 110: * .. Parameters .. 111: DOUBLE PRECISION ZERO, ONE 112: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 113: DOUBLE PRECISION EIGHT, SEVTEN 114: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 115: * .. 116: * .. Local Scalars .. 117: LOGICAL UPPER 118: INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC, 119: $ KSTEP, KX, NPP 120: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1, 121: $ ROWMAX, T, WK, WKM1, WKP1 122: * .. 123: * .. External Functions .. 124: LOGICAL LSAME 125: INTEGER IDAMAX 126: EXTERNAL LSAME, IDAMAX 127: * .. 128: * .. External Subroutines .. 129: EXTERNAL DSCAL, DSPR, DSWAP, XERBLA 130: * .. 131: * .. Intrinsic Functions .. 132: INTRINSIC ABS, MAX, SQRT 133: * .. 134: * .. Executable Statements .. 135: * 136: * Test the input parameters. 137: * 138: INFO = 0 139: UPPER = LSAME( UPLO, 'U' ) 140: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 141: INFO = -1 142: ELSE IF( N.LT.0 ) THEN 143: INFO = -2 144: END IF 145: IF( INFO.NE.0 ) THEN 146: CALL XERBLA( 'DSPTRF', -INFO ) 147: RETURN 148: END IF 149: * 150: * Initialize ALPHA for use in choosing pivot block size. 151: * 152: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 153: * 154: IF( UPPER ) THEN 155: * 156: * Factorize A as U*D*U' using the upper triangle of A 157: * 158: * K is the main loop index, decreasing from N to 1 in steps of 159: * 1 or 2 160: * 161: K = N 162: KC = ( N-1 )*N / 2 + 1 163: 10 CONTINUE 164: KNC = KC 165: * 166: * If K < 1, exit from loop 167: * 168: IF( K.LT.1 ) 169: $ GO TO 110 170: KSTEP = 1 171: * 172: * Determine rows and columns to be interchanged and whether 173: * a 1-by-1 or 2-by-2 pivot block will be used 174: * 175: ABSAKK = ABS( AP( KC+K-1 ) ) 176: * 177: * IMAX is the row-index of the largest off-diagonal element in 178: * column K, and COLMAX is its absolute value 179: * 180: IF( K.GT.1 ) THEN 181: IMAX = IDAMAX( K-1, AP( KC ), 1 ) 182: COLMAX = ABS( AP( KC+IMAX-1 ) ) 183: ELSE 184: COLMAX = ZERO 185: END IF 186: * 187: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 188: * 189: * Column K is zero: set INFO and continue 190: * 191: IF( INFO.EQ.0 ) 192: $ INFO = K 193: KP = K 194: ELSE 195: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 196: * 197: * no interchange, use 1-by-1 pivot block 198: * 199: KP = K 200: ELSE 201: * 202: * JMAX is the column-index of the largest off-diagonal 203: * element in row IMAX, and ROWMAX is its absolute value 204: * 205: ROWMAX = ZERO 206: JMAX = IMAX 207: KX = IMAX*( IMAX+1 ) / 2 + IMAX 208: DO 20 J = IMAX + 1, K 209: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN 210: ROWMAX = ABS( AP( KX ) ) 211: JMAX = J 212: END IF 213: KX = KX + J 214: 20 CONTINUE 215: KPC = ( IMAX-1 )*IMAX / 2 + 1 216: IF( IMAX.GT.1 ) THEN 217: JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 ) 218: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) ) 219: END IF 220: * 221: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 222: * 223: * no interchange, use 1-by-1 pivot block 224: * 225: KP = K 226: ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN 227: * 228: * interchange rows and columns K and IMAX, use 1-by-1 229: * pivot block 230: * 231: KP = IMAX 232: ELSE 233: * 234: * interchange rows and columns K-1 and IMAX, use 2-by-2 235: * pivot block 236: * 237: KP = IMAX 238: KSTEP = 2 239: END IF 240: END IF 241: * 242: KK = K - KSTEP + 1 243: IF( KSTEP.EQ.2 ) 244: $ KNC = KNC - K + 1 245: IF( KP.NE.KK ) THEN 246: * 247: * Interchange rows and columns KK and KP in the leading 248: * submatrix A(1:k,1:k) 249: * 250: CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 ) 251: KX = KPC + KP - 1 252: DO 30 J = KP + 1, KK - 1 253: KX = KX + J - 1 254: T = AP( KNC+J-1 ) 255: AP( KNC+J-1 ) = AP( KX ) 256: AP( KX ) = T 257: 30 CONTINUE 258: T = AP( KNC+KK-1 ) 259: AP( KNC+KK-1 ) = AP( KPC+KP-1 ) 260: AP( KPC+KP-1 ) = T 261: IF( KSTEP.EQ.2 ) THEN 262: T = AP( KC+K-2 ) 263: AP( KC+K-2 ) = AP( KC+KP-1 ) 264: AP( KC+KP-1 ) = T 265: END IF 266: END IF 267: * 268: * Update the leading submatrix 269: * 270: IF( KSTEP.EQ.1 ) THEN 271: * 272: * 1-by-1 pivot block D(k): column k now holds 273: * 274: * W(k) = U(k)*D(k) 275: * 276: * where U(k) is the k-th column of U 277: * 278: * Perform a rank-1 update of A(1:k-1,1:k-1) as 279: * 280: * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' 281: * 282: R1 = ONE / AP( KC+K-1 ) 283: CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP ) 284: * 285: * Store U(k) in column k 286: * 287: CALL DSCAL( K-1, R1, AP( KC ), 1 ) 288: ELSE 289: * 290: * 2-by-2 pivot block D(k): columns k and k-1 now hold 291: * 292: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 293: * 294: * where U(k) and U(k-1) are the k-th and (k-1)-th columns 295: * of U 296: * 297: * Perform a rank-2 update of A(1:k-2,1:k-2) as 298: * 299: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )' 300: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' 301: * 302: IF( K.GT.2 ) THEN 303: * 304: D12 = AP( K-1+( K-1 )*K / 2 ) 305: D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 306: D11 = AP( K+( K-1 )*K / 2 ) / D12 307: T = ONE / ( D11*D22-ONE ) 308: D12 = T / D12 309: * 310: DO 50 J = K - 2, 1, -1 311: WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- 312: $ AP( J+( K-1 )*K / 2 ) ) 313: WK = D12*( D22*AP( J+( K-1 )*K / 2 )- 314: $ AP( J+( K-2 )*( K-1 ) / 2 ) ) 315: DO 40 I = J, 1, -1 316: AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - 317: $ AP( I+( K-1 )*K / 2 )*WK - 318: $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1 319: 40 CONTINUE 320: AP( J+( K-1 )*K / 2 ) = WK 321: AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1 322: 50 CONTINUE 323: * 324: END IF 325: * 326: END IF 327: END IF 328: * 329: * Store details of the interchanges in IPIV 330: * 331: IF( KSTEP.EQ.1 ) THEN 332: IPIV( K ) = KP 333: ELSE 334: IPIV( K ) = -KP 335: IPIV( K-1 ) = -KP 336: END IF 337: * 338: * Decrease K and return to the start of the main loop 339: * 340: K = K - KSTEP 341: KC = KNC - K 342: GO TO 10 343: * 344: ELSE 345: * 346: * Factorize A as L*D*L' using the lower triangle of A 347: * 348: * K is the main loop index, increasing from 1 to N in steps of 349: * 1 or 2 350: * 351: K = 1 352: KC = 1 353: NPP = N*( N+1 ) / 2 354: 60 CONTINUE 355: KNC = KC 356: * 357: * If K > N, exit from loop 358: * 359: IF( K.GT.N ) 360: $ GO TO 110 361: KSTEP = 1 362: * 363: * Determine rows and columns to be interchanged and whether 364: * a 1-by-1 or 2-by-2 pivot block will be used 365: * 366: ABSAKK = ABS( AP( KC ) ) 367: * 368: * IMAX is the row-index of the largest off-diagonal element in 369: * column K, and COLMAX is its absolute value 370: * 371: IF( K.LT.N ) THEN 372: IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 ) 373: COLMAX = ABS( AP( KC+IMAX-K ) ) 374: ELSE 375: COLMAX = ZERO 376: END IF 377: * 378: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 379: * 380: * Column K is zero: set INFO and continue 381: * 382: IF( INFO.EQ.0 ) 383: $ INFO = K 384: KP = K 385: ELSE 386: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 387: * 388: * no interchange, use 1-by-1 pivot block 389: * 390: KP = K 391: ELSE 392: * 393: * JMAX is the column-index of the largest off-diagonal 394: * element in row IMAX, and ROWMAX is its absolute value 395: * 396: ROWMAX = ZERO 397: KX = KC + IMAX - K 398: DO 70 J = K, IMAX - 1 399: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN 400: ROWMAX = ABS( AP( KX ) ) 401: JMAX = J 402: END IF 403: KX = KX + N - J 404: 70 CONTINUE 405: KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1 406: IF( IMAX.LT.N ) THEN 407: JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 ) 408: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) ) 409: END IF 410: * 411: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 412: * 413: * no interchange, use 1-by-1 pivot block 414: * 415: KP = K 416: ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN 417: * 418: * interchange rows and columns K and IMAX, use 1-by-1 419: * pivot block 420: * 421: KP = IMAX 422: ELSE 423: * 424: * interchange rows and columns K+1 and IMAX, use 2-by-2 425: * pivot block 426: * 427: KP = IMAX 428: KSTEP = 2 429: END IF 430: END IF 431: * 432: KK = K + KSTEP - 1 433: IF( KSTEP.EQ.2 ) 434: $ KNC = KNC + N - K + 1 435: IF( KP.NE.KK ) THEN 436: * 437: * Interchange rows and columns KK and KP in the trailing 438: * submatrix A(k:n,k:n) 439: * 440: IF( KP.LT.N ) 441: $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), 442: $ 1 ) 443: KX = KNC + KP - KK 444: DO 80 J = KK + 1, KP - 1 445: KX = KX + N - J + 1 446: T = AP( KNC+J-KK ) 447: AP( KNC+J-KK ) = AP( KX ) 448: AP( KX ) = T 449: 80 CONTINUE 450: T = AP( KNC ) 451: AP( KNC ) = AP( KPC ) 452: AP( KPC ) = T 453: IF( KSTEP.EQ.2 ) THEN 454: T = AP( KC+1 ) 455: AP( KC+1 ) = AP( KC+KP-K ) 456: AP( KC+KP-K ) = T 457: END IF 458: END IF 459: * 460: * Update the trailing submatrix 461: * 462: IF( KSTEP.EQ.1 ) THEN 463: * 464: * 1-by-1 pivot block D(k): column k now holds 465: * 466: * W(k) = L(k)*D(k) 467: * 468: * where L(k) is the k-th column of L 469: * 470: IF( K.LT.N ) THEN 471: * 472: * Perform a rank-1 update of A(k+1:n,k+1:n) as 473: * 474: * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' 475: * 476: R1 = ONE / AP( KC ) 477: CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1, 478: $ AP( KC+N-K+1 ) ) 479: * 480: * Store L(k) in column K 481: * 482: CALL DSCAL( N-K, R1, AP( KC+1 ), 1 ) 483: END IF 484: ELSE 485: * 486: * 2-by-2 pivot block D(k): columns K and K+1 now hold 487: * 488: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 489: * 490: * where L(k) and L(k+1) are the k-th and (k+1)-th columns 491: * of L 492: * 493: IF( K.LT.N-1 ) THEN 494: * 495: * Perform a rank-2 update of A(k+2:n,k+2:n) as 496: * 497: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )' 498: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )' 499: * 500: D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) 501: D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 502: D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 503: T = ONE / ( D11*D22-ONE ) 504: D21 = T / D21 505: * 506: DO 100 J = K + 2, N 507: WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- 508: $ AP( J+K*( 2*N-K-1 ) / 2 ) ) 509: WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- 510: $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) 511: * 512: DO 90 I = J, N 513: AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* 514: $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / 515: $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1 516: 90 CONTINUE 517: * 518: AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK 519: AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1 520: * 521: 100 CONTINUE 522: END IF 523: END IF 524: END IF 525: * 526: * Store details of the interchanges in IPIV 527: * 528: IF( KSTEP.EQ.1 ) THEN 529: IPIV( K ) = KP 530: ELSE 531: IPIV( K ) = -KP 532: IPIV( K+1 ) = -KP 533: END IF 534: * 535: * Increase K and return to the start of the main loop 536: * 537: K = K + KSTEP 538: KC = KNC + N - K + 2 539: GO TO 60 540: * 541: END IF 542: * 543: 110 CONTINUE 544: RETURN 545: * 546: * End of DSPTRF 547: * 548: END