Annotation of rpl/lapack/lapack/dsytri_3x.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DSYTRI_3X
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSYTRI_3X + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri_3x.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri_3x.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri_3x.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSYTRI_3X( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * CHARACTER UPLO
! 25: * INTEGER INFO, LDA, N, NB
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IPIV( * )
! 29: * DOUBLE PRECISION A( LDA, * ), E( * ), WORK( N+NB+1, * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *> DSYTRI_3X computes the inverse of a real symmetric indefinite
! 38: *> matrix A using the factorization computed by DSYTRF_RK or DSYTRF_BK:
! 39: *>
! 40: *> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
! 41: *>
! 42: *> where U (or L) is unit upper (or lower) triangular matrix,
! 43: *> U**T (or L**T) is the transpose of U (or L), P is a permutation
! 44: *> matrix, P**T is the transpose of P, and D is symmetric and block
! 45: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
! 46: *>
! 47: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
! 48: *> \endverbatim
! 49: *
! 50: * Arguments:
! 51: * ==========
! 52: *
! 53: *> \param[in] UPLO
! 54: *> \verbatim
! 55: *> UPLO is CHARACTER*1
! 56: *> Specifies whether the details of the factorization are
! 57: *> stored as an upper or lower triangular matrix.
! 58: *> = 'U': Upper triangle of A is stored;
! 59: *> = 'L': Lower triangle of A is stored.
! 60: *> \endverbatim
! 61: *>
! 62: *> \param[in] N
! 63: *> \verbatim
! 64: *> N is INTEGER
! 65: *> The order of the matrix A. N >= 0.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in,out] A
! 69: *> \verbatim
! 70: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 71: *> On entry, diagonal of the block diagonal matrix D and
! 72: *> factors U or L as computed by DSYTRF_RK and DSYTRF_BK:
! 73: *> a) ONLY diagonal elements of the symmetric block diagonal
! 74: *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
! 75: *> (superdiagonal (or subdiagonal) elements of D
! 76: *> should be provided on entry in array E), and
! 77: *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
! 78: *> If UPLO = 'L': factor L in the subdiagonal part of A.
! 79: *>
! 80: *> On exit, if INFO = 0, the symmetric inverse of the original
! 81: *> matrix.
! 82: *> If UPLO = 'U': the upper triangular part of the inverse
! 83: *> is formed and the part of A below the diagonal is not
! 84: *> referenced;
! 85: *> If UPLO = 'L': the lower triangular part of the inverse
! 86: *> is formed and the part of A above the diagonal is not
! 87: *> referenced.
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[in] LDA
! 91: *> \verbatim
! 92: *> LDA is INTEGER
! 93: *> The leading dimension of the array A. LDA >= max(1,N).
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] E
! 97: *> \verbatim
! 98: *> E is DOUBLE PRECISION array, dimension (N)
! 99: *> On entry, contains the superdiagonal (or subdiagonal)
! 100: *> elements of the symmetric block diagonal matrix D
! 101: *> with 1-by-1 or 2-by-2 diagonal blocks, where
! 102: *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) not refernced;
! 103: *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) not referenced.
! 104: *>
! 105: *> NOTE: For 1-by-1 diagonal block D(k), where
! 106: *> 1 <= k <= N, the element E(k) is not referenced in both
! 107: *> UPLO = 'U' or UPLO = 'L' cases.
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in] IPIV
! 111: *> \verbatim
! 112: *> IPIV is INTEGER array, dimension (N)
! 113: *> Details of the interchanges and the block structure of D
! 114: *> as determined by DSYTRF_RK or DSYTRF_BK.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[out] WORK
! 118: *> \verbatim
! 119: *> WORK is DOUBLE PRECISION array, dimension (N+NB+1,NB+3).
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[in] NB
! 123: *> \verbatim
! 124: *> NB is INTEGER
! 125: *> Block size.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[out] INFO
! 129: *> \verbatim
! 130: *> INFO is INTEGER
! 131: *> = 0: successful exit
! 132: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 133: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
! 134: *> inverse could not be computed.
! 135: *> \endverbatim
! 136: *
! 137: * Authors:
! 138: * ========
! 139: *
! 140: *> \author Univ. of Tennessee
! 141: *> \author Univ. of California Berkeley
! 142: *> \author Univ. of Colorado Denver
! 143: *> \author NAG Ltd.
! 144: *
! 145: *> \date December 2016
! 146: *
! 147: *> \ingroup doubleSYcomputational
! 148: *
! 149: *> \par Contributors:
! 150: * ==================
! 151: *> \verbatim
! 152: *>
! 153: *> December 2016, Igor Kozachenko,
! 154: *> Computer Science Division,
! 155: *> University of California, Berkeley
! 156: *>
! 157: *> \endverbatim
! 158: *
! 159: * =====================================================================
! 160: SUBROUTINE DSYTRI_3X( UPLO, N, A, LDA, E, IPIV, WORK, NB, INFO )
! 161: *
! 162: * -- LAPACK computational routine (version 3.7.0) --
! 163: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 164: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 165: * December 2016
! 166: *
! 167: * .. Scalar Arguments ..
! 168: CHARACTER UPLO
! 169: INTEGER INFO, LDA, N, NB
! 170: * ..
! 171: * .. Array Arguments ..
! 172: INTEGER IPIV( * )
! 173: DOUBLE PRECISION A( LDA, * ), E( * ), WORK( N+NB+1, * )
! 174: * ..
! 175: *
! 176: * =====================================================================
! 177: *
! 178: * .. Parameters ..
! 179: DOUBLE PRECISION ONE, ZERO
! 180: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 181: * ..
! 182: * .. Local Scalars ..
! 183: LOGICAL UPPER
! 184: INTEGER CUT, I, ICOUNT, INVD, IP, K, NNB, J, U11
! 185: DOUBLE PRECISION AK, AKKP1, AKP1, D, T, U01_I_J, U01_IP1_J,
! 186: $ U11_I_J, U11_IP1_J
! 187: * ..
! 188: * .. External Functions ..
! 189: LOGICAL LSAME
! 190: EXTERNAL LSAME
! 191: * ..
! 192: * .. External Subroutines ..
! 193: EXTERNAL DGEMM, DSYSWAPR, DTRTRI, DTRMM, XERBLA
! 194: * ..
! 195: * .. Intrinsic Functions ..
! 196: INTRINSIC ABS, MAX, MOD
! 197: * ..
! 198: * .. Executable Statements ..
! 199: *
! 200: * Test the input parameters.
! 201: *
! 202: INFO = 0
! 203: UPPER = LSAME( UPLO, 'U' )
! 204: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 205: INFO = -1
! 206: ELSE IF( N.LT.0 ) THEN
! 207: INFO = -2
! 208: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 209: INFO = -4
! 210: END IF
! 211: *
! 212: * Quick return if possible
! 213: *
! 214: IF( INFO.NE.0 ) THEN
! 215: CALL XERBLA( 'DSYTRI_3X', -INFO )
! 216: RETURN
! 217: END IF
! 218: IF( N.EQ.0 )
! 219: $ RETURN
! 220: *
! 221: * Workspace got Non-diag elements of D
! 222: *
! 223: DO K = 1, N
! 224: WORK( K, 1 ) = E( K )
! 225: END DO
! 226: *
! 227: * Check that the diagonal matrix D is nonsingular.
! 228: *
! 229: IF( UPPER ) THEN
! 230: *
! 231: * Upper triangular storage: examine D from bottom to top
! 232: *
! 233: DO INFO = N, 1, -1
! 234: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
! 235: $ RETURN
! 236: END DO
! 237: ELSE
! 238: *
! 239: * Lower triangular storage: examine D from top to bottom.
! 240: *
! 241: DO INFO = 1, N
! 242: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
! 243: $ RETURN
! 244: END DO
! 245: END IF
! 246: *
! 247: INFO = 0
! 248: *
! 249: * Splitting Workspace
! 250: * U01 is a block ( N, NB+1 )
! 251: * The first element of U01 is in WORK( 1, 1 )
! 252: * U11 is a block ( NB+1, NB+1 )
! 253: * The first element of U11 is in WORK( N+1, 1 )
! 254: *
! 255: U11 = N
! 256: *
! 257: * INVD is a block ( N, 2 )
! 258: * The first element of INVD is in WORK( 1, INVD )
! 259: *
! 260: INVD = NB + 2
! 261:
! 262: IF( UPPER ) THEN
! 263: *
! 264: * Begin Upper
! 265: *
! 266: * invA = P * inv(U**T) * inv(D) * inv(U) * P**T.
! 267: *
! 268: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
! 269: *
! 270: * inv(D) and inv(D) * inv(U)
! 271: *
! 272: K = 1
! 273: DO WHILE( K.LE.N )
! 274: IF( IPIV( K ).GT.0 ) THEN
! 275: * 1 x 1 diagonal NNB
! 276: WORK( K, INVD ) = ONE / A( K, K )
! 277: WORK( K, INVD+1 ) = ZERO
! 278: ELSE
! 279: * 2 x 2 diagonal NNB
! 280: T = WORK( K+1, 1 )
! 281: AK = A( K, K ) / T
! 282: AKP1 = A( K+1, K+1 ) / T
! 283: AKKP1 = WORK( K+1, 1 ) / T
! 284: D = T*( AK*AKP1-ONE )
! 285: WORK( K, INVD ) = AKP1 / D
! 286: WORK( K+1, INVD+1 ) = AK / D
! 287: WORK( K, INVD+1 ) = -AKKP1 / D
! 288: WORK( K+1, INVD ) = WORK( K, INVD+1 )
! 289: K = K + 1
! 290: END IF
! 291: K = K + 1
! 292: END DO
! 293: *
! 294: * inv(U**T) = (inv(U))**T
! 295: *
! 296: * inv(U**T) * inv(D) * inv(U)
! 297: *
! 298: CUT = N
! 299: DO WHILE( CUT.GT.0 )
! 300: NNB = NB
! 301: IF( CUT.LE.NNB ) THEN
! 302: NNB = CUT
! 303: ELSE
! 304: ICOUNT = 0
! 305: * count negative elements,
! 306: DO I = CUT+1-NNB, CUT
! 307: IF( IPIV( I ).LT.0 ) ICOUNT = ICOUNT + 1
! 308: END DO
! 309: * need a even number for a clear cut
! 310: IF( MOD( ICOUNT, 2 ).EQ.1 ) NNB = NNB + 1
! 311: END IF
! 312:
! 313: CUT = CUT - NNB
! 314: *
! 315: * U01 Block
! 316: *
! 317: DO I = 1, CUT
! 318: DO J = 1, NNB
! 319: WORK( I, J ) = A( I, CUT+J )
! 320: END DO
! 321: END DO
! 322: *
! 323: * U11 Block
! 324: *
! 325: DO I = 1, NNB
! 326: WORK( U11+I, I ) = ONE
! 327: DO J = 1, I-1
! 328: WORK( U11+I, J ) = ZERO
! 329: END DO
! 330: DO J = I+1, NNB
! 331: WORK( U11+I, J ) = A( CUT+I, CUT+J )
! 332: END DO
! 333: END DO
! 334: *
! 335: * invD * U01
! 336: *
! 337: I = 1
! 338: DO WHILE( I.LE.CUT )
! 339: IF( IPIV( I ).GT.0 ) THEN
! 340: DO J = 1, NNB
! 341: WORK( I, J ) = WORK( I, INVD ) * WORK( I, J )
! 342: END DO
! 343: ELSE
! 344: DO J = 1, NNB
! 345: U01_I_J = WORK( I, J )
! 346: U01_IP1_J = WORK( I+1, J )
! 347: WORK( I, J ) = WORK( I, INVD ) * U01_I_J
! 348: $ + WORK( I, INVD+1 ) * U01_IP1_J
! 349: WORK( I+1, J ) = WORK( I+1, INVD ) * U01_I_J
! 350: $ + WORK( I+1, INVD+1 ) * U01_IP1_J
! 351: END DO
! 352: I = I + 1
! 353: END IF
! 354: I = I + 1
! 355: END DO
! 356: *
! 357: * invD1 * U11
! 358: *
! 359: I = 1
! 360: DO WHILE ( I.LE.NNB )
! 361: IF( IPIV( CUT+I ).GT.0 ) THEN
! 362: DO J = I, NNB
! 363: WORK( U11+I, J ) = WORK(CUT+I,INVD) * WORK(U11+I,J)
! 364: END DO
! 365: ELSE
! 366: DO J = I, NNB
! 367: U11_I_J = WORK(U11+I,J)
! 368: U11_IP1_J = WORK(U11+I+1,J)
! 369: WORK( U11+I, J ) = WORK(CUT+I,INVD) * WORK(U11+I,J)
! 370: $ + WORK(CUT+I,INVD+1) * WORK(U11+I+1,J)
! 371: WORK( U11+I+1, J ) = WORK(CUT+I+1,INVD) * U11_I_J
! 372: $ + WORK(CUT+I+1,INVD+1) * U11_IP1_J
! 373: END DO
! 374: I = I + 1
! 375: END IF
! 376: I = I + 1
! 377: END DO
! 378: *
! 379: * U11**T * invD1 * U11 -> U11
! 380: *
! 381: CALL DTRMM( 'L', 'U', 'T', 'U', NNB, NNB,
! 382: $ ONE, A( CUT+1, CUT+1 ), LDA, WORK( U11+1, 1 ),
! 383: $ N+NB+1 )
! 384: *
! 385: DO I = 1, NNB
! 386: DO J = I, NNB
! 387: A( CUT+I, CUT+J ) = WORK( U11+I, J )
! 388: END DO
! 389: END DO
! 390: *
! 391: * U01**T * invD * U01 -> A( CUT+I, CUT+J )
! 392: *
! 393: CALL DGEMM( 'T', 'N', NNB, NNB, CUT, ONE, A( 1, CUT+1 ),
! 394: $ LDA, WORK, N+NB+1, ZERO, WORK(U11+1,1), N+NB+1 )
! 395:
! 396: *
! 397: * U11 = U11**T * invD1 * U11 + U01**T * invD * U01
! 398: *
! 399: DO I = 1, NNB
! 400: DO J = I, NNB
! 401: A( CUT+I, CUT+J ) = A( CUT+I, CUT+J ) + WORK(U11+I,J)
! 402: END DO
! 403: END DO
! 404: *
! 405: * U01 = U00**T * invD0 * U01
! 406: *
! 407: CALL DTRMM( 'L', UPLO, 'T', 'U', CUT, NNB,
! 408: $ ONE, A, LDA, WORK, N+NB+1 )
! 409:
! 410: *
! 411: * Update U01
! 412: *
! 413: DO I = 1, CUT
! 414: DO J = 1, NNB
! 415: A( I, CUT+J ) = WORK( I, J )
! 416: END DO
! 417: END DO
! 418: *
! 419: * Next Block
! 420: *
! 421: END DO
! 422: *
! 423: * Apply PERMUTATIONS P and P**T:
! 424: * P * inv(U**T) * inv(D) * inv(U) * P**T.
! 425: * Interchange rows and columns I and IPIV(I) in reverse order
! 426: * from the formation order of IPIV vector for Upper case.
! 427: *
! 428: * ( We can use a loop over IPIV with increment 1,
! 429: * since the ABS value of IPIV(I) represents the row (column)
! 430: * index of the interchange with row (column) i in both 1x1
! 431: * and 2x2 pivot cases, i.e. we don't need separate code branches
! 432: * for 1x1 and 2x2 pivot cases )
! 433: *
! 434: DO I = 1, N
! 435: IP = ABS( IPIV( I ) )
! 436: IF( IP.NE.I ) THEN
! 437: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
! 438: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
! 439: END IF
! 440: END DO
! 441: *
! 442: ELSE
! 443: *
! 444: * Begin Lower
! 445: *
! 446: * inv A = P * inv(L**T) * inv(D) * inv(L) * P**T.
! 447: *
! 448: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
! 449: *
! 450: * inv(D) and inv(D) * inv(L)
! 451: *
! 452: K = N
! 453: DO WHILE ( K .GE. 1 )
! 454: IF( IPIV( K ).GT.0 ) THEN
! 455: * 1 x 1 diagonal NNB
! 456: WORK( K, INVD ) = ONE / A( K, K )
! 457: WORK( K, INVD+1 ) = ZERO
! 458: ELSE
! 459: * 2 x 2 diagonal NNB
! 460: T = WORK( K-1, 1 )
! 461: AK = A( K-1, K-1 ) / T
! 462: AKP1 = A( K, K ) / T
! 463: AKKP1 = WORK( K-1, 1 ) / T
! 464: D = T*( AK*AKP1-ONE )
! 465: WORK( K-1, INVD ) = AKP1 / D
! 466: WORK( K, INVD ) = AK / D
! 467: WORK( K, INVD+1 ) = -AKKP1 / D
! 468: WORK( K-1, INVD+1 ) = WORK( K, INVD+1 )
! 469: K = K - 1
! 470: END IF
! 471: K = K - 1
! 472: END DO
! 473: *
! 474: * inv(L**T) = (inv(L))**T
! 475: *
! 476: * inv(L**T) * inv(D) * inv(L)
! 477: *
! 478: CUT = 0
! 479: DO WHILE( CUT.LT.N )
! 480: NNB = NB
! 481: IF( (CUT + NNB).GT.N ) THEN
! 482: NNB = N - CUT
! 483: ELSE
! 484: ICOUNT = 0
! 485: * count negative elements,
! 486: DO I = CUT + 1, CUT+NNB
! 487: IF ( IPIV( I ).LT.0 ) ICOUNT = ICOUNT + 1
! 488: END DO
! 489: * need a even number for a clear cut
! 490: IF( MOD( ICOUNT, 2 ).EQ.1 ) NNB = NNB + 1
! 491: END IF
! 492: *
! 493: * L21 Block
! 494: *
! 495: DO I = 1, N-CUT-NNB
! 496: DO J = 1, NNB
! 497: WORK( I, J ) = A( CUT+NNB+I, CUT+J )
! 498: END DO
! 499: END DO
! 500: *
! 501: * L11 Block
! 502: *
! 503: DO I = 1, NNB
! 504: WORK( U11+I, I) = ONE
! 505: DO J = I+1, NNB
! 506: WORK( U11+I, J ) = ZERO
! 507: END DO
! 508: DO J = 1, I-1
! 509: WORK( U11+I, J ) = A( CUT+I, CUT+J )
! 510: END DO
! 511: END DO
! 512: *
! 513: * invD*L21
! 514: *
! 515: I = N-CUT-NNB
! 516: DO WHILE( I.GE.1 )
! 517: IF( IPIV( CUT+NNB+I ).GT.0 ) THEN
! 518: DO J = 1, NNB
! 519: WORK( I, J ) = WORK( CUT+NNB+I, INVD) * WORK( I, J)
! 520: END DO
! 521: ELSE
! 522: DO J = 1, NNB
! 523: U01_I_J = WORK(I,J)
! 524: U01_IP1_J = WORK(I-1,J)
! 525: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
! 526: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
! 527: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
! 528: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
! 529: END DO
! 530: I = I - 1
! 531: END IF
! 532: I = I - 1
! 533: END DO
! 534: *
! 535: * invD1*L11
! 536: *
! 537: I = NNB
! 538: DO WHILE( I.GE.1 )
! 539: IF( IPIV( CUT+I ).GT.0 ) THEN
! 540: DO J = 1, NNB
! 541: WORK( U11+I, J ) = WORK( CUT+I, INVD)*WORK(U11+I,J)
! 542: END DO
! 543:
! 544: ELSE
! 545: DO J = 1, NNB
! 546: U11_I_J = WORK( U11+I, J )
! 547: U11_IP1_J = WORK( U11+I-1, J )
! 548: WORK( U11+I, J ) = WORK(CUT+I,INVD) * WORK(U11+I,J)
! 549: $ + WORK(CUT+I,INVD+1) * U11_IP1_J
! 550: WORK( U11+I-1, J ) = WORK(CUT+I-1,INVD+1) * U11_I_J
! 551: $ + WORK(CUT+I-1,INVD) * U11_IP1_J
! 552: END DO
! 553: I = I - 1
! 554: END IF
! 555: I = I - 1
! 556: END DO
! 557: *
! 558: * L11**T * invD1 * L11 -> L11
! 559: *
! 560: CALL DTRMM( 'L', UPLO, 'T', 'U', NNB, NNB, ONE,
! 561: $ A( CUT+1, CUT+1 ), LDA, WORK( U11+1, 1 ),
! 562: $ N+NB+1 )
! 563:
! 564: *
! 565: DO I = 1, NNB
! 566: DO J = 1, I
! 567: A( CUT+I, CUT+J ) = WORK( U11+I, J )
! 568: END DO
! 569: END DO
! 570: *
! 571: IF( (CUT+NNB).LT.N ) THEN
! 572: *
! 573: * L21**T * invD2*L21 -> A( CUT+I, CUT+J )
! 574: *
! 575: CALL DGEMM( 'T', 'N', NNB, NNB, N-NNB-CUT, ONE,
! 576: $ A( CUT+NNB+1, CUT+1 ), LDA, WORK, N+NB+1,
! 577: $ ZERO, WORK( U11+1, 1 ), N+NB+1 )
! 578:
! 579: *
! 580: * L11 = L11**T * invD1 * L11 + U01**T * invD * U01
! 581: *
! 582: DO I = 1, NNB
! 583: DO J = 1, I
! 584: A( CUT+I, CUT+J ) = A( CUT+I, CUT+J )+WORK(U11+I,J)
! 585: END DO
! 586: END DO
! 587: *
! 588: * L01 = L22**T * invD2 * L21
! 589: *
! 590: CALL DTRMM( 'L', UPLO, 'T', 'U', N-NNB-CUT, NNB, ONE,
! 591: $ A( CUT+NNB+1, CUT+NNB+1 ), LDA, WORK,
! 592: $ N+NB+1 )
! 593: *
! 594: * Update L21
! 595: *
! 596: DO I = 1, N-CUT-NNB
! 597: DO J = 1, NNB
! 598: A( CUT+NNB+I, CUT+J ) = WORK( I, J )
! 599: END DO
! 600: END DO
! 601: *
! 602: ELSE
! 603: *
! 604: * L11 = L11**T * invD1 * L11
! 605: *
! 606: DO I = 1, NNB
! 607: DO J = 1, I
! 608: A( CUT+I, CUT+J ) = WORK( U11+I, J )
! 609: END DO
! 610: END DO
! 611: END IF
! 612: *
! 613: * Next Block
! 614: *
! 615: CUT = CUT + NNB
! 616: *
! 617: END DO
! 618: *
! 619: * Apply PERMUTATIONS P and P**T:
! 620: * P * inv(L**T) * inv(D) * inv(L) * P**T.
! 621: * Interchange rows and columns I and IPIV(I) in reverse order
! 622: * from the formation order of IPIV vector for Lower case.
! 623: *
! 624: * ( We can use a loop over IPIV with increment -1,
! 625: * since the ABS value of IPIV(I) represents the row (column)
! 626: * index of the interchange with row (column) i in both 1x1
! 627: * and 2x2 pivot cases, i.e. we don't need separate code branches
! 628: * for 1x1 and 2x2 pivot cases )
! 629: *
! 630: DO I = N, 1, -1
! 631: IP = ABS( IPIV( I ) )
! 632: IF( IP.NE.I ) THEN
! 633: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
! 634: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
! 635: END IF
! 636: END DO
! 637: *
! 638: END IF
! 639: *
! 640: RETURN
! 641: *
! 642: * End of DSYTRI_3X
! 643: *
! 644: END
! 645:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytri_3x.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack