1: SUBROUTINE DTRTTF( TRANSR, UPLO, N, A, LDA, ARF, 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: * DTRTTF copies a triangular matrix A from standard full format (TR)
23: * to rectangular full packed format (TF) .
24: *
25: * Arguments
26: * =========
27: *
28: * TRANSR (input) CHARACTER
29: * = 'N': ARF in Normal form is wanted;
30: * = 'T': ARF in Transpose form is wanted.
31: *
32: * UPLO (input) CHARACTER
33: * = 'U': Upper triangle of A is stored;
34: * = 'L': Lower triangle of A is stored.
35: *
36: * N (input) INTEGER
37: * The order of the matrix A. N >= 0.
38: *
39: * A (input) DOUBLE PRECISION array, dimension (LDA,N).
40: * On entry, the triangular matrix A. If UPLO = 'U', the
41: * leading N-by-N upper triangular part of the array A contains
42: * the upper triangular matrix, and the strictly lower
43: * triangular part of A is not referenced. If UPLO = 'L', the
44: * leading N-by-N lower triangular part of the array A contains
45: * the lower triangular matrix, and the strictly upper
46: * triangular part of A is not referenced.
47: *
48: * LDA (input) INTEGER
49: * The leading dimension of the matrix A. LDA >= max(1,N).
50: *
51: * ARF (output) DOUBLE PRECISION array, dimension (NT).
52: * NT=N*(N+1)/2. On exit, the triangular matrix A in RFP format.
53: *
54: * INFO (output) INTEGER
55: * = 0: successful exit
56: * < 0: if INFO = -i, the i-th argument had an illegal value
57: *
58: * Further Details
59: * ===============
60: *
61: * We first consider Rectangular Full Packed (RFP) Format when N is
62: * even. We give an example where N = 6.
63: *
64: * AP is Upper AP is Lower
65: *
66: * 00 01 02 03 04 05 00
67: * 11 12 13 14 15 10 11
68: * 22 23 24 25 20 21 22
69: * 33 34 35 30 31 32 33
70: * 44 45 40 41 42 43 44
71: * 55 50 51 52 53 54 55
72: *
73: *
74: * Let TRANSR = 'N'. RFP holds AP as follows:
75: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
76: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
77: * the transpose of the first three columns of AP upper.
78: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
79: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
80: * the transpose of the last three columns of AP lower.
81: * This covers the case N even and TRANSR = 'N'.
82: *
83: * RFP A RFP A
84: *
85: * 03 04 05 33 43 53
86: * 13 14 15 00 44 54
87: * 23 24 25 10 11 55
88: * 33 34 35 20 21 22
89: * 00 44 45 30 31 32
90: * 01 11 55 40 41 42
91: * 02 12 22 50 51 52
92: *
93: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
94: * transpose of RFP A above. One therefore gets:
95: *
96: *
97: * RFP A RFP A
98: *
99: * 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 then 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrttf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack