Annotation of rpl/lapack/lapack/dsytrf_aa_2stage.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DSYTRF_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 DSYTRF_AA_2STAGE + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSYTRF_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: * DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DSYTRF_AA_2STAGE computes the factorization of a real 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 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 DOUBLE PRECISION array, dimension (LDA,N)
! 70: *> On entry, the symmetric 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 DOUBLE PRECISION 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] WORK
! 106: *> \verbatim
! 107: *> WORK is DOUBLE PRECISION workspace of size LWORK
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in] LWORK
! 111: *> \verbatim
! 112: *> The size of WORK. LWORK >= N, internally used to select NB
! 113: *> such that LWORK >= N*NB.
! 114: *>
! 115: *> If LWORK = -1, then a workspace query is assumed; the
! 116: *> routine only calculates the optimal size of the WORK array,
! 117: *> returns this value as the first entry of the WORK array, and
! 118: *> no error message related to LWORK is issued by XERBLA.
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[out] IPIV
! 122: *> \verbatim
! 123: *> IPIV is INTEGER array, dimension (N)
! 124: *> On exit, it contains the details of the interchanges, i.e.,
! 125: *> the row and column k of A were interchanged with the
! 126: *> row and column IPIV(k).
! 127: *> \endverbatim
! 128: *>
! 129: *> \param[out] IPIV2
! 130: *> \verbatim
! 131: *> IPIV is INTEGER array, dimension (N)
! 132: *> On exit, it contains the details of the interchanges, i.e.,
! 133: *> the row and column k of T were interchanged with the
! 134: *> row and column IPIV(k).
! 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 doubleSYcomputational
! 156: *
! 157: * =====================================================================
! 158: SUBROUTINE DSYTRF_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: DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
! 175: * ..
! 176: *
! 177: * =====================================================================
! 178: * .. Parameters ..
! 179: DOUBLE PRECISION ZERO, ONE
! 180: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 181: *
! 182: * .. Local Scalars ..
! 183: LOGICAL UPPER, TQUERY, WQUERY
! 184: INTEGER I, J, K, I1, I2, TD
! 185: INTEGER LDTB, NB, KB, JB, NT, IINFO
! 186: DOUBLE PRECISION PIV
! 187: * ..
! 188: * .. External Functions ..
! 189: LOGICAL LSAME
! 190: INTEGER ILAENV
! 191: EXTERNAL LSAME, ILAENV
! 192: * ..
! 193: * .. External Subroutines ..
! 194: EXTERNAL XERBLA, DCOPY, DLACGV, DLACPY,
! 195: $ DLASET, DGBTRF, DGEMM, DGETRF,
! 196: $ DSYGST, DSWAP, DTRSM
! 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( 'DSYTRF_AA_2STAGE', -INFO )
! 223: RETURN
! 224: END IF
! 225: *
! 226: * Answer the query
! 227: *
! 228: NB = ILAENV( 1, 'DSYTRF_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,I+1)*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 DGEMM( 'NoTranspose', 'NoTranspose',
! 293: $ NB, KB, JB,
! 294: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 295: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
! 296: $ ZERO, 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 DGEMM( 'NoTranspose', 'NoTranspose',
! 305: $ NB, KB, JB,
! 306: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 307: $ LDTB-1,
! 308: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
! 309: $ ZERO, WORK( I*NB+1 ), N )
! 310: END IF
! 311: END DO
! 312: *
! 313: * Compute T(J,J)
! 314: *
! 315: CALL DLACPY( '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 DGEMM( 'Transpose', 'NoTranspose',
! 320: $ KB, KB, (J-1)*NB,
! 321: $ -ONE, A( 1, J*NB+1 ), LDA,
! 322: $ WORK( NB+1 ), N,
! 323: $ ONE, 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 DGEMM( 'Transpose', 'NoTranspose',
! 326: $ KB, NB, KB,
! 327: $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
! 328: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 329: $ ZERO, WORK( 1 ), N )
! 330: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 331: $ KB, KB, NB,
! 332: $ -ONE, WORK( 1 ), N,
! 333: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 334: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 335: END IF
! 336: IF( J.GT.0 ) THEN
! 337: CALL DSYGST( 1, 'Upper', KB,
! 338: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 339: $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
! 340: END IF
! 341: *
! 342: * Expand T(J,J) into full format
! 343: *
! 344: DO I = 1, KB
! 345: DO K = I+1, KB
! 346: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 347: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 348: END DO
! 349: END DO
! 350: *
! 351: IF( J.LT.NT-1 ) THEN
! 352: IF( J.GT.0 ) THEN
! 353: *
! 354: * Compute H(J,J)
! 355: *
! 356: IF( J.EQ.1 ) THEN
! 357: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 358: $ KB, KB, KB,
! 359: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 360: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 361: $ ZERO, WORK( J*NB+1 ), N )
! 362: ELSE
! 363: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 364: $ KB, KB, NB+KB,
! 365: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 366: $ LDTB-1,
! 367: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
! 368: $ ZERO, WORK( J*NB+1 ), N )
! 369: END IF
! 370: *
! 371: * Update with the previous column
! 372: *
! 373: CALL DGEMM( 'Transpose', 'NoTranspose',
! 374: $ NB, N-(J+1)*NB, J*NB,
! 375: $ -ONE, WORK( NB+1 ), N,
! 376: $ A( 1, (J+1)*NB+1 ), LDA,
! 377: $ ONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
! 378: END IF
! 379: *
! 380: * Copy panel to workspace to call DGETRF
! 381: *
! 382: DO K = 1, NB
! 383: CALL DCOPY( N-(J+1)*NB,
! 384: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
! 385: $ WORK( 1+(K-1)*N ), 1 )
! 386: END DO
! 387: *
! 388: * Factorize panel
! 389: *
! 390: CALL DGETRF( N-(J+1)*NB, NB,
! 391: $ WORK, N,
! 392: $ IPIV( (J+1)*NB+1 ), IINFO )
! 393: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 394: c INFO = IINFO+(J+1)*NB
! 395: c END IF
! 396: *
! 397: * Copy panel back
! 398: *
! 399: DO K = 1, NB
! 400: CALL DCOPY( N-(J+1)*NB,
! 401: $ WORK( 1+(K-1)*N ), 1,
! 402: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
! 403: END DO
! 404: *
! 405: * Compute T(J+1, J), zero out for GEMM update
! 406: *
! 407: KB = MIN(NB, N-(J+1)*NB)
! 408: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
! 409: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
! 410: CALL DLACPY( 'Upper', KB, NB,
! 411: $ WORK, N,
! 412: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 413: IF( J.GT.0 ) THEN
! 414: CALL DTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE,
! 415: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
! 416: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 417: END IF
! 418: *
! 419: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
! 420: * updates
! 421: *
! 422: DO K = 1, NB
! 423: DO I = 1, KB
! 424: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
! 425: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
! 426: END DO
! 427: END DO
! 428: CALL DLASET( 'Lower', KB, NB, ZERO, ONE,
! 429: $ A( J*NB+1, (J+1)*NB+1), LDA )
! 430: *
! 431: * Apply pivots to trailing submatrix of A
! 432: *
! 433: DO K = 1, KB
! 434: * > Adjust ipiv
! 435: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 436: *
! 437: I1 = (J+1)*NB+K
! 438: I2 = IPIV( (J+1)*NB+K )
! 439: IF( I1.NE.I2 ) THEN
! 440: * > Apply pivots to previous columns of L
! 441: CALL DSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
! 442: $ A( (J+1)*NB+1, I2 ), 1 )
! 443: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 444: CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
! 445: $ A( I1+1, I2 ), 1 )
! 446: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 447: CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
! 448: $ A( I2, I2+1 ), LDA )
! 449: * > Swap A(I1, I1) with A(I2, I2)
! 450: PIV = A( I1, I1 )
! 451: A( I1, I1 ) = A( I2, I2 )
! 452: A( I2, I2 ) = PIV
! 453: * > Apply pivots to previous columns of L
! 454: IF( J.GT.0 ) THEN
! 455: CALL DSWAP( J*NB, A( 1, I1 ), 1,
! 456: $ A( 1, I2 ), 1 )
! 457: END IF
! 458: ENDIF
! 459: END DO
! 460: END IF
! 461: END DO
! 462: ELSE
! 463: *
! 464: * .....................................................
! 465: * Factorize A as L*D*L**T using the lower triangle of A
! 466: * .....................................................
! 467: *
! 468: DO J = 0, NT-1
! 469: *
! 470: * Generate Jth column of W and H
! 471: *
! 472: KB = MIN(NB, N-J*NB)
! 473: DO I = 1, J-1
! 474: IF( I.EQ.1 ) THEN
! 475: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
! 476: IF( I .EQ. J-1) THEN
! 477: JB = NB+KB
! 478: ELSE
! 479: JB = 2*NB
! 480: END IF
! 481: CALL DGEMM( 'NoTranspose', 'Transpose',
! 482: $ NB, KB, JB,
! 483: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
! 484: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
! 485: $ ZERO, WORK( I*NB+1 ), N )
! 486: ELSE
! 487: * 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)'
! 488: IF( I .EQ. J-1) THEN
! 489: JB = 2*NB+KB
! 490: ELSE
! 491: JB = 3*NB
! 492: END IF
! 493: CALL DGEMM( 'NoTranspose', 'Transpose',
! 494: $ NB, KB, JB,
! 495: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
! 496: $ LDTB-1,
! 497: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
! 498: $ ZERO, WORK( I*NB+1 ), N )
! 499: END IF
! 500: END DO
! 501: *
! 502: * Compute T(J,J)
! 503: *
! 504: CALL DLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
! 505: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 506: IF( J.GT.1 ) THEN
! 507: * T(J,J) = L(J,1:J)*H(1:J)
! 508: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 509: $ KB, KB, (J-1)*NB,
! 510: $ -ONE, A( J*NB+1, 1 ), LDA,
! 511: $ WORK( NB+1 ), N,
! 512: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 513: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
! 514: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 515: $ KB, NB, KB,
! 516: $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
! 517: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
! 518: $ ZERO, WORK( 1 ), N )
! 519: CALL DGEMM( 'NoTranspose', 'Transpose',
! 520: $ KB, KB, NB,
! 521: $ -ONE, WORK( 1 ), N,
! 522: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 523: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
! 524: END IF
! 525: IF( J.GT.0 ) THEN
! 526: CALL DSYGST( 1, 'Lower', KB,
! 527: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 528: $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
! 529: END IF
! 530: *
! 531: * Expand T(J,J) into full format
! 532: *
! 533: DO I = 1, KB
! 534: DO K = I+1, KB
! 535: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
! 536: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
! 537: END DO
! 538: END DO
! 539: *
! 540: IF( J.LT.NT-1 ) THEN
! 541: IF( J.GT.0 ) THEN
! 542: *
! 543: * Compute H(J,J)
! 544: *
! 545: IF( J.EQ.1 ) THEN
! 546: CALL DGEMM( 'NoTranspose', 'Transpose',
! 547: $ KB, KB, KB,
! 548: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
! 549: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 550: $ ZERO, WORK( J*NB+1 ), N )
! 551: ELSE
! 552: CALL DGEMM( 'NoTranspose', 'Transpose',
! 553: $ KB, KB, NB+KB,
! 554: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
! 555: $ LDTB-1,
! 556: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
! 557: $ ZERO, WORK( J*NB+1 ), N )
! 558: END IF
! 559: *
! 560: * Update with the previous column
! 561: *
! 562: CALL DGEMM( 'NoTranspose', 'NoTranspose',
! 563: $ N-(J+1)*NB, NB, J*NB,
! 564: $ -ONE, A( (J+1)*NB+1, 1 ), LDA,
! 565: $ WORK( NB+1 ), N,
! 566: $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
! 567: END IF
! 568: *
! 569: * Factorize panel
! 570: *
! 571: CALL DGETRF( N-(J+1)*NB, NB,
! 572: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 573: $ IPIV( (J+1)*NB+1 ), IINFO )
! 574: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
! 575: c INFO = IINFO+(J+1)*NB
! 576: c END IF
! 577: *
! 578: * Compute T(J+1, J), zero out for GEMM update
! 579: *
! 580: KB = MIN(NB, N-(J+1)*NB)
! 581: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
! 582: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
! 583: CALL DLACPY( 'Upper', KB, NB,
! 584: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
! 585: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 586: IF( J.GT.0 ) THEN
! 587: CALL DTRSM( 'R', 'L', 'T', 'U', KB, NB, ONE,
! 588: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
! 589: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
! 590: END IF
! 591: *
! 592: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
! 593: * updates
! 594: *
! 595: DO K = 1, NB
! 596: DO I = 1, KB
! 597: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
! 598: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
! 599: END DO
! 600: END DO
! 601: CALL DLASET( 'Upper', KB, NB, ZERO, ONE,
! 602: $ A( (J+1)*NB+1, J*NB+1), LDA )
! 603: *
! 604: * Apply pivots to trailing submatrix of A
! 605: *
! 606: DO K = 1, KB
! 607: * > Adjust ipiv
! 608: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
! 609: *
! 610: I1 = (J+1)*NB+K
! 611: I2 = IPIV( (J+1)*NB+K )
! 612: IF( I1.NE.I2 ) THEN
! 613: * > Apply pivots to previous columns of L
! 614: CALL DSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
! 615: $ A( I2, (J+1)*NB+1 ), LDA )
! 616: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 617: CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
! 618: $ A( I2, I1+1 ), LDA )
! 619: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
! 620: CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
! 621: $ A( I2+1, I2 ), 1 )
! 622: * > Swap A(I1, I1) with A(I2, I2)
! 623: PIV = A( I1, I1 )
! 624: A( I1, I1 ) = A( I2, I2 )
! 625: A( I2, I2 ) = PIV
! 626: * > Apply pivots to previous columns of L
! 627: IF( J.GT.0 ) THEN
! 628: CALL DSWAP( J*NB, A( I1, 1 ), LDA,
! 629: $ A( I2, 1 ), LDA )
! 630: END IF
! 631: ENDIF
! 632: END DO
! 633: *
! 634: * Apply pivots to previous columns of L
! 635: *
! 636: c CALL DLASWP( J*NB, A( 1, 1 ), LDA,
! 637: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
! 638: END IF
! 639: END DO
! 640: END IF
! 641: *
! 642: * Factor the band matrix
! 643: CALL DGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
! 644: *
! 645: * End of DSYTRF_AA_2STAGE
! 646: *
! 647: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrf_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