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