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