![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTRTTF( TRANSR, UPLO, N, A, LDA, ARF, 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: * ZTRTTF copies a triangular matrix A from standard full format (TR) 23: * to rectangular full packed format (TF) . 24: * 25: * Arguments 26: * ========= 27: * 28: * TRANSR (input) CHARACTER*1 29: * = 'N': ARF in Normal mode is wanted; 30: * = 'C': ARF in Conjugate Transpose mode is wanted; 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: * A (input) COMPLEX*16 array, dimension ( LDA, N ) 40: * On entry, the triangular matrix A. If UPLO = 'U', the 41: * leading N-by-N upper triangular part of the array A contains 42: * the upper triangular matrix, and the strictly lower 43: * triangular part of A is not referenced. If UPLO = 'L', the 44: * leading N-by-N lower triangular part of the array A contains 45: * the lower triangular matrix, and the strictly upper 46: * triangular part of A is not referenced. 47: * 48: * LDA (input) INTEGER 49: * The leading dimension of the matrix A. LDA >= max(1,N). 50: * 51: * ARF (output) COMPLEX*16 array, dimension ( N*(N+1)/2 ), 52: * On exit, the upper or lower triangular matrix A stored in 53: * RFP format. For a further discussion see Notes below. 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 I, IJ, J, K, L, N1, N2, NT, NX2, NP1X2 171: * .. 172: * .. External Functions .. 173: LOGICAL LSAME 174: EXTERNAL LSAME 175: * .. 176: * .. External Subroutines .. 177: EXTERNAL XERBLA 178: * .. 179: * .. Intrinsic Functions .. 180: INTRINSIC DCONJG, MAX, MOD 181: * .. 182: * .. Executable Statements .. 183: * 184: * Test the input parameters. 185: * 186: INFO = 0 187: NORMALTRANSR = LSAME( TRANSR, 'N' ) 188: LOWER = LSAME( UPLO, 'L' ) 189: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 190: INFO = -1 191: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 192: INFO = -2 193: ELSE IF( N.LT.0 ) THEN 194: INFO = -3 195: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 196: INFO = -5 197: END IF 198: IF( INFO.NE.0 ) THEN 199: CALL XERBLA( 'ZTRTTF', -INFO ) 200: RETURN 201: END IF 202: * 203: * Quick return if possible 204: * 205: IF( N.LE.1 ) THEN 206: IF( N.EQ.1 ) THEN 207: IF( NORMALTRANSR ) THEN 208: ARF( 0 ) = A( 0, 0 ) 209: ELSE 210: ARF( 0 ) = DCONJG( A( 0, 0 ) ) 211: END IF 212: END IF 213: RETURN 214: END IF 215: * 216: * Size of array ARF(1:2,0:nt-1) 217: * 218: NT = N*( N+1 ) / 2 219: * 220: * set N1 and N2 depending on LOWER: for N even N1=N2=K 221: * 222: IF( LOWER ) THEN 223: N2 = N / 2 224: N1 = N - N2 225: ELSE 226: N1 = N / 2 227: N2 = N - N1 228: END IF 229: * 230: * If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2. 231: * If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is 232: * N--by--(N+1)/2. 233: * 234: IF( MOD( N, 2 ).EQ.0 ) THEN 235: K = N / 2 236: NISODD = .FALSE. 237: IF( .NOT.LOWER ) 238: + NP1X2 = N + N + 2 239: ELSE 240: NISODD = .TRUE. 241: IF( .NOT.LOWER ) 242: + NX2 = N + N 243: END IF 244: * 245: IF( NISODD ) THEN 246: * 247: * N is odd 248: * 249: IF( NORMALTRANSR ) THEN 250: * 251: * N is odd and TRANSR = 'N' 252: * 253: IF( LOWER ) THEN 254: * 255: * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 256: * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 257: * T1 -> a(0), T2 -> a(n), S -> a(n1); lda=n 258: * 259: IJ = 0 260: DO J = 0, N2 261: DO I = N1, N2 + J 262: ARF( IJ ) = DCONJG( A( N2+J, I ) ) 263: IJ = IJ + 1 264: END DO 265: DO I = J, N - 1 266: ARF( IJ ) = A( I, J ) 267: IJ = IJ + 1 268: END DO 269: END DO 270: * 271: ELSE 272: * 273: * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 274: * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 275: * T1 -> a(n2), T2 -> a(n1), S -> a(0); lda=n 276: * 277: IJ = NT - N 278: DO J = N - 1, N1, -1 279: DO I = 0, J 280: ARF( IJ ) = A( I, J ) 281: IJ = IJ + 1 282: END DO 283: DO L = J - N1, N1 - 1 284: ARF( IJ ) = DCONJG( A( J-N1, L ) ) 285: IJ = IJ + 1 286: END DO 287: IJ = IJ - NX2 288: END DO 289: * 290: END IF 291: * 292: ELSE 293: * 294: * N is odd and TRANSR = 'C' 295: * 296: IF( LOWER ) THEN 297: * 298: * SRPA for LOWER, TRANSPOSE and N is odd 299: * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 300: * T1 -> A(0+0) , T2 -> A(1+0) , S -> A(0+n1*n1); lda=n1 301: * 302: IJ = 0 303: DO J = 0, N2 - 1 304: DO I = 0, J 305: ARF( IJ ) = DCONJG( A( J, I ) ) 306: IJ = IJ + 1 307: END DO 308: DO I = N1 + J, N - 1 309: ARF( IJ ) = A( I, N1+J ) 310: IJ = IJ + 1 311: END DO 312: END DO 313: DO J = N2, N - 1 314: DO I = 0, N1 - 1 315: ARF( IJ ) = DCONJG( A( J, I ) ) 316: IJ = IJ + 1 317: END DO 318: END DO 319: * 320: ELSE 321: * 322: * SRPA for UPPER, TRANSPOSE and N is odd 323: * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 324: * T1 -> A(n2*n2), T2 -> A(n1*n2), S -> A(0); lda=n2 325: * 326: IJ = 0 327: DO J = 0, N1 328: DO I = N1, N - 1 329: ARF( IJ ) = DCONJG( A( J, I ) ) 330: IJ = IJ + 1 331: END DO 332: END DO 333: DO J = 0, N1 - 1 334: DO I = 0, J 335: ARF( IJ ) = A( I, J ) 336: IJ = IJ + 1 337: END DO 338: DO L = N2 + J, N - 1 339: ARF( IJ ) = DCONJG( A( N2+J, L ) ) 340: IJ = IJ + 1 341: END DO 342: END DO 343: * 344: END IF 345: * 346: END IF 347: * 348: ELSE 349: * 350: * N is even 351: * 352: IF( NORMALTRANSR ) THEN 353: * 354: * N is even and TRANSR = 'N' 355: * 356: IF( LOWER ) THEN 357: * 358: * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 359: * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 360: * T1 -> a(1), T2 -> a(0), S -> a(k+1); lda=n+1 361: * 362: IJ = 0 363: DO J = 0, K - 1 364: DO I = K, K + J 365: ARF( IJ ) = DCONJG( A( K+J, I ) ) 366: IJ = IJ + 1 367: END DO 368: DO I = J, N - 1 369: ARF( IJ ) = A( I, J ) 370: IJ = IJ + 1 371: END DO 372: END DO 373: * 374: ELSE 375: * 376: * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 377: * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 378: * T1 -> a(k+1), T2 -> a(k), S -> a(0); lda=n+1 379: * 380: IJ = NT - N - 1 381: DO J = N - 1, K, -1 382: DO I = 0, J 383: ARF( IJ ) = A( I, J ) 384: IJ = IJ + 1 385: END DO 386: DO L = J - K, K - 1 387: ARF( IJ ) = DCONJG( A( J-K, L ) ) 388: IJ = IJ + 1 389: END DO 390: IJ = IJ - NP1X2 391: END DO 392: * 393: END IF 394: * 395: ELSE 396: * 397: * N is even and TRANSR = 'C' 398: * 399: IF( LOWER ) THEN 400: * 401: * SRPA for LOWER, TRANSPOSE and N is even (see paper, A=B) 402: * T1 -> A(0,1) , T2 -> A(0,0) , S -> A(0,k+1) : 403: * T1 -> A(0+k) , T2 -> A(0+0) , S -> A(0+k*(k+1)); lda=k 404: * 405: IJ = 0 406: J = K 407: DO I = K, N - 1 408: ARF( IJ ) = A( I, J ) 409: IJ = IJ + 1 410: END DO 411: DO J = 0, K - 2 412: DO I = 0, J 413: ARF( IJ ) = DCONJG( A( J, I ) ) 414: IJ = IJ + 1 415: END DO 416: DO I = K + 1 + J, N - 1 417: ARF( IJ ) = A( I, K+1+J ) 418: IJ = IJ + 1 419: END DO 420: END DO 421: DO J = K - 1, N - 1 422: DO I = 0, K - 1 423: ARF( IJ ) = DCONJG( A( J, I ) ) 424: IJ = IJ + 1 425: END DO 426: END DO 427: * 428: ELSE 429: * 430: * SRPA for UPPER, TRANSPOSE and N is even (see paper, A=B) 431: * T1 -> A(0,k+1) , T2 -> A(0,k) , S -> A(0,0) 432: * T1 -> A(0+k*(k+1)) , T2 -> A(0+k*k) , S -> A(0+0)); lda=k 433: * 434: IJ = 0 435: DO J = 0, K 436: DO I = K, N - 1 437: ARF( IJ ) = DCONJG( A( J, I ) ) 438: IJ = IJ + 1 439: END DO 440: END DO 441: DO J = 0, K - 2 442: DO I = 0, J 443: ARF( IJ ) = A( I, J ) 444: IJ = IJ + 1 445: END DO 446: DO L = K + 1 + J, N - 1 447: ARF( IJ ) = DCONJG( A( K+1+J, L ) ) 448: IJ = IJ + 1 449: END DO 450: END DO 451: * 452: * Note that here J = K-1 453: * 454: DO I = 0, J 455: ARF( IJ ) = A( I, J ) 456: IJ = IJ + 1 457: END DO 458: * 459: END IF 460: * 461: END IF 462: * 463: END IF 464: * 465: RETURN 466: * 467: * End of ZTRTTF 468: * 469: END