Annotation of rpl/lapack/lapack/zhetrf_aa.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETRF_AA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHETRF_AA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * CHARACTER UPLO
! 25: * INTEGER N, LDA, LWORK, INFO
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IPIV( * )
! 29: * COMPLEX*16 A( LDA, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *> \par Purpose:
! 33: * =============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *> ZHETRF_AA computes the factorization of a complex hermitian matrix A
! 38: *> using the Aasen's algorithm. The form of the factorization is
! 39: *>
! 40: *> A = U*T*U**H or A = L*T*L**H
! 41: *>
! 42: *> where U (or L) is a product of permutation and unit upper (lower)
! 43: *> triangular matrices, and T is a hermitian tridiagonal matrix.
! 44: *>
! 45: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
! 46: *> \endverbatim
! 47: *
! 48: * Arguments:
! 49: * ==========
! 50: *
! 51: *> \param[in] UPLO
! 52: *> \verbatim
! 53: *> UPLO is CHARACTER*1
! 54: *> = 'U': Upper triangle of A is stored;
! 55: *> = 'L': Lower triangle of A is stored.
! 56: *> \endverbatim
! 57: *>
! 58: *> \param[in] N
! 59: *> \verbatim
! 60: *> N is INTEGER
! 61: *> The order of the matrix A. N >= 0.
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in,out] A
! 65: *> \verbatim
! 66: *> A is COMPLEX*16 array, dimension (LDA,N)
! 67: *> On entry, the hermitian matrix A. If UPLO = 'U', the leading
! 68: *> N-by-N upper triangular part of A contains the upper
! 69: *> triangular part of the matrix A, and the strictly lower
! 70: *> triangular part of A is not referenced. If UPLO = 'L', the
! 71: *> leading N-by-N lower triangular part of A contains the lower
! 72: *> triangular part of the matrix A, and the strictly upper
! 73: *> triangular part of A is not referenced.
! 74: *>
! 75: *> On exit, the tridiagonal matrix is stored in the diagonals
! 76: *> and the subdiagonals of A just below (or above) the diagonals,
! 77: *> and L is stored below (or above) the subdiaonals, when UPLO
! 78: *> is 'L' (or 'U').
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] LDA
! 82: *> \verbatim
! 83: *> LDA is INTEGER
! 84: *> The leading dimension of the array A. LDA >= max(1,N).
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[out] IPIV
! 88: *> \verbatim
! 89: *> IPIV is INTEGER array, dimension (N)
! 90: *> On exit, it contains the details of the interchanges, i.e.,
! 91: *> the row and column k of A were interchanged with the
! 92: *> row and column IPIV(k).
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[out] WORK
! 96: *> \verbatim
! 97: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 98: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 99: *> \endverbatim
! 100: *>
! 101: *> \param[in] LWORK
! 102: *> \verbatim
! 103: *> LWORK is INTEGER
! 104: *> The length of WORK. LWORK >= MAX(1,2*N). For optimum performance
! 105: *> LWORK >= N*(1+NB), where NB is the optimal blocksize.
! 106: *>
! 107: *> If LWORK = -1, then a workspace query is assumed; the routine
! 108: *> only calculates the optimal size of the WORK array, returns
! 109: *> this value as the first entry of the WORK array, and no error
! 110: *> message related to LWORK is issued by XERBLA.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[out] INFO
! 114: *> \verbatim
! 115: *> INFO is INTEGER
! 116: *> = 0: successful exit
! 117: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 118: *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
! 119: *> has been completed, but the block diagonal matrix D is
! 120: *> exactly singular, and division by zero will occur if it
! 121: *> is used to solve a system of equations.
! 122: *> \endverbatim
! 123: *
! 124: * Authors:
! 125: * ========
! 126: *
! 127: *> \author Univ. of Tennessee
! 128: *> \author Univ. of California Berkeley
! 129: *> \author Univ. of Colorado Denver
! 130: *> \author NAG Ltd.
! 131: *
! 132: *> \date December 2016
! 133: *
! 134: *> \ingroup complex16HEcomputational
! 135: *
! 136: * =====================================================================
! 137: SUBROUTINE ZHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
! 138: *
! 139: * -- LAPACK computational routine (version 3.7.0) --
! 140: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 141: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 142: * December 2016
! 143: *
! 144: IMPLICIT NONE
! 145: *
! 146: * .. Scalar Arguments ..
! 147: CHARACTER UPLO
! 148: INTEGER N, LDA, LWORK, INFO
! 149: * ..
! 150: * .. Array Arguments ..
! 151: INTEGER IPIV( * )
! 152: COMPLEX*16 A( LDA, * ), WORK( * )
! 153: * ..
! 154: *
! 155: * =====================================================================
! 156: * .. Parameters ..
! 157: COMPLEX*16 ZERO, ONE
! 158: PARAMETER ( ZERO = (0.0D+0, 0.0D+0), ONE = (1.0D+0, 0.0D+0) )
! 159: *
! 160: * .. Local Scalars ..
! 161: LOGICAL LQUERY, UPPER
! 162: INTEGER J, LWKOPT, IINFO
! 163: INTEGER NB, MJ, NJ, K1, K2, J1, J2, J3, JB
! 164: COMPLEX*16 ALPHA
! 165: * ..
! 166: * .. External Functions ..
! 167: LOGICAL LSAME
! 168: INTEGER ILAENV
! 169: EXTERNAL LSAME, ILAENV
! 170: * ..
! 171: * .. External Subroutines ..
! 172: EXTERNAL XERBLA
! 173: * ..
! 174: * .. Intrinsic Functions ..
! 175: INTRINSIC DBLE, DCONJG, MAX
! 176: * ..
! 177: * .. Executable Statements ..
! 178: *
! 179: * Determine the block size
! 180: *
! 181: NB = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 )
! 182: *
! 183: * Test the input parameters.
! 184: *
! 185: INFO = 0
! 186: UPPER = LSAME( UPLO, 'U' )
! 187: LQUERY = ( LWORK.EQ.-1 )
! 188: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 189: INFO = -1
! 190: ELSE IF( N.LT.0 ) THEN
! 191: INFO = -2
! 192: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 193: INFO = -4
! 194: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
! 195: INFO = -7
! 196: END IF
! 197: *
! 198: IF( INFO.EQ.0 ) THEN
! 199: LWKOPT = (NB+1)*N
! 200: WORK( 1 ) = LWKOPT
! 201: END IF
! 202: *
! 203: IF( INFO.NE.0 ) THEN
! 204: CALL XERBLA( 'ZHETRF_AA', -INFO )
! 205: RETURN
! 206: ELSE IF( LQUERY ) THEN
! 207: RETURN
! 208: END IF
! 209: *
! 210: * Quick return
! 211: *
! 212: IF ( N.EQ.0 ) THEN
! 213: RETURN
! 214: ENDIF
! 215: IPIV( 1 ) = 1
! 216: IF ( N.EQ.1 ) THEN
! 217: A( 1, 1 ) = DBLE( A( 1, 1 ) )
! 218: IF ( A( 1, 1 ).EQ.ZERO ) THEN
! 219: INFO = 1
! 220: END IF
! 221: RETURN
! 222: END IF
! 223: *
! 224: * Adjubst block size based on the workspace size
! 225: *
! 226: IF( LWORK.LT.((1+NB)*N) ) THEN
! 227: NB = ( LWORK-N ) / N
! 228: END IF
! 229: *
! 230: IF( UPPER ) THEN
! 231: *
! 232: * .....................................................
! 233: * Factorize A as L*D*L**H using the upper triangle of A
! 234: * .....................................................
! 235: *
! 236: * copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
! 237: *
! 238: CALL ZCOPY( N, A( 1, 1 ), LDA, WORK( 1 ), 1 )
! 239: *
! 240: * J is the main loop index, increasing from 1 to N in steps of
! 241: * JB, where JB is the number of columns factorized by ZLAHEF;
! 242: * JB is either NB, or N-J+1 for the last block
! 243: *
! 244: J = 0
! 245: 10 CONTINUE
! 246: IF( J.GE.N )
! 247: $ GO TO 20
! 248: *
! 249: * each step of the main loop
! 250: * J is the last column of the previous panel
! 251: * J1 is the first column of the current panel
! 252: * K1 identifies if the previous column of the panel has been
! 253: * explicitly stored, e.g., K1=1 for the first panel, and
! 254: * K1=0 for the rest
! 255: *
! 256: J1 = J + 1
! 257: JB = MIN( N-J1+1, NB )
! 258: K1 = MAX(1, J)-J
! 259: *
! 260: * Panel factorization
! 261: *
! 262: CALL ZLAHEF_AA( UPLO, 2-K1, N-J, JB,
! 263: $ A( MAX(1, J), J+1 ), LDA,
! 264: $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ),
! 265: $ IINFO )
! 266: IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN
! 267: INFO = IINFO+J
! 268: ENDIF
! 269: *
! 270: * Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
! 271: *
! 272: DO J2 = J+2, MIN(N, J+JB+1)
! 273: IPIV( J2 ) = IPIV( J2 ) + J
! 274: IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN
! 275: CALL ZSWAP( J1-K1-2, A( 1, J2 ), 1,
! 276: $ A( 1, IPIV(J2) ), 1 )
! 277: END IF
! 278: END DO
! 279: J = J + JB
! 280: *
! 281: * Trailing submatrix update, where
! 282: * the row A(J1-1, J2-1:N) stores U(J1, J2+1:N) and
! 283: * WORK stores the current block of the auxiriarly matrix H
! 284: *
! 285: IF( J.LT.N ) THEN
! 286: *
! 287: * if the first panel and JB=1 (NB=1), then nothing to do
! 288: *
! 289: IF( J1.GT.1 .OR. JB.GT.1 ) THEN
! 290: *
! 291: * Merge rank-1 update with BLAS-3 update
! 292: *
! 293: ALPHA = DCONJG( A( J, J+1 ) )
! 294: A( J, J+1 ) = ONE
! 295: CALL ZCOPY( N-J, A( J-1, J+1 ), LDA,
! 296: $ WORK( (J+1-J1+1)+JB*N ), 1 )
! 297: CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 )
! 298: *
! 299: * K1 identifies if the previous column of the panel has been
! 300: * explicitly stored, e.g., K1=0 and K2=1 for the first panel,
! 301: * and K1=1 and K2=0 for the rest
! 302: *
! 303: IF( J1.GT.1 ) THEN
! 304: *
! 305: * Not first panel
! 306: *
! 307: K2 = 1
! 308: ELSE
! 309: *
! 310: * First panel
! 311: *
! 312: K2 = 0
! 313: *
! 314: * First update skips the first column
! 315: *
! 316: JB = JB - 1
! 317: END IF
! 318: *
! 319: DO J2 = J+1, N, NB
! 320: NJ = MIN( NB, N-J2+1 )
! 321: *
! 322: * Update (J2, J2) diagonal block with ZGEMV
! 323: *
! 324: J3 = J2
! 325: DO MJ = NJ-1, 1, -1
! 326: CALL ZGEMM( 'Conjugate transpose', 'Transpose',
! 327: $ 1, MJ, JB+1,
! 328: $ -ONE, A( J1-K2, J3 ), LDA,
! 329: $ WORK( (J3-J1+1)+K1*N ), N,
! 330: $ ONE, A( J3, J3 ), LDA )
! 331: J3 = J3 + 1
! 332: END DO
! 333: *
! 334: * Update off-diagonal block of J2-th block row with ZGEMM
! 335: *
! 336: CALL ZGEMM( 'Conjugate transpose', 'Transpose',
! 337: $ NJ, N-J3+1, JB+1,
! 338: $ -ONE, A( J1-K2, J2 ), LDA,
! 339: $ WORK( (J3-J1+1)+K1*N ), N,
! 340: $ ONE, A( J2, J3 ), LDA )
! 341: END DO
! 342: *
! 343: * Recover T( J, J+1 )
! 344: *
! 345: A( J, J+1 ) = DCONJG( ALPHA )
! 346: END IF
! 347: *
! 348: * WORK(J+1, 1) stores H(J+1, 1)
! 349: *
! 350: CALL ZCOPY( N-J, A( J+1, J+1 ), LDA, WORK( 1 ), 1 )
! 351: END IF
! 352: GO TO 10
! 353: ELSE
! 354: *
! 355: * .....................................................
! 356: * Factorize A as L*D*L**H using the lower triangle of A
! 357: * .....................................................
! 358: *
! 359: * copy first column A(1:N, 1) into H(1:N, 1)
! 360: * (stored in WORK(1:N))
! 361: *
! 362: CALL ZCOPY( N, A( 1, 1 ), 1, WORK( 1 ), 1 )
! 363: *
! 364: * J is the main loop index, increasing from 1 to N in steps of
! 365: * JB, where JB is the number of columns factorized by ZLAHEF;
! 366: * JB is either NB, or N-J+1 for the last block
! 367: *
! 368: J = 0
! 369: 11 CONTINUE
! 370: IF( J.GE.N )
! 371: $ GO TO 20
! 372: *
! 373: * each step of the main loop
! 374: * J is the last column of the previous panel
! 375: * J1 is the first column of the current panel
! 376: * K1 identifies if the previous column of the panel has been
! 377: * explicitly stored, e.g., K1=1 for the first panel, and
! 378: * K1=0 for the rest
! 379: *
! 380: J1 = J+1
! 381: JB = MIN( N-J1+1, NB )
! 382: K1 = MAX(1, J)-J
! 383: *
! 384: * Panel factorization
! 385: *
! 386: CALL ZLAHEF_AA( UPLO, 2-K1, N-J, JB,
! 387: $ A( J+1, MAX(1, J) ), LDA,
! 388: $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ), IINFO)
! 389: IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN
! 390: INFO = IINFO+J
! 391: ENDIF
! 392: *
! 393: * Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
! 394: *
! 395: DO J2 = J+2, MIN(N, J+JB+1)
! 396: IPIV( J2 ) = IPIV( J2 ) + J
! 397: IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN
! 398: CALL ZSWAP( J1-K1-2, A( J2, 1 ), LDA,
! 399: $ A( IPIV(J2), 1 ), LDA )
! 400: END IF
! 401: END DO
! 402: J = J + JB
! 403: *
! 404: * Trailing submatrix update, where
! 405: * A(J2+1, J1-1) stores L(J2+1, J1) and
! 406: * WORK(J2+1, 1) stores H(J2+1, 1)
! 407: *
! 408: IF( J.LT.N ) THEN
! 409: *
! 410: * if the first panel and JB=1 (NB=1), then nothing to do
! 411: *
! 412: IF( J1.GT.1 .OR. JB.GT.1 ) THEN
! 413: *
! 414: * Merge rank-1 update with BLAS-3 update
! 415: *
! 416: ALPHA = DCONJG( A( J+1, J ) )
! 417: A( J+1, J ) = ONE
! 418: CALL ZCOPY( N-J, A( J+1, J-1 ), 1,
! 419: $ WORK( (J+1-J1+1)+JB*N ), 1 )
! 420: CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 )
! 421: *
! 422: * K1 identifies if the previous column of the panel has been
! 423: * explicitly stored, e.g., K1=0 and K2=1 for the first panel,
! 424: * and K1=1 and K2=0 for the rest
! 425: *
! 426: IF( J1.GT.1 ) THEN
! 427: *
! 428: * Not first panel
! 429: *
! 430: K2 = 1
! 431: ELSE
! 432: *
! 433: * First panel
! 434: *
! 435: K2 = 0
! 436: *
! 437: * First update skips the first column
! 438: *
! 439: JB = JB - 1
! 440: END IF
! 441: *
! 442: DO J2 = J+1, N, NB
! 443: NJ = MIN( NB, N-J2+1 )
! 444: *
! 445: * Update (J2, J2) diagonal block with ZGEMV
! 446: *
! 447: J3 = J2
! 448: DO MJ = NJ-1, 1, -1
! 449: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
! 450: $ MJ, 1, JB+1,
! 451: $ -ONE, WORK( (J3-J1+1)+K1*N ), N,
! 452: $ A( J3, J1-K2 ), LDA,
! 453: $ ONE, A( J3, J3 ), LDA )
! 454: J3 = J3 + 1
! 455: END DO
! 456: *
! 457: * Update off-diagonal block of J2-th block column with ZGEMM
! 458: *
! 459: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
! 460: $ N-J3+1, NJ, JB+1,
! 461: $ -ONE, WORK( (J3-J1+1)+K1*N ), N,
! 462: $ A( J2, J1-K2 ), LDA,
! 463: $ ONE, A( J3, J2 ), LDA )
! 464: END DO
! 465: *
! 466: * Recover T( J+1, J )
! 467: *
! 468: A( J+1, J ) = DCONJG( ALPHA )
! 469: END IF
! 470: *
! 471: * WORK(J+1, 1) stores H(J+1, 1)
! 472: *
! 473: CALL ZCOPY( N-J, A( J+1, J+1 ), 1, WORK( 1 ), 1 )
! 474: END IF
! 475: GO TO 11
! 476: END IF
! 477: *
! 478: 20 CONTINUE
! 479: RETURN
! 480: *
! 481: * End of ZHETRF_AA
! 482: *
! 483: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrf_aa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack