Annotation of rpl/lapack/lapack/zhetrf_aa_2stage.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETRF_AA_2STAGE
! 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_2STAGE + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
! 22: * IPIV2, WORK, LWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER UPLO
! 26: * INTEGER N, LDA, LTB, LWORK, INFO
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IPIV( * ), IPIV2( * )
! 30: * COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
! 39: *> using the Aasen's algorithm. The form of the factorization is
! 40: *>
! 41: *> A = U*T*U**T or A = L*T*L**T
! 42: *>
! 43: *> where U (or L) is a product of permutation and unit upper (lower)
! 44: *> triangular matrices, and T is a hermitian band matrix with the
! 45: *> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is
! 46: *> LU factorized with partial pivoting).
! 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 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, L is stored below (or above) the subdiaonal blocks,
! 79: *> when UPLO is 'L' (or 'U').
! 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] TB
! 89: *> \verbatim
! 90: *> TB is COMPLEX array, dimension (LTB)
! 91: *> On exit, details of the LU factorization of the band matrix.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] LTB
! 95: *> \verbatim
! 96: *> The size of the array TB. LTB >= 4*N, internally
! 97: *> used to select NB such that LTB >= (3*NB+1)*N.
! 98: *>
! 99: *> If LTB = -1, then a workspace query is assumed; the
! 100: *> routine only calculates the optimal size of LTB,
! 101: *> returns this value as the first entry of TB, and
! 102: *> no error message related to LTB is issued by XERBLA.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[out] IPIV
! 106: *> \verbatim
! 107: *> IPIV is INTEGER array, dimension (N)
! 108: *> On exit, it contains the details of the interchanges, i.e.,
! 109: *> the row and column k of A were interchanged with the
! 110: *> row and column IPIV(k).
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[out] IPIV2
! 114: *> \verbatim
! 115: *> IPIV is INTEGER array, dimension (N)
! 116: *> On exit, it contains the details of the interchanges, i.e.,
! 117: *> the row and column k of T were interchanged with the
! 118: *> row and column IPIV(k).
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[out] WORK
! 122: *> \verbatim
! 123: *> WORK is COMPLEX workspace of size LWORK
! 124: *> \endverbatim
! 125: *>
! 126: *> \param[in] LWORK
! 127: *> \verbatim
! 128: *> The size of WORK. LWORK >= N, internally used to select NB
! 129: *> such that LWORK >= N*NB.
! 130: *>
! 131: *> If LWORK = -1, then a workspace query is assumed; the
! 132: *> routine only calculates the optimal size of the WORK array,
! 133: *> returns this value as the first entry of the WORK array, and
! 134: *> no error message related to LWORK is issued by XERBLA.
! 135: *> \endverbatim
! 136: *>
! 137: *> \param[out] INFO
! 138: *> \verbatim
! 139: *> INFO is INTEGER
! 140: *> = 0: successful exit
! 141: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 142: *> > 0: if INFO = i, band LU factorization failed on i-th column
! 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 2017
! 154: *
! 155: *> \ingroup complex16SYcomputational
! 156: *
! 157: * =====================================================================
! 158: SUBROUTINE ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
! 159: $ IPIV2, WORK, LWORK, INFO )
! 160: *
! 161: * -- LAPACK computational routine (version 3.8.0) --
! 162: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 163: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 164: * November 2017
! 165: *
! 166: IMPLICIT NONE
! 167: *
! 168: * .. Scalar Arguments ..
! 169: CHARACTER UPLO
! 170: INTEGER N, LDA, LTB, LWORK, INFO
! 171: * ..
! 172: * .. Array Arguments ..
! 173: INTEGER IPIV( * ), IPIV2( * )
! 174: COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
! 175: * ..
! 176: *
! 177: * =====================================================================
! 178: * .. Parameters ..
! 179: COMPLEX*16 ZERO, ONE
! 180: PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
! 181: $ ONE = ( 1.0E+0, 0.0E+0 ) )
! 182: *
! 183: * .. Local Scalars ..
! 184: LOGICAL UPPER, TQUERY, WQUERY
! 185: INTEGER I, J, K, I1, I2, TD
! 186: INTEGER LDTB, NB, KB, JB, NT, IINFO
! 187: COMPLEX*16 PIV
! 188: * ..
! 189: * .. External Functions ..
! 190: LOGICAL LSAME
! 191: INTEGER ILAENV
! 192: EXTERNAL LSAME, ILAENV
! 193: * ..
! 194: * .. External Subroutines ..
! 195: EXTERNAL XERBLA, ZCOPY, ZLACGV, ZLACPY,
! 196: $ ZLASET, ZGBTRF, ZGEMM, ZGETRF,
! 197: $ ZHEGST, ZSWAP, ZTRSM
! 198: * ..
! 199: * .. Intrinsic Functions ..
! 200: INTRINSIC DCONJG, MIN, MAX
! 201: * ..
! 202: * .. Executable Statements ..
! 203: *
! 204: * Test the input parameters.
! 205: *
! 206: INFO = 0
! 207: UPPER = LSAME( UPLO, 'U' )
! 208: WQUERY = ( LWORK.EQ.-1 )
! 209: TQUERY = ( LTB.EQ.-1 )
! 210: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 211: INFO = -1
! 212: ELSE IF( N.LT.0 ) THEN
! 213: INFO = -2
! 214: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 215: INFO = -4
! 216: ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN
! 217: INFO = -6
! 218: ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN
! 219: INFO = -10
! 220: END IF
! 221: *
! 222: IF( INFO.NE.0 ) THEN
! 223: CALL XERBLA( 'ZHETRF_AA_2STAGE', -INFO )
! 224: RETURN
! 225: END IF
! 226: *
! 227: * Answer the query
! 228: *
! 229: NB = ILAENV( 1, 'ZHETRF_AA_2STAGE', UPLO, N, -1, -1, -1 )
! 230: IF( INFO.EQ.0 ) THEN
! 231: IF( TQUERY ) THEN
! 232: TB( 1 ) = (3*NB+1)*N
! 233: END IF
! 234: IF( WQUERY ) THEN
! 235: WORK( 1 ) = N*NB
! 236: END IF
! 237: END IF
! 238: IF( TQUERY .OR. WQUERY ) THEN
! 239: RETURN
! 240: END IF
! 241: *
! 242: * Quick return
! 243: *
! 244: IF ( N.EQ.0 ) THEN
! 245: RETURN
! 246: ENDIF
! 247: *
! 248: * Determine the number of the block size
! 249: *
! 250: LDTB = LTB/N
! 251: IF( LDTB .LT. 3*NB+1 ) THEN
! 252: NB = (LDTB-1)/3
! 253: END IF
! 254: IF( LWORK .LT. NB*N ) THEN
! 255: NB = LWORK/N
! 256: END IF
! 257: *
! 258: * Determine the number of the block columns
! 259: *
! 260: NT = (N+NB-1)/NB
! 261: TD = 2*NB
! 262: KB = MIN(NB, N)
! 263: *
! 264: * Initialize vectors/matrices
! 265: *
! 266: DO J = 1, KB
! 267: IPIV( J ) = J
! 268: END DO
! 269: *
! 270: * Save NB
! 271: *
! 272: TB( 1 ) = NB
! 273: *
! 274: IF( UPPER ) THEN
! 275: *
! 276: * .....................................................
! 277: * Factorize A as L*D*L**T using the upper triangle of A
! 278: * .....................................................
! 279: *
! 280: DO J = 0, NT-1
! 281: *
! 282: * Generate Jth column of W and H
! 283: *
! 284: KB = MIN(NB, N-J*NB)
! 285: DO I = 1, J-1
! 286: IF( I.EQ.1 ) THEN
! 287: * H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
! 288: IF( I .EQ. (J-1) ) THEN
! 289: JB = NB+KB
! 290: ELSE
! 291: JB = 2*NB
! 292: END IF
! 293: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 294: $ NB, KB, JB,
! 295: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 296: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
! 297: $ ZERO, WORK( I*NB+1 ), N )
! 298: ELSE
! 299: * H(I,J) = T(I,I-1)*U(I-1,J) + T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
! 300: IF( I .EQ. (J-1) ) THEN
! 301: JB = 2*NB+KB
! 302: ELSE
! 303: JB = 3*NB
! 304: END IF
! 305: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 306: $ NB, KB, JB,
! 307: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 308: $ LDTB-1,
! 309: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
! 310: $ ZERO, WORK( I*NB+1 ), N )
! 311: END IF
! 312: END DO
! 313: *
! 314: * Compute T(J,J)
! 315: *
! 316: CALL ZLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
! 317: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 318: IF( J.GT.1 ) THEN
! 319: * T(J,J) = U(1:J,J)'*H(1:J)
! 320: CALL ZGEMM( 'Conjugate transpose', 'NoTranspose',
! 321: $ KB, KB, (J-1)*NB,
! 322: $ -ONE, A( 1, J*NB+1 ), LDA,
! 323: $ WORK( NB+1 ), N,
! 324: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 325: * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
! 326: CALL ZGEMM( 'Conjugate transpose', 'NoTranspose',
! 327: $ KB, NB, KB,
! 328: $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
! 329: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 330: $ ZERO, WORK( 1 ), N )
! 331: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 332: $ KB, KB, NB,
! 333: $ -ONE, WORK( 1 ), N,
! 334: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 335: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 336: END IF
! 337: IF( J.GT.0 ) THEN
! 338: CALL ZHEGST( 1, 'Upper', KB,
! 339: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 340: $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
! 341: END IF
! 342: *
! 343: * Expand T(J,J) into full format
! 344: *
! 345: DO I = 1, KB
! 346: TB( TD+1 + (J*NB+I-1)*LDTB )
! 347: $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) )
! 348: DO K = I+1, KB
! 349: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 350: $ = DCONJG( TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB ) )
! 351: END DO
! 352: END DO
! 353: *
! 354: IF( J.LT.NT-1 ) THEN
! 355: IF( J.GT.0 ) THEN
! 356: *
! 357: * Compute H(J,J)
! 358: *
! 359: IF( J.EQ.1 ) THEN
! 360: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 361: $ KB, KB, KB,
! 362: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 363: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 364: $ ZERO, WORK( J*NB+1 ), N )
! 365: ELSE
! 366: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 367: $ KB, KB, NB+KB,
! 368: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 369: $ LDTB-1,
! 370: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 371: $ ZERO, WORK( J*NB+1 ), N )
! 372: END IF
! 373: *
! 374: * Update with the previous column
! 375: *
! 376: CALL ZGEMM( 'Conjugate transpose', 'NoTranspose',
! 377: $ NB, N-(J+1)*NB, J*NB,
! 378: $ -ONE, WORK( NB+1 ), N,
! 379: $ A( 1, (J+1)*NB+1 ), LDA,
! 380: $ ONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
! 381: END IF
! 382: *
! 383: * Copy panel to workspace to call ZGETRF
! 384: *
! 385: DO K = 1, NB
! 386: CALL ZCOPY( N-(J+1)*NB,
! 387: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
! 388: $ WORK( 1+(K-1)*N ), 1 )
! 389: END DO
! 390: *
! 391: * Factorize panel
! 392: *
! 393: CALL ZGETRF( N-(J+1)*NB, NB,
! 394: $ WORK, N,
! 395: $ IPIV( (J+1)*NB+1 ), IINFO )
! 396: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 397: c INFO = IINFO+(J+1)*NB
! 398: c END IF
! 399: *
! 400: * Copy panel back
! 401: *
! 402: DO K = 1, NB
! 403: *
! 404: * Copy only L-factor
! 405: *
! 406: CALL ZCOPY( N-K-(J+1)*NB,
! 407: $ WORK( K+1+(K-1)*N ), 1,
! 408: $ A( J*NB+K, (J+1)*NB+K+1 ), LDA )
! 409: *
! 410: * Transpose U-factor to be copied back into T(J+1, J)
! 411: *
! 412: CALL ZLACGV( K, WORK( 1+(K-1)*N ), 1 )
! 413: END DO
! 414: *
! 415: * Compute T(J+1, J), zero out for GEMM update
! 416: *
! 417: KB = MIN(NB, N-(J+1)*NB)
! 418: CALL ZLASET( 'Full', KB, NB, ZERO, ZERO,
! 419: $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 )
! 420: CALL ZLACPY( 'Upper', KB, NB,
! 421: $ WORK, N,
! 422: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 423: IF( J.GT.0 ) THEN
! 424: CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE,
! 425: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 426: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 427: END IF
! 428: *
! 429: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
! 430: * updates
! 431: *
! 432: DO K = 1, NB
! 433: DO I = 1, KB
! 434: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
! 435: $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) )
! 436: END DO
! 437: END DO
! 438: CALL ZLASET( 'Lower', KB, NB, ZERO, ONE,
! 439: $ A( J*NB+1, (J+1)*NB+1), LDA )
! 440: *
! 441: * Apply pivots to trailing submatrix of A
! 442: *
! 443: DO K = 1, KB
! 444: * > Adjust ipiv
! 445: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 446: *
! 447: I1 = (J+1)*NB+K
! 448: I2 = IPIV( (J+1)*NB+K )
! 449: IF( I1.NE.I2 ) THEN
! 450: * > Apply pivots to previous columns of L
! 451: CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
! 452: $ A( (J+1)*NB+1, I2 ), 1 )
! 453: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 454: CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
! 455: $ A( I1+1, I2 ), 1 )
! 456: CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA )
! 457: CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 )
! 458: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 459: CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
! 460: $ A( I2, I2+1 ), LDA )
! 461: * > Swap A(I1, I1) with A(I2, I2)
! 462: PIV = A( I1, I1 )
! 463: A( I1, I1 ) = A( I2, I2 )
! 464: A( I2, I2 ) = PIV
! 465: * > Apply pivots to previous columns of L
! 466: IF( J.GT.0 ) THEN
! 467: CALL ZSWAP( J*NB, A( 1, I1 ), 1,
! 468: $ A( 1, I2 ), 1 )
! 469: END IF
! 470: ENDIF
! 471: END DO
! 472: END IF
! 473: END DO
! 474: ELSE
! 475: *
! 476: * .....................................................
! 477: * Factorize A as L*D*L**T using the lower triangle of A
! 478: * .....................................................
! 479: *
! 480: DO J = 0, NT-1
! 481: *
! 482: * Generate Jth column of W and H
! 483: *
! 484: KB = MIN(NB, N-J*NB)
! 485: DO I = 1, J-1
! 486: IF( I.EQ.1 ) THEN
! 487: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
! 488: IF( I .EQ. (J-1) ) THEN
! 489: JB = NB+KB
! 490: ELSE
! 491: JB = 2*NB
! 492: END IF
! 493: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
! 494: $ NB, KB, JB,
! 495: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 496: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
! 497: $ ZERO, WORK( I*NB+1 ), N )
! 498: ELSE
! 499: * H(I,J) = T(I,I-1)*L(J,I-1)' + T(I,I)*L(J,I)' + T(I,I+1)*L(J,I+1)'
! 500: IF( I .EQ. (J-1) ) THEN
! 501: JB = 2*NB+KB
! 502: ELSE
! 503: JB = 3*NB
! 504: END IF
! 505: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
! 506: $ NB, KB, JB,
! 507: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 508: $ LDTB-1,
! 509: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
! 510: $ ZERO, WORK( I*NB+1 ), N )
! 511: END IF
! 512: END DO
! 513: *
! 514: * Compute T(J,J)
! 515: *
! 516: CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
! 517: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 518: IF( J.GT.1 ) THEN
! 519: * T(J,J) = L(J,1:J)*H(1:J)
! 520: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 521: $ KB, KB, (J-1)*NB,
! 522: $ -ONE, A( J*NB+1, 1 ), LDA,
! 523: $ WORK( NB+1 ), N,
! 524: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 525: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
! 526: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 527: $ KB, NB, KB,
! 528: $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
! 529: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 530: $ ZERO, WORK( 1 ), N )
! 531: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
! 532: $ KB, KB, NB,
! 533: $ -ONE, WORK( 1 ), N,
! 534: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 535: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 536: END IF
! 537: IF( J.GT.0 ) THEN
! 538: CALL ZHEGST( 1, 'Lower', KB,
! 539: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 540: $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
! 541: END IF
! 542: *
! 543: * Expand T(J,J) into full format
! 544: *
! 545: DO I = 1, KB
! 546: TB( TD+1 + (J*NB+I-1)*LDTB )
! 547: $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) )
! 548: DO K = I+1, KB
! 549: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 550: $ = DCONJG( TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) )
! 551: END DO
! 552: END DO
! 553: *
! 554: IF( J.LT.NT-1 ) THEN
! 555: IF( J.GT.0 ) THEN
! 556: *
! 557: * Compute H(J,J)
! 558: *
! 559: IF( J.EQ.1 ) THEN
! 560: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
! 561: $ KB, KB, KB,
! 562: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 563: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 564: $ ZERO, WORK( J*NB+1 ), N )
! 565: ELSE
! 566: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
! 567: $ KB, KB, NB+KB,
! 568: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 569: $ LDTB-1,
! 570: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 571: $ ZERO, WORK( J*NB+1 ), N )
! 572: END IF
! 573: *
! 574: * Update with the previous column
! 575: *
! 576: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 577: $ N-(J+1)*NB, NB, J*NB,
! 578: $ -ONE, A( (J+1)*NB+1, 1 ), LDA,
! 579: $ WORK( NB+1 ), N,
! 580: $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
! 581: END IF
! 582: *
! 583: * Factorize panel
! 584: *
! 585: CALL ZGETRF( N-(J+1)*NB, NB,
! 586: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 587: $ IPIV( (J+1)*NB+1 ), IINFO )
! 588: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 589: c INFO = IINFO+(J+1)*NB
! 590: c END IF
! 591: *
! 592: * Compute T(J+1, J), zero out for GEMM update
! 593: *
! 594: KB = MIN(NB, N-(J+1)*NB)
! 595: CALL ZLASET( 'Full', KB, NB, ZERO, ZERO,
! 596: $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 )
! 597: CALL ZLACPY( 'Upper', KB, NB,
! 598: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 599: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 600: IF( J.GT.0 ) THEN
! 601: CALL ZTRSM( 'R', 'L', 'C', 'U', KB, NB, ONE,
! 602: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 603: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 604: END IF
! 605: *
! 606: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
! 607: * updates
! 608: *
! 609: DO K = 1, NB
! 610: DO I = 1, KB
! 611: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
! 612: $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) )
! 613: END DO
! 614: END DO
! 615: CALL ZLASET( 'Upper', KB, NB, ZERO, ONE,
! 616: $ A( (J+1)*NB+1, J*NB+1), LDA )
! 617: *
! 618: * Apply pivots to trailing submatrix of A
! 619: *
! 620: DO K = 1, KB
! 621: * > Adjust ipiv
! 622: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 623: *
! 624: I1 = (J+1)*NB+K
! 625: I2 = IPIV( (J+1)*NB+K )
! 626: IF( I1.NE.I2 ) THEN
! 627: * > Apply pivots to previous columns of L
! 628: CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
! 629: $ A( I2, (J+1)*NB+1 ), LDA )
! 630: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 631: CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
! 632: $ A( I2, I1+1 ), LDA )
! 633: CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 )
! 634: CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA )
! 635: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 636: CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
! 637: $ A( I2+1, I2 ), 1 )
! 638: * > Swap A(I1, I1) with A(I2, I2)
! 639: PIV = A( I1, I1 )
! 640: A( I1, I1 ) = A( I2, I2 )
! 641: A( I2, I2 ) = PIV
! 642: * > Apply pivots to previous columns of L
! 643: IF( J.GT.0 ) THEN
! 644: CALL ZSWAP( J*NB, A( I1, 1 ), LDA,
! 645: $ A( I2, 1 ), LDA )
! 646: END IF
! 647: ENDIF
! 648: END DO
! 649: *
! 650: * Apply pivots to previous columns of L
! 651: *
! 652: c CALL ZLASWP( J*NB, A( 1, 1 ), LDA,
! 653: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
! 654: END IF
! 655: END DO
! 656: END IF
! 657: *
! 658: * Factor the band matrix
! 659: CALL ZGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
! 660: *
! 661: * End of ZHETRF_AA_2STAGE
! 662: *
! 663: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrf_aa_2stage.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack