![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPFTRI( 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: DOUBLE PRECISION A( 0: * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DPFTRI computes the inverse of a (real) symmetric positive definite 22: * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T 23: * computed by DPFTRF. 24: * 25: * Arguments 26: * ========= 27: * 28: * TRANSR (input) CHARACTER*1 29: * = 'N': The Normal TRANSR of RFP A is stored; 30: * = 'T': The 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) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ) 40: * On entry, the symmetric 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 = 'T' then RFP is 44: * the 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: * 'T'. 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 symmetric 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 Rectangular Full Packed (RFP) Format when N is 65: * even. 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: * the 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: * the transpose of the last three columns of AP lower. 84: * This covers the case N even and TRANSR = 'N'. 85: * 86: * RFP A RFP A 87: * 88: * 03 04 05 33 43 53 89: * 13 14 15 00 44 54 90: * 23 24 25 10 11 55 91: * 33 34 35 20 21 22 92: * 00 44 45 30 31 32 93: * 01 11 55 40 41 42 94: * 02 12 22 50 51 52 95: * 96: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 97: * transpose of RFP A above. One therefore gets: 98: * 99: * 100: * RFP A RFP A 101: * 102: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 103: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 104: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 105: * 106: * 107: * We then consider Rectangular Full Packed (RFP) Format when N is 108: * odd. We give an example where N = 5. 109: * 110: * AP is Upper AP is Lower 111: * 112: * 00 01 02 03 04 00 113: * 11 12 13 14 10 11 114: * 22 23 24 20 21 22 115: * 33 34 30 31 32 33 116: * 44 40 41 42 43 44 117: * 118: * 119: * Let TRANSR = 'N'. RFP holds AP as follows: 120: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 121: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 122: * the transpose of the first two columns of AP upper. 123: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 124: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 125: * the transpose of the last two columns of AP lower. 126: * This covers the case N odd and TRANSR = 'N'. 127: * 128: * RFP A RFP A 129: * 130: * 02 03 04 00 33 43 131: * 12 13 14 10 11 44 132: * 22 23 24 20 21 22 133: * 00 33 34 30 31 32 134: * 01 11 44 40 41 42 135: * 136: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 137: * transpose of RFP A above. One therefore gets: 138: * 139: * RFP A RFP A 140: * 141: * 02 12 22 00 01 00 10 20 30 40 50 142: * 03 13 23 33 11 33 11 21 31 41 51 143: * 04 14 24 34 44 43 44 22 32 42 52 144: * 145: * ===================================================================== 146: * 147: * .. Parameters .. 148: DOUBLE PRECISION ONE 149: PARAMETER ( ONE = 1.0D+0 ) 150: * .. 151: * .. Local Scalars .. 152: LOGICAL LOWER, NISODD, NORMALTRANSR 153: INTEGER N1, N2, K 154: * .. 155: * .. External Functions .. 156: LOGICAL LSAME 157: EXTERNAL LSAME 158: * .. 159: * .. External Subroutines .. 160: EXTERNAL XERBLA, DTFTRI, DLAUUM, DTRMM, DSYRK 161: * .. 162: * .. Intrinsic Functions .. 163: INTRINSIC MOD 164: * .. 165: * .. Executable Statements .. 166: * 167: * Test the input parameters. 168: * 169: INFO = 0 170: NORMALTRANSR = LSAME( TRANSR, 'N' ) 171: LOWER = LSAME( UPLO, 'L' ) 172: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 173: INFO = -1 174: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 175: INFO = -2 176: ELSE IF( N.LT.0 ) THEN 177: INFO = -3 178: END IF 179: IF( INFO.NE.0 ) THEN 180: CALL XERBLA( 'DPFTRI', -INFO ) 181: RETURN 182: END IF 183: * 184: * Quick return if possible 185: * 186: IF( N.EQ.0 ) 187: + RETURN 188: * 189: * Invert the triangular Cholesky factor U or L. 190: * 191: CALL DTFTRI( TRANSR, UPLO, 'N', N, A, INFO ) 192: IF( INFO.GT.0 ) 193: + RETURN 194: * 195: * If N is odd, set NISODD = .TRUE. 196: * If N is even, set K = N/2 and NISODD = .FALSE. 197: * 198: IF( MOD( N, 2 ).EQ.0 ) THEN 199: K = N / 2 200: NISODD = .FALSE. 201: ELSE 202: NISODD = .TRUE. 203: END IF 204: * 205: * Set N1 and N2 depending on LOWER 206: * 207: IF( LOWER ) THEN 208: N2 = N / 2 209: N1 = N - N2 210: ELSE 211: N1 = N / 2 212: N2 = N - N1 213: END IF 214: * 215: * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or 216: * inv(L)^C*inv(L). There are eight cases. 217: * 218: IF( NISODD ) THEN 219: * 220: * N is odd 221: * 222: IF( NORMALTRANSR ) THEN 223: * 224: * N is odd and TRANSR = 'N' 225: * 226: IF( LOWER ) THEN 227: * 228: * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) ) 229: * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0) 230: * T1 -> a(0), T2 -> a(n), S -> a(N1) 231: * 232: CALL DLAUUM( 'L', N1, A( 0 ), N, INFO ) 233: CALL DSYRK( 'L', 'T', N1, N2, ONE, A( N1 ), N, ONE, 234: + A( 0 ), N ) 235: CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N, 236: + A( N1 ), N ) 237: CALL DLAUUM( 'U', N2, A( N ), N, INFO ) 238: * 239: ELSE 240: * 241: * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1) 242: * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0) 243: * T1 -> a(N2), T2 -> a(N1), S -> a(0) 244: * 245: CALL DLAUUM( 'L', N1, A( N2 ), N, INFO ) 246: CALL DSYRK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE, 247: + A( N2 ), N ) 248: CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N, 249: + A( 0 ), N ) 250: CALL DLAUUM( 'U', N2, A( N1 ), N, INFO ) 251: * 252: END IF 253: * 254: ELSE 255: * 256: * N is odd and TRANSR = 'T' 257: * 258: IF( LOWER ) THEN 259: * 260: * SRPA for LOWER, TRANSPOSE, and N is odd 261: * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1) 262: * 263: CALL DLAUUM( 'U', N1, A( 0 ), N1, INFO ) 264: CALL DSYRK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE, 265: + A( 0 ), N1 ) 266: CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1, 267: + A( N1*N1 ), N1 ) 268: CALL DLAUUM( 'L', N2, A( 1 ), N1, INFO ) 269: * 270: ELSE 271: * 272: * SRPA for UPPER, TRANSPOSE, and N is odd 273: * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0) 274: * 275: CALL DLAUUM( 'U', N1, A( N2*N2 ), N2, INFO ) 276: CALL DSYRK( 'U', 'T', N1, N2, ONE, A( 0 ), N2, ONE, 277: + A( N2*N2 ), N2 ) 278: CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ), 279: + N2, A( 0 ), N2 ) 280: CALL DLAUUM( 'L', N2, A( N1*N2 ), N2, INFO ) 281: * 282: END IF 283: * 284: END IF 285: * 286: ELSE 287: * 288: * N is even 289: * 290: IF( NORMALTRANSR ) THEN 291: * 292: * N is even and TRANSR = 'N' 293: * 294: IF( LOWER ) THEN 295: * 296: * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 297: * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 298: * T1 -> a(1), T2 -> a(0), S -> a(k+1) 299: * 300: CALL DLAUUM( 'L', K, A( 1 ), N+1, INFO ) 301: CALL DSYRK( 'L', 'T', K, K, ONE, A( K+1 ), N+1, ONE, 302: + A( 1 ), N+1 ) 303: CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1, 304: + A( K+1 ), N+1 ) 305: CALL DLAUUM( 'U', K, A( 0 ), N+1, INFO ) 306: * 307: ELSE 308: * 309: * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 310: * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 311: * T1 -> a(k+1), T2 -> a(k), S -> a(0) 312: * 313: CALL DLAUUM( 'L', K, A( K+1 ), N+1, INFO ) 314: CALL DSYRK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE, 315: + A( K+1 ), N+1 ) 316: CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1, 317: + A( 0 ), N+1 ) 318: CALL DLAUUM( 'U', K, A( K ), N+1, INFO ) 319: * 320: END IF 321: * 322: ELSE 323: * 324: * N is even and TRANSR = 'T' 325: * 326: IF( LOWER ) THEN 327: * 328: * SRPA for LOWER, TRANSPOSE, and N is even (see paper) 329: * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1), 330: * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 331: * 332: CALL DLAUUM( 'U', K, A( K ), K, INFO ) 333: CALL DSYRK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE, 334: + A( K ), K ) 335: CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K, 336: + A( K*( K+1 ) ), K ) 337: CALL DLAUUM( 'L', K, A( 0 ), K, INFO ) 338: * 339: ELSE 340: * 341: * SRPA for UPPER, TRANSPOSE, and N is even (see paper) 342: * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0), 343: * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 344: * 345: CALL DLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO ) 346: CALL DSYRK( 'U', 'T', K, K, ONE, A( 0 ), K, ONE, 347: + A( K*( K+1 ) ), K ) 348: CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K, 349: + A( 0 ), K ) 350: CALL DLAUUM( 'L', K, A( K*K ), K, INFO ) 351: * 352: END IF 353: * 354: END IF 355: * 356: END IF 357: * 358: RETURN 359: * 360: * End of DPFTRI 361: * 362: END