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