Annotation of rpl/lapack/lapack/zhetrf_rook.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETRF_ROOK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (blocked algorithm, calling Level 3 BLAS).
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHETRF_ROOK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_rook.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_rook.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_rook.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * CHARACTER UPLO
! 25: * INTEGER INFO, LDA, LWORK, N
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IPIV( * )
! 29: * COMPLEX*16 A( LDA, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> ZHETRF_ROOK computes the factorization of a complex Hermitian matrix A
! 39: *> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
! 40: *> The form of the factorization is
! 41: *>
! 42: *> A = U*D*U**T or A = L*D*L**T
! 43: *>
! 44: *> where U (or L) is a product of permutation and unit upper (lower)
! 45: *> triangular matrices, and D is Hermitian and block diagonal with
! 46: *> 1-by-1 and 2-by-2 diagonal blocks.
! 47: *>
! 48: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
! 49: *> \endverbatim
! 50: *
! 51: * Arguments:
! 52: * ==========
! 53: *
! 54: *> \param[in] UPLO
! 55: *> \verbatim
! 56: *> UPLO is CHARACTER*1
! 57: *> = 'U': Upper triangle of A is stored;
! 58: *> = 'L': Lower triangle of A is stored.
! 59: *> \endverbatim
! 60: *>
! 61: *> \param[in] N
! 62: *> \verbatim
! 63: *> N is INTEGER
! 64: *> The order of the matrix A. N >= 0.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in,out] A
! 68: *> \verbatim
! 69: *> A is COMPLEX*16 array, dimension (LDA,N)
! 70: *> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
! 71: *> N-by-N upper triangular part of A contains the upper
! 72: *> triangular part of the matrix A, and the strictly lower
! 73: *> triangular part of A is not referenced. If UPLO = 'L', the
! 74: *> leading N-by-N lower triangular part of A contains the lower
! 75: *> triangular part of the matrix A, and the strictly upper
! 76: *> triangular part of A is not referenced.
! 77: *>
! 78: *> On exit, the block diagonal matrix D and the multipliers used
! 79: *> to obtain the factor U or L (see below for further details).
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] LDA
! 83: *> \verbatim
! 84: *> LDA is INTEGER
! 85: *> The leading dimension of the array A. LDA >= max(1,N).
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[out] IPIV
! 89: *> \verbatim
! 90: *> IPIV is INTEGER array, dimension (N)
! 91: *> Details of the interchanges and the block structure of D.
! 92: *>
! 93: *> If UPLO = 'U':
! 94: *> Only the last KB elements of IPIV are set.
! 95: *>
! 96: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
! 97: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
! 98: *>
! 99: *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
! 100: *> columns k and -IPIV(k) were interchanged and rows and
! 101: *> columns k-1 and -IPIV(k-1) were inerchaged,
! 102: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
! 103: *>
! 104: *> If UPLO = 'L':
! 105: *> Only the first KB elements of IPIV are set.
! 106: *>
! 107: *> If IPIV(k) > 0, then rows and columns k and IPIV(k)
! 108: *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
! 109: *>
! 110: *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
! 111: *> columns k and -IPIV(k) were interchanged and rows and
! 112: *> columns k+1 and -IPIV(k+1) were inerchaged,
! 113: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[out] WORK
! 117: *> \verbatim
! 118: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)).
! 119: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[in] LWORK
! 123: *> \verbatim
! 124: *> LWORK is INTEGER
! 125: *> The length of WORK. LWORK >=1. For best performance
! 126: *> LWORK >= N*NB, where NB is the block size returned by ILAENV.
! 127: *>
! 128: *> If LWORK = -1, then a workspace query is assumed; the routine
! 129: *> only calculates the optimal size of the WORK array, returns
! 130: *> this value as the first entry of the WORK array, and no error
! 131: *> message related to LWORK is issued by XERBLA.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[out] INFO
! 135: *> \verbatim
! 136: *> INFO is INTEGER
! 137: *> = 0: successful exit
! 138: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 139: *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
! 140: *> has been completed, but the block diagonal matrix D is
! 141: *> exactly singular, and division by zero will occur if it
! 142: *> is used to solve a system of equations.
! 143: *> \endverbatim
! 144: *
! 145: * Authors:
! 146: * ========
! 147: *
! 148: *> \author Univ. of Tennessee
! 149: *> \author Univ. of California Berkeley
! 150: *> \author Univ. of Colorado Denver
! 151: *> \author NAG Ltd.
! 152: *
! 153: *> \date November 2013
! 154: *
! 155: *> \ingroup complex16HEcomputational
! 156: *
! 157: *> \par Further Details:
! 158: * =====================
! 159: *>
! 160: *> \verbatim
! 161: *>
! 162: *> If UPLO = 'U', then A = U*D*U**T, where
! 163: *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
! 164: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
! 165: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
! 166: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
! 167: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
! 168: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
! 169: *>
! 170: *> ( I v 0 ) k-s
! 171: *> U(k) = ( 0 I 0 ) s
! 172: *> ( 0 0 I ) n-k
! 173: *> k-s s n-k
! 174: *>
! 175: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
! 176: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
! 177: *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
! 178: *>
! 179: *> If UPLO = 'L', then A = L*D*L**T, where
! 180: *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
! 181: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
! 182: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
! 183: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
! 184: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
! 185: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
! 186: *>
! 187: *> ( I 0 0 ) k-1
! 188: *> L(k) = ( 0 I 0 ) s
! 189: *> ( 0 v I ) n-k-s+1
! 190: *> k-1 s n-k-s+1
! 191: *>
! 192: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
! 193: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
! 194: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
! 195: *> \endverbatim
! 196: *
! 197: *> \par Contributors:
! 198: * ==================
! 199: *>
! 200: *> \verbatim
! 201: *>
! 202: *> November 2013, Igor Kozachenko,
! 203: *> Computer Science Division,
! 204: *> University of California, Berkeley
! 205: *>
! 206: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
! 207: *> School of Mathematics,
! 208: *> University of Manchester
! 209: *>
! 210: *> \endverbatim
! 211: *
! 212: * =====================================================================
! 213: SUBROUTINE ZHETRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
! 214: *
! 215: * -- LAPACK computational routine (version 3.5.0) --
! 216: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 217: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 218: * November 2013
! 219: *
! 220: * .. Scalar Arguments ..
! 221: CHARACTER UPLO
! 222: INTEGER INFO, LDA, LWORK, N
! 223: * ..
! 224: * .. Array Arguments ..
! 225: INTEGER IPIV( * )
! 226: COMPLEX*16 A( LDA, * ), WORK( * )
! 227: * ..
! 228: *
! 229: * =====================================================================
! 230: *
! 231: * .. Local Scalars ..
! 232: LOGICAL LQUERY, UPPER
! 233: INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
! 234: * ..
! 235: * .. External Functions ..
! 236: LOGICAL LSAME
! 237: INTEGER ILAENV
! 238: EXTERNAL LSAME, ILAENV
! 239: * ..
! 240: * .. External Subroutines ..
! 241: EXTERNAL ZLAHEF_ROOK, ZHETF2_ROOK, XERBLA
! 242: * ..
! 243: * .. Intrinsic Functions ..
! 244: INTRINSIC MAX
! 245: * ..
! 246: * .. Executable Statements ..
! 247: *
! 248: * Test the input parameters.
! 249: *
! 250: INFO = 0
! 251: UPPER = LSAME( UPLO, 'U' )
! 252: LQUERY = ( LWORK.EQ.-1 )
! 253: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 254: INFO = -1
! 255: ELSE IF( N.LT.0 ) THEN
! 256: INFO = -2
! 257: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 258: INFO = -4
! 259: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
! 260: INFO = -7
! 261: END IF
! 262: *
! 263: IF( INFO.EQ.0 ) THEN
! 264: *
! 265: * Determine the block size
! 266: *
! 267: NB = ILAENV( 1, 'ZHETRF_ROOK', UPLO, N, -1, -1, -1 )
! 268: LWKOPT = N*NB
! 269: WORK( 1 ) = LWKOPT
! 270: END IF
! 271: *
! 272: IF( INFO.NE.0 ) THEN
! 273: CALL XERBLA( 'ZHETRF_ROOK', -INFO )
! 274: RETURN
! 275: ELSE IF( LQUERY ) THEN
! 276: RETURN
! 277: END IF
! 278: *
! 279: NBMIN = 2
! 280: LDWORK = N
! 281: IF( NB.GT.1 .AND. NB.LT.N ) THEN
! 282: IWS = LDWORK*NB
! 283: IF( LWORK.LT.IWS ) THEN
! 284: NB = MAX( LWORK / LDWORK, 1 )
! 285: NBMIN = MAX( 2, ILAENV( 2, 'ZHETRF_ROOK',
! 286: $ UPLO, N, -1, -1, -1 ) )
! 287: END IF
! 288: ELSE
! 289: IWS = 1
! 290: END IF
! 291: IF( NB.LT.NBMIN )
! 292: $ NB = N
! 293: *
! 294: IF( UPPER ) THEN
! 295: *
! 296: * Factorize A as U*D*U**T using the upper triangle of A
! 297: *
! 298: * K is the main loop index, decreasing from N to 1 in steps of
! 299: * KB, where KB is the number of columns factorized by ZLAHEF_ROOK;
! 300: * KB is either NB or NB-1, or K for the last block
! 301: *
! 302: K = N
! 303: 10 CONTINUE
! 304: *
! 305: * If K < 1, exit from loop
! 306: *
! 307: IF( K.LT.1 )
! 308: $ GO TO 40
! 309: *
! 310: IF( K.GT.NB ) THEN
! 311: *
! 312: * Factorize columns k-kb+1:k of A and use blocked code to
! 313: * update columns 1:k-kb
! 314: *
! 315: CALL ZLAHEF_ROOK( UPLO, K, NB, KB, A, LDA,
! 316: $ IPIV, WORK, LDWORK, IINFO )
! 317: ELSE
! 318: *
! 319: * Use unblocked code to factorize columns 1:k of A
! 320: *
! 321: CALL ZHETF2_ROOK( UPLO, K, A, LDA, IPIV, IINFO )
! 322: KB = K
! 323: END IF
! 324: *
! 325: * Set INFO on the first occurrence of a zero pivot
! 326: *
! 327: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
! 328: $ INFO = IINFO
! 329: *
! 330: * No need to adjust IPIV
! 331: *
! 332: * Decrease K and return to the start of the main loop
! 333: *
! 334: K = K - KB
! 335: GO TO 10
! 336: *
! 337: ELSE
! 338: *
! 339: * Factorize A as L*D*L**T using the lower triangle of A
! 340: *
! 341: * K is the main loop index, increasing from 1 to N in steps of
! 342: * KB, where KB is the number of columns factorized by ZLAHEF_ROOK;
! 343: * KB is either NB or NB-1, or N-K+1 for the last block
! 344: *
! 345: K = 1
! 346: 20 CONTINUE
! 347: *
! 348: * If K > N, exit from loop
! 349: *
! 350: IF( K.GT.N )
! 351: $ GO TO 40
! 352: *
! 353: IF( K.LE.N-NB ) THEN
! 354: *
! 355: * Factorize columns k:k+kb-1 of A and use blocked code to
! 356: * update columns k+kb:n
! 357: *
! 358: CALL ZLAHEF_ROOK( UPLO, N-K+1, NB, KB, A( K, K ), LDA,
! 359: $ IPIV( K ), WORK, LDWORK, IINFO )
! 360: ELSE
! 361: *
! 362: * Use unblocked code to factorize columns k:n of A
! 363: *
! 364: CALL ZHETF2_ROOK( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ),
! 365: $ IINFO )
! 366: KB = N - K + 1
! 367: END IF
! 368: *
! 369: * Set INFO on the first occurrence of a zero pivot
! 370: *
! 371: IF( INFO.EQ.0 .AND. IINFO.GT.0 )
! 372: $ INFO = IINFO + K - 1
! 373: *
! 374: * Adjust IPIV
! 375: *
! 376: DO 30 J = K, K + KB - 1
! 377: IF( IPIV( J ).GT.0 ) THEN
! 378: IPIV( J ) = IPIV( J ) + K - 1
! 379: ELSE
! 380: IPIV( J ) = IPIV( J ) - K + 1
! 381: END IF
! 382: 30 CONTINUE
! 383: *
! 384: * Increase K and return to the start of the main loop
! 385: *
! 386: K = K + KB
! 387: GO TO 20
! 388: *
! 389: END IF
! 390: *
! 391: 40 CONTINUE
! 392: WORK( 1 ) = LWKOPT
! 393: RETURN
! 394: *
! 395: * End of ZHETRF_ROOK
! 396: *
! 397: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrf_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack