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