![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZSYTF2( 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: COMPLEX*16 A( LDA, * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZSYTF2 computes the factorization of a complex symmetric matrix A 21: * using 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) COMPLEX*16 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.209 and l.377 83: * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 84: * by 85: * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN 86: * 87: * 1-96 - Based on modifications by J. Lewis, Boeing Computer Services 88: * Company 89: * 90: * If UPLO = 'U', then A = U*D*U', where 91: * U = P(n)*U(n)* ... *P(k)U(k)* ..., 92: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to 93: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 94: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 95: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 96: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 97: * 98: * ( I v 0 ) k-s 99: * U(k) = ( 0 I 0 ) s 100: * ( 0 0 I ) n-k 101: * k-s s n-k 102: * 103: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 104: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 105: * and A(k,k), and v overwrites A(1:k-2,k-1:k). 106: * 107: * If UPLO = 'L', then A = L*D*L', where 108: * L = P(1)*L(1)* ... *P(k)*L(k)* ..., 109: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 110: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 111: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 112: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 113: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 114: * 115: * ( I 0 0 ) k-1 116: * L(k) = ( 0 I 0 ) s 117: * ( 0 v I ) n-k-s+1 118: * k-1 s n-k-s+1 119: * 120: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 121: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 122: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 123: * 124: * ===================================================================== 125: * 126: * .. Parameters .. 127: DOUBLE PRECISION ZERO, ONE 128: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 129: DOUBLE PRECISION EIGHT, SEVTEN 130: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 131: COMPLEX*16 CONE 132: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 133: * .. 134: * .. Local Scalars .. 135: LOGICAL UPPER 136: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP 137: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX 138: COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z 139: * .. 140: * .. External Functions .. 141: LOGICAL DISNAN, LSAME 142: INTEGER IZAMAX 143: EXTERNAL DISNAN, LSAME, IZAMAX 144: * .. 145: * .. External Subroutines .. 146: EXTERNAL XERBLA, ZSCAL, ZSWAP, ZSYR 147: * .. 148: * .. Intrinsic Functions .. 149: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 150: * .. 151: * .. Statement Functions .. 152: DOUBLE PRECISION CABS1 153: * .. 154: * .. Statement Function definitions .. 155: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) ) 156: * .. 157: * .. Executable Statements .. 158: * 159: * Test the input parameters. 160: * 161: INFO = 0 162: UPPER = LSAME( UPLO, 'U' ) 163: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 164: INFO = -1 165: ELSE IF( N.LT.0 ) THEN 166: INFO = -2 167: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 168: INFO = -4 169: END IF 170: IF( INFO.NE.0 ) THEN 171: CALL XERBLA( 'ZSYTF2', -INFO ) 172: RETURN 173: END IF 174: * 175: * Initialize ALPHA for use in choosing pivot block size. 176: * 177: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 178: * 179: IF( UPPER ) THEN 180: * 181: * Factorize A as U*D*U' using the upper triangle of A 182: * 183: * K is the main loop index, decreasing from N to 1 in steps of 184: * 1 or 2 185: * 186: K = N 187: 10 CONTINUE 188: * 189: * If K < 1, exit from loop 190: * 191: IF( K.LT.1 ) 192: $ GO TO 70 193: KSTEP = 1 194: * 195: * Determine rows and columns to be interchanged and whether 196: * a 1-by-1 or 2-by-2 pivot block will be used 197: * 198: ABSAKK = CABS1( A( K, K ) ) 199: * 200: * IMAX is the row-index of the largest off-diagonal element in 201: * column K, and COLMAX is its absolute value 202: * 203: IF( K.GT.1 ) THEN 204: IMAX = IZAMAX( K-1, A( 1, K ), 1 ) 205: COLMAX = CABS1( A( IMAX, K ) ) 206: ELSE 207: COLMAX = ZERO 208: END IF 209: * 210: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN 211: * 212: * Column K is zero or contains a NaN: set INFO and continue 213: * 214: IF( INFO.EQ.0 ) 215: $ INFO = K 216: KP = K 217: ELSE 218: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 219: * 220: * no interchange, use 1-by-1 pivot block 221: * 222: KP = K 223: ELSE 224: * 225: * JMAX is the column-index of the largest off-diagonal 226: * element in row IMAX, and ROWMAX is its absolute value 227: * 228: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA ) 229: ROWMAX = CABS1( A( IMAX, JMAX ) ) 230: IF( IMAX.GT.1 ) THEN 231: JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 ) 232: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 233: END IF 234: * 235: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 236: * 237: * no interchange, use 1-by-1 pivot block 238: * 239: KP = K 240: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 241: * 242: * interchange rows and columns K and IMAX, use 1-by-1 243: * pivot block 244: * 245: KP = IMAX 246: ELSE 247: * 248: * interchange rows and columns K-1 and IMAX, use 2-by-2 249: * pivot block 250: * 251: KP = IMAX 252: KSTEP = 2 253: END IF 254: END IF 255: * 256: KK = K - KSTEP + 1 257: IF( KP.NE.KK ) THEN 258: * 259: * Interchange rows and columns KK and KP in the leading 260: * submatrix A(1:k,1:k) 261: * 262: CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 263: CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), 264: $ LDA ) 265: T = A( KK, KK ) 266: A( KK, KK ) = A( KP, KP ) 267: A( KP, KP ) = T 268: IF( KSTEP.EQ.2 ) THEN 269: T = A( K-1, K ) 270: A( K-1, K ) = A( KP, K ) 271: A( KP, K ) = T 272: END IF 273: END IF 274: * 275: * Update the leading submatrix 276: * 277: IF( KSTEP.EQ.1 ) THEN 278: * 279: * 1-by-1 pivot block D(k): column k now holds 280: * 281: * W(k) = U(k)*D(k) 282: * 283: * where U(k) is the k-th column of U 284: * 285: * Perform a rank-1 update of A(1:k-1,1:k-1) as 286: * 287: * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' 288: * 289: R1 = CONE / A( K, K ) 290: CALL ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA ) 291: * 292: * Store U(k) in column k 293: * 294: CALL ZSCAL( K-1, R1, A( 1, K ), 1 ) 295: ELSE 296: * 297: * 2-by-2 pivot block D(k): columns k and k-1 now hold 298: * 299: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 300: * 301: * where U(k) and U(k-1) are the k-th and (k-1)-th columns 302: * of U 303: * 304: * Perform a rank-2 update of A(1:k-2,1:k-2) as 305: * 306: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )' 307: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' 308: * 309: IF( K.GT.2 ) THEN 310: * 311: D12 = A( K-1, K ) 312: D22 = A( K-1, K-1 ) / D12 313: D11 = A( K, K ) / D12 314: T = CONE / ( D11*D22-CONE ) 315: D12 = T / D12 316: * 317: DO 30 J = K - 2, 1, -1 318: WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) 319: WK = D12*( D22*A( J, K )-A( J, K-1 ) ) 320: DO 20 I = J, 1, -1 321: A( I, J ) = A( I, J ) - A( I, K )*WK - 322: $ A( I, K-1 )*WKM1 323: 20 CONTINUE 324: A( J, K ) = WK 325: A( J, K-1 ) = WKM1 326: 30 CONTINUE 327: * 328: END IF 329: * 330: END IF 331: END IF 332: * 333: * Store details of the interchanges in IPIV 334: * 335: IF( KSTEP.EQ.1 ) THEN 336: IPIV( K ) = KP 337: ELSE 338: IPIV( K ) = -KP 339: IPIV( K-1 ) = -KP 340: END IF 341: * 342: * Decrease K and return to the start of the main loop 343: * 344: K = K - KSTEP 345: GO TO 10 346: * 347: ELSE 348: * 349: * Factorize A as L*D*L' using the lower triangle of A 350: * 351: * K is the main loop index, increasing from 1 to N in steps of 352: * 1 or 2 353: * 354: K = 1 355: 40 CONTINUE 356: * 357: * If K > N, exit from loop 358: * 359: IF( K.GT.N ) 360: $ GO TO 70 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 = CABS1( A( K, K ) ) 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 + IZAMAX( N-K, A( K+1, K ), 1 ) 373: COLMAX = CABS1( A( IMAX, K ) ) 374: ELSE 375: COLMAX = ZERO 376: END IF 377: * 378: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN 379: * 380: * Column K is zero or contains a NaN: 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: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA ) 397: ROWMAX = CABS1( A( IMAX, JMAX ) ) 398: IF( IMAX.LT.N ) THEN 399: JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 ) 400: ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 401: END IF 402: * 403: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 404: * 405: * no interchange, use 1-by-1 pivot block 406: * 407: KP = K 408: ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 409: * 410: * interchange rows and columns K and IMAX, use 1-by-1 411: * pivot block 412: * 413: KP = IMAX 414: ELSE 415: * 416: * interchange rows and columns K+1 and IMAX, use 2-by-2 417: * pivot block 418: * 419: KP = IMAX 420: KSTEP = 2 421: END IF 422: END IF 423: * 424: KK = K + KSTEP - 1 425: IF( KP.NE.KK ) THEN 426: * 427: * Interchange rows and columns KK and KP in the trailing 428: * submatrix A(k:n,k:n) 429: * 430: IF( KP.LT.N ) 431: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 432: CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 433: $ LDA ) 434: T = A( KK, KK ) 435: A( KK, KK ) = A( KP, KP ) 436: A( KP, KP ) = T 437: IF( KSTEP.EQ.2 ) THEN 438: T = A( K+1, K ) 439: A( K+1, K ) = A( KP, K ) 440: A( KP, K ) = T 441: END IF 442: END IF 443: * 444: * Update the trailing submatrix 445: * 446: IF( KSTEP.EQ.1 ) THEN 447: * 448: * 1-by-1 pivot block D(k): column k now holds 449: * 450: * W(k) = L(k)*D(k) 451: * 452: * where L(k) is the k-th column of L 453: * 454: IF( K.LT.N ) THEN 455: * 456: * Perform a rank-1 update of A(k+1:n,k+1:n) as 457: * 458: * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' 459: * 460: R1 = CONE / A( K, K ) 461: CALL ZSYR( UPLO, N-K, -R1, A( K+1, K ), 1, 462: $ A( K+1, K+1 ), LDA ) 463: * 464: * Store L(k) in column K 465: * 466: CALL ZSCAL( N-K, R1, A( K+1, K ), 1 ) 467: END IF 468: ELSE 469: * 470: * 2-by-2 pivot block D(k) 471: * 472: IF( K.LT.N-1 ) THEN 473: * 474: * Perform a rank-2 update of A(k+2:n,k+2:n) as 475: * 476: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )' 477: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )' 478: * 479: * where L(k) and L(k+1) are the k-th and (k+1)-th 480: * columns of L 481: * 482: D21 = A( K+1, K ) 483: D11 = A( K+1, K+1 ) / D21 484: D22 = A( K, K ) / D21 485: T = CONE / ( D11*D22-CONE ) 486: D21 = T / D21 487: * 488: DO 60 J = K + 2, N 489: WK = D21*( D11*A( J, K )-A( J, K+1 ) ) 490: WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) 491: DO 50 I = J, N 492: A( I, J ) = A( I, J ) - A( I, K )*WK - 493: $ A( I, K+1 )*WKP1 494: 50 CONTINUE 495: A( J, K ) = WK 496: A( J, K+1 ) = WKP1 497: 60 CONTINUE 498: END IF 499: END IF 500: END IF 501: * 502: * Store details of the interchanges in IPIV 503: * 504: IF( KSTEP.EQ.1 ) THEN 505: IPIV( K ) = KP 506: ELSE 507: IPIV( K ) = -KP 508: IPIV( K+1 ) = -KP 509: END IF 510: * 511: * Increase K and return to the start of the main loop 512: * 513: K = K + KSTEP 514: GO TO 40 515: * 516: END IF 517: * 518: 70 CONTINUE 519: RETURN 520: * 521: * End of ZSYTF2 522: * 523: END