001:       SUBROUTINE CTFTTP( TRANSR, UPLO, N, ARF, AP, 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:       COMPLEX            AP( 0: * ), ARF( 0: * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  CTFTTP 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: *          = 'C':  ARF is in Conjugate-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) COMPLEX 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) COMPLEX 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: *  Further Details
056: *  ===============
057: *
058: *  We first consider Standard Packed Format when N is even.
059: *  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: *  conjugate-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: *  conjugate-transpose of the last three columns of AP lower.
078: *  To denote conjugate we place -- above the element. This covers the
079: *  case N even and TRANSR = 'N'.
080: *
081: *         RFP A                   RFP A
082: *
083: *                                -- -- --
084: *        03 04 05                33 43 53
085: *                                   -- --
086: *        13 14 15                00 44 54
087: *                                      --
088: *        23 24 25                10 11 55
089: *
090: *        33 34 35                20 21 22
091: *        --
092: *        00 44 45                30 31 32
093: *        -- --
094: *        01 11 55                40 41 42
095: *        -- -- --
096: *        02 12 22                50 51 52
097: *
098: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
099: *  transpose of RFP A above. One therefore gets:
100: *
101: *
102: *           RFP A                   RFP A
103: *
104: *     -- -- -- --                -- -- -- -- -- --
105: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
106: *     -- -- -- -- --                -- -- -- -- --
107: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
108: *     -- -- -- -- -- --                -- -- -- --
109: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
110: *
111: *
112: *  We next  consider Standard Packed Format when N is odd.
113: *  We give an example where N = 5.
114: *
115: *     AP is Upper                 AP is Lower
116: *
117: *   00 01 02 03 04              00
118: *      11 12 13 14              10 11
119: *         22 23 24              20 21 22
120: *            33 34              30 31 32 33
121: *               44              40 41 42 43 44
122: *
123: *
124: *  Let TRANSR = 'N'. RFP holds AP as follows:
125: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
126: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
127: *  conjugate-transpose of the first two   columns of AP upper.
128: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
129: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
130: *  conjugate-transpose of the last two   columns of AP lower.
131: *  To denote conjugate we place -- above the element. This covers the
132: *  case N odd  and TRANSR = 'N'.
133: *
134: *         RFP A                   RFP A
135: *
136: *                                   -- --
137: *        02 03 04                00 33 43
138: *                                      --
139: *        12 13 14                10 11 44
140: *
141: *        22 23 24                20 21 22
142: *        --
143: *        00 33 34                30 31 32
144: *        -- --
145: *        01 11 44                40 41 42
146: *
147: *  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
148: *  transpose of RFP A above. One therefore gets:
149: *
150: *
151: *           RFP A                   RFP A
152: *
153: *     -- -- --                   -- -- -- -- -- --
154: *     02 12 22 00 01             00 10 20 30 40 50
155: *     -- -- -- --                   -- -- -- -- --
156: *     03 13 23 33 11             33 11 21 31 41 51
157: *     -- -- -- -- --                   -- -- -- --
158: *     04 14 24 34 44             43 44 22 32 42 52
159: *
160: *  =====================================================================
161: *
162: *     .. Parameters ..
163: *     ..
164: *     .. Local Scalars ..
165:       LOGICAL            LOWER, NISODD, NORMALTRANSR
166:       INTEGER            N1, N2, K, NT
167:       INTEGER            I, J, IJ
168:       INTEGER            IJP, JP, LDA, JS
169: *     ..
170: *     .. External Functions ..
171:       LOGICAL            LSAME
172:       EXTERNAL           LSAME
173: *     ..
174: *     .. External Subroutines ..
175:       EXTERNAL           XERBLA
176: *     ..
177: *     .. Intrinsic Functions ..
178:       INTRINSIC          CONJG
179: *     ..
180: *     .. Intrinsic Functions ..
181: *     ..
182: *     .. Executable Statements ..
183: *
184: *     Test the input parameters.
185: *
186:       INFO = 0
187:       NORMALTRANSR = LSAME( TRANSR, 'N' )
188:       LOWER = LSAME( UPLO, 'L' )
189:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
190:          INFO = -1
191:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
192:          INFO = -2
193:       ELSE IF( N.LT.0 ) THEN
194:          INFO = -3
195:       END IF
196:       IF( INFO.NE.0 ) THEN
197:          CALL XERBLA( 'CTFTTP', -INFO )
198:          RETURN
199:       END IF
200: *
201: *     Quick return if possible
202: *
203:       IF( N.EQ.0 )
204:      +   RETURN
205: *
206:       IF( N.EQ.1 ) THEN
207:          IF( NORMALTRANSR ) THEN
208:             AP( 0 ) = ARF( 0 )
209:          ELSE
210:             AP( 0 ) = CONJG( ARF( 0 ) )
211:          END IF
212:          RETURN
213:       END IF
214: *
215: *     Size of array ARF(0:NT-1)
216: *
217:       NT = N*( N+1 ) / 2
218: *
219: *     Set N1 and N2 depending on LOWER
220: *
221:       IF( LOWER ) THEN
222:          N2 = N / 2
223:          N1 = N - N2
224:       ELSE
225:          N1 = N / 2
226:          N2 = N - N1
227:       END IF
228: *
229: *     If N is odd, set NISODD = .TRUE.
230: *     If N is even, set K = N/2 and NISODD = .FALSE.
231: *
232: *     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
233: *     where noe = 0 if n is even, noe = 1 if n is odd
234: *
235:       IF( MOD( N, 2 ).EQ.0 ) THEN
236:          K = N / 2
237:          NISODD = .FALSE.
238:          LDA = N + 1
239:       ELSE
240:          NISODD = .TRUE.
241:          LDA = N
242:       END IF
243: *
244: *     ARF^C has lda rows and n+1-noe cols
245: *
246:       IF( .NOT.NORMALTRANSR )
247:      +   LDA = ( N+1 ) / 2
248: *
249: *     start execution: there are eight cases
250: *
251:       IF( NISODD ) THEN
252: *
253: *        N is odd
254: *
255:          IF( NORMALTRANSR ) THEN
256: *
257: *           N is odd and TRANSR = 'N'
258: *
259:             IF( LOWER ) THEN
260: *
261: *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
262: *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
263: *             T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
264: *
265:                IJP = 0
266:                JP = 0
267:                DO J = 0, N2
268:                   DO I = J, N - 1
269:                      IJ = I + JP
270:                      AP( IJP ) = ARF( IJ )
271:                      IJP = IJP + 1
272:                   END DO
273:                   JP = JP + LDA
274:                END DO
275:                DO I = 0, N2 - 1
276:                   DO J = 1 + I, N2
277:                      IJ = I + J*LDA
278:                      AP( IJP ) = CONJG( ARF( IJ ) )
279:                      IJP = IJP + 1
280:                   END DO
281:                END DO
282: *
283:             ELSE
284: *
285: *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
286: *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
287: *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
288: *
289:                IJP = 0
290:                DO J = 0, N1 - 1
291:                   IJ = N2 + J
292:                   DO I = 0, J
293:                      AP( IJP ) = CONJG( ARF( IJ ) )
294:                      IJP = IJP + 1
295:                      IJ = IJ + LDA
296:                   END DO
297:                END DO
298:                JS = 0
299:                DO J = N1, N - 1
300:                   IJ = JS
301:                   DO IJ = JS, JS + J
302:                      AP( IJP ) = ARF( IJ )
303:                      IJP = IJP + 1
304:                   END DO
305:                   JS = JS + LDA
306:                END DO
307: *
308:             END IF
309: *
310:          ELSE
311: *
312: *           N is odd and TRANSR = 'C'
313: *
314:             IF( LOWER ) THEN
315: *
316: *              SRPA for LOWER, TRANSPOSE and N is odd
317: *              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
318: *              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
319: *
320:                IJP = 0
321:                DO I = 0, N2
322:                   DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
323:                      AP( IJP ) = CONJG( ARF( IJ ) )
324:                      IJP = IJP + 1
325:                   END DO
326:                END DO
327:                JS = 1
328:                DO J = 0, N2 - 1
329:                   DO IJ = JS, JS + N2 - J - 1
330:                      AP( IJP ) = ARF( IJ )
331:                      IJP = IJP + 1
332:                   END DO
333:                   JS = JS + LDA + 1
334:                END DO
335: *
336:             ELSE
337: *
338: *              SRPA for UPPER, TRANSPOSE and N is odd
339: *              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
340: *              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
341: *
342:                IJP = 0
343:                JS = N2*LDA
344:                DO J = 0, N1 - 1
345:                   DO IJ = JS, JS + J
346:                      AP( IJP ) = ARF( IJ )
347:                      IJP = IJP + 1
348:                   END DO
349:                   JS = JS + LDA
350:                END DO
351:                DO I = 0, N1
352:                   DO IJ = I, I + ( N1+I )*LDA, LDA
353:                      AP( IJP ) = CONJG( ARF( IJ ) )
354:                      IJP = IJP + 1
355:                   END DO
356:                END DO
357: *
358:             END IF
359: *
360:          END IF
361: *
362:       ELSE
363: *
364: *        N is even
365: *
366:          IF( NORMALTRANSR ) THEN
367: *
368: *           N is even and TRANSR = 'N'
369: *
370:             IF( LOWER ) THEN
371: *
372: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
373: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
374: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
375: *
376:                IJP = 0
377:                JP = 0
378:                DO J = 0, K - 1
379:                   DO I = J, N - 1
380:                      IJ = 1 + I + JP
381:                      AP( IJP ) = ARF( IJ )
382:                      IJP = IJP + 1
383:                   END DO
384:                   JP = JP + LDA
385:                END DO
386:                DO I = 0, K - 1
387:                   DO J = I, K - 1
388:                      IJ = I + J*LDA
389:                      AP( IJP ) = CONJG( ARF( IJ ) )
390:                      IJP = IJP + 1
391:                   END DO
392:                END DO
393: *
394:             ELSE
395: *
396: *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
397: *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
398: *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
399: *
400:                IJP = 0
401:                DO J = 0, K - 1
402:                   IJ = K + 1 + J
403:                   DO I = 0, J
404:                      AP( IJP ) = CONJG( ARF( IJ ) )
405:                      IJP = IJP + 1
406:                      IJ = IJ + LDA
407:                   END DO
408:                END DO
409:                JS = 0
410:                DO J = K, N - 1
411:                   IJ = JS
412:                   DO IJ = JS, JS + J
413:                      AP( IJP ) = ARF( IJ )
414:                      IJP = IJP + 1
415:                   END DO
416:                   JS = JS + LDA
417:                END DO
418: *
419:             END IF
420: *
421:          ELSE
422: *
423: *           N is even and TRANSR = 'C'
424: *
425:             IF( LOWER ) THEN
426: *
427: *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
428: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
429: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
430: *
431:                IJP = 0
432:                DO I = 0, K - 1
433:                   DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
434:                      AP( IJP ) = CONJG( ARF( IJ ) )
435:                      IJP = IJP + 1
436:                   END DO
437:                END DO
438:                JS = 0
439:                DO J = 0, K - 1
440:                   DO IJ = JS, JS + K - J - 1
441:                      AP( IJP ) = ARF( IJ )
442:                      IJP = IJP + 1
443:                   END DO
444:                   JS = JS + LDA + 1
445:                END DO
446: *
447:             ELSE
448: *
449: *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
450: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
451: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
452: *
453:                IJP = 0
454:                JS = ( K+1 )*LDA
455:                DO J = 0, K - 1
456:                   DO IJ = JS, JS + J
457:                      AP( IJP ) = ARF( IJ )
458:                      IJP = IJP + 1
459:                   END DO
460:                   JS = JS + LDA
461:                END DO
462:                DO I = 0, K - 1
463:                   DO IJ = I, I + ( K+I )*LDA, LDA
464:                      AP( IJP ) = CONJG( ARF( IJ ) )
465:                      IJP = IJP + 1
466:                   END DO
467:                END DO
468: *
469:             END IF
470: *
471:          END IF
472: *
473:       END IF
474: *
475:       RETURN
476: *
477: *     End of CTFTTP
478: *
479:       END
480: