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