![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO ) 2: * 3: * -- LAPACK routine (version 3.3.0) -- 4: * 5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 6: * November 2010 7: * 8: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 10: * 11: * .. Scalar Arguments .. 12: CHARACTER TRANSR, UPLO 13: INTEGER INFO, N 14: * .. Array Arguments .. 15: COMPLEX*16 A( 0: * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * ZPFTRI computes the inverse of a complex Hermitian positive definite 22: * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 23: * computed by ZPFTRF. 24: * 25: * Arguments 26: * ========= 27: * 28: * TRANSR (input) CHARACTER*1 29: * = 'N': The Normal TRANSR of RFP A is stored; 30: * = 'C': The Conjugate-transpose TRANSR of RFP A is stored. 31: * 32: * UPLO (input) CHARACTER*1 33: * = 'U': Upper triangle of A is stored; 34: * = 'L': Lower triangle of A is stored. 35: * 36: * N (input) INTEGER 37: * The order of the matrix A. N >= 0. 38: * 39: * A (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 ); 40: * On entry, the Hermitian matrix A in RFP format. RFP format is 41: * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' 42: * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is 43: * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is 44: * the Conjugate-transpose of RFP A as defined when 45: * TRANSR = 'N'. The contents of RFP A are defined by UPLO as 46: * follows: If UPLO = 'U' the RFP A contains the nt elements of 47: * upper packed A. If UPLO = 'L' the RFP A contains the elements 48: * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 49: * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N 50: * is odd. See the Note below for more details. 51: * 52: * On exit, the Hermitian inverse of the original matrix, in the 53: * same storage format. 54: * 55: * INFO (output) INTEGER 56: * = 0: successful exit 57: * < 0: if INFO = -i, the i-th argument had an illegal value 58: * > 0: if INFO = i, the (i,i) element of the factor U or L is 59: * zero, and the inverse could not be computed. 60: * 61: * Further Details 62: * =============== 63: * 64: * We first consider Standard Packed Format when N is even. 65: * We give an example where N = 6. 66: * 67: * AP is Upper AP is Lower 68: * 69: * 00 01 02 03 04 05 00 70: * 11 12 13 14 15 10 11 71: * 22 23 24 25 20 21 22 72: * 33 34 35 30 31 32 33 73: * 44 45 40 41 42 43 44 74: * 55 50 51 52 53 54 55 75: * 76: * 77: * Let TRANSR = 'N'. RFP holds AP as follows: 78: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 79: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 80: * conjugate-transpose of the first three columns of AP upper. 81: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 82: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 83: * conjugate-transpose of the last three columns of AP lower. 84: * To denote conjugate we place -- above the element. This covers the 85: * case N even and TRANSR = 'N'. 86: * 87: * RFP A RFP A 88: * 89: * -- -- -- 90: * 03 04 05 33 43 53 91: * -- -- 92: * 13 14 15 00 44 54 93: * -- 94: * 23 24 25 10 11 55 95: * 96: * 33 34 35 20 21 22 97: * -- 98: * 00 44 45 30 31 32 99: * -- -- 100: * 01 11 55 40 41 42 101: * -- -- -- 102: * 02 12 22 50 51 52 103: * 104: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 105: * transpose of RFP A above. One therefore gets: 106: * 107: * 108: * RFP A RFP A 109: * 110: * -- -- -- -- -- -- -- -- -- -- 111: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 112: * -- -- -- -- -- -- -- -- -- -- 113: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 114: * -- -- -- -- -- -- -- -- -- -- 115: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 116: * 117: * 118: * We next consider Standard Packed Format when N is odd. 119: * We give an example where N = 5. 120: * 121: * AP is Upper AP is Lower 122: * 123: * 00 01 02 03 04 00 124: * 11 12 13 14 10 11 125: * 22 23 24 20 21 22 126: * 33 34 30 31 32 33 127: * 44 40 41 42 43 44 128: * 129: * 130: * Let TRANSR = 'N'. RFP holds AP as follows: 131: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 132: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 133: * conjugate-transpose of the first two columns of AP upper. 134: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 135: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 136: * conjugate-transpose of the last two columns of AP lower. 137: * To denote conjugate we place -- above the element. This covers the 138: * case N odd and TRANSR = 'N'. 139: * 140: * RFP A RFP A 141: * 142: * -- -- 143: * 02 03 04 00 33 43 144: * -- 145: * 12 13 14 10 11 44 146: * 147: * 22 23 24 20 21 22 148: * -- 149: * 00 33 34 30 31 32 150: * -- -- 151: * 01 11 44 40 41 42 152: * 153: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 154: * transpose of RFP A above. One therefore gets: 155: * 156: * 157: * RFP A RFP A 158: * 159: * -- -- -- -- -- -- -- -- -- 160: * 02 12 22 00 01 00 10 20 30 40 50 161: * -- -- -- -- -- -- -- -- -- 162: * 03 13 23 33 11 33 11 21 31 41 51 163: * -- -- -- -- -- -- -- -- -- 164: * 04 14 24 34 44 43 44 22 32 42 52 165: * 166: * ===================================================================== 167: * 168: * .. Parameters .. 169: DOUBLE PRECISION ONE 170: COMPLEX*16 CONE 171: PARAMETER ( ONE = 1.D0, CONE = ( 1.D0, 0.D0 ) ) 172: * .. 173: * .. Local Scalars .. 174: LOGICAL LOWER, NISODD, NORMALTRANSR 175: INTEGER N1, N2, K 176: * .. 177: * .. External Functions .. 178: LOGICAL LSAME 179: EXTERNAL LSAME 180: * .. 181: * .. External Subroutines .. 182: EXTERNAL XERBLA, ZTFTRI, ZLAUUM, ZTRMM, ZHERK 183: * .. 184: * .. Intrinsic Functions .. 185: INTRINSIC MOD 186: * .. 187: * .. Executable Statements .. 188: * 189: * Test the input parameters. 190: * 191: INFO = 0 192: NORMALTRANSR = LSAME( TRANSR, 'N' ) 193: LOWER = LSAME( UPLO, 'L' ) 194: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 195: INFO = -1 196: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 197: INFO = -2 198: ELSE IF( N.LT.0 ) THEN 199: INFO = -3 200: END IF 201: IF( INFO.NE.0 ) THEN 202: CALL XERBLA( 'ZPFTRI', -INFO ) 203: RETURN 204: END IF 205: * 206: * Quick return if possible 207: * 208: IF( N.EQ.0 ) 209: + RETURN 210: * 211: * Invert the triangular Cholesky factor U or L. 212: * 213: CALL ZTFTRI( TRANSR, UPLO, 'N', N, A, INFO ) 214: IF( INFO.GT.0 ) 215: + RETURN 216: * 217: * If N is odd, set NISODD = .TRUE. 218: * If N is even, set K = N/2 and NISODD = .FALSE. 219: * 220: IF( MOD( N, 2 ).EQ.0 ) THEN 221: K = N / 2 222: NISODD = .FALSE. 223: ELSE 224: NISODD = .TRUE. 225: END IF 226: * 227: * Set N1 and N2 depending on LOWER 228: * 229: IF( LOWER ) THEN 230: N2 = N / 2 231: N1 = N - N2 232: ELSE 233: N1 = N / 2 234: N2 = N - N1 235: END IF 236: * 237: * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or 238: * inv(L)^C*inv(L). There are eight cases. 239: * 240: IF( NISODD ) THEN 241: * 242: * N is odd 243: * 244: IF( NORMALTRANSR ) THEN 245: * 246: * N is odd and TRANSR = 'N' 247: * 248: IF( LOWER ) THEN 249: * 250: * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) ) 251: * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0) 252: * T1 -> a(0), T2 -> a(n), S -> a(N1) 253: * 254: CALL ZLAUUM( 'L', N1, A( 0 ), N, INFO ) 255: CALL ZHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE, 256: + A( 0 ), N ) 257: CALL ZTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N, 258: + A( N1 ), N ) 259: CALL ZLAUUM( 'U', N2, A( N ), N, INFO ) 260: * 261: ELSE 262: * 263: * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1) 264: * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0) 265: * T1 -> a(N2), T2 -> a(N1), S -> a(0) 266: * 267: CALL ZLAUUM( 'L', N1, A( N2 ), N, INFO ) 268: CALL ZHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE, 269: + A( N2 ), N ) 270: CALL ZTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N, 271: + A( 0 ), N ) 272: CALL ZLAUUM( 'U', N2, A( N1 ), N, INFO ) 273: * 274: END IF 275: * 276: ELSE 277: * 278: * N is odd and TRANSR = 'C' 279: * 280: IF( LOWER ) THEN 281: * 282: * SRPA for LOWER, TRANSPOSE, and N is odd 283: * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1) 284: * 285: CALL ZLAUUM( 'U', N1, A( 0 ), N1, INFO ) 286: CALL ZHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE, 287: + A( 0 ), N1 ) 288: CALL ZTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1, 289: + A( N1*N1 ), N1 ) 290: CALL ZLAUUM( 'L', N2, A( 1 ), N1, INFO ) 291: * 292: ELSE 293: * 294: * SRPA for UPPER, TRANSPOSE, and N is odd 295: * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0) 296: * 297: CALL ZLAUUM( 'U', N1, A( N2*N2 ), N2, INFO ) 298: CALL ZHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE, 299: + A( N2*N2 ), N2 ) 300: CALL ZTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ), 301: + N2, A( 0 ), N2 ) 302: CALL ZLAUUM( 'L', N2, A( N1*N2 ), N2, INFO ) 303: * 304: END IF 305: * 306: END IF 307: * 308: ELSE 309: * 310: * N is even 311: * 312: IF( NORMALTRANSR ) THEN 313: * 314: * N is even and TRANSR = 'N' 315: * 316: IF( LOWER ) THEN 317: * 318: * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 319: * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 320: * T1 -> a(1), T2 -> a(0), S -> a(k+1) 321: * 322: CALL ZLAUUM( 'L', K, A( 1 ), N+1, INFO ) 323: CALL ZHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE, 324: + A( 1 ), N+1 ) 325: CALL ZTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1, 326: + A( K+1 ), N+1 ) 327: CALL ZLAUUM( 'U', K, A( 0 ), N+1, INFO ) 328: * 329: ELSE 330: * 331: * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 332: * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 333: * T1 -> a(k+1), T2 -> a(k), S -> a(0) 334: * 335: CALL ZLAUUM( 'L', K, A( K+1 ), N+1, INFO ) 336: CALL ZHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE, 337: + A( K+1 ), N+1 ) 338: CALL ZTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1, 339: + A( 0 ), N+1 ) 340: CALL ZLAUUM( 'U', K, A( K ), N+1, INFO ) 341: * 342: END IF 343: * 344: ELSE 345: * 346: * N is even and TRANSR = 'C' 347: * 348: IF( LOWER ) THEN 349: * 350: * SRPA for LOWER, TRANSPOSE, and N is even (see paper) 351: * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1), 352: * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 353: * 354: CALL ZLAUUM( 'U', K, A( K ), K, INFO ) 355: CALL ZHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE, 356: + A( K ), K ) 357: CALL ZTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K, 358: + A( K*( K+1 ) ), K ) 359: CALL ZLAUUM( 'L', K, A( 0 ), K, INFO ) 360: * 361: ELSE 362: * 363: * SRPA for UPPER, TRANSPOSE, and N is even (see paper) 364: * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0), 365: * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 366: * 367: CALL ZLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO ) 368: CALL ZHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE, 369: + A( K*( K+1 ) ), K ) 370: CALL ZTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K, 371: + A( 0 ), K ) 372: CALL ZLAUUM( 'L', K, A( K*K ), K, INFO ) 373: * 374: END IF 375: * 376: END IF 377: * 378: END IF 379: * 380: RETURN 381: * 382: * End of ZPFTRI 383: * 384: END