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