1: SUBROUTINE DTFTTR( TRANSR, UPLO, N, ARF, A, LDA, INFO )
2: *
3: * -- LAPACK routine (version 3.2.2) --
4: *
5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
6: * -- June 2010 --
7: *
8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10: *
11: * .. Scalar Arguments ..
12: CHARACTER TRANSR, UPLO
13: INTEGER INFO, N, LDA
14: * ..
15: * .. Array Arguments ..
16: DOUBLE PRECISION A( 0: LDA-1, 0: * ), ARF( 0: * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * DTFTTR copies a triangular matrix A from rectangular full packed
23: * format (TF) to standard full format (TR).
24: *
25: * Arguments
26: * =========
27: *
28: * TRANSR (input) CHARACTER
29: * = 'N': ARF is in Normal format;
30: * = 'T': ARF is in Transpose format.
31: *
32: * UPLO (input) CHARACTER
33: * = 'U': A is upper triangular;
34: * = 'L': A is lower triangular.
35: *
36: * N (input) INTEGER
37: * The order of the matrices ARF and A. N >= 0.
38: *
39: * ARF (input) DOUBLE PRECISION array, dimension (N*(N+1)/2).
40: * On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
41: * matrix A in RFP format. See the "Notes" below for more
42: * details.
43: *
44: * A (output) DOUBLE PRECISION array, dimension (LDA,N)
45: * On exit, the triangular matrix A. If UPLO = 'U', the
46: * leading N-by-N upper triangular part of the array A contains
47: * the upper triangular matrix, and the strictly lower
48: * triangular part of A is not referenced. If UPLO = 'L', the
49: * leading N-by-N lower triangular part of the array A contains
50: * the lower triangular matrix, and the strictly upper
51: * triangular part of A is not referenced.
52: *
53: * LDA (input) INTEGER
54: * The leading dimension of the array A. LDA >= max(1,N).
55: *
56: * INFO (output) INTEGER
57: * = 0: successful exit
58: * < 0: if INFO = -i, the i-th argument had an illegal value
59: *
60: * Further Details
61: * ===============
62: *
63: * We first consider Rectangular Full Packed (RFP) Format when N is
64: * even. We give an example where N = 6.
65: *
66: * AP is Upper AP is Lower
67: *
68: * 00 01 02 03 04 05 00
69: * 11 12 13 14 15 10 11
70: * 22 23 24 25 20 21 22
71: * 33 34 35 30 31 32 33
72: * 44 45 40 41 42 43 44
73: * 55 50 51 52 53 54 55
74: *
75: *
76: * Let TRANSR = 'N'. RFP holds AP as follows:
77: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
78: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
79: * the transpose of the first three columns of AP upper.
80: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
81: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
82: * the transpose of the last three columns of AP lower.
83: * This covers the case N even and TRANSR = 'N'.
84: *
85: * RFP A RFP A
86: *
87: * 03 04 05 33 43 53
88: * 13 14 15 00 44 54
89: * 23 24 25 10 11 55
90: * 33 34 35 20 21 22
91: * 00 44 45 30 31 32
92: * 01 11 55 40 41 42
93: * 02 12 22 50 51 52
94: *
95: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
96: * transpose of RFP A above. One therefore gets:
97: *
98: *
99: * RFP A RFP A
100: *
101: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
102: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
103: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
104: *
105: *
106: * We then consider Rectangular Full Packed (RFP) Format when N is
107: * odd. We give an example where N = 5.
108: *
109: * AP is Upper AP is Lower
110: *
111: * 00 01 02 03 04 00
112: * 11 12 13 14 10 11
113: * 22 23 24 20 21 22
114: * 33 34 30 31 32 33
115: * 44 40 41 42 43 44
116: *
117: *
118: * Let TRANSR = 'N'. RFP holds AP as follows:
119: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
120: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
121: * the transpose of the first two columns of AP upper.
122: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
123: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
124: * the transpose of the last two columns of AP lower.
125: * This covers the case N odd and TRANSR = 'N'.
126: *
127: * RFP A RFP A
128: *
129: * 02 03 04 00 33 43
130: * 12 13 14 10 11 44
131: * 22 23 24 20 21 22
132: * 00 33 34 30 31 32
133: * 01 11 44 40 41 42
134: *
135: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
136: * transpose of RFP A above. One therefore gets:
137: *
138: * RFP A RFP A
139: *
140: * 02 12 22 00 01 00 10 20 30 40 50
141: * 03 13 23 33 11 33 11 21 31 41 51
142: * 04 14 24 34 44 43 44 22 32 42 52
143: *
144: * Reference
145: * =========
146: *
147: * =====================================================================
148: *
149: * ..
150: * .. Local Scalars ..
151: LOGICAL LOWER, NISODD, NORMALTRANSR
152: INTEGER N1, N2, K, NT, NX2, NP1X2
153: INTEGER I, J, L, IJ
154: * ..
155: * .. External Functions ..
156: LOGICAL LSAME
157: EXTERNAL LSAME
158: * ..
159: * .. External Subroutines ..
160: EXTERNAL XERBLA
161: * ..
162: * .. Intrinsic Functions ..
163: INTRINSIC MAX, MOD
164: * ..
165: * .. Executable Statements ..
166: *
167: * Test the input parameters.
168: *
169: INFO = 0
170: NORMALTRANSR = LSAME( TRANSR, 'N' )
171: LOWER = LSAME( UPLO, 'L' )
172: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
173: INFO = -1
174: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
175: INFO = -2
176: ELSE IF( N.LT.0 ) THEN
177: INFO = -3
178: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
179: INFO = -6
180: END IF
181: IF( INFO.NE.0 ) THEN
182: CALL XERBLA( 'DTFTTR', -INFO )
183: RETURN
184: END IF
185: *
186: * Quick return if possible
187: *
188: IF( N.LE.1 ) THEN
189: IF( N.EQ.1 ) THEN
190: A( 0, 0 ) = ARF( 0 )
191: END IF
192: RETURN
193: END IF
194: *
195: * Size of array ARF(0:nt-1)
196: *
197: NT = N*( N+1 ) / 2
198: *
199: * set N1 and N2 depending on LOWER: for N even N1=N2=K
200: *
201: IF( LOWER ) THEN
202: N2 = N / 2
203: N1 = N - N2
204: ELSE
205: N1 = N / 2
206: N2 = N - N1
207: END IF
208: *
209: * If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2.
210: * If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is
211: * N--by--(N+1)/2.
212: *
213: IF( MOD( N, 2 ).EQ.0 ) THEN
214: K = N / 2
215: NISODD = .FALSE.
216: IF( .NOT.LOWER )
217: + NP1X2 = N + N + 2
218: ELSE
219: NISODD = .TRUE.
220: IF( .NOT.LOWER )
221: + NX2 = N + N
222: END IF
223: *
224: IF( NISODD ) THEN
225: *
226: * N is odd
227: *
228: IF( NORMALTRANSR ) THEN
229: *
230: * N is odd and TRANSR = 'N'
231: *
232: IF( LOWER ) THEN
233: *
234: * N is odd, TRANSR = 'N', and UPLO = 'L'
235: *
236: IJ = 0
237: DO J = 0, N2
238: DO I = N1, N2 + J
239: A( N2+J, I ) = ARF( IJ )
240: IJ = IJ + 1
241: END DO
242: DO I = J, N - 1
243: A( I, J ) = ARF( IJ )
244: IJ = IJ + 1
245: END DO
246: END DO
247: *
248: ELSE
249: *
250: * N is odd, TRANSR = 'N', and UPLO = 'U'
251: *
252: IJ = NT - N
253: DO J = N - 1, N1, -1
254: DO I = 0, J
255: A( I, J ) = ARF( IJ )
256: IJ = IJ + 1
257: END DO
258: DO L = J - N1, N1 - 1
259: A( J-N1, L ) = ARF( IJ )
260: IJ = IJ + 1
261: END DO
262: IJ = IJ - NX2
263: END DO
264: *
265: END IF
266: *
267: ELSE
268: *
269: * N is odd and TRANSR = 'T'
270: *
271: IF( LOWER ) THEN
272: *
273: * N is odd, TRANSR = 'T', and UPLO = 'L'
274: *
275: IJ = 0
276: DO J = 0, N2 - 1
277: DO I = 0, J
278: A( J, I ) = ARF( IJ )
279: IJ = IJ + 1
280: END DO
281: DO I = N1 + J, N - 1
282: A( I, N1+J ) = ARF( IJ )
283: IJ = IJ + 1
284: END DO
285: END DO
286: DO J = N2, N - 1
287: DO I = 0, N1 - 1
288: A( J, I ) = ARF( IJ )
289: IJ = IJ + 1
290: END DO
291: END DO
292: *
293: ELSE
294: *
295: * N is odd, TRANSR = 'T', and UPLO = 'U'
296: *
297: IJ = 0
298: DO J = 0, N1
299: DO I = N1, N - 1
300: A( J, I ) = ARF( IJ )
301: IJ = IJ + 1
302: END DO
303: END DO
304: DO J = 0, N1 - 1
305: DO I = 0, J
306: A( I, J ) = ARF( IJ )
307: IJ = IJ + 1
308: END DO
309: DO L = N2 + J, N - 1
310: A( N2+J, L ) = ARF( IJ )
311: IJ = IJ + 1
312: END DO
313: END DO
314: *
315: END IF
316: *
317: END IF
318: *
319: ELSE
320: *
321: * N is even
322: *
323: IF( NORMALTRANSR ) THEN
324: *
325: * N is even and TRANSR = 'N'
326: *
327: IF( LOWER ) THEN
328: *
329: * N is even, TRANSR = 'N', and UPLO = 'L'
330: *
331: IJ = 0
332: DO J = 0, K - 1
333: DO I = K, K + J
334: A( K+J, I ) = ARF( IJ )
335: IJ = IJ + 1
336: END DO
337: DO I = J, N - 1
338: A( I, J ) = ARF( IJ )
339: IJ = IJ + 1
340: END DO
341: END DO
342: *
343: ELSE
344: *
345: * N is even, TRANSR = 'N', and UPLO = 'U'
346: *
347: IJ = NT - N - 1
348: DO J = N - 1, K, -1
349: DO I = 0, J
350: A( I, J ) = ARF( IJ )
351: IJ = IJ + 1
352: END DO
353: DO L = J - K, K - 1
354: A( J-K, L ) = ARF( IJ )
355: IJ = IJ + 1
356: END DO
357: IJ = IJ - NP1X2
358: END DO
359: *
360: END IF
361: *
362: ELSE
363: *
364: * N is even and TRANSR = 'T'
365: *
366: IF( LOWER ) THEN
367: *
368: * N is even, TRANSR = 'T', and UPLO = 'L'
369: *
370: IJ = 0
371: J = K
372: DO I = K, N - 1
373: A( I, J ) = ARF( IJ )
374: IJ = IJ + 1
375: END DO
376: DO J = 0, K - 2
377: DO I = 0, J
378: A( J, I ) = ARF( IJ )
379: IJ = IJ + 1
380: END DO
381: DO I = K + 1 + J, N - 1
382: A( I, K+1+J ) = ARF( IJ )
383: IJ = IJ + 1
384: END DO
385: END DO
386: DO J = K - 1, N - 1
387: DO I = 0, K - 1
388: A( J, I ) = ARF( IJ )
389: IJ = IJ + 1
390: END DO
391: END DO
392: *
393: ELSE
394: *
395: * N is even, TRANSR = 'T', and UPLO = 'U'
396: *
397: IJ = 0
398: DO J = 0, K
399: DO I = K, N - 1
400: A( J, I ) = ARF( IJ )
401: IJ = IJ + 1
402: END DO
403: END DO
404: DO J = 0, K - 2
405: DO I = 0, J
406: A( I, J ) = ARF( IJ )
407: IJ = IJ + 1
408: END DO
409: DO L = K + 1 + J, N - 1
410: A( K+1+J, L ) = ARF( IJ )
411: IJ = IJ + 1
412: END DO
413: END DO
414: * Note that here, on exit of the loop, J = K-1
415: DO I = 0, J
416: A( I, J ) = ARF( IJ )
417: IJ = IJ + 1
418: END DO
419: *
420: END IF
421: *
422: END IF
423: *
424: END IF
425: *
426: RETURN
427: *
428: * End of DTFTTR
429: *
430: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtfttr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack