![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.5.0.
1: *> \brief \b DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm). 2: * 3: * =========== DOCUMENTATION =========== 4: * 5: * Online html documentation available at 6: * http://www.netlib.org/lapack/explore-html/ 7: * 8: *> \htmlonly 9: *> Download DSYTF2_ROOK + dependencies 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2_rook.f"> 11: *> [TGZ]</a> 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rook.f"> 13: *> [ZIP]</a> 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rook.f"> 15: *> [TXT]</a> 16: *> \endhtmlonly 17: * 18: * Definition: 19: * =========== 20: * 21: * SUBROUTINE DSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) 22: * 23: * .. Scalar Arguments .. 24: * CHARACTER UPLO 25: * INTEGER INFO, LDA, N 26: * .. 27: * .. Array Arguments .. 28: * INTEGER IPIV( * ) 29: * DOUBLE PRECISION A( LDA, * ) 30: * .. 31: * 32: * 33: *> \par Purpose: 34: * ============= 35: *> 36: *> \verbatim 37: *> 38: *> DSYTF2_ROOK computes the factorization of a real symmetric matrix A 39: *> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method: 40: *> 41: *> A = U*D*U**T or A = L*D*L**T 42: *> 43: *> where U (or L) is a product of permutation and unit upper (lower) 44: *> triangular matrices, U**T is the transpose of U, and D is symmetric and 45: *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks. 46: *> 47: *> This is the unblocked version of the algorithm, calling Level 2 BLAS. 48: *> \endverbatim 49: * 50: * Arguments: 51: * ========== 52: * 53: *> \param[in] UPLO 54: *> \verbatim 55: *> UPLO is CHARACTER*1 56: *> Specifies whether the upper or lower triangular part of the 57: *> symmetric matrix A is stored: 58: *> = 'U': Upper triangular 59: *> = 'L': Lower triangular 60: *> \endverbatim 61: *> 62: *> \param[in] N 63: *> \verbatim 64: *> N is INTEGER 65: *> The order of the matrix A. N >= 0. 66: *> \endverbatim 67: *> 68: *> \param[in,out] A 69: *> \verbatim 70: *> A is DOUBLE PRECISION array, dimension (LDA,N) 71: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading 72: *> n-by-n upper triangular part of A contains the upper 73: *> triangular part of the matrix A, and the strictly lower 74: *> triangular part of A is not referenced. If UPLO = 'L', the 75: *> leading n-by-n lower triangular part of A contains the lower 76: *> triangular part of the matrix A, and the strictly upper 77: *> triangular part of A is not referenced. 78: *> 79: *> On exit, the block diagonal matrix D and the multipliers used 80: *> to obtain the factor U or L (see below for further details). 81: *> \endverbatim 82: *> 83: *> \param[in] LDA 84: *> \verbatim 85: *> LDA is INTEGER 86: *> The leading dimension of the array A. LDA >= max(1,N). 87: *> \endverbatim 88: *> 89: *> \param[out] IPIV 90: *> \verbatim 91: *> IPIV is INTEGER array, dimension (N) 92: *> Details of the interchanges and the block structure of D. 93: *> 94: *> If UPLO = 'U': 95: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) 96: *> were interchanged and D(k,k) is a 1-by-1 diagonal block. 97: *> 98: *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and 99: *> columns k and -IPIV(k) were interchanged and rows and 100: *> columns k-1 and -IPIV(k-1) were inerchaged, 101: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 102: *> 103: *> If UPLO = 'L': 104: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) 105: *> were interchanged and D(k,k) is a 1-by-1 diagonal block. 106: *> 107: *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and 108: *> columns k and -IPIV(k) were interchanged and rows and 109: *> columns k+1 and -IPIV(k+1) were inerchaged, 110: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 111: *> \endverbatim 112: *> 113: *> \param[out] INFO 114: *> \verbatim 115: *> INFO is INTEGER 116: *> = 0: successful exit 117: *> < 0: if INFO = -k, the k-th argument had an illegal value 118: *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization 119: *> has been completed, but the block diagonal matrix D is 120: *> exactly singular, and division by zero will occur if it 121: *> is used to solve a system of equations. 122: *> \endverbatim 123: * 124: * Authors: 125: * ======== 126: * 127: *> \author Univ. of Tennessee 128: *> \author Univ. of California Berkeley 129: *> \author Univ. of Colorado Denver 130: *> \author NAG Ltd. 131: * 132: *> \date November 2013 133: * 134: *> \ingroup doubleSYcomputational 135: * 136: *> \par Further Details: 137: * ===================== 138: *> 139: *> \verbatim 140: *> 141: *> If UPLO = 'U', then A = U*D*U**T, where 142: *> U = P(n)*U(n)* ... *P(k)U(k)* ..., 143: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to 144: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 145: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 146: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 147: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 148: *> 149: *> ( I v 0 ) k-s 150: *> U(k) = ( 0 I 0 ) s 151: *> ( 0 0 I ) n-k 152: *> k-s s n-k 153: *> 154: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 155: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 156: *> and A(k,k), and v overwrites A(1:k-2,k-1:k). 157: *> 158: *> If UPLO = 'L', then A = L*D*L**T, where 159: *> L = P(1)*L(1)* ... *P(k)*L(k)* ..., 160: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 161: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 162: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 163: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 164: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 165: *> 166: *> ( I 0 0 ) k-1 167: *> L(k) = ( 0 I 0 ) s 168: *> ( 0 v I ) n-k-s+1 169: *> k-1 s n-k-s+1 170: *> 171: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 172: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 173: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 174: *> \endverbatim 175: * 176: *> \par Contributors: 177: * ================== 178: *> 179: *> \verbatim 180: *> 181: *> November 2013, Igor Kozachenko, 182: *> Computer Science Division, 183: *> University of California, Berkeley 184: *> 185: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 186: *> School of Mathematics, 187: *> University of Manchester 188: *> 189: *> 01-01-96 - Based on modifications by 190: *> J. Lewis, Boeing Computer Services Company 191: *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA 192: *> \endverbatim 193: * 194: * ===================================================================== 195: SUBROUTINE DSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) 196: * 197: * -- LAPACK computational routine (version 3.5.0) -- 198: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 199: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 200: * November 2013 201: * 202: * .. Scalar Arguments .. 203: CHARACTER UPLO 204: INTEGER INFO, LDA, N 205: * .. 206: * .. Array Arguments .. 207: INTEGER IPIV( * ) 208: DOUBLE PRECISION A( LDA, * ) 209: * .. 210: * 211: * ===================================================================== 212: * 213: * .. Parameters .. 214: DOUBLE PRECISION ZERO, ONE 215: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 216: DOUBLE PRECISION EIGHT, SEVTEN 217: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 218: * .. 219: * .. Local Scalars .. 220: LOGICAL UPPER, DONE 221: INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP, 222: $ P, II 223: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, 224: $ ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN 225: * .. 226: * .. External Functions .. 227: LOGICAL LSAME 228: INTEGER IDAMAX 229: DOUBLE PRECISION DLAMCH 230: EXTERNAL LSAME, IDAMAX, DLAMCH 231: * .. 232: * .. External Subroutines .. 233: EXTERNAL DSCAL, DSWAP, DSYR, XERBLA 234: * .. 235: * .. Intrinsic Functions .. 236: INTRINSIC ABS, MAX, SQRT 237: * .. 238: * .. Executable Statements .. 239: * 240: * Test the input parameters. 241: * 242: INFO = 0 243: UPPER = LSAME( UPLO, 'U' ) 244: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 245: INFO = -1 246: ELSE IF( N.LT.0 ) THEN 247: INFO = -2 248: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 249: INFO = -4 250: END IF 251: IF( INFO.NE.0 ) THEN 252: CALL XERBLA( 'DSYTF2_ROOK', -INFO ) 253: RETURN 254: END IF 255: * 256: * Initialize ALPHA for use in choosing pivot block size. 257: * 258: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 259: * 260: * Compute machine safe minimum 261: * 262: SFMIN = DLAMCH( 'S' ) 263: * 264: IF( UPPER ) THEN 265: * 266: * Factorize A as U*D*U**T using the upper triangle of A 267: * 268: * K is the main loop index, decreasing from N to 1 in steps of 269: * 1 or 2 270: * 271: K = N 272: 10 CONTINUE 273: * 274: * If K < 1, exit from loop 275: * 276: IF( K.LT.1 ) 277: $ GO TO 70 278: KSTEP = 1 279: P = K 280: * 281: * Determine rows and columns to be interchanged and whether 282: * a 1-by-1 or 2-by-2 pivot block will be used 283: * 284: ABSAKK = ABS( A( K, K ) ) 285: * 286: * IMAX is the row-index of the largest off-diagonal element in 287: * column K, and COLMAX is its absolute value. 288: * Determine both COLMAX and IMAX. 289: * 290: IF( K.GT.1 ) THEN 291: IMAX = IDAMAX( K-1, A( 1, K ), 1 ) 292: COLMAX = ABS( A( IMAX, K ) ) 293: ELSE 294: COLMAX = ZERO 295: END IF 296: * 297: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN 298: * 299: * Column K is zero or underflow: set INFO and continue 300: * 301: IF( INFO.EQ.0 ) 302: $ INFO = K 303: KP = K 304: ELSE 305: * 306: * Test for interchange 307: * 308: * Equivalent to testing for (used to handle NaN and Inf) 309: * ABSAKK.GE.ALPHA*COLMAX 310: * 311: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 312: * 313: * no interchange, 314: * use 1-by-1 pivot block 315: * 316: KP = K 317: ELSE 318: * 319: DONE = .FALSE. 320: * 321: * Loop until pivot found 322: * 323: 12 CONTINUE 324: * 325: * Begin pivot search loop body 326: * 327: * JMAX is the column-index of the largest off-diagonal 328: * element in row IMAX, and ROWMAX is its absolute value. 329: * Determine both ROWMAX and JMAX. 330: * 331: IF( IMAX.NE.K ) THEN 332: JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), 333: $ LDA ) 334: ROWMAX = ABS( A( IMAX, JMAX ) ) 335: ELSE 336: ROWMAX = ZERO 337: END IF 338: * 339: IF( IMAX.GT.1 ) THEN 340: ITEMP = IDAMAX( IMAX-1, A( 1, IMAX ), 1 ) 341: DTEMP = ABS( A( ITEMP, IMAX ) ) 342: IF( DTEMP.GT.ROWMAX ) THEN 343: ROWMAX = DTEMP 344: JMAX = ITEMP 345: END IF 346: END IF 347: * 348: * Equivalent to testing for (used to handle NaN and Inf) 349: * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX 350: * 351: IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) ) 352: $ THEN 353: * 354: * interchange rows and columns K and IMAX, 355: * use 1-by-1 pivot block 356: * 357: KP = IMAX 358: DONE = .TRUE. 359: * 360: * Equivalent to testing for ROWMAX .EQ. COLMAX, 361: * used to handle NaN and Inf 362: * 363: ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN 364: * 365: * interchange rows and columns K+1 and IMAX, 366: * use 2-by-2 pivot block 367: * 368: KP = IMAX 369: KSTEP = 2 370: DONE = .TRUE. 371: ELSE 372: * 373: * Pivot NOT found, set variables and repeat 374: * 375: P = IMAX 376: COLMAX = ROWMAX 377: IMAX = JMAX 378: END IF 379: * 380: * End pivot search loop body 381: * 382: IF( .NOT. DONE ) GOTO 12 383: * 384: END IF 385: * 386: * Swap TWO rows and TWO columns 387: * 388: * First swap 389: * 390: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 391: * 392: * Interchange rows and column K and P in the leading 393: * submatrix A(1:k,1:k) if we have a 2-by-2 pivot 394: * 395: IF( P.GT.1 ) 396: $ CALL DSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 ) 397: IF( P.LT.(K-1) ) 398: $ CALL DSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ), 399: $ LDA ) 400: T = A( K, K ) 401: A( K, K ) = A( P, P ) 402: A( P, P ) = T 403: END IF 404: * 405: * Second swap 406: * 407: KK = K - KSTEP + 1 408: IF( KP.NE.KK ) THEN 409: * 410: * Interchange rows and columns KK and KP in the leading 411: * submatrix A(1:k,1:k) 412: * 413: IF( KP.GT.1 ) 414: $ CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 415: IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) ) 416: $ CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), 417: $ LDA ) 418: T = A( KK, KK ) 419: A( KK, KK ) = A( KP, KP ) 420: A( KP, KP ) = T 421: IF( KSTEP.EQ.2 ) THEN 422: T = A( K-1, K ) 423: A( K-1, K ) = A( KP, K ) 424: A( KP, K ) = T 425: END IF 426: END IF 427: * 428: * Update the leading submatrix 429: * 430: IF( KSTEP.EQ.1 ) THEN 431: * 432: * 1-by-1 pivot block D(k): column k now holds 433: * 434: * W(k) = U(k)*D(k) 435: * 436: * where U(k) is the k-th column of U 437: * 438: IF( K.GT.1 ) THEN 439: * 440: * Perform a rank-1 update of A(1:k-1,1:k-1) and 441: * store U(k) in column k 442: * 443: IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 444: * 445: * Perform a rank-1 update of A(1:k-1,1:k-1) as 446: * A := A - U(k)*D(k)*U(k)**T 447: * = A - W(k)*1/D(k)*W(k)**T 448: * 449: D11 = ONE / A( K, K ) 450: CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 451: * 452: * Store U(k) in column k 453: * 454: CALL DSCAL( K-1, D11, A( 1, K ), 1 ) 455: ELSE 456: * 457: * Store L(k) in column K 458: * 459: D11 = A( K, K ) 460: DO 16 II = 1, K - 1 461: A( II, K ) = A( II, K ) / D11 462: 16 CONTINUE 463: * 464: * Perform a rank-1 update of A(k+1:n,k+1:n) as 465: * A := A - U(k)*D(k)*U(k)**T 466: * = A - W(k)*(1/D(k))*W(k)**T 467: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 468: * 469: CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA ) 470: END IF 471: END IF 472: * 473: ELSE 474: * 475: * 2-by-2 pivot block D(k): columns k and k-1 now hold 476: * 477: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 478: * 479: * where U(k) and U(k-1) are the k-th and (k-1)-th columns 480: * of U 481: * 482: * Perform a rank-2 update of A(1:k-2,1:k-2) as 483: * 484: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T 485: * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T 486: * 487: * and store L(k) and L(k+1) in columns k and k+1 488: * 489: IF( K.GT.2 ) THEN 490: * 491: D12 = A( K-1, K ) 492: D22 = A( K-1, K-1 ) / D12 493: D11 = A( K, K ) / D12 494: T = ONE / ( D11*D22-ONE ) 495: * 496: DO 30 J = K - 2, 1, -1 497: * 498: WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) 499: WK = T*( D22*A( J, K )-A( J, K-1 ) ) 500: * 501: DO 20 I = J, 1, -1 502: A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - 503: $ ( A( I, K-1 ) / D12 )*WKM1 504: 20 CONTINUE 505: * 506: * Store U(k) and U(k-1) in cols k and k-1 for row J 507: * 508: A( J, K ) = WK / D12 509: A( J, K-1 ) = WKM1 / D12 510: * 511: 30 CONTINUE 512: * 513: END IF 514: * 515: END IF 516: END IF 517: * 518: * Store details of the interchanges in IPIV 519: * 520: IF( KSTEP.EQ.1 ) THEN 521: IPIV( K ) = KP 522: ELSE 523: IPIV( K ) = -P 524: IPIV( K-1 ) = -KP 525: END IF 526: * 527: * Decrease K and return to the start of the main loop 528: * 529: K = K - KSTEP 530: GO TO 10 531: * 532: ELSE 533: * 534: * Factorize A as L*D*L**T using the lower triangle of A 535: * 536: * K is the main loop index, increasing from 1 to N in steps of 537: * 1 or 2 538: * 539: K = 1 540: 40 CONTINUE 541: * 542: * If K > N, exit from loop 543: * 544: IF( K.GT.N ) 545: $ GO TO 70 546: KSTEP = 1 547: P = K 548: * 549: * Determine rows and columns to be interchanged and whether 550: * a 1-by-1 or 2-by-2 pivot block will be used 551: * 552: ABSAKK = ABS( A( K, K ) ) 553: * 554: * IMAX is the row-index of the largest off-diagonal element in 555: * column K, and COLMAX is its absolute value. 556: * Determine both COLMAX and IMAX. 557: * 558: IF( K.LT.N ) THEN 559: IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 ) 560: COLMAX = ABS( A( IMAX, K ) ) 561: ELSE 562: COLMAX = ZERO 563: END IF 564: * 565: IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN 566: * 567: * Column K is zero or underflow: set INFO and continue 568: * 569: IF( INFO.EQ.0 ) 570: $ INFO = K 571: KP = K 572: ELSE 573: * 574: * Test for interchange 575: * 576: * Equivalent to testing for (used to handle NaN and Inf) 577: * ABSAKK.GE.ALPHA*COLMAX 578: * 579: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN 580: * 581: * no interchange, use 1-by-1 pivot block 582: * 583: KP = K 584: ELSE 585: * 586: DONE = .FALSE. 587: * 588: * Loop until pivot found 589: * 590: 42 CONTINUE 591: * 592: * Begin pivot search loop body 593: * 594: * JMAX is the column-index of the largest off-diagonal 595: * element in row IMAX, and ROWMAX is its absolute value. 596: * Determine both ROWMAX and JMAX. 597: * 598: IF( IMAX.NE.K ) THEN 599: JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA ) 600: ROWMAX = ABS( A( IMAX, JMAX ) ) 601: ELSE 602: ROWMAX = ZERO 603: END IF 604: * 605: IF( IMAX.LT.N ) THEN 606: ITEMP = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 607: $ 1 ) 608: DTEMP = ABS( A( ITEMP, IMAX ) ) 609: IF( DTEMP.GT.ROWMAX ) THEN 610: ROWMAX = DTEMP 611: JMAX = ITEMP 612: END IF 613: END IF 614: * 615: * Equivalent to testing for (used to handle NaN and Inf) 616: * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX 617: * 618: IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) ) 619: $ THEN 620: * 621: * interchange rows and columns K and IMAX, 622: * use 1-by-1 pivot block 623: * 624: KP = IMAX 625: DONE = .TRUE. 626: * 627: * Equivalent to testing for ROWMAX .EQ. COLMAX, 628: * used to handle NaN and Inf 629: * 630: ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN 631: * 632: * interchange rows and columns K+1 and IMAX, 633: * use 2-by-2 pivot block 634: * 635: KP = IMAX 636: KSTEP = 2 637: DONE = .TRUE. 638: ELSE 639: * 640: * Pivot NOT found, set variables and repeat 641: * 642: P = IMAX 643: COLMAX = ROWMAX 644: IMAX = JMAX 645: END IF 646: * 647: * End pivot search loop body 648: * 649: IF( .NOT. DONE ) GOTO 42 650: * 651: END IF 652: * 653: * Swap TWO rows and TWO columns 654: * 655: * First swap 656: * 657: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN 658: * 659: * Interchange rows and column K and P in the trailing 660: * submatrix A(k:n,k:n) if we have a 2-by-2 pivot 661: * 662: IF( P.LT.N ) 663: $ CALL DSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 ) 664: IF( P.GT.(K+1) ) 665: $ CALL DSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA ) 666: T = A( K, K ) 667: A( K, K ) = A( P, P ) 668: A( P, P ) = T 669: END IF 670: * 671: * Second swap 672: * 673: KK = K + KSTEP - 1 674: IF( KP.NE.KK ) THEN 675: * 676: * Interchange rows and columns KK and KP in the trailing 677: * submatrix A(k:n,k:n) 678: * 679: IF( KP.LT.N ) 680: $ CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 681: IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) ) 682: $ CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 683: $ LDA ) 684: T = A( KK, KK ) 685: A( KK, KK ) = A( KP, KP ) 686: A( KP, KP ) = T 687: IF( KSTEP.EQ.2 ) THEN 688: T = A( K+1, K ) 689: A( K+1, K ) = A( KP, K ) 690: A( KP, K ) = T 691: END IF 692: END IF 693: * 694: * Update the trailing submatrix 695: * 696: IF( KSTEP.EQ.1 ) THEN 697: * 698: * 1-by-1 pivot block D(k): column k now holds 699: * 700: * W(k) = L(k)*D(k) 701: * 702: * where L(k) is the k-th column of L 703: * 704: IF( K.LT.N ) THEN 705: * 706: * Perform a rank-1 update of A(k+1:n,k+1:n) and 707: * store L(k) in column k 708: * 709: IF( ABS( A( K, K ) ).GE.SFMIN ) THEN 710: * 711: * Perform a rank-1 update of A(k+1:n,k+1:n) as 712: * A := A - L(k)*D(k)*L(k)**T 713: * = A - W(k)*(1/D(k))*W(k)**T 714: * 715: D11 = ONE / A( K, K ) 716: CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1, 717: $ A( K+1, K+1 ), LDA ) 718: * 719: * Store L(k) in column k 720: * 721: CALL DSCAL( N-K, D11, A( K+1, K ), 1 ) 722: ELSE 723: * 724: * Store L(k) in column k 725: * 726: D11 = A( K, K ) 727: DO 46 II = K + 1, N 728: A( II, K ) = A( II, K ) / D11 729: 46 CONTINUE 730: * 731: * Perform a rank-1 update of A(k+1:n,k+1:n) as 732: * A := A - L(k)*D(k)*L(k)**T 733: * = A - W(k)*(1/D(k))*W(k)**T 734: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T 735: * 736: CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1, 737: $ A( K+1, K+1 ), LDA ) 738: END IF 739: END IF 740: * 741: ELSE 742: * 743: * 2-by-2 pivot block D(k): columns k and k+1 now hold 744: * 745: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 746: * 747: * where L(k) and L(k+1) are the k-th and (k+1)-th columns 748: * of L 749: * 750: * 751: * Perform a rank-2 update of A(k+2:n,k+2:n) as 752: * 753: * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T 754: * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T 755: * 756: * and store L(k) and L(k+1) in columns k and k+1 757: * 758: IF( K.LT.N-1 ) THEN 759: * 760: D21 = A( K+1, K ) 761: D11 = A( K+1, K+1 ) / D21 762: D22 = A( K, K ) / D21 763: T = ONE / ( D11*D22-ONE ) 764: * 765: DO 60 J = K + 2, N 766: * 767: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J 768: * 769: WK = T*( D11*A( J, K )-A( J, K+1 ) ) 770: WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) 771: * 772: * Perform a rank-2 update of A(k+2:n,k+2:n) 773: * 774: DO 50 I = J, N 775: A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - 776: $ ( A( I, K+1 ) / D21 )*WKP1 777: 50 CONTINUE 778: * 779: * Store L(k) and L(k+1) in cols k and k+1 for row J 780: * 781: A( J, K ) = WK / D21 782: A( J, K+1 ) = WKP1 / D21 783: * 784: 60 CONTINUE 785: * 786: END IF 787: * 788: END IF 789: END IF 790: * 791: * Store details of the interchanges in IPIV 792: * 793: IF( KSTEP.EQ.1 ) THEN 794: IPIV( K ) = KP 795: ELSE 796: IPIV( K ) = -P 797: IPIV( K+1 ) = -KP 798: END IF 799: * 800: * Increase K and return to the start of the main loop 801: * 802: K = K + KSTEP 803: GO TO 40 804: * 805: END IF 806: * 807: 70 CONTINUE 808: * 809: RETURN 810: * 811: * End of DSYTF2_ROOK 812: * 813: END