Annotation of rpl/lapack/lapack/zsytrf_aa_2stage.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZSYTRF_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 ZSYTRF_AA_2STAGE + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZSYTRF_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: *> ZSYTRF_AA_2STAGE computes the factorization of a complex symmetric 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 complex symmetric 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*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, 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*16 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*16 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 ZSYTRF_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 CZERO, CONE
! 180: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 181: $ CONE = ( 1.0D+0, 0.0D+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, ZGBTRF, ZGEMM, ZGETRF,
! 196: $ ZLACPY, ZLASET, ZLASWP, ZTRSM, ZSWAP
! 197: * ..
! 198: * .. Intrinsic Functions ..
! 199: INTRINSIC MIN, MAX
! 200: * ..
! 201: * .. Executable Statements ..
! 202: *
! 203: * Test the input parameters.
! 204: *
! 205: INFO = 0
! 206: UPPER = LSAME( UPLO, 'U' )
! 207: WQUERY = ( LWORK.EQ.-1 )
! 208: TQUERY = ( LTB.EQ.-1 )
! 209: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 210: INFO = -1
! 211: ELSE IF( N.LT.0 ) THEN
! 212: INFO = -2
! 213: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 214: INFO = -4
! 215: ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN
! 216: INFO = -6
! 217: ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN
! 218: INFO = -10
! 219: END IF
! 220: *
! 221: IF( INFO.NE.0 ) THEN
! 222: CALL XERBLA( 'ZSYTRF_AA_2STAGE', -INFO )
! 223: RETURN
! 224: END IF
! 225: *
! 226: * Answer the query
! 227: *
! 228: NB = ILAENV( 1, 'ZSYTRF_AA_2STAGE', UPLO, N, -1, -1, -1 )
! 229: IF( INFO.EQ.0 ) THEN
! 230: IF( TQUERY ) THEN
! 231: TB( 1 ) = (3*NB+1)*N
! 232: END IF
! 233: IF( WQUERY ) THEN
! 234: WORK( 1 ) = N*NB
! 235: END IF
! 236: END IF
! 237: IF( TQUERY .OR. WQUERY ) THEN
! 238: RETURN
! 239: END IF
! 240: *
! 241: * Quick return
! 242: *
! 243: IF ( N.EQ.0 ) THEN
! 244: RETURN
! 245: ENDIF
! 246: *
! 247: * Determine the number of the block size
! 248: *
! 249: LDTB = LTB/N
! 250: IF( LDTB .LT. 3*NB+1 ) THEN
! 251: NB = (LDTB-1)/3
! 252: END IF
! 253: IF( LWORK .LT. NB*N ) THEN
! 254: NB = LWORK/N
! 255: END IF
! 256: *
! 257: * Determine the number of the block columns
! 258: *
! 259: NT = (N+NB-1)/NB
! 260: TD = 2*NB
! 261: KB = MIN(NB, N)
! 262: *
! 263: * Initialize vectors/matrices
! 264: *
! 265: DO J = 1, KB
! 266: IPIV( J ) = J
! 267: END DO
! 268: *
! 269: * Save NB
! 270: *
! 271: TB( 1 ) = NB
! 272: *
! 273: IF( UPPER ) THEN
! 274: *
! 275: * .....................................................
! 276: * Factorize A as L*D*L**T using the upper triangle of A
! 277: * .....................................................
! 278: *
! 279: DO J = 0, NT-1
! 280: *
! 281: * Generate Jth column of W and H
! 282: *
! 283: KB = MIN(NB, N-J*NB)
! 284: DO I = 1, J-1
! 285: IF( I.EQ.1 ) THEN
! 286: * H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
! 287: IF( I .EQ. (J-1) ) THEN
! 288: JB = NB+KB
! 289: ELSE
! 290: JB = 2*NB
! 291: END IF
! 292: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 293: $ NB, KB, JB,
! 294: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 295: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
! 296: $ CZERO, WORK( I*NB+1 ), N )
! 297: ELSE
! 298: * 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)
! 299: IF( I .EQ. (J-1) ) THEN
! 300: JB = 2*NB+KB
! 301: ELSE
! 302: JB = 3*NB
! 303: END IF
! 304: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 305: $ NB, KB, JB,
! 306: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 307: $ LDTB-1,
! 308: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
! 309: $ CZERO, WORK( I*NB+1 ), N )
! 310: END IF
! 311: END DO
! 312: *
! 313: * Compute T(J,J)
! 314: *
! 315: CALL ZLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
! 316: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 317: IF( J.GT.1 ) THEN
! 318: * T(J,J) = U(1:J,J)'*H(1:J)
! 319: CALL ZGEMM( 'Transpose', 'NoTranspose',
! 320: $ KB, KB, (J-1)*NB,
! 321: $ -CONE, A( 1, J*NB+1 ), LDA,
! 322: $ WORK( NB+1 ), N,
! 323: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 324: * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
! 325: CALL ZGEMM( 'Transpose', 'NoTranspose',
! 326: $ KB, NB, KB,
! 327: $ CONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
! 328: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 329: $ CZERO, WORK( 1 ), N )
! 330: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 331: $ KB, KB, NB,
! 332: $ -CONE, WORK( 1 ), N,
! 333: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 334: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 335: END IF
! 336: *
! 337: * Expand T(J,J) into full format
! 338: *
! 339: DO I = 1, KB
! 340: DO K = I+1, KB
! 341: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 342: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 343: END DO
! 344: END DO
! 345: IF( J.GT.0 ) THEN
! 346: c CALL CHEGST( 1, 'Upper', KB,
! 347: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 348: c $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
! 349: CALL ZTRSM( 'L', 'U', 'T', 'N', KB, KB, CONE,
! 350: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 351: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 352: CALL ZTRSM( 'R', 'U', 'N', 'N', KB, KB, CONE,
! 353: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 354: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 355: END IF
! 356: *
! 357: IF( J.LT.NT-1 ) THEN
! 358: IF( J.GT.0 ) THEN
! 359: *
! 360: * Compute H(J,J)
! 361: *
! 362: IF( J.EQ.1 ) THEN
! 363: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 364: $ KB, KB, KB,
! 365: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 366: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 367: $ CZERO, WORK( J*NB+1 ), N )
! 368: ELSE
! 369: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 370: $ KB, KB, NB+KB,
! 371: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 372: $ LDTB-1,
! 373: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 374: $ CZERO, WORK( J*NB+1 ), N )
! 375: END IF
! 376: *
! 377: * Update with the previous column
! 378: *
! 379: CALL ZGEMM( 'Transpose', 'NoTranspose',
! 380: $ NB, N-(J+1)*NB, J*NB,
! 381: $ -CONE, WORK( NB+1 ), N,
! 382: $ A( 1, (J+1)*NB+1 ), LDA,
! 383: $ CONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
! 384: END IF
! 385: *
! 386: * Copy panel to workspace to call ZGETRF
! 387: *
! 388: DO K = 1, NB
! 389: CALL ZCOPY( N-(J+1)*NB,
! 390: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
! 391: $ WORK( 1+(K-1)*N ), 1 )
! 392: END DO
! 393: *
! 394: * Factorize panel
! 395: *
! 396: CALL ZGETRF( N-(J+1)*NB, NB,
! 397: $ WORK, N,
! 398: $ IPIV( (J+1)*NB+1 ), IINFO )
! 399: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 400: c INFO = IINFO+(J+1)*NB
! 401: c END IF
! 402: *
! 403: * Copy panel back
! 404: *
! 405: DO K = 1, NB
! 406: CALL ZCOPY( N-(J+1)*NB,
! 407: $ WORK( 1+(K-1)*N ), 1,
! 408: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
! 409: END DO
! 410: *
! 411: * Compute T(J+1, J), zero out for GEMM update
! 412: *
! 413: KB = MIN(NB, N-(J+1)*NB)
! 414: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
! 415: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
! 416: CALL ZLACPY( 'Upper', KB, NB,
! 417: $ WORK, N,
! 418: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 419: IF( J.GT.0 ) THEN
! 420: CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, CONE,
! 421: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 422: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 423: END IF
! 424: *
! 425: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
! 426: * updates
! 427: *
! 428: DO K = 1, NB
! 429: DO I = 1, KB
! 430: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
! 431: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
! 432: END DO
! 433: END DO
! 434: CALL ZLASET( 'Lower', KB, NB, CZERO, CONE,
! 435: $ A( J*NB+1, (J+1)*NB+1), LDA )
! 436: *
! 437: * Apply pivots to trailing submatrix of A
! 438: *
! 439: DO K = 1, KB
! 440: * > Adjust ipiv
! 441: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 442: *
! 443: I1 = (J+1)*NB+K
! 444: I2 = IPIV( (J+1)*NB+K )
! 445: IF( I1.NE.I2 ) THEN
! 446: * > Apply pivots to previous columns of L
! 447: CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
! 448: $ A( (J+1)*NB+1, I2 ), 1 )
! 449: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 450: CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
! 451: $ A( I1+1, I2 ), 1 )
! 452: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 453: CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
! 454: $ A( I2, I2+1 ), LDA )
! 455: * > Swap A(I1, I1) with A(I2, I2)
! 456: PIV = A( I1, I1 )
! 457: A( I1, I1 ) = A( I2, I2 )
! 458: A( I2, I2 ) = PIV
! 459: * > Apply pivots to previous columns of L
! 460: IF( J.GT.0 ) THEN
! 461: CALL ZSWAP( J*NB, A( 1, I1 ), 1,
! 462: $ A( 1, I2 ), 1 )
! 463: END IF
! 464: ENDIF
! 465: END DO
! 466: END IF
! 467: END DO
! 468: ELSE
! 469: *
! 470: * .....................................................
! 471: * Factorize A as L*D*L**T using the lower triangle of A
! 472: * .....................................................
! 473: *
! 474: DO J = 0, NT-1
! 475: *
! 476: * Generate Jth column of W and H
! 477: *
! 478: KB = MIN(NB, N-J*NB)
! 479: DO I = 1, J-1
! 480: IF( I.EQ.1 ) THEN
! 481: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
! 482: IF( I .EQ. (J-1) ) THEN
! 483: JB = NB+KB
! 484: ELSE
! 485: JB = 2*NB
! 486: END IF
! 487: CALL ZGEMM( 'NoTranspose', 'Transpose',
! 488: $ NB, KB, JB,
! 489: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 490: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
! 491: $ CZERO, WORK( I*NB+1 ), N )
! 492: ELSE
! 493: * 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)'
! 494: IF( I .EQ. (J-1) ) THEN
! 495: JB = 2*NB+KB
! 496: ELSE
! 497: JB = 3*NB
! 498: END IF
! 499: CALL ZGEMM( 'NoTranspose', 'Transpose',
! 500: $ NB, KB, JB,
! 501: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 502: $ LDTB-1,
! 503: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
! 504: $ CZERO, WORK( I*NB+1 ), N )
! 505: END IF
! 506: END DO
! 507: *
! 508: * Compute T(J,J)
! 509: *
! 510: CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
! 511: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 512: IF( J.GT.1 ) THEN
! 513: * T(J,J) = L(J,1:J)*H(1:J)
! 514: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 515: $ KB, KB, (J-1)*NB,
! 516: $ -CONE, A( J*NB+1, 1 ), LDA,
! 517: $ WORK( NB+1 ), N,
! 518: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 519: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
! 520: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 521: $ KB, NB, KB,
! 522: $ CONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
! 523: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 524: $ CZERO, WORK( 1 ), N )
! 525: CALL ZGEMM( 'NoTranspose', 'Transpose',
! 526: $ KB, KB, NB,
! 527: $ -CONE, WORK( 1 ), N,
! 528: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 529: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 530: END IF
! 531: *
! 532: * Expand T(J,J) into full format
! 533: *
! 534: DO I = 1, KB
! 535: DO K = I+1, KB
! 536: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 537: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 538: END DO
! 539: END DO
! 540: IF( J.GT.0 ) THEN
! 541: c CALL CHEGST( 1, 'Lower', KB,
! 542: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 543: c $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
! 544: CALL ZTRSM( 'L', 'L', 'N', 'N', KB, KB, CONE,
! 545: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 546: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 547: CALL ZTRSM( 'R', 'L', 'T', 'N', KB, KB, CONE,
! 548: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 549: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 550: END IF
! 551: *
! 552: * Symmetrize T(J,J)
! 553: *
! 554: DO I = 1, KB
! 555: DO K = I+1, KB
! 556: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 557: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 558: END DO
! 559: END DO
! 560: *
! 561: IF( J.LT.NT-1 ) THEN
! 562: IF( J.GT.0 ) THEN
! 563: *
! 564: * Compute H(J,J)
! 565: *
! 566: IF( J.EQ.1 ) THEN
! 567: CALL ZGEMM( 'NoTranspose', 'Transpose',
! 568: $ KB, KB, KB,
! 569: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 570: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 571: $ CZERO, WORK( J*NB+1 ), N )
! 572: ELSE
! 573: CALL ZGEMM( 'NoTranspose', 'Transpose',
! 574: $ KB, KB, NB+KB,
! 575: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 576: $ LDTB-1,
! 577: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 578: $ CZERO, WORK( J*NB+1 ), N )
! 579: END IF
! 580: *
! 581: * Update with the previous column
! 582: *
! 583: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
! 584: $ N-(J+1)*NB, NB, J*NB,
! 585: $ -CONE, A( (J+1)*NB+1, 1 ), LDA,
! 586: $ WORK( NB+1 ), N,
! 587: $ CONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
! 588: END IF
! 589: *
! 590: * Factorize panel
! 591: *
! 592: CALL ZGETRF( N-(J+1)*NB, NB,
! 593: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 594: $ IPIV( (J+1)*NB+1 ), IINFO )
! 595: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 596: c INFO = IINFO+(J+1)*NB
! 597: c END IF
! 598: *
! 599: * Compute T(J+1, J), zero out for GEMM update
! 600: *
! 601: KB = MIN(NB, N-(J+1)*NB)
! 602: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
! 603: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
! 604: CALL ZLACPY( 'Upper', KB, NB,
! 605: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 606: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 607: IF( J.GT.0 ) THEN
! 608: CALL ZTRSM( 'R', 'L', 'T', 'U', KB, NB, CONE,
! 609: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 610: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 611: END IF
! 612: *
! 613: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
! 614: * updates
! 615: *
! 616: DO K = 1, NB
! 617: DO I = 1, KB
! 618: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) =
! 619: $ TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
! 620: END DO
! 621: END DO
! 622: CALL ZLASET( 'Upper', KB, NB, CZERO, CONE,
! 623: $ A( (J+1)*NB+1, J*NB+1 ), LDA )
! 624: *
! 625: * Apply pivots to trailing submatrix of A
! 626: *
! 627: DO K = 1, KB
! 628: * > Adjust ipiv
! 629: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 630: *
! 631: I1 = (J+1)*NB+K
! 632: I2 = IPIV( (J+1)*NB+K )
! 633: IF( I1.NE.I2 ) THEN
! 634: * > Apply pivots to previous columns of L
! 635: CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
! 636: $ A( I2, (J+1)*NB+1 ), LDA )
! 637: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 638: CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
! 639: $ A( I2, I1+1 ), LDA )
! 640: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 641: CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
! 642: $ A( I2+1, I2 ), 1 )
! 643: * > Swap A(I1, I1) with A(I2, I2)
! 644: PIV = A( I1, I1 )
! 645: A( I1, I1 ) = A( I2, I2 )
! 646: A( I2, I2 ) = PIV
! 647: * > Apply pivots to previous columns of L
! 648: IF( J.GT.0 ) THEN
! 649: CALL ZSWAP( J*NB, A( I1, 1 ), LDA,
! 650: $ A( I2, 1 ), LDA )
! 651: END IF
! 652: ENDIF
! 653: END DO
! 654: *
! 655: * Apply pivots to previous columns of L
! 656: *
! 657: c CALL ZLASWP( J*NB, A( 1, 1 ), LDA,
! 658: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
! 659: END IF
! 660: END DO
! 661: END IF
! 662: *
! 663: * Factor the band matrix
! 664: CALL ZGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
! 665: *
! 666: * End of ZSYTRF_AA_2STAGE
! 667: *
! 668: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytrf_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