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