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