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