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