![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.5.0.
1: *> \brief \b ZSYTF2 computes the factorization of a real symmetric indefinite matrix, using the 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 ZSYTF2 + dependencies 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2.f"> 11: *> [TGZ]</a> 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2.f"> 13: *> [ZIP]</a> 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2.f"> 15: *> [TXT]</a> 16: *> \endhtmlonly 17: * 18: * Definition: 19: * =========== 20: * 21: * SUBROUTINE ZSYTF2( 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: * COMPLEX*16 A( LDA, * ) 30: * .. 31: * 32: * 33: *> \par Purpose: 34: * ============= 35: *> 36: *> \verbatim 37: *> 38: *> ZSYTF2 computes the factorization of a complex symmetric matrix A 39: *> using the Bunch-Kaufman 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 COMPLEX*16 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) were 96: *> interchanged and D(k,k) is a 1-by-1 diagonal block. 97: *> 98: *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns 99: *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 100: *> is a 2-by-2 diagonal block. 101: *> 102: *> If UPLO = 'L': 103: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 104: *> interchanged and D(k,k) is a 1-by-1 diagonal block. 105: *> 106: *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns 107: *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) 108: *> is a 2-by-2 diagonal block. 109: *> \endverbatim 110: *> 111: *> \param[out] INFO 112: *> \verbatim 113: *> INFO is INTEGER 114: *> = 0: successful exit 115: *> < 0: if INFO = -k, the k-th argument had an illegal value 116: *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization 117: *> has been completed, but the block diagonal matrix D is 118: *> exactly singular, and division by zero will occur if it 119: *> is used to solve a system of equations. 120: *> \endverbatim 121: * 122: * Authors: 123: * ======== 124: * 125: *> \author Univ. of Tennessee 126: *> \author Univ. of California Berkeley 127: *> \author Univ. of Colorado Denver 128: *> \author NAG Ltd. 129: * 130: *> \date November 2013 131: * 132: *> \ingroup complex16SYcomputational 133: * 134: *> \par Further Details: 135: * ===================== 136: *> 137: *> \verbatim 138: *> 139: *> If UPLO = 'U', then A = U*D*U**T, where 140: *> U = P(n)*U(n)* ... *P(k)U(k)* ..., 141: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to 142: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 143: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 144: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 145: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 146: *> 147: *> ( I v 0 ) k-s 148: *> U(k) = ( 0 I 0 ) s 149: *> ( 0 0 I ) n-k 150: *> k-s s n-k 151: *> 152: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 153: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 154: *> and A(k,k), and v overwrites A(1:k-2,k-1:k). 155: *> 156: *> If UPLO = 'L', then A = L*D*L**T, where 157: *> L = P(1)*L(1)* ... *P(k)*L(k)* ..., 158: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 159: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 160: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 161: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 162: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 163: *> 164: *> ( I 0 0 ) k-1 165: *> L(k) = ( 0 I 0 ) s 166: *> ( 0 v I ) n-k-s+1 167: *> k-1 s n-k-s+1 168: *> 169: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 170: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 171: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 172: *> \endverbatim 173: * 174: *> \par Contributors: 175: * ================== 176: *> 177: *> \verbatim 178: *> 179: *> 09-29-06 - patch from 180: *> Bobby Cheng, MathWorks 181: *> 182: *> Replace l.209 and l.377 183: *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 184: *> by 185: *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN 186: *> 187: *> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services 188: *> Company 189: *> \endverbatim 190: * 191: * ===================================================================== 192: SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO ) 193: * 194: * -- LAPACK computational routine (version 3.5.0) -- 195: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 196: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 197: * November 2013 198: * 199: * .. Scalar Arguments .. 200: CHARACTER UPLO 201: INTEGER INFO, LDA, N 202: * .. 203: * .. Array Arguments .. 204: INTEGER IPIV( * ) 205: COMPLEX*16 A( LDA, * ) 206: * .. 207: * 208: * ===================================================================== 209: * 210: * .. Parameters .. 211: DOUBLE PRECISION ZERO, ONE 212: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 213: DOUBLE PRECISION EIGHT, SEVTEN 214: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 215: COMPLEX*16 CONE 216: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 217: * .. 218: * .. Local Scalars .. 219: LOGICAL UPPER 220: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP 221: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX 222: COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z 223: * .. 224: * .. External Functions .. 225: LOGICAL DISNAN, LSAME 226: INTEGER IZAMAX 227: EXTERNAL DISNAN, LSAME, IZAMAX 228: * .. 229: * .. External Subroutines .. 230: EXTERNAL XERBLA, ZSCAL, ZSWAP, ZSYR 231: * .. 232: * .. Intrinsic Functions .. 233: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 234: * .. 235: * .. Statement Functions .. 236: DOUBLE PRECISION CABS1 237: * .. 238: * .. Statement Function definitions .. 239: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) ) 240: * .. 241: * .. Executable Statements .. 242: * 243: * Test the input parameters. 244: * 245: INFO = 0 246: UPPER = LSAME( UPLO, 'U' ) 247: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 248: INFO = -1 249: ELSE IF( N.LT.0 ) THEN 250: INFO = -2 251: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 252: INFO = -4 253: END IF 254: IF( INFO.NE.0 ) THEN 255: CALL XERBLA( 'ZSYTF2', -INFO ) 256: RETURN 257: END IF 258: * 259: * Initialize ALPHA for use in choosing pivot block size. 260: * 261: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 262: * 263: IF( UPPER ) THEN 264: * 265: * Factorize A as U*D*U**T using the upper triangle of A 266: * 267: * K is the main loop index, decreasing from N to 1 in steps of 268: * 1 or 2 269: * 270: K = N 271: 10 CONTINUE 272: * 273: * If K < 1, exit from loop 274: * 275: IF( K.LT.1 ) 276: $ GO TO 70 277: KSTEP = 1 278: * 279: * Determine rows and columns to be interchanged and whether 280: * a 1-by-1 or 2-by-2 pivot block will be used 281: * 282: ABSAKK = CABS1( A( K, K ) ) 283: * 284: * IMAX is the row-index of the largest off-diagonal element in 285: * column K, and COLMAX is its absolute value. 286: * Determine both COLMAX and IMAX. 287: * 288: IF( K.GT.1 ) THEN 289: IMAX = IZAMAX( K-1, A( 1, K ), 1 ) 290: COLMAX = CABS1( A( IMAX, K ) ) 291: ELSE 292: COLMAX = ZERO 293: END IF 294: * 295: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN 296: * 297: * Column K is zero or underflow, or contains a NaN: 298: * set INFO and continue 299: * 300: IF( INFO.EQ.0 ) 301: $ INFO = K 302: KP = K 303: ELSE 304: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 305: * 306: * no interchange, use 1-by-1 pivot block 307: * 308: KP = K 309: ELSE 310: * 311: * JMAX is the column-index of the largest off-diagonal 312: * element in row IMAX, and ROWMAX is its absolute value 313: * 314: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA ) 315: ROWMAX = CABS1( A( IMAX, JMAX ) ) 316: IF( IMAX.GT.1 ) THEN 317: JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 ) 318: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 319: END IF 320: * 321: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 322: * 323: * no interchange, use 1-by-1 pivot block 324: * 325: KP = K 326: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 327: * 328: * interchange rows and columns K and IMAX, use 1-by-1 329: * pivot block 330: * 331: KP = IMAX 332: ELSE 333: * 334: * interchange rows and columns K-1 and IMAX, use 2-by-2 335: * pivot block 336: * 337: KP = IMAX 338: KSTEP = 2 339: END IF 340: END IF 341: * 342: KK = K - KSTEP + 1 343: IF( KP.NE.KK ) THEN 344: * 345: * Interchange rows and columns KK and KP in the leading 346: * submatrix A(1:k,1:k) 347: * 348: CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 349: CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), 350: $ LDA ) 351: T = A( KK, KK ) 352: A( KK, KK ) = A( KP, KP ) 353: A( KP, KP ) = T 354: IF( KSTEP.EQ.2 ) THEN 355: T = A( K-1, K ) 356: A( K-1, K ) = A( KP, K ) 357: A( KP, K ) = T 358: END IF 359: END IF 360: * 361: * Update the leading submatrix 362: * 363: IF( KSTEP.EQ.1 ) THEN 364: * 365: * 1-by-1 pivot block D(k): column k now holds 366: * 367: * W(k) = U(k)*D(k) 368: * 369: * where U(k) is the k-th column of U 370: * 371: * Perform a rank-1 update of A(1:k-1,1:k-1) as 372: * 373: * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T 374: * 375: R1 = CONE / A( K, K ) 376: CALL ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA ) 377: * 378: * Store U(k) in column k 379: * 380: CALL ZSCAL( K-1, R1, A( 1, K ), 1 ) 381: ELSE 382: * 383: * 2-by-2 pivot block D(k): columns k and k-1 now hold 384: * 385: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 386: * 387: * where U(k) and U(k-1) are the k-th and (k-1)-th columns 388: * of U 389: * 390: * Perform a rank-2 update of A(1:k-2,1:k-2) as 391: * 392: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T 393: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T 394: * 395: IF( K.GT.2 ) THEN 396: * 397: D12 = A( K-1, K ) 398: D22 = A( K-1, K-1 ) / D12 399: D11 = A( K, K ) / D12 400: T = CONE / ( D11*D22-CONE ) 401: D12 = T / D12 402: * 403: DO 30 J = K - 2, 1, -1 404: WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) 405: WK = D12*( D22*A( J, K )-A( J, K-1 ) ) 406: DO 20 I = J, 1, -1 407: A( I, J ) = A( I, J ) - A( I, K )*WK - 408: $ A( I, K-1 )*WKM1 409: 20 CONTINUE 410: A( J, K ) = WK 411: A( J, K-1 ) = WKM1 412: 30 CONTINUE 413: * 414: END IF 415: * 416: END IF 417: END IF 418: * 419: * Store details of the interchanges in IPIV 420: * 421: IF( KSTEP.EQ.1 ) THEN 422: IPIV( K ) = KP 423: ELSE 424: IPIV( K ) = -KP 425: IPIV( K-1 ) = -KP 426: END IF 427: * 428: * Decrease K and return to the start of the main loop 429: * 430: K = K - KSTEP 431: GO TO 10 432: * 433: ELSE 434: * 435: * Factorize A as L*D*L**T using the lower triangle of A 436: * 437: * K is the main loop index, increasing from 1 to N in steps of 438: * 1 or 2 439: * 440: K = 1 441: 40 CONTINUE 442: * 443: * If K > N, exit from loop 444: * 445: IF( K.GT.N ) 446: $ GO TO 70 447: KSTEP = 1 448: * 449: * Determine rows and columns to be interchanged and whether 450: * a 1-by-1 or 2-by-2 pivot block will be used 451: * 452: ABSAKK = CABS1( A( K, K ) ) 453: * 454: * IMAX is the row-index of the largest off-diagonal element in 455: * column K, and COLMAX is its absolute value. 456: * Determine both COLMAX and IMAX. 457: * 458: IF( K.LT.N ) THEN 459: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 ) 460: COLMAX = CABS1( A( IMAX, K ) ) 461: ELSE 462: COLMAX = ZERO 463: END IF 464: * 465: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN 466: * 467: * Column K is zero or underflow, or contains a NaN: 468: * set INFO and continue 469: * 470: IF( INFO.EQ.0 ) 471: $ INFO = K 472: KP = K 473: ELSE 474: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 475: * 476: * no interchange, use 1-by-1 pivot block 477: * 478: KP = K 479: ELSE 480: * 481: * JMAX is the column-index of the largest off-diagonal 482: * element in row IMAX, and ROWMAX is its absolute value 483: * 484: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA ) 485: ROWMAX = CABS1( A( IMAX, JMAX ) ) 486: IF( IMAX.LT.N ) THEN 487: JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 ) 488: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 489: END IF 490: * 491: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 492: * 493: * no interchange, use 1-by-1 pivot block 494: * 495: KP = K 496: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 497: * 498: * interchange rows and columns K and IMAX, use 1-by-1 499: * pivot block 500: * 501: KP = IMAX 502: ELSE 503: * 504: * interchange rows and columns K+1 and IMAX, use 2-by-2 505: * pivot block 506: * 507: KP = IMAX 508: KSTEP = 2 509: END IF 510: END IF 511: * 512: KK = K + KSTEP - 1 513: IF( KP.NE.KK ) THEN 514: * 515: * Interchange rows and columns KK and KP in the trailing 516: * submatrix A(k:n,k:n) 517: * 518: IF( KP.LT.N ) 519: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 520: CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 521: $ LDA ) 522: T = A( KK, KK ) 523: A( KK, KK ) = A( KP, KP ) 524: A( KP, KP ) = T 525: IF( KSTEP.EQ.2 ) THEN 526: T = A( K+1, K ) 527: A( K+1, K ) = A( KP, K ) 528: A( KP, K ) = T 529: END IF 530: END IF 531: * 532: * Update the trailing submatrix 533: * 534: IF( KSTEP.EQ.1 ) THEN 535: * 536: * 1-by-1 pivot block D(k): column k now holds 537: * 538: * W(k) = L(k)*D(k) 539: * 540: * where L(k) is the k-th column of L 541: * 542: IF( K.LT.N ) THEN 543: * 544: * Perform a rank-1 update of A(k+1:n,k+1:n) as 545: * 546: * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T 547: * 548: R1 = CONE / A( K, K ) 549: CALL ZSYR( UPLO, N-K, -R1, A( K+1, K ), 1, 550: $ A( K+1, K+1 ), LDA ) 551: * 552: * Store L(k) in column K 553: * 554: CALL ZSCAL( N-K, R1, A( K+1, K ), 1 ) 555: END IF 556: ELSE 557: * 558: * 2-by-2 pivot block D(k) 559: * 560: IF( K.LT.N-1 ) THEN 561: * 562: * Perform a rank-2 update of A(k+2:n,k+2:n) as 563: * 564: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T 565: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T 566: * 567: * where L(k) and L(k+1) are the k-th and (k+1)-th 568: * columns of L 569: * 570: D21 = A( K+1, K ) 571: D11 = A( K+1, K+1 ) / D21 572: D22 = A( K, K ) / D21 573: T = CONE / ( D11*D22-CONE ) 574: D21 = T / D21 575: * 576: DO 60 J = K + 2, N 577: WK = D21*( D11*A( J, K )-A( J, K+1 ) ) 578: WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) 579: DO 50 I = J, N 580: A( I, J ) = A( I, J ) - A( I, K )*WK - 581: $ A( I, K+1 )*WKP1 582: 50 CONTINUE 583: A( J, K ) = WK 584: A( J, K+1 ) = WKP1 585: 60 CONTINUE 586: END IF 587: END IF 588: END IF 589: * 590: * Store details of the interchanges in IPIV 591: * 592: IF( KSTEP.EQ.1 ) THEN 593: IPIV( K ) = KP 594: ELSE 595: IPIV( K ) = -KP 596: IPIV( K+1 ) = -KP 597: END IF 598: * 599: * Increase K and return to the start of the main loop 600: * 601: K = K + KSTEP 602: GO TO 40 603: * 604: END IF 605: * 606: 70 CONTINUE 607: RETURN 608: * 609: * End of ZSYTF2 610: * 611: END