Annotation of rpl/lapack/lapack/dsytri2x.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.3.0) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2010
! 7: *
! 8: * -- Written by Julie Langou of the Univ. of TN --
! 9: *
! 10: * .. Scalar Arguments ..
! 11: CHARACTER UPLO
! 12: INTEGER INFO, LDA, N, NB
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER IPIV( * )
! 16: DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * DSYTRI2X computes the inverse of a real symmetric indefinite matrix
! 23: * A using the factorization A = U*D*U**T or A = L*D*L**T computed by
! 24: * DSYTRF.
! 25: *
! 26: * Arguments
! 27: * =========
! 28: *
! 29: * UPLO (input) CHARACTER*1
! 30: * Specifies whether the details of the factorization are stored
! 31: * as an upper or lower triangular matrix.
! 32: * = 'U': Upper triangular, form is A = U*D*U**T;
! 33: * = 'L': Lower triangular, form is A = L*D*L**T.
! 34: *
! 35: * N (input) INTEGER
! 36: * The order of the matrix A. N >= 0.
! 37: *
! 38: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! 39: * On entry, the NNB diagonal matrix D and the multipliers
! 40: * used to obtain the factor U or L as computed by DSYTRF.
! 41: *
! 42: * On exit, if INFO = 0, the (symmetric) inverse of the original
! 43: * matrix. If UPLO = 'U', the upper triangular part of the
! 44: * inverse is formed and the part of A below the diagonal is not
! 45: * referenced; if UPLO = 'L' the lower triangular part of the
! 46: * inverse is formed and the part of A above the diagonal is
! 47: * not referenced.
! 48: *
! 49: * LDA (input) INTEGER
! 50: * The leading dimension of the array A. LDA >= max(1,N).
! 51: *
! 52: * IPIV (input) INTEGER array, dimension (N)
! 53: * Details of the interchanges and the NNB structure of D
! 54: * as determined by DSYTRF.
! 55: *
! 56: * WORK (workspace) DOUBLE PRECISION array, dimension (N+NNB+1,NNB+3)
! 57: *
! 58: * NB (input) INTEGER
! 59: * Block size
! 60: *
! 61: * INFO (output) INTEGER
! 62: * = 0: successful exit
! 63: * < 0: if INFO = -i, the i-th argument had an illegal value
! 64: * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
! 65: * inverse could not be computed.
! 66: *
! 67: * =====================================================================
! 68: *
! 69: * .. Parameters ..
! 70: DOUBLE PRECISION ONE, ZERO
! 71: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 72: * ..
! 73: * .. Local Scalars ..
! 74: LOGICAL UPPER
! 75: INTEGER I, IINFO, IP, K, CUT, NNB
! 76: INTEGER COUNT
! 77: INTEGER J, U11, INVD
! 78:
! 79: DOUBLE PRECISION AK, AKKP1, AKP1, D, T
! 80: DOUBLE PRECISION U01_I_J, U01_IP1_J
! 81: DOUBLE PRECISION U11_I_J, U11_IP1_J
! 82: * ..
! 83: * .. External Functions ..
! 84: LOGICAL LSAME
! 85: EXTERNAL LSAME
! 86: * ..
! 87: * .. External Subroutines ..
! 88: EXTERNAL DSYCONV, XERBLA, DTRTRI
! 89: EXTERNAL DGEMM, DTRMM, DSYSWAPR
! 90: * ..
! 91: * .. Intrinsic Functions ..
! 92: INTRINSIC MAX
! 93: * ..
! 94: * .. Executable Statements ..
! 95: *
! 96: * Test the input parameters.
! 97: *
! 98: INFO = 0
! 99: UPPER = LSAME( UPLO, 'U' )
! 100: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 101: INFO = -1
! 102: ELSE IF( N.LT.0 ) THEN
! 103: INFO = -2
! 104: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 105: INFO = -4
! 106: END IF
! 107: *
! 108: * Quick return if possible
! 109: *
! 110: *
! 111: IF( INFO.NE.0 ) THEN
! 112: CALL XERBLA( 'DSYTRI2X', -INFO )
! 113: RETURN
! 114: END IF
! 115: IF( N.EQ.0 )
! 116: $ RETURN
! 117: *
! 118: * Convert A
! 119: * Workspace got Non-diag elements of D
! 120: *
! 121: CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
! 122: *
! 123: * Check that the diagonal matrix D is nonsingular.
! 124: *
! 125: IF( UPPER ) THEN
! 126: *
! 127: * Upper triangular storage: examine D from bottom to top
! 128: *
! 129: DO INFO = N, 1, -1
! 130: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
! 131: $ RETURN
! 132: END DO
! 133: ELSE
! 134: *
! 135: * Lower triangular storage: examine D from top to bottom.
! 136: *
! 137: DO INFO = 1, N
! 138: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
! 139: $ RETURN
! 140: END DO
! 141: END IF
! 142: INFO = 0
! 143: *
! 144: * Splitting Workspace
! 145: * U01 is a block (N,NB+1)
! 146: * The first element of U01 is in WORK(1,1)
! 147: * U11 is a block (NB+1,NB+1)
! 148: * The first element of U11 is in WORK(N+1,1)
! 149: U11 = N
! 150: * INVD is a block (N,2)
! 151: * The first element of INVD is in WORK(1,INVD)
! 152: INVD = NB+2
! 153:
! 154: IF( UPPER ) THEN
! 155: *
! 156: * invA = P * inv(U')*inv(D)*inv(U)*P'.
! 157: *
! 158: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
! 159: *
! 160: * inv(D) and inv(D)*inv(U)
! 161: *
! 162: K=1
! 163: DO WHILE ( K .LE. N )
! 164: IF( IPIV( K ).GT.0 ) THEN
! 165: * 1 x 1 diagonal NNB
! 166: WORK(K,INVD) = ONE / A( K, K )
! 167: WORK(K,INVD+1) = 0
! 168: K=K+1
! 169: ELSE
! 170: * 2 x 2 diagonal NNB
! 171: T = WORK(K+1,1)
! 172: AK = A( K, K ) / T
! 173: AKP1 = A( K+1, K+1 ) / T
! 174: AKKP1 = WORK(K+1,1) / T
! 175: D = T*( AK*AKP1-ONE )
! 176: WORK(K,INVD) = AKP1 / D
! 177: WORK(K+1,INVD+1) = AK / D
! 178: WORK(K,INVD+1) = -AKKP1 / D
! 179: WORK(K+1,INVD) = -AKKP1 / D
! 180: K=K+2
! 181: END IF
! 182: END DO
! 183: *
! 184: * inv(U') = (inv(U))'
! 185: *
! 186: * inv(U')*inv(D)*inv(U)
! 187: *
! 188: CUT=N
! 189: DO WHILE (CUT .GT. 0)
! 190: NNB=NB
! 191: IF (CUT .LE. NNB) THEN
! 192: NNB=CUT
! 193: ELSE
! 194: COUNT = 0
! 195: * count negative elements,
! 196: DO I=CUT+1-NNB,CUT
! 197: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
! 198: END DO
! 199: * need a even number for a clear cut
! 200: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
! 201: END IF
! 202:
! 203: CUT=CUT-NNB
! 204: *
! 205: * U01 Block
! 206: *
! 207: DO I=1,CUT
! 208: DO J=1,NNB
! 209: WORK(I,J)=A(I,CUT+J)
! 210: END DO
! 211: END DO
! 212: *
! 213: * U11 Block
! 214: *
! 215: DO I=1,NNB
! 216: WORK(U11+I,I)=ONE
! 217: DO J=1,I-1
! 218: WORK(U11+I,J)=ZERO
! 219: END DO
! 220: DO J=I+1,NNB
! 221: WORK(U11+I,J)=A(CUT+I,CUT+J)
! 222: END DO
! 223: END DO
! 224: *
! 225: * invD*U01
! 226: *
! 227: I=1
! 228: DO WHILE (I .LE. CUT)
! 229: IF (IPIV(I) > 0) THEN
! 230: DO J=1,NNB
! 231: WORK(I,J)=WORK(I,INVD)*WORK(I,J)
! 232: END DO
! 233: I=I+1
! 234: ELSE
! 235: DO J=1,NNB
! 236: U01_I_J = WORK(I,J)
! 237: U01_IP1_J = WORK(I+1,J)
! 238: WORK(I,J)=WORK(I,INVD)*U01_I_J+
! 239: $ WORK(I,INVD+1)*U01_IP1_J
! 240: WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
! 241: $ WORK(I+1,INVD+1)*U01_IP1_J
! 242: END DO
! 243: I=I+2
! 244: END IF
! 245: END DO
! 246: *
! 247: * invD1*U11
! 248: *
! 249: I=1
! 250: DO WHILE (I .LE. NNB)
! 251: IF (IPIV(CUT+I) > 0) THEN
! 252: DO J=I,NNB
! 253: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
! 254: END DO
! 255: I=I+1
! 256: ELSE
! 257: DO J=I,NNB
! 258: U11_I_J = WORK(U11+I,J)
! 259: U11_IP1_J = WORK(U11+I+1,J)
! 260: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
! 261: $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
! 262: WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
! 263: $ WORK(CUT+I+1,INVD+1)*U11_IP1_J
! 264: END DO
! 265: I=I+2
! 266: END IF
! 267: END DO
! 268: *
! 269: * U11T*invD1*U11->U11
! 270: *
! 271: CALL DTRMM('L','U','T','U',NNB, NNB,
! 272: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
! 273: *
! 274: * U01'invD*U01->A(CUT+I,CUT+J)
! 275: *
! 276: CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
! 277: $ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
! 278: *
! 279: * U11 = U11T*invD1*U11 + U01'invD*U01
! 280: *
! 281: DO I=1,NNB
! 282: DO J=I,NNB
! 283: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
! 284: END DO
! 285: END DO
! 286: *
! 287: * U01 = U00T*invD0*U01
! 288: *
! 289: CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
! 290: $ ONE,A,LDA,WORK,N+NB+1)
! 291:
! 292: *
! 293: * Update U01
! 294: *
! 295: DO I=1,CUT
! 296: DO J=1,NNB
! 297: A(I,CUT+J)=WORK(I,J)
! 298: END DO
! 299: END DO
! 300: *
! 301: * Next Block
! 302: *
! 303: END DO
! 304: *
! 305: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
! 306: *
! 307: I=1
! 308: DO WHILE ( I .LE. N )
! 309: IF( IPIV(I) .GT. 0 ) THEN
! 310: IP=IPIV(I)
! 311: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
! 312: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
! 313: ELSE
! 314: IP=-IPIV(I)
! 315: I=I+1
! 316: IF ( (I-1) .LT. IP)
! 317: $ CALL DSYSWAPR( UPLO, N, A, I-1 ,IP )
! 318: IF ( (I-1) .GT. IP)
! 319: $ CALL DSYSWAPR( UPLO, N, A, IP ,I-1 )
! 320: ENDIF
! 321: I=I+1
! 322: END DO
! 323: ELSE
! 324: *
! 325: * LOWER...
! 326: *
! 327: * invA = P * inv(U')*inv(D)*inv(U)*P'.
! 328: *
! 329: CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
! 330: *
! 331: * inv(D) and inv(D)*inv(U)
! 332: *
! 333: K=N
! 334: DO WHILE ( K .GE. 1 )
! 335: IF( IPIV( K ).GT.0 ) THEN
! 336: * 1 x 1 diagonal NNB
! 337: WORK(K,INVD) = ONE / A( K, K )
! 338: WORK(K,INVD+1) = 0
! 339: K=K-1
! 340: ELSE
! 341: * 2 x 2 diagonal NNB
! 342: T = WORK(K-1,1)
! 343: AK = A( K-1, K-1 ) / T
! 344: AKP1 = A( K, K ) / T
! 345: AKKP1 = WORK(K-1,1) / T
! 346: D = T*( AK*AKP1-ONE )
! 347: WORK(K-1,INVD) = AKP1 / D
! 348: WORK(K,INVD) = AK / D
! 349: WORK(K,INVD+1) = -AKKP1 / D
! 350: WORK(K-1,INVD+1) = -AKKP1 / D
! 351: K=K-2
! 352: END IF
! 353: END DO
! 354: *
! 355: * inv(U') = (inv(U))'
! 356: *
! 357: * inv(U')*inv(D)*inv(U)
! 358: *
! 359: CUT=0
! 360: DO WHILE (CUT .LT. N)
! 361: NNB=NB
! 362: IF (CUT + NNB .GT. N) THEN
! 363: NNB=N-CUT
! 364: ELSE
! 365: COUNT = 0
! 366: * count negative elements,
! 367: DO I=CUT+1,CUT+NNB
! 368: IF (IPIV(I) .LT. 0) COUNT=COUNT+1
! 369: END DO
! 370: * need a even number for a clear cut
! 371: IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
! 372: END IF
! 373: * L21 Block
! 374: DO I=1,N-CUT-NNB
! 375: DO J=1,NNB
! 376: WORK(I,J)=A(CUT+NNB+I,CUT+J)
! 377: END DO
! 378: END DO
! 379: * L11 Block
! 380: DO I=1,NNB
! 381: WORK(U11+I,I)=ONE
! 382: DO J=I+1,NNB
! 383: WORK(U11+I,J)=ZERO
! 384: END DO
! 385: DO J=1,I-1
! 386: WORK(U11+I,J)=A(CUT+I,CUT+J)
! 387: END DO
! 388: END DO
! 389: *
! 390: * invD*L21
! 391: *
! 392: I=N-CUT-NNB
! 393: DO WHILE (I .GE. 1)
! 394: IF (IPIV(CUT+NNB+I) > 0) THEN
! 395: DO J=1,NNB
! 396: WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
! 397: END DO
! 398: I=I-1
! 399: ELSE
! 400: DO J=1,NNB
! 401: U01_I_J = WORK(I,J)
! 402: U01_IP1_J = WORK(I-1,J)
! 403: WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
! 404: $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
! 405: WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
! 406: $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
! 407: END DO
! 408: I=I-2
! 409: END IF
! 410: END DO
! 411: *
! 412: * invD1*L11
! 413: *
! 414: I=NNB
! 415: DO WHILE (I .GE. 1)
! 416: IF (IPIV(CUT+I) > 0) THEN
! 417: DO J=1,NNB
! 418: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
! 419: END DO
! 420: I=I-1
! 421: ELSE
! 422: DO J=1,NNB
! 423: U11_I_J = WORK(U11+I,J)
! 424: U11_IP1_J = WORK(U11+I-1,J)
! 425: WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
! 426: $ WORK(CUT+I,INVD+1)*U11_IP1_J
! 427: WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
! 428: $ WORK(CUT+I-1,INVD)*U11_IP1_J
! 429: END DO
! 430: I=I-2
! 431: END IF
! 432: END DO
! 433: *
! 434: * L11T*invD1*L11->L11
! 435: *
! 436: CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
! 437: $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
! 438:
! 439: IF ( (CUT+NNB) .LT. N ) THEN
! 440: *
! 441: * L21T*invD2*L21->A(CUT+I,CUT+J)
! 442: *
! 443: CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
! 444: $ ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
! 445:
! 446: *
! 447: * L11 = L11T*invD1*L11 + U01'invD*U01
! 448: *
! 449: DO I=1,NNB
! 450: DO J=1,I
! 451: A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
! 452: END DO
! 453: END DO
! 454: *
! 455: * L01 = L22T*invD2*L21
! 456: *
! 457: CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
! 458: $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
! 459: *
! 460: * Update L21
! 461: *
! 462: DO I=1,N-CUT-NNB
! 463: DO J=1,NNB
! 464: A(CUT+NNB+I,CUT+J)=WORK(I,J)
! 465: END DO
! 466: END DO
! 467:
! 468: ELSE
! 469: *
! 470: * L11 = L11T*invD1*L11
! 471: *
! 472: DO I=1,NNB
! 473: DO J=1,I
! 474: A(CUT+I,CUT+J)=WORK(U11+I,J)
! 475: END DO
! 476: END DO
! 477: END IF
! 478: *
! 479: * Next Block
! 480: *
! 481: CUT=CUT+NNB
! 482: END DO
! 483: *
! 484: * Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
! 485: *
! 486: I=N
! 487: DO WHILE ( I .GE. 1 )
! 488: IF( IPIV(I) .GT. 0 ) THEN
! 489: IP=IPIV(I)
! 490: IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
! 491: IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
! 492: ELSE
! 493: IP=-IPIV(I)
! 494: IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
! 495: IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
! 496: I=I-1
! 497: ENDIF
! 498: I=I-1
! 499: END DO
! 500: END IF
! 501: *
! 502: RETURN
! 503: *
! 504: * End of DSYTRI2X
! 505: *
! 506: END
! 507:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytri2x.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack