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