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