001:       SUBROUTINE STPTTF( TRANSR, UPLO, N, AP, ARF, INFO )
002: *
003: *  -- LAPACK routine (version 3.2.1)                                    --
004: *
005: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
006: *  -- April 2009                                                      --
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: *  Further Details
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: