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