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