![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, 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, KB, LDA, LDW, N, NB 11: * .. 12: * .. Array Arguments .. 13: INTEGER IPIV( * ) 14: DOUBLE PRECISION A( LDA, * ), W( LDW, * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLASYF computes a partial factorization of a real symmetric matrix A 21: * using the Bunch-Kaufman diagonal pivoting method. The partial 22: * factorization has the form: 23: * 24: * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or: 25: * ( 0 U22 ) ( 0 D ) ( U12' U22' ) 26: * 27: * A = ( L11 0 ) ( D 0 ) ( L11' L21' ) if UPLO = 'L' 28: * ( L21 I ) ( 0 A22 ) ( 0 I ) 29: * 30: * where the order of D is at most NB. The actual order is returned in 31: * the argument KB, and is either NB or NB-1, or N if N <= NB. 32: * 33: * DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code 34: * (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or 35: * A22 (if UPLO = 'L'). 36: * 37: * Arguments 38: * ========= 39: * 40: * UPLO (input) CHARACTER*1 41: * Specifies whether the upper or lower triangular part of the 42: * symmetric matrix A is stored: 43: * = 'U': Upper triangular 44: * = 'L': Lower triangular 45: * 46: * N (input) INTEGER 47: * The order of the matrix A. N >= 0. 48: * 49: * NB (input) INTEGER 50: * The maximum number of columns of the matrix A that should be 51: * factored. NB should be at least 2 to allow for 2-by-2 pivot 52: * blocks. 53: * 54: * KB (output) INTEGER 55: * The number of columns of A that were actually factored. 56: * KB is either NB-1 or NB, or N if N <= NB. 57: * 58: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 59: * On entry, the symmetric matrix A. If UPLO = 'U', the leading 60: * n-by-n upper triangular part of A contains the upper 61: * triangular part of the matrix A, and the strictly lower 62: * triangular part of A is not referenced. If UPLO = 'L', the 63: * leading n-by-n lower triangular part of A contains the lower 64: * triangular part of the matrix A, and the strictly upper 65: * triangular part of A is not referenced. 66: * On exit, A contains details of the partial factorization. 67: * 68: * LDA (input) INTEGER 69: * The leading dimension of the array A. LDA >= max(1,N). 70: * 71: * IPIV (output) INTEGER array, dimension (N) 72: * Details of the interchanges and the block structure of D. 73: * If UPLO = 'U', only the last KB elements of IPIV are set; 74: * if UPLO = 'L', only the first KB elements are set. 75: * 76: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 77: * interchanged and D(k,k) is a 1-by-1 diagonal block. 78: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 79: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 80: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 81: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 82: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 83: * 84: * W (workspace) DOUBLE PRECISION array, dimension (LDW,NB) 85: * 86: * LDW (input) INTEGER 87: * The leading dimension of the array W. LDW >= max(1,N). 88: * 89: * INFO (output) INTEGER 90: * = 0: successful exit 91: * > 0: if INFO = k, D(k,k) is exactly zero. The factorization 92: * has been completed, but the block diagonal matrix D is 93: * exactly singular. 94: * 95: * ===================================================================== 96: * 97: * .. Parameters .. 98: DOUBLE PRECISION ZERO, ONE 99: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 100: DOUBLE PRECISION EIGHT, SEVTEN 101: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 102: * .. 103: * .. Local Scalars .. 104: INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP, 105: $ KSTEP, KW 106: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1, 107: $ ROWMAX, T 108: * .. 109: * .. External Functions .. 110: LOGICAL LSAME 111: INTEGER IDAMAX 112: EXTERNAL LSAME, IDAMAX 113: * .. 114: * .. External Subroutines .. 115: EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP 116: * .. 117: * .. Intrinsic Functions .. 118: INTRINSIC ABS, MAX, MIN, SQRT 119: * .. 120: * .. Executable Statements .. 121: * 122: INFO = 0 123: * 124: * Initialize ALPHA for use in choosing pivot block size. 125: * 126: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 127: * 128: IF( LSAME( UPLO, 'U' ) ) THEN 129: * 130: * Factorize the trailing columns of A using the upper triangle 131: * of A and working backwards, and compute the matrix W = U12*D 132: * for use in updating A11 133: * 134: * K is the main loop index, decreasing from N in steps of 1 or 2 135: * 136: * KW is the column of W which corresponds to column K of A 137: * 138: K = N 139: 10 CONTINUE 140: KW = NB + K - N 141: * 142: * Exit from loop 143: * 144: IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 ) 145: $ GO TO 30 146: * 147: * Copy column K of A to column KW of W and update it 148: * 149: CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 ) 150: IF( K.LT.N ) 151: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA, 152: $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 ) 153: * 154: KSTEP = 1 155: * 156: * Determine rows and columns to be interchanged and whether 157: * a 1-by-1 or 2-by-2 pivot block will be used 158: * 159: ABSAKK = ABS( W( K, KW ) ) 160: * 161: * IMAX is the row-index of the largest off-diagonal element in 162: * column K, and COLMAX is its absolute value 163: * 164: IF( K.GT.1 ) THEN 165: IMAX = IDAMAX( K-1, W( 1, KW ), 1 ) 166: COLMAX = ABS( W( IMAX, KW ) ) 167: ELSE 168: COLMAX = ZERO 169: END IF 170: * 171: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 172: * 173: * Column K is zero: set INFO and continue 174: * 175: IF( INFO.EQ.0 ) 176: $ INFO = K 177: KP = K 178: ELSE 179: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 180: * 181: * no interchange, use 1-by-1 pivot block 182: * 183: KP = K 184: ELSE 185: * 186: * Copy column IMAX to column KW-1 of W and update it 187: * 188: CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) 189: CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, 190: $ W( IMAX+1, KW-1 ), 1 ) 191: IF( K.LT.N ) 192: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), 193: $ LDA, W( IMAX, KW+1 ), LDW, ONE, 194: $ W( 1, KW-1 ), 1 ) 195: * 196: * JMAX is the column-index of the largest off-diagonal 197: * element in row IMAX, and ROWMAX is its absolute value 198: * 199: JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 ) 200: ROWMAX = ABS( W( JMAX, KW-1 ) ) 201: IF( IMAX.GT.1 ) THEN 202: JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 ) 203: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) ) 204: END IF 205: * 206: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 207: * 208: * no interchange, use 1-by-1 pivot block 209: * 210: KP = K 211: ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN 212: * 213: * interchange rows and columns K and IMAX, use 1-by-1 214: * pivot block 215: * 216: KP = IMAX 217: * 218: * copy column KW-1 of W to column KW 219: * 220: CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) 221: ELSE 222: * 223: * interchange rows and columns K-1 and IMAX, use 2-by-2 224: * pivot block 225: * 226: KP = IMAX 227: KSTEP = 2 228: END IF 229: END IF 230: * 231: KK = K - KSTEP + 1 232: KKW = NB + KK - N 233: * 234: * Updated column KP is already stored in column KKW of W 235: * 236: IF( KP.NE.KK ) THEN 237: * 238: * Copy non-updated column KK to column KP 239: * 240: A( KP, K ) = A( KK, K ) 241: CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), 242: $ LDA ) 243: CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 ) 244: * 245: * Interchange rows KK and KP in last KK columns of A and W 246: * 247: CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA ) 248: CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), 249: $ LDW ) 250: END IF 251: * 252: IF( KSTEP.EQ.1 ) THEN 253: * 254: * 1-by-1 pivot block D(k): column KW of W now holds 255: * 256: * W(k) = U(k)*D(k) 257: * 258: * where U(k) is the k-th column of U 259: * 260: * Store U(k) in column k of A 261: * 262: CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) 263: R1 = ONE / A( K, K ) 264: CALL DSCAL( K-1, R1, A( 1, K ), 1 ) 265: ELSE 266: * 267: * 2-by-2 pivot block D(k): columns KW and KW-1 of W now 268: * hold 269: * 270: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 271: * 272: * where U(k) and U(k-1) are the k-th and (k-1)-th columns 273: * of U 274: * 275: IF( K.GT.2 ) THEN 276: * 277: * Store U(k) and U(k-1) in columns k and k-1 of A 278: * 279: D21 = W( K-1, KW ) 280: D11 = W( K, KW ) / D21 281: D22 = W( K-1, KW-1 ) / D21 282: T = ONE / ( D11*D22-ONE ) 283: D21 = T / D21 284: DO 20 J = 1, K - 2 285: A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) 286: A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) 287: 20 CONTINUE 288: END IF 289: * 290: * Copy D(k) to A 291: * 292: A( K-1, K-1 ) = W( K-1, KW-1 ) 293: A( K-1, K ) = W( K-1, KW ) 294: A( K, K ) = W( K, KW ) 295: END IF 296: END IF 297: * 298: * Store details of the interchanges in IPIV 299: * 300: IF( KSTEP.EQ.1 ) THEN 301: IPIV( K ) = KP 302: ELSE 303: IPIV( K ) = -KP 304: IPIV( K-1 ) = -KP 305: END IF 306: * 307: * Decrease K and return to the start of the main loop 308: * 309: K = K - KSTEP 310: GO TO 10 311: * 312: 30 CONTINUE 313: * 314: * Update the upper triangle of A11 (= A(1:k,1:k)) as 315: * 316: * A11 := A11 - U12*D*U12' = A11 - U12*W' 317: * 318: * computing blocks of NB columns at a time 319: * 320: DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB 321: JB = MIN( NB, K-J+1 ) 322: * 323: * Update the upper triangle of the diagonal block 324: * 325: DO 40 JJ = J, J + JB - 1 326: CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE, 327: $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE, 328: $ A( J, JJ ), 1 ) 329: 40 CONTINUE 330: * 331: * Update the rectangular superdiagonal block 332: * 333: CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE, 334: $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE, 335: $ A( 1, J ), LDA ) 336: 50 CONTINUE 337: * 338: * Put U12 in standard form by partially undoing the interchanges 339: * in columns k+1:n 340: * 341: J = K + 1 342: 60 CONTINUE 343: JJ = J 344: JP = IPIV( J ) 345: IF( JP.LT.0 ) THEN 346: JP = -JP 347: J = J + 1 348: END IF 349: J = J + 1 350: IF( JP.NE.JJ .AND. J.LE.N ) 351: $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA ) 352: IF( J.LE.N ) 353: $ GO TO 60 354: * 355: * Set KB to the number of columns factorized 356: * 357: KB = N - K 358: * 359: ELSE 360: * 361: * Factorize the leading columns of A using the lower triangle 362: * of A and working forwards, and compute the matrix W = L21*D 363: * for use in updating A22 364: * 365: * K is the main loop index, increasing from 1 in steps of 1 or 2 366: * 367: K = 1 368: 70 CONTINUE 369: * 370: * Exit from loop 371: * 372: IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) 373: $ GO TO 90 374: * 375: * Copy column K of A to column K of W and update it 376: * 377: CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 ) 378: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA, 379: $ W( K, 1 ), LDW, ONE, W( K, K ), 1 ) 380: * 381: KSTEP = 1 382: * 383: * Determine rows and columns to be interchanged and whether 384: * a 1-by-1 or 2-by-2 pivot block will be used 385: * 386: ABSAKK = ABS( W( K, K ) ) 387: * 388: * IMAX is the row-index of the largest off-diagonal element in 389: * column K, and COLMAX is its absolute value 390: * 391: IF( K.LT.N ) THEN 392: IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 ) 393: COLMAX = ABS( W( IMAX, K ) ) 394: ELSE 395: COLMAX = ZERO 396: END IF 397: * 398: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 399: * 400: * Column K is zero: set INFO and continue 401: * 402: IF( INFO.EQ.0 ) 403: $ INFO = K 404: KP = K 405: ELSE 406: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 407: * 408: * no interchange, use 1-by-1 pivot block 409: * 410: KP = K 411: ELSE 412: * 413: * Copy column IMAX to column K+1 of W and update it 414: * 415: CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 ) 416: CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ), 417: $ 1 ) 418: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), 419: $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 ) 420: * 421: * JMAX is the column-index of the largest off-diagonal 422: * element in row IMAX, and ROWMAX is its absolute value 423: * 424: JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 ) 425: ROWMAX = ABS( W( JMAX, K+1 ) ) 426: IF( IMAX.LT.N ) THEN 427: JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 ) 428: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) ) 429: END IF 430: * 431: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 432: * 433: * no interchange, use 1-by-1 pivot block 434: * 435: KP = K 436: ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN 437: * 438: * interchange rows and columns K and IMAX, use 1-by-1 439: * pivot block 440: * 441: KP = IMAX 442: * 443: * copy column K+1 of W to column K 444: * 445: CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) 446: ELSE 447: * 448: * interchange rows and columns K+1 and IMAX, use 2-by-2 449: * pivot block 450: * 451: KP = IMAX 452: KSTEP = 2 453: END IF 454: END IF 455: * 456: KK = K + KSTEP - 1 457: * 458: * Updated column KP is already stored in column KK of W 459: * 460: IF( KP.NE.KK ) THEN 461: * 462: * Copy non-updated column KK to column KP 463: * 464: A( KP, K ) = A( KK, K ) 465: CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA ) 466: CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 ) 467: * 468: * Interchange rows KK and KP in first KK columns of A and W 469: * 470: CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) 471: CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) 472: END IF 473: * 474: IF( KSTEP.EQ.1 ) THEN 475: * 476: * 1-by-1 pivot block D(k): column k of W now holds 477: * 478: * W(k) = L(k)*D(k) 479: * 480: * where L(k) is the k-th column of L 481: * 482: * Store L(k) in column k of A 483: * 484: CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) 485: IF( K.LT.N ) THEN 486: R1 = ONE / A( K, K ) 487: CALL DSCAL( N-K, R1, A( K+1, K ), 1 ) 488: END IF 489: ELSE 490: * 491: * 2-by-2 pivot block D(k): columns k and k+1 of W now hold 492: * 493: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 494: * 495: * where L(k) and L(k+1) are the k-th and (k+1)-th columns 496: * of L 497: * 498: IF( K.LT.N-1 ) THEN 499: * 500: * Store L(k) and L(k+1) in columns k and k+1 of A 501: * 502: D21 = W( K+1, K ) 503: D11 = W( K+1, K+1 ) / D21 504: D22 = W( K, K ) / D21 505: T = ONE / ( D11*D22-ONE ) 506: D21 = T / D21 507: DO 80 J = K + 2, N 508: A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) 509: A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) 510: 80 CONTINUE 511: END IF 512: * 513: * Copy D(k) to A 514: * 515: A( K, K ) = W( K, K ) 516: A( K+1, K ) = W( K+1, K ) 517: A( K+1, K+1 ) = W( K+1, K+1 ) 518: END IF 519: END IF 520: * 521: * Store details of the interchanges in IPIV 522: * 523: IF( KSTEP.EQ.1 ) THEN 524: IPIV( K ) = KP 525: ELSE 526: IPIV( K ) = -KP 527: IPIV( K+1 ) = -KP 528: END IF 529: * 530: * Increase K and return to the start of the main loop 531: * 532: K = K + KSTEP 533: GO TO 70 534: * 535: 90 CONTINUE 536: * 537: * Update the lower triangle of A22 (= A(k:n,k:n)) as 538: * 539: * A22 := A22 - L21*D*L21' = A22 - L21*W' 540: * 541: * computing blocks of NB columns at a time 542: * 543: DO 110 J = K, N, NB 544: JB = MIN( NB, N-J+1 ) 545: * 546: * Update the lower triangle of the diagonal block 547: * 548: DO 100 JJ = J, J + JB - 1 549: CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE, 550: $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE, 551: $ A( JJ, JJ ), 1 ) 552: 100 CONTINUE 553: * 554: * Update the rectangular subdiagonal block 555: * 556: IF( J+JB.LE.N ) 557: $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 558: $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW, 559: $ ONE, A( J+JB, J ), LDA ) 560: 110 CONTINUE 561: * 562: * Put L21 in standard form by partially undoing the interchanges 563: * in columns 1:k-1 564: * 565: J = K - 1 566: 120 CONTINUE 567: JJ = J 568: JP = IPIV( J ) 569: IF( JP.LT.0 ) THEN 570: JP = -JP 571: J = J - 1 572: END IF 573: J = J - 1 574: IF( JP.NE.JJ .AND. J.GE.1 ) 575: $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA ) 576: IF( J.GE.1 ) 577: $ GO TO 120 578: * 579: * Set KB to the number of columns factorized 580: * 581: KB = K - 1 582: * 583: END IF 584: RETURN 585: * 586: * End of DLASYF 587: * 588: END