Annotation of rpl/lapack/lapack/zhetrf_rk.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETRF_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHETRF_RK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_rk.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_rk.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_rk.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
! 22: * INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER UPLO
! 26: * INTEGER INFO, LDA, LWORK, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IPIV( * )
! 30: * COMPLEX*16 A( LDA, * ), E ( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *> ZHETRF_RK computes the factorization of a complex Hermitian matrix A
! 39: *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
! 40: *>
! 41: *> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
! 42: *>
! 43: *> where U (or L) is unit upper (or lower) triangular matrix,
! 44: *> U**H (or L**H) is the conjugate of U (or L), P is a permutation
! 45: *> matrix, P**T is the transpose of P, and D is Hermitian and block
! 46: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
! 47: *>
! 48: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
! 49: *> For more information see Further Details section.
! 50: *> \endverbatim
! 51: *
! 52: * Arguments:
! 53: * ==========
! 54: *
! 55: *> \param[in] UPLO
! 56: *> \verbatim
! 57: *> UPLO is CHARACTER*1
! 58: *> Specifies whether the upper or lower triangular part of the
! 59: *> Hermitian matrix A is stored:
! 60: *> = 'U': Upper triangular
! 61: *> = 'L': Lower triangular
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in] N
! 65: *> \verbatim
! 66: *> N is INTEGER
! 67: *> The order of the matrix A. N >= 0.
! 68: *> \endverbatim
! 69: *>
! 70: *> \param[in,out] A
! 71: *> \verbatim
! 72: *> A is COMPLEX*16 array, dimension (LDA,N)
! 73: *> On entry, the Hermitian matrix A.
! 74: *> If UPLO = 'U': the leading N-by-N upper triangular part
! 75: *> of A contains the upper triangular part of the matrix A,
! 76: *> and the strictly lower triangular part of A is not
! 77: *> referenced.
! 78: *>
! 79: *> If UPLO = 'L': the leading N-by-N lower triangular part
! 80: *> of A contains the lower triangular part of the matrix A,
! 81: *> and the strictly upper triangular part of A is not
! 82: *> referenced.
! 83: *>
! 84: *> On exit, contains:
! 85: *> a) ONLY diagonal elements of the Hermitian block diagonal
! 86: *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
! 87: *> (superdiagonal (or subdiagonal) elements of D
! 88: *> are stored on exit in array E), and
! 89: *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
! 90: *> If UPLO = 'L': factor L in the subdiagonal part of A.
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in] LDA
! 94: *> \verbatim
! 95: *> LDA is INTEGER
! 96: *> The leading dimension of the array A. LDA >= max(1,N).
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[out] E
! 100: *> \verbatim
! 101: *> E is COMPLEX*16 array, dimension (N)
! 102: *> On exit, contains the superdiagonal (or subdiagonal)
! 103: *> elements of the Hermitian block diagonal matrix D
! 104: *> with 1-by-1 or 2-by-2 diagonal blocks, where
! 105: *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
! 106: *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
! 107: *>
! 108: *> NOTE: For 1-by-1 diagonal block D(k), where
! 109: *> 1 <= k <= N, the element E(k) is set to 0 in both
! 110: *> UPLO = 'U' or UPLO = 'L' cases.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[out] IPIV
! 114: *> \verbatim
! 115: *> IPIV is INTEGER array, dimension (N)
! 116: *> IPIV describes the permutation matrix P in the factorization
! 117: *> of matrix A as follows. The absolute value of IPIV(k)
! 118: *> represents the index of row and column that were
! 119: *> interchanged with the k-th row and column. The value of UPLO
! 120: *> describes the order in which the interchanges were applied.
! 121: *> Also, the sign of IPIV represents the block structure of
! 122: *> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
! 123: *> diagonal blocks which correspond to 1 or 2 interchanges
! 124: *> at each factorization step. For more info see Further
! 125: *> Details section.
! 126: *>
! 127: *> If UPLO = 'U',
! 128: *> ( in factorization order, k decreases from N to 1 ):
! 129: *> a) A single positive entry IPIV(k) > 0 means:
! 130: *> D(k,k) is a 1-by-1 diagonal block.
! 131: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
! 132: *> interchanged in the matrix A(1:N,1:N);
! 133: *> If IPIV(k) = k, no interchange occurred.
! 134: *>
! 135: *> b) A pair of consecutive negative entries
! 136: *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
! 137: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
! 138: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
! 139: *> 1) If -IPIV(k) != k, rows and columns
! 140: *> k and -IPIV(k) were interchanged
! 141: *> in the matrix A(1:N,1:N).
! 142: *> If -IPIV(k) = k, no interchange occurred.
! 143: *> 2) If -IPIV(k-1) != k-1, rows and columns
! 144: *> k-1 and -IPIV(k-1) were interchanged
! 145: *> in the matrix A(1:N,1:N).
! 146: *> If -IPIV(k-1) = k-1, no interchange occurred.
! 147: *>
! 148: *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
! 149: *>
! 150: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
! 151: *>
! 152: *> If UPLO = 'L',
! 153: *> ( in factorization order, k increases from 1 to N ):
! 154: *> a) A single positive entry IPIV(k) > 0 means:
! 155: *> D(k,k) is a 1-by-1 diagonal block.
! 156: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
! 157: *> interchanged in the matrix A(1:N,1:N).
! 158: *> If IPIV(k) = k, no interchange occurred.
! 159: *>
! 160: *> b) A pair of consecutive negative entries
! 161: *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
! 162: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
! 163: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
! 164: *> 1) If -IPIV(k) != k, rows and columns
! 165: *> k and -IPIV(k) were interchanged
! 166: *> in the matrix A(1:N,1:N).
! 167: *> If -IPIV(k) = k, no interchange occurred.
! 168: *> 2) If -IPIV(k+1) != k+1, rows and columns
! 169: *> k-1 and -IPIV(k-1) were interchanged
! 170: *> in the matrix A(1:N,1:N).
! 171: *> If -IPIV(k+1) = k+1, no interchange occurred.
! 172: *>
! 173: *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
! 174: *>
! 175: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
! 176: *> \endverbatim
! 177: *>
! 178: *> \param[out] WORK
! 179: *> \verbatim
! 180: *> WORK is COMPLEX*16 array, dimension ( MAX(1,LWORK) ).
! 181: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 182: *> \endverbatim
! 183: *>
! 184: *> \param[in] LWORK
! 185: *> \verbatim
! 186: *> LWORK is INTEGER
! 187: *> The length of WORK. LWORK >=1. For best performance
! 188: *> LWORK >= N*NB, where NB is the block size returned
! 189: *> by ILAENV.
! 190: *>
! 191: *> If LWORK = -1, then a workspace query is assumed;
! 192: *> the routine only calculates the optimal size of the WORK
! 193: *> array, returns this value as the first entry of the WORK
! 194: *> array, and no error message related to LWORK is issued
! 195: *> by XERBLA.
! 196: *> \endverbatim
! 197: *>
! 198: *> \param[out] INFO
! 199: *> \verbatim
! 200: *> INFO is INTEGER
! 201: *> = 0: successful exit
! 202: *>
! 203: *> < 0: If INFO = -k, the k-th argument had an illegal value
! 204: *>
! 205: *> > 0: If INFO = k, the matrix A is singular, because:
! 206: *> If UPLO = 'U': column k in the upper
! 207: *> triangular part of A contains all zeros.
! 208: *> If UPLO = 'L': column k in the lower
! 209: *> triangular part of A contains all zeros.
! 210: *>
! 211: *> Therefore D(k,k) is exactly zero, and superdiagonal
! 212: *> elements of column k of U (or subdiagonal elements of
! 213: *> column k of L ) are all zeros. The factorization has
! 214: *> been completed, but the block diagonal matrix D is
! 215: *> exactly singular, and division by zero will occur if
! 216: *> it is used to solve a system of equations.
! 217: *>
! 218: *> NOTE: INFO only stores the first occurrence of
! 219: *> a singularity, any subsequent occurrence of singularity
! 220: *> is not stored in INFO even though the factorization
! 221: *> always completes.
! 222: *> \endverbatim
! 223: *
! 224: * Authors:
! 225: * ========
! 226: *
! 227: *> \author Univ. of Tennessee
! 228: *> \author Univ. of California Berkeley
! 229: *> \author Univ. of Colorado Denver
! 230: *> \author NAG Ltd.
! 231: *
! 232: *> \date December 2016
! 233: *
! 234: *> \ingroup complex16HEcomputational
! 235: *
! 236: *> \par Further Details:
! 237: * =====================
! 238: *>
! 239: *> \verbatim
! 240: *> TODO: put correct description
! 241: *> \endverbatim
! 242: *
! 243: *> \par Contributors:
! 244: * ==================
! 245: *>
! 246: *> \verbatim
! 247: *>
! 248: *> December 2016, Igor Kozachenko,
! 249: *> Computer Science Division,
! 250: *> University of California, Berkeley
! 251: *>
! 252: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
! 253: *> School of Mathematics,
! 254: *> University of Manchester
! 255: *>
! 256: *> \endverbatim
! 257: *
! 258: * =====================================================================
! 259: SUBROUTINE ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
! 260: $ INFO )
! 261: *
! 262: * -- LAPACK computational routine (version 3.7.0) --
! 263: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 264: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 265: * December 2016
! 266: *
! 267: * .. Scalar Arguments ..
! 268: CHARACTER UPLO
! 269: INTEGER INFO, LDA, LWORK, N
! 270: * ..
! 271: * .. Array Arguments ..
! 272: INTEGER IPIV( * )
! 273: COMPLEX*16 A( LDA, * ), E( * ), WORK( * )
! 274: * ..
! 275: *
! 276: * =====================================================================
! 277: *
! 278: * .. Local Scalars ..
! 279: LOGICAL LQUERY, UPPER
! 280: INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
! 281: $ NB, NBMIN
! 282: * ..
! 283: * .. External Functions ..
! 284: LOGICAL LSAME
! 285: INTEGER ILAENV
! 286: EXTERNAL LSAME, ILAENV
! 287: * ..
! 288: * .. External Subroutines ..
! 289: EXTERNAL ZLAHEF_RK, ZHETF2_RK, ZSWAP, XERBLA
! 290: * ..
! 291: * .. Intrinsic Functions ..
! 292: INTRINSIC ABS, MAX
! 293: * ..
! 294: * .. Executable Statements ..
! 295: *
! 296: * Test the input parameters.
! 297: *
! 298: INFO = 0
! 299: UPPER = LSAME( UPLO, 'U' )
! 300: LQUERY = ( LWORK.EQ.-1 )
! 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: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
! 308: INFO = -8
! 309: END IF
! 310: *
! 311: IF( INFO.EQ.0 ) THEN
! 312: *
! 313: * Determine the block size
! 314: *
! 315: NB = ILAENV( 1, 'ZHETRF_RK', UPLO, N, -1, -1, -1 )
! 316: LWKOPT = N*NB
! 317: WORK( 1 ) = LWKOPT
! 318: END IF
! 319: *
! 320: IF( INFO.NE.0 ) THEN
! 321: CALL XERBLA( 'ZHETRF_RK', -INFO )
! 322: RETURN
! 323: ELSE IF( LQUERY ) THEN
! 324: RETURN
! 325: END IF
! 326: *
! 327: NBMIN = 2
! 328: LDWORK = N
! 329: IF( NB.GT.1 .AND. NB.LT.N ) THEN
! 330: IWS = LDWORK*NB
! 331: IF( LWORK.LT.IWS ) THEN
! 332: NB = MAX( LWORK / LDWORK, 1 )
! 333: NBMIN = MAX( 2, ILAENV( 2, 'ZHETRF_RK',
! 334: $ UPLO, N, -1, -1, -1 ) )
! 335: END IF
! 336: ELSE
! 337: IWS = 1
! 338: END IF
! 339: IF( NB.LT.NBMIN )
! 340: $ NB = N
! 341: *
! 342: IF( UPPER ) THEN
! 343: *
! 344: * Factorize A as U*D*U**T using the upper triangle of A
! 345: *
! 346: * K is the main loop index, decreasing from N to 1 in steps of
! 347: * KB, where KB is the number of columns factorized by ZLAHEF_RK;
! 348: * KB is either NB or NB-1, or K for the last block
! 349: *
! 350: K = N
! 351: 10 CONTINUE
! 352: *
! 353: * If K < 1, exit from loop
! 354: *
! 355: IF( K.LT.1 )
! 356: $ GO TO 15
! 357: *
! 358: IF( K.GT.NB ) THEN
! 359: *
! 360: * Factorize columns k-kb+1:k of A and use blocked code to
! 361: * update columns 1:k-kb
! 362: *
! 363: CALL ZLAHEF_RK( UPLO, K, NB, KB, A, LDA, E,
! 364: $ IPIV, WORK, LDWORK, IINFO )
! 365: ELSE
! 366: *
! 367: * Use unblocked code to factorize columns 1:k of A
! 368: *
! 369: CALL ZHETF2_RK( UPLO, K, A, LDA, E, IPIV, IINFO )
! 370: KB = K
! 371: END IF
! 372: *
! 373: * Set INFO on the first occurrence of a zero pivot
! 374: *
! 375: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
! 376: $ INFO = IINFO
! 377: *
! 378: * No need to adjust IPIV
! 379: *
! 380: *
! 381: * Apply permutations to the leading panel 1:k-1
! 382: *
! 383: * Read IPIV from the last block factored, i.e.
! 384: * indices k-kb+1:k and apply row permutations to the
! 385: * last k+1 colunms k+1:N after that block
! 386: * (We can do the simple loop over IPIV with decrement -1,
! 387: * since the ABS value of IPIV( I ) represents the row index
! 388: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
! 389: *
! 390: IF( K.LT.N ) THEN
! 391: DO I = K, ( K - KB + 1 ), -1
! 392: IP = ABS( IPIV( I ) )
! 393: IF( IP.NE.I ) THEN
! 394: CALL ZSWAP( N-K, A( I, K+1 ), LDA,
! 395: $ A( IP, K+1 ), LDA )
! 396: END IF
! 397: END DO
! 398: END IF
! 399: *
! 400: * Decrease K and return to the start of the main loop
! 401: *
! 402: K = K - KB
! 403: GO TO 10
! 404: *
! 405: * This label is the exit from main loop over K decreasing
! 406: * from N to 1 in steps of KB
! 407: *
! 408: 15 CONTINUE
! 409: *
! 410: ELSE
! 411: *
! 412: * Factorize A as L*D*L**T using the lower triangle of A
! 413: *
! 414: * K is the main loop index, increasing from 1 to N in steps of
! 415: * KB, where KB is the number of columns factorized by ZLAHEF_RK;
! 416: * KB is either NB or NB-1, or N-K+1 for the last block
! 417: *
! 418: K = 1
! 419: 20 CONTINUE
! 420: *
! 421: * If K > N, exit from loop
! 422: *
! 423: IF( K.GT.N )
! 424: $ GO TO 35
! 425: *
! 426: IF( K.LE.N-NB ) THEN
! 427: *
! 428: * Factorize columns k:k+kb-1 of A and use blocked code to
! 429: * update columns k+kb:n
! 430: *
! 431: CALL ZLAHEF_RK( UPLO, N-K+1, NB, KB, A( K, K ), LDA, E( K ),
! 432: $ IPIV( K ), WORK, LDWORK, IINFO )
! 433:
! 434:
! 435: ELSE
! 436: *
! 437: * Use unblocked code to factorize columns k:n of A
! 438: *
! 439: CALL ZHETF2_RK( UPLO, N-K+1, A( K, K ), LDA, E( K ),
! 440: $ IPIV( K ), IINFO )
! 441: KB = N - K + 1
! 442: *
! 443: END IF
! 444: *
! 445: * Set INFO on the first occurrence of a zero pivot
! 446: *
! 447: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
! 448: $ INFO = IINFO + K - 1
! 449: *
! 450: * Adjust IPIV
! 451: *
! 452: DO I = K, K + KB - 1
! 453: IF( IPIV( I ).GT.0 ) THEN
! 454: IPIV( I ) = IPIV( I ) + K - 1
! 455: ELSE
! 456: IPIV( I ) = IPIV( I ) - K + 1
! 457: END IF
! 458: END DO
! 459: *
! 460: * Apply permutations to the leading panel 1:k-1
! 461: *
! 462: * Read IPIV from the last block factored, i.e.
! 463: * indices k:k+kb-1 and apply row permutations to the
! 464: * first k-1 colunms 1:k-1 before that block
! 465: * (We can do the simple loop over IPIV with increment 1,
! 466: * since the ABS value of IPIV( I ) represents the row index
! 467: * of the interchange with row i in both 1x1 and 2x2 pivot cases)
! 468: *
! 469: IF( K.GT.1 ) THEN
! 470: DO I = K, ( K + KB - 1 ), 1
! 471: IP = ABS( IPIV( I ) )
! 472: IF( IP.NE.I ) THEN
! 473: CALL ZSWAP( K-1, A( I, 1 ), LDA,
! 474: $ A( IP, 1 ), LDA )
! 475: END IF
! 476: END DO
! 477: END IF
! 478: *
! 479: * Increase K and return to the start of the main loop
! 480: *
! 481: K = K + KB
! 482: GO TO 20
! 483: *
! 484: * This label is the exit from main loop over K increasing
! 485: * from 1 to N in steps of KB
! 486: *
! 487: 35 CONTINUE
! 488: *
! 489: * End Lower
! 490: *
! 491: END IF
! 492: *
! 493: WORK( 1 ) = LWKOPT
! 494: RETURN
! 495: *
! 496: * End of ZHETRF_RK
! 497: *
! 498: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrf_rk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack