![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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: