001: SUBROUTINE STPTTF( TRANSR, UPLO, N, AP, ARF, INFO ) 002: * 003: * -- LAPACK routine (version 3.2) -- 004: * 005: * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 006: * -- November 2008 -- 007: * 008: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 009: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 010: * 011: * .. 012: * .. Scalar Arguments .. 013: CHARACTER TRANSR, UPLO 014: INTEGER INFO, N 015: * .. 016: * .. Array Arguments .. 017: REAL AP( 0: * ), ARF( 0: * ) 018: * 019: * Purpose 020: * ======= 021: * 022: * STPTTF copies a triangular matrix A from standard packed format (TP) 023: * to rectangular full packed format (TF). 024: * 025: * Arguments 026: * ========= 027: * 028: * TRANSR (input) CHARACTER 029: * = 'N': ARF in Normal format is wanted; 030: * = 'T': ARF in Conjugate-transpose format is wanted. 031: * 032: * UPLO (input) CHARACTER 033: * = 'U': A is upper triangular; 034: * = 'L': A is lower triangular. 035: * 036: * N (input) INTEGER 037: * The order of the matrix A. N >= 0. 038: * 039: * AP (input) REAL array, dimension ( N*(N+1)/2 ), 040: * On entry, the upper or lower triangular matrix A, packed 041: * columnwise in a linear array. The j-th column of A is stored 042: * in the array AP as follows: 043: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 044: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 045: * 046: * ARF (output) REAL array, dimension ( N*(N+1)/2 ), 047: * On exit, the upper or lower triangular matrix A stored in 048: * RFP format. For a further discussion see Notes below. 049: * 050: * INFO (output) INTEGER 051: * = 0: successful exit 052: * < 0: if INFO = -i, the i-th argument had an illegal value 053: * 054: * Notes 055: * ===== 056: * 057: * We first consider Rectangular Full Packed (RFP) Format when N is 058: * even. We give an example where N = 6. 059: * 060: * AP is Upper AP is Lower 061: * 062: * 00 01 02 03 04 05 00 063: * 11 12 13 14 15 10 11 064: * 22 23 24 25 20 21 22 065: * 33 34 35 30 31 32 33 066: * 44 45 40 41 42 43 44 067: * 55 50 51 52 53 54 55 068: * 069: * 070: * Let TRANSR = 'N'. RFP holds AP as follows: 071: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 072: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 073: * the transpose of the first three columns of AP upper. 074: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 075: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 076: * the transpose of the last three columns of AP lower. 077: * This covers the case N even and TRANSR = 'N'. 078: * 079: * RFP A RFP A 080: * 081: * 03 04 05 33 43 53 082: * 13 14 15 00 44 54 083: * 23 24 25 10 11 55 084: * 33 34 35 20 21 22 085: * 00 44 45 30 31 32 086: * 01 11 55 40 41 42 087: * 02 12 22 50 51 52 088: * 089: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 090: * transpose of RFP A above. One therefore gets: 091: * 092: * 093: * RFP A RFP A 094: * 095: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 096: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 097: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 098: * 099: * 100: * We first 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( 'STPTTF', -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 STPTTF 438: * 439: END 440: