Annotation of rpl/lapack/lapack/zhetf2_rk.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 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 ZHETF2_RK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2_rk.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2_rk.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2_rk.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, 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, * ), E ( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *> ZHETF2_RK computes the factorization of a complex Hermitian matrix A
! 38: *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
! 39: *>
! 40: *> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
! 41: *>
! 42: *> where U (or L) is unit upper (or lower) triangular matrix,
! 43: *> U**H (or L**H) is the conjugate of U (or L), P is a permutation
! 44: *> matrix, P**T is the transpose of P, and D is Hermitian and block
! 45: *> 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: *> For more information see Further Details section.
! 49: *> \endverbatim
! 50: *
! 51: * Arguments:
! 52: * ==========
! 53: *
! 54: *> \param[in] UPLO
! 55: *> \verbatim
! 56: *> UPLO is CHARACTER*1
! 57: *> Specifies whether the upper or lower triangular part of the
! 58: *> Hermitian matrix A is stored:
! 59: *> = 'U': Upper triangular
! 60: *> = 'L': Lower triangular
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] N
! 64: *> \verbatim
! 65: *> N is INTEGER
! 66: *> The order of the matrix A. N >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in,out] A
! 70: *> \verbatim
! 71: *> A is COMPLEX*16 array, dimension (LDA,N)
! 72: *> On entry, the Hermitian matrix A.
! 73: *> If UPLO = 'U': the leading N-by-N upper triangular part
! 74: *> of A contains the upper triangular part of the matrix A,
! 75: *> and the strictly lower triangular part of A is not
! 76: *> referenced.
! 77: *>
! 78: *> If UPLO = 'L': the leading N-by-N lower triangular part
! 79: *> of A contains the lower triangular part of the matrix A,
! 80: *> and the strictly upper triangular part of A is not
! 81: *> referenced.
! 82: *>
! 83: *> On exit, contains:
! 84: *> a) ONLY diagonal elements of the Hermitian block diagonal
! 85: *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
! 86: *> (superdiagonal (or subdiagonal) elements of D
! 87: *> are stored on exit in array E), and
! 88: *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
! 89: *> If UPLO = 'L': factor L in the subdiagonal part of A.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] LDA
! 93: *> \verbatim
! 94: *> LDA is INTEGER
! 95: *> The leading dimension of the array A. LDA >= max(1,N).
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[out] E
! 99: *> \verbatim
! 100: *> E is COMPLEX*16 array, dimension (N)
! 101: *> On exit, contains the superdiagonal (or subdiagonal)
! 102: *> elements of the Hermitian block diagonal matrix D
! 103: *> with 1-by-1 or 2-by-2 diagonal blocks, where
! 104: *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
! 105: *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
! 106: *>
! 107: *> NOTE: For 1-by-1 diagonal block D(k), where
! 108: *> 1 <= k <= N, the element E(k) is set to 0 in both
! 109: *> UPLO = 'U' or UPLO = 'L' cases.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[out] IPIV
! 113: *> \verbatim
! 114: *> IPIV is INTEGER array, dimension (N)
! 115: *> IPIV describes the permutation matrix P in the factorization
! 116: *> of matrix A as follows. The absolute value of IPIV(k)
! 117: *> represents the index of row and column that were
! 118: *> interchanged with the k-th row and column. The value of UPLO
! 119: *> describes the order in which the interchanges were applied.
! 120: *> Also, the sign of IPIV represents the block structure of
! 121: *> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
! 122: *> diagonal blocks which correspond to 1 or 2 interchanges
! 123: *> at each factorization step. For more info see Further
! 124: *> Details section.
! 125: *>
! 126: *> If UPLO = 'U',
! 127: *> ( in factorization order, k decreases from N to 1 ):
! 128: *> a) A single positive entry IPIV(k) > 0 means:
! 129: *> D(k,k) is a 1-by-1 diagonal block.
! 130: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
! 131: *> interchanged in the matrix A(1:N,1:N);
! 132: *> If IPIV(k) = k, no interchange occurred.
! 133: *>
! 134: *> b) A pair of consecutive negative entries
! 135: *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
! 136: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
! 137: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
! 138: *> 1) If -IPIV(k) != k, rows and columns
! 139: *> k and -IPIV(k) were interchanged
! 140: *> in the matrix A(1:N,1:N).
! 141: *> If -IPIV(k) = k, no interchange occurred.
! 142: *> 2) If -IPIV(k-1) != k-1, rows and columns
! 143: *> k-1 and -IPIV(k-1) were interchanged
! 144: *> in the matrix A(1:N,1:N).
! 145: *> If -IPIV(k-1) = k-1, no interchange occurred.
! 146: *>
! 147: *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
! 148: *>
! 149: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
! 150: *>
! 151: *> If UPLO = 'L',
! 152: *> ( in factorization order, k increases from 1 to N ):
! 153: *> a) A single positive entry IPIV(k) > 0 means:
! 154: *> D(k,k) is a 1-by-1 diagonal block.
! 155: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
! 156: *> interchanged in the matrix A(1:N,1:N).
! 157: *> If IPIV(k) = k, no interchange occurred.
! 158: *>
! 159: *> b) A pair of consecutive negative entries
! 160: *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
! 161: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
! 162: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
! 163: *> 1) If -IPIV(k) != k, rows and columns
! 164: *> k and -IPIV(k) were interchanged
! 165: *> in the matrix A(1:N,1:N).
! 166: *> If -IPIV(k) = k, no interchange occurred.
! 167: *> 2) If -IPIV(k+1) != k+1, rows and columns
! 168: *> k-1 and -IPIV(k-1) were interchanged
! 169: *> in the matrix A(1:N,1:N).
! 170: *> If -IPIV(k+1) = k+1, no interchange occurred.
! 171: *>
! 172: *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
! 173: *>
! 174: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
! 175: *> \endverbatim
! 176: *>
! 177: *> \param[out] INFO
! 178: *> \verbatim
! 179: *> INFO is INTEGER
! 180: *> = 0: successful exit
! 181: *>
! 182: *> < 0: If INFO = -k, the k-th argument had an illegal value
! 183: *>
! 184: *> > 0: If INFO = k, the matrix A is singular, because:
! 185: *> If UPLO = 'U': column k in the upper
! 186: *> triangular part of A contains all zeros.
! 187: *> If UPLO = 'L': column k in the lower
! 188: *> triangular part of A contains all zeros.
! 189: *>
! 190: *> Therefore D(k,k) is exactly zero, and superdiagonal
! 191: *> elements of column k of U (or subdiagonal elements of
! 192: *> column k of L ) are all zeros. The factorization has
! 193: *> been completed, but the block diagonal matrix D is
! 194: *> exactly singular, and division by zero will occur if
! 195: *> it is used to solve a system of equations.
! 196: *>
! 197: *> NOTE: INFO only stores the first occurrence of
! 198: *> a singularity, any subsequent occurrence of singularity
! 199: *> is not stored in INFO even though the factorization
! 200: *> always completes.
! 201: *> \endverbatim
! 202: *
! 203: * Authors:
! 204: * ========
! 205: *
! 206: *> \author Univ. of Tennessee
! 207: *> \author Univ. of California Berkeley
! 208: *> \author Univ. of Colorado Denver
! 209: *> \author NAG Ltd.
! 210: *
! 211: *> \date December 2016
! 212: *
! 213: *> \ingroup complex16HEcomputational
! 214: *
! 215: *> \par Further Details:
! 216: * =====================
! 217: *>
! 218: *> \verbatim
! 219: *> TODO: put further details
! 220: *> \endverbatim
! 221: *
! 222: *> \par Contributors:
! 223: * ==================
! 224: *>
! 225: *> \verbatim
! 226: *>
! 227: *> December 2016, Igor Kozachenko,
! 228: *> Computer Science Division,
! 229: *> University of California, Berkeley
! 230: *>
! 231: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
! 232: *> School of Mathematics,
! 233: *> University of Manchester
! 234: *>
! 235: *> 01-01-96 - Based on modifications by
! 236: *> J. Lewis, Boeing Computer Services Company
! 237: *> A. Petitet, Computer Science Dept.,
! 238: *> Univ. of Tenn., Knoxville abd , USA
! 239: *> \endverbatim
! 240: *
! 241: * =====================================================================
! 242: SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
! 243: *
! 244: * -- LAPACK computational routine (version 3.7.0) --
! 245: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 246: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 247: * December 2016
! 248: *
! 249: * .. Scalar Arguments ..
! 250: CHARACTER UPLO
! 251: INTEGER INFO, LDA, N
! 252: * ..
! 253: * .. Array Arguments ..
! 254: INTEGER IPIV( * )
! 255: COMPLEX*16 A( LDA, * ), E( * )
! 256: * ..
! 257: *
! 258: * ======================================================================
! 259: *
! 260: * .. Parameters ..
! 261: DOUBLE PRECISION ZERO, ONE
! 262: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 263: DOUBLE PRECISION EIGHT, SEVTEN
! 264: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
! 265: COMPLEX*16 CZERO
! 266: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
! 267: * ..
! 268: * .. Local Scalars ..
! 269: LOGICAL DONE, UPPER
! 270: INTEGER I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP,
! 271: $ P
! 272: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, DTEMP,
! 273: $ ROWMAX, TT, SFMIN
! 274: COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, Z
! 275: * ..
! 276: * .. External Functions ..
! 277: *
! 278: LOGICAL LSAME
! 279: INTEGER IZAMAX
! 280: DOUBLE PRECISION DLAMCH, DLAPY2
! 281: EXTERNAL LSAME, IZAMAX, DLAMCH, DLAPY2
! 282: * ..
! 283: * .. External Subroutines ..
! 284: EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
! 285: * ..
! 286: * .. Intrinsic Functions ..
! 287: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
! 288: * ..
! 289: * .. Statement Functions ..
! 290: DOUBLE PRECISION CABS1
! 291: * ..
! 292: * .. Statement Function definitions ..
! 293: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
! 294: * ..
! 295: * .. Executable Statements ..
! 296: *
! 297: * Test the input parameters.
! 298: *
! 299: INFO = 0
! 300: UPPER = LSAME( UPLO, 'U' )
! 301: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 302: INFO = -1
! 303: ELSE IF( N.LT.0 ) THEN
! 304: INFO = -2
! 305: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 306: INFO = -4
! 307: END IF
! 308: IF( INFO.NE.0 ) THEN
! 309: CALL XERBLA( 'ZHETF2_RK', -INFO )
! 310: RETURN
! 311: END IF
! 312: *
! 313: * Initialize ALPHA for use in choosing pivot block size.
! 314: *
! 315: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
! 316: *
! 317: * Compute machine safe minimum
! 318: *
! 319: SFMIN = DLAMCH( 'S' )
! 320: *
! 321: IF( UPPER ) THEN
! 322: *
! 323: * Factorize A as U*D*U**H using the upper triangle of A
! 324: *
! 325: * Initilize the first entry of array E, where superdiagonal
! 326: * elements of D are stored
! 327: *
! 328: E( 1 ) = CZERO
! 329: *
! 330: * K is the main loop index, decreasing from N to 1 in steps of
! 331: * 1 or 2
! 332: *
! 333: K = N
! 334: 10 CONTINUE
! 335: *
! 336: * If K < 1, exit from loop
! 337: *
! 338: IF( K.LT.1 )
! 339: $ GO TO 34
! 340: KSTEP = 1
! 341: P = K
! 342: *
! 343: * Determine rows and columns to be interchanged and whether
! 344: * a 1-by-1 or 2-by-2 pivot block will be used
! 345: *
! 346: ABSAKK = ABS( DBLE( A( K, K ) ) )
! 347: *
! 348: * IMAX is the row-index of the largest off-diagonal element in
! 349: * column K, and COLMAX is its absolute value.
! 350: * Determine both COLMAX and IMAX.
! 351: *
! 352: IF( K.GT.1 ) THEN
! 353: IMAX = IZAMAX( K-1, A( 1, K ), 1 )
! 354: COLMAX = CABS1( A( IMAX, K ) )
! 355: ELSE
! 356: COLMAX = ZERO
! 357: END IF
! 358: *
! 359: IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN
! 360: *
! 361: * Column K is zero or underflow: set INFO and continue
! 362: *
! 363: IF( INFO.EQ.0 )
! 364: $ INFO = K
! 365: KP = K
! 366: A( K, K ) = DBLE( A( K, K ) )
! 367: *
! 368: * Set E( K ) to zero
! 369: *
! 370: IF( K.GT.1 )
! 371: $ E( K ) = CZERO
! 372: *
! 373: ELSE
! 374: *
! 375: * ============================================================
! 376: *
! 377: * BEGIN pivot search
! 378: *
! 379: * Case(1)
! 380: * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
! 381: * (used to handle NaN and Inf)
! 382: *
! 383: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
! 384: *
! 385: * no interchange, use 1-by-1 pivot block
! 386: *
! 387: KP = K
! 388: *
! 389: ELSE
! 390: *
! 391: DONE = .FALSE.
! 392: *
! 393: * Loop until pivot found
! 394: *
! 395: 12 CONTINUE
! 396: *
! 397: * BEGIN pivot search loop body
! 398: *
! 399: *
! 400: * JMAX is the column-index of the largest off-diagonal
! 401: * element in row IMAX, and ROWMAX is its absolute value.
! 402: * Determine both ROWMAX and JMAX.
! 403: *
! 404: IF( IMAX.NE.K ) THEN
! 405: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ),
! 406: $ LDA )
! 407: ROWMAX = CABS1( A( IMAX, JMAX ) )
! 408: ELSE
! 409: ROWMAX = ZERO
! 410: END IF
! 411: *
! 412: IF( IMAX.GT.1 ) THEN
! 413: ITEMP = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
! 414: DTEMP = CABS1( A( ITEMP, IMAX ) )
! 415: IF( DTEMP.GT.ROWMAX ) THEN
! 416: ROWMAX = DTEMP
! 417: JMAX = ITEMP
! 418: END IF
! 419: END IF
! 420: *
! 421: * Case(2)
! 422: * Equivalent to testing for
! 423: * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
! 424: * (used to handle NaN and Inf)
! 425: *
! 426: IF( .NOT.( ABS( DBLE( A( IMAX, IMAX ) ) )
! 427: $ .LT.ALPHA*ROWMAX ) ) THEN
! 428: *
! 429: * interchange rows and columns K and IMAX,
! 430: * use 1-by-1 pivot block
! 431: *
! 432: KP = IMAX
! 433: DONE = .TRUE.
! 434: *
! 435: * Case(3)
! 436: * Equivalent to testing for ROWMAX.EQ.COLMAX,
! 437: * (used to handle NaN and Inf)
! 438: *
! 439: ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
! 440: $ THEN
! 441: *
! 442: * interchange rows and columns K-1 and IMAX,
! 443: * use 2-by-2 pivot block
! 444: *
! 445: KP = IMAX
! 446: KSTEP = 2
! 447: DONE = .TRUE.
! 448: *
! 449: * Case(4)
! 450: ELSE
! 451: *
! 452: * Pivot not found: set params and repeat
! 453: *
! 454: P = IMAX
! 455: COLMAX = ROWMAX
! 456: IMAX = JMAX
! 457: END IF
! 458: *
! 459: * END pivot search loop body
! 460: *
! 461: IF( .NOT.DONE ) GOTO 12
! 462: *
! 463: END IF
! 464: *
! 465: * END pivot search
! 466: *
! 467: * ============================================================
! 468: *
! 469: * KK is the column of A where pivoting step stopped
! 470: *
! 471: KK = K - KSTEP + 1
! 472: *
! 473: * For only a 2x2 pivot, interchange rows and columns K and P
! 474: * in the leading submatrix A(1:k,1:k)
! 475: *
! 476: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
! 477: * (1) Swap columnar parts
! 478: IF( P.GT.1 )
! 479: $ CALL ZSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
! 480: * (2) Swap and conjugate middle parts
! 481: DO 14 J = P + 1, K - 1
! 482: T = DCONJG( A( J, K ) )
! 483: A( J, K ) = DCONJG( A( P, J ) )
! 484: A( P, J ) = T
! 485: 14 CONTINUE
! 486: * (3) Swap and conjugate corner elements at row-col interserction
! 487: A( P, K ) = DCONJG( A( P, K ) )
! 488: * (4) Swap diagonal elements at row-col intersection
! 489: R1 = DBLE( A( K, K ) )
! 490: A( K, K ) = DBLE( A( P, P ) )
! 491: A( P, P ) = R1
! 492: *
! 493: * Convert upper triangle of A into U form by applying
! 494: * the interchanges in columns k+1:N.
! 495: *
! 496: IF( K.LT.N )
! 497: $ CALL ZSWAP( N-K, A( K, K+1 ), LDA, A( P, K+1 ), LDA )
! 498: *
! 499: END IF
! 500: *
! 501: * For both 1x1 and 2x2 pivots, interchange rows and
! 502: * columns KK and KP in the leading submatrix A(1:k,1:k)
! 503: *
! 504: IF( KP.NE.KK ) THEN
! 505: * (1) Swap columnar parts
! 506: IF( KP.GT.1 )
! 507: $ CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
! 508: * (2) Swap and conjugate middle parts
! 509: DO 15 J = KP + 1, KK - 1
! 510: T = DCONJG( A( J, KK ) )
! 511: A( J, KK ) = DCONJG( A( KP, J ) )
! 512: A( KP, J ) = T
! 513: 15 CONTINUE
! 514: * (3) Swap and conjugate corner elements at row-col interserction
! 515: A( KP, KK ) = DCONJG( A( KP, KK ) )
! 516: * (4) Swap diagonal elements at row-col intersection
! 517: R1 = DBLE( A( KK, KK ) )
! 518: A( KK, KK ) = DBLE( A( KP, KP ) )
! 519: A( KP, KP ) = R1
! 520: *
! 521: IF( KSTEP.EQ.2 ) THEN
! 522: * (*) Make sure that diagonal element of pivot is real
! 523: A( K, K ) = DBLE( A( K, K ) )
! 524: * (5) Swap row elements
! 525: T = A( K-1, K )
! 526: A( K-1, K ) = A( KP, K )
! 527: A( KP, K ) = T
! 528: END IF
! 529: *
! 530: * Convert upper triangle of A into U form by applying
! 531: * the interchanges in columns k+1:N.
! 532: *
! 533: IF( K.LT.N )
! 534: $ CALL ZSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
! 535: $ LDA )
! 536: *
! 537: ELSE
! 538: * (*) Make sure that diagonal element of pivot is real
! 539: A( K, K ) = DBLE( A( K, K ) )
! 540: IF( KSTEP.EQ.2 )
! 541: $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
! 542: END IF
! 543: *
! 544: * Update the leading submatrix
! 545: *
! 546: IF( KSTEP.EQ.1 ) THEN
! 547: *
! 548: * 1-by-1 pivot block D(k): column k now holds
! 549: *
! 550: * W(k) = U(k)*D(k)
! 551: *
! 552: * where U(k) is the k-th column of U
! 553: *
! 554: IF( K.GT.1 ) THEN
! 555: *
! 556: * Perform a rank-1 update of A(1:k-1,1:k-1) and
! 557: * store U(k) in column k
! 558: *
! 559: IF( ABS( DBLE( A( K, K ) ) ).GE.SFMIN ) THEN
! 560: *
! 561: * Perform a rank-1 update of A(1:k-1,1:k-1) as
! 562: * A := A - U(k)*D(k)*U(k)**T
! 563: * = A - W(k)*1/D(k)*W(k)**T
! 564: *
! 565: D11 = ONE / DBLE( A( K, K ) )
! 566: CALL ZHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
! 567: *
! 568: * Store U(k) in column k
! 569: *
! 570: CALL ZDSCAL( K-1, D11, A( 1, K ), 1 )
! 571: ELSE
! 572: *
! 573: * Store L(k) in column K
! 574: *
! 575: D11 = DBLE( A( K, K ) )
! 576: DO 16 II = 1, K - 1
! 577: A( II, K ) = A( II, K ) / D11
! 578: 16 CONTINUE
! 579: *
! 580: * Perform a rank-1 update of A(k+1:n,k+1:n) as
! 581: * A := A - U(k)*D(k)*U(k)**T
! 582: * = A - W(k)*(1/D(k))*W(k)**T
! 583: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
! 584: *
! 585: CALL ZHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
! 586: END IF
! 587: *
! 588: * Store the superdiagonal element of D in array E
! 589: *
! 590: E( K ) = CZERO
! 591: *
! 592: END IF
! 593: *
! 594: ELSE
! 595: *
! 596: * 2-by-2 pivot block D(k): columns k and k-1 now hold
! 597: *
! 598: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
! 599: *
! 600: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
! 601: * of U
! 602: *
! 603: * Perform a rank-2 update of A(1:k-2,1:k-2) as
! 604: *
! 605: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
! 606: * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
! 607: *
! 608: * and store L(k) and L(k+1) in columns k and k+1
! 609: *
! 610: IF( K.GT.2 ) THEN
! 611: * D = |A12|
! 612: D = DLAPY2( DBLE( A( K-1, K ) ),
! 613: $ DIMAG( A( K-1, K ) ) )
! 614: D11 = A( K, K ) / D
! 615: D22 = A( K-1, K-1 ) / D
! 616: D12 = A( K-1, K ) / D
! 617: TT = ONE / ( D11*D22-ONE )
! 618: *
! 619: DO 30 J = K - 2, 1, -1
! 620: *
! 621: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
! 622: *
! 623: WKM1 = TT*( D11*A( J, K-1 )-DCONJG( D12 )*
! 624: $ A( J, K ) )
! 625: WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) )
! 626: *
! 627: * Perform a rank-2 update of A(1:k-2,1:k-2)
! 628: *
! 629: DO 20 I = J, 1, -1
! 630: A( I, J ) = A( I, J ) -
! 631: $ ( A( I, K ) / D )*DCONJG( WK ) -
! 632: $ ( A( I, K-1 ) / D )*DCONJG( WKM1 )
! 633: 20 CONTINUE
! 634: *
! 635: * Store U(k) and U(k-1) in cols k and k-1 for row J
! 636: *
! 637: A( J, K ) = WK / D
! 638: A( J, K-1 ) = WKM1 / D
! 639: * (*) Make sure that diagonal element of pivot is real
! 640: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO )
! 641: *
! 642: 30 CONTINUE
! 643: *
! 644: END IF
! 645: *
! 646: * Copy superdiagonal elements of D(K) to E(K) and
! 647: * ZERO out superdiagonal entry of A
! 648: *
! 649: E( K ) = A( K-1, K )
! 650: E( K-1 ) = CZERO
! 651: A( K-1, K ) = CZERO
! 652: *
! 653: END IF
! 654: *
! 655: * End column K is nonsingular
! 656: *
! 657: END IF
! 658: *
! 659: * Store details of the interchanges in IPIV
! 660: *
! 661: IF( KSTEP.EQ.1 ) THEN
! 662: IPIV( K ) = KP
! 663: ELSE
! 664: IPIV( K ) = -P
! 665: IPIV( K-1 ) = -KP
! 666: END IF
! 667: *
! 668: * Decrease K and return to the start of the main loop
! 669: *
! 670: K = K - KSTEP
! 671: GO TO 10
! 672: *
! 673: 34 CONTINUE
! 674: *
! 675: ELSE
! 676: *
! 677: * Factorize A as L*D*L**H using the lower triangle of A
! 678: *
! 679: * Initilize the unused last entry of the subdiagonal array E.
! 680: *
! 681: E( N ) = CZERO
! 682: *
! 683: * K is the main loop index, increasing from 1 to N in steps of
! 684: * 1 or 2
! 685: *
! 686: K = 1
! 687: 40 CONTINUE
! 688: *
! 689: * If K > N, exit from loop
! 690: *
! 691: IF( K.GT.N )
! 692: $ GO TO 64
! 693: KSTEP = 1
! 694: P = K
! 695: *
! 696: * Determine rows and columns to be interchanged and whether
! 697: * a 1-by-1 or 2-by-2 pivot block will be used
! 698: *
! 699: ABSAKK = ABS( DBLE( A( K, K ) ) )
! 700: *
! 701: * IMAX is the row-index of the largest off-diagonal element in
! 702: * column K, and COLMAX is its absolute value.
! 703: * Determine both COLMAX and IMAX.
! 704: *
! 705: IF( K.LT.N ) THEN
! 706: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
! 707: COLMAX = CABS1( A( IMAX, K ) )
! 708: ELSE
! 709: COLMAX = ZERO
! 710: END IF
! 711: *
! 712: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
! 713: *
! 714: * Column K is zero or underflow: set INFO and continue
! 715: *
! 716: IF( INFO.EQ.0 )
! 717: $ INFO = K
! 718: KP = K
! 719: A( K, K ) = DBLE( A( K, K ) )
! 720: *
! 721: * Set E( K ) to zero
! 722: *
! 723: IF( K.LT.N )
! 724: $ E( K ) = CZERO
! 725: *
! 726: ELSE
! 727: *
! 728: * ============================================================
! 729: *
! 730: * BEGIN pivot search
! 731: *
! 732: * Case(1)
! 733: * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
! 734: * (used to handle NaN and Inf)
! 735: *
! 736: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
! 737: *
! 738: * no interchange, use 1-by-1 pivot block
! 739: *
! 740: KP = K
! 741: *
! 742: ELSE
! 743: *
! 744: DONE = .FALSE.
! 745: *
! 746: * Loop until pivot found
! 747: *
! 748: 42 CONTINUE
! 749: *
! 750: * BEGIN pivot search loop body
! 751: *
! 752: *
! 753: * JMAX is the column-index of the largest off-diagonal
! 754: * element in row IMAX, and ROWMAX is its absolute value.
! 755: * Determine both ROWMAX and JMAX.
! 756: *
! 757: IF( IMAX.NE.K ) THEN
! 758: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
! 759: ROWMAX = CABS1( A( IMAX, JMAX ) )
! 760: ELSE
! 761: ROWMAX = ZERO
! 762: END IF
! 763: *
! 764: IF( IMAX.LT.N ) THEN
! 765: ITEMP = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ),
! 766: $ 1 )
! 767: DTEMP = CABS1( A( ITEMP, IMAX ) )
! 768: IF( DTEMP.GT.ROWMAX ) THEN
! 769: ROWMAX = DTEMP
! 770: JMAX = ITEMP
! 771: END IF
! 772: END IF
! 773: *
! 774: * Case(2)
! 775: * Equivalent to testing for
! 776: * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
! 777: * (used to handle NaN and Inf)
! 778: *
! 779: IF( .NOT.( ABS( DBLE( A( IMAX, IMAX ) ) )
! 780: $ .LT.ALPHA*ROWMAX ) ) THEN
! 781: *
! 782: * interchange rows and columns K and IMAX,
! 783: * use 1-by-1 pivot block
! 784: *
! 785: KP = IMAX
! 786: DONE = .TRUE.
! 787: *
! 788: * Case(3)
! 789: * Equivalent to testing for ROWMAX.EQ.COLMAX,
! 790: * (used to handle NaN and Inf)
! 791: *
! 792: ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
! 793: $ THEN
! 794: *
! 795: * interchange rows and columns K+1 and IMAX,
! 796: * use 2-by-2 pivot block
! 797: *
! 798: KP = IMAX
! 799: KSTEP = 2
! 800: DONE = .TRUE.
! 801: *
! 802: * Case(4)
! 803: ELSE
! 804: *
! 805: * Pivot not found: set params and repeat
! 806: *
! 807: P = IMAX
! 808: COLMAX = ROWMAX
! 809: IMAX = JMAX
! 810: END IF
! 811: *
! 812: *
! 813: * END pivot search loop body
! 814: *
! 815: IF( .NOT.DONE ) GOTO 42
! 816: *
! 817: END IF
! 818: *
! 819: * END pivot search
! 820: *
! 821: * ============================================================
! 822: *
! 823: * KK is the column of A where pivoting step stopped
! 824: *
! 825: KK = K + KSTEP - 1
! 826: *
! 827: * For only a 2x2 pivot, interchange rows and columns K and P
! 828: * in the trailing submatrix A(k:n,k:n)
! 829: *
! 830: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
! 831: * (1) Swap columnar parts
! 832: IF( P.LT.N )
! 833: $ CALL ZSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
! 834: * (2) Swap and conjugate middle parts
! 835: DO 44 J = K + 1, P - 1
! 836: T = DCONJG( A( J, K ) )
! 837: A( J, K ) = DCONJG( A( P, J ) )
! 838: A( P, J ) = T
! 839: 44 CONTINUE
! 840: * (3) Swap and conjugate corner elements at row-col interserction
! 841: A( P, K ) = DCONJG( A( P, K ) )
! 842: * (4) Swap diagonal elements at row-col intersection
! 843: R1 = DBLE( A( K, K ) )
! 844: A( K, K ) = DBLE( A( P, P ) )
! 845: A( P, P ) = R1
! 846: *
! 847: * Convert lower triangle of A into L form by applying
! 848: * the interchanges in columns 1:k-1.
! 849: *
! 850: IF ( K.GT.1 )
! 851: $ CALL ZSWAP( K-1, A( K, 1 ), LDA, A( P, 1 ), LDA )
! 852: *
! 853: END IF
! 854: *
! 855: * For both 1x1 and 2x2 pivots, interchange rows and
! 856: * columns KK and KP in the trailing submatrix A(k:n,k:n)
! 857: *
! 858: IF( KP.NE.KK ) THEN
! 859: * (1) Swap columnar parts
! 860: IF( KP.LT.N )
! 861: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
! 862: * (2) Swap and conjugate middle parts
! 863: DO 45 J = KK + 1, KP - 1
! 864: T = DCONJG( A( J, KK ) )
! 865: A( J, KK ) = DCONJG( A( KP, J ) )
! 866: A( KP, J ) = T
! 867: 45 CONTINUE
! 868: * (3) Swap and conjugate corner elements at row-col interserction
! 869: A( KP, KK ) = DCONJG( A( KP, KK ) )
! 870: * (4) Swap diagonal elements at row-col intersection
! 871: R1 = DBLE( A( KK, KK ) )
! 872: A( KK, KK ) = DBLE( A( KP, KP ) )
! 873: A( KP, KP ) = R1
! 874: *
! 875: IF( KSTEP.EQ.2 ) THEN
! 876: * (*) Make sure that diagonal element of pivot is real
! 877: A( K, K ) = DBLE( A( K, K ) )
! 878: * (5) Swap row elements
! 879: T = A( K+1, K )
! 880: A( K+1, K ) = A( KP, K )
! 881: A( KP, K ) = T
! 882: END IF
! 883: *
! 884: * Convert lower triangle of A into L form by applying
! 885: * the interchanges in columns 1:k-1.
! 886: *
! 887: IF ( K.GT.1 )
! 888: $ CALL ZSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
! 889: *
! 890: ELSE
! 891: * (*) Make sure that diagonal element of pivot is real
! 892: A( K, K ) = DBLE( A( K, K ) )
! 893: IF( KSTEP.EQ.2 )
! 894: $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
! 895: END IF
! 896: *
! 897: * Update the trailing submatrix
! 898: *
! 899: IF( KSTEP.EQ.1 ) THEN
! 900: *
! 901: * 1-by-1 pivot block D(k): column k of A now holds
! 902: *
! 903: * W(k) = L(k)*D(k),
! 904: *
! 905: * where L(k) is the k-th column of L
! 906: *
! 907: IF( K.LT.N ) THEN
! 908: *
! 909: * Perform a rank-1 update of A(k+1:n,k+1:n) and
! 910: * store L(k) in column k
! 911: *
! 912: * Handle division by a small number
! 913: *
! 914: IF( ABS( DBLE( A( K, K ) ) ).GE.SFMIN ) THEN
! 915: *
! 916: * Perform a rank-1 update of A(k+1:n,k+1:n) as
! 917: * A := A - L(k)*D(k)*L(k)**T
! 918: * = A - W(k)*(1/D(k))*W(k)**T
! 919: *
! 920: D11 = ONE / DBLE( A( K, K ) )
! 921: CALL ZHER( UPLO, N-K, -D11, A( K+1, K ), 1,
! 922: $ A( K+1, K+1 ), LDA )
! 923: *
! 924: * Store L(k) in column k
! 925: *
! 926: CALL ZDSCAL( N-K, D11, A( K+1, K ), 1 )
! 927: ELSE
! 928: *
! 929: * Store L(k) in column k
! 930: *
! 931: D11 = DBLE( A( K, K ) )
! 932: DO 46 II = K + 1, N
! 933: A( II, K ) = A( II, K ) / D11
! 934: 46 CONTINUE
! 935: *
! 936: * Perform a rank-1 update of A(k+1:n,k+1:n) as
! 937: * A := A - L(k)*D(k)*L(k)**T
! 938: * = A - W(k)*(1/D(k))*W(k)**T
! 939: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
! 940: *
! 941: CALL ZHER( UPLO, N-K, -D11, A( K+1, K ), 1,
! 942: $ A( K+1, K+1 ), LDA )
! 943: END IF
! 944: *
! 945: * Store the subdiagonal element of D in array E
! 946: *
! 947: E( K ) = CZERO
! 948: *
! 949: END IF
! 950: *
! 951: ELSE
! 952: *
! 953: * 2-by-2 pivot block D(k): columns k and k+1 now hold
! 954: *
! 955: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
! 956: *
! 957: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
! 958: * of L
! 959: *
! 960: *
! 961: * Perform a rank-2 update of A(k+2:n,k+2:n) as
! 962: *
! 963: * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
! 964: * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
! 965: *
! 966: * and store L(k) and L(k+1) in columns k and k+1
! 967: *
! 968: IF( K.LT.N-1 ) THEN
! 969: * D = |A21|
! 970: D = DLAPY2( DBLE( A( K+1, K ) ),
! 971: $ DIMAG( A( K+1, K ) ) )
! 972: D11 = DBLE( A( K+1, K+1 ) ) / D
! 973: D22 = DBLE( A( K, K ) ) / D
! 974: D21 = A( K+1, K ) / D
! 975: TT = ONE / ( D11*D22-ONE )
! 976: *
! 977: DO 60 J = K + 2, N
! 978: *
! 979: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
! 980: *
! 981: WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) )
! 982: WKP1 = TT*( D22*A( J, K+1 )-DCONJG( D21 )*
! 983: $ A( J, K ) )
! 984: *
! 985: * Perform a rank-2 update of A(k+2:n,k+2:n)
! 986: *
! 987: DO 50 I = J, N
! 988: A( I, J ) = A( I, J ) -
! 989: $ ( A( I, K ) / D )*DCONJG( WK ) -
! 990: $ ( A( I, K+1 ) / D )*DCONJG( WKP1 )
! 991: 50 CONTINUE
! 992: *
! 993: * Store L(k) and L(k+1) in cols k and k+1 for row J
! 994: *
! 995: A( J, K ) = WK / D
! 996: A( J, K+1 ) = WKP1 / D
! 997: * (*) Make sure that diagonal element of pivot is real
! 998: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO )
! 999: *
! 1000: 60 CONTINUE
! 1001: *
! 1002: END IF
! 1003: *
! 1004: * Copy subdiagonal elements of D(K) to E(K) and
! 1005: * ZERO out subdiagonal entry of A
! 1006: *
! 1007: E( K ) = A( K+1, K )
! 1008: E( K+1 ) = CZERO
! 1009: A( K+1, K ) = CZERO
! 1010: *
! 1011: END IF
! 1012: *
! 1013: * End column K is nonsingular
! 1014: *
! 1015: END IF
! 1016: *
! 1017: * Store details of the interchanges in IPIV
! 1018: *
! 1019: IF( KSTEP.EQ.1 ) THEN
! 1020: IPIV( K ) = KP
! 1021: ELSE
! 1022: IPIV( K ) = -P
! 1023: IPIV( K+1 ) = -KP
! 1024: END IF
! 1025: *
! 1026: * Increase K and return to the start of the main loop
! 1027: *
! 1028: K = K + KSTEP
! 1029: GO TO 40
! 1030: *
! 1031: 64 CONTINUE
! 1032: *
! 1033: END IF
! 1034: *
! 1035: RETURN
! 1036: *
! 1037: * End of ZHETF2_RK
! 1038: *
! 1039: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetf2_rk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack