![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrttf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to ztrttf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZTRTTF( TRANSR, UPLO, N, A, LDA, ARF, INFO )
2: *
1.6 ! bertrand 3: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 4: *
5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
1.6 ! bertrand 6: * -- April 2011 --
1.1 bertrand 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: COMPLEX*16 A( 0: LDA-1, 0: * ), ARF( 0: * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZTRTTF copies a triangular matrix A from standard full format (TR)
23: * to rectangular full packed format (TF) .
24: *
25: * Arguments
26: * =========
27: *
1.4 bertrand 28: * TRANSR (input) CHARACTER*1
1.1 bertrand 29: * = 'N': ARF in Normal mode is wanted;
30: * = 'C': ARF in Conjugate Transpose mode is wanted;
31: *
1.4 bertrand 32: * UPLO (input) CHARACTER*1
1.1 bertrand 33: * = 'U': A is upper triangular;
34: * = 'L': A is lower triangular.
35: *
36: * N (input) INTEGER
37: * The order of the matrix A. N >= 0.
38: *
39: * A (input) COMPLEX*16 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) COMPLEX*16 array, dimension ( N*(N+1)/2 ),
52: * On exit, the upper or lower triangular matrix A stored in
53: * RFP format. For a further discussion see Notes below.
54: *
55: * INFO (output) INTEGER
56: * = 0: successful exit
57: * < 0: if INFO = -i, the i-th argument had an illegal value
58: *
59: * Further Details
60: * ===============
61: *
62: * We first consider Standard Packed Format when N is even.
63: * We give an example where N = 6.
64: *
65: * AP is Upper AP is Lower
66: *
67: * 00 01 02 03 04 05 00
68: * 11 12 13 14 15 10 11
69: * 22 23 24 25 20 21 22
70: * 33 34 35 30 31 32 33
71: * 44 45 40 41 42 43 44
72: * 55 50 51 52 53 54 55
73: *
74: *
1.6 ! bertrand 75: * Let TRANSR = 'N'. RFP holds AP as follows:
! 76: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
1.1 bertrand 77: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
78: * conjugate-transpose of the first three columns of AP upper.
1.6 ! bertrand 79: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
1.1 bertrand 80: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
81: * conjugate-transpose of the last three columns of AP lower.
82: * To denote conjugate we place -- above the element. This covers the
1.6 ! bertrand 83: * case N even and TRANSR = 'N'.
1.1 bertrand 84: *
85: * RFP A RFP A
86: *
87: * -- -- --
88: * 03 04 05 33 43 53
89: * -- --
90: * 13 14 15 00 44 54
91: * --
92: * 23 24 25 10 11 55
93: *
94: * 33 34 35 20 21 22
95: * --
96: * 00 44 45 30 31 32
97: * -- --
98: * 01 11 55 40 41 42
99: * -- -- --
100: * 02 12 22 50 51 52
101: *
1.6 ! bertrand 102: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
1.1 bertrand 103: * transpose of RFP A above. One therefore gets:
104: *
105: *
106: * RFP A RFP A
107: *
108: * -- -- -- -- -- -- -- -- -- --
109: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
110: * -- -- -- -- -- -- -- -- -- --
111: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
112: * -- -- -- -- -- -- -- -- -- --
113: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
114: *
115: *
116: * We next consider Standard Packed Format when N is odd.
117: * We give an example where N = 5.
118: *
119: * AP is Upper AP is Lower
120: *
121: * 00 01 02 03 04 00
122: * 11 12 13 14 10 11
123: * 22 23 24 20 21 22
124: * 33 34 30 31 32 33
125: * 44 40 41 42 43 44
126: *
127: *
1.6 ! bertrand 128: * Let TRANSR = 'N'. RFP holds AP as follows:
! 129: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
1.1 bertrand 130: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
131: * conjugate-transpose of the first two columns of AP upper.
1.6 ! bertrand 132: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
1.1 bertrand 133: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
134: * conjugate-transpose of the last two columns of AP lower.
135: * To denote conjugate we place -- above the element. This covers the
1.6 ! bertrand 136: * case N odd and TRANSR = 'N'.
1.1 bertrand 137: *
138: * RFP A RFP A
139: *
140: * -- --
141: * 02 03 04 00 33 43
142: * --
143: * 12 13 14 10 11 44
144: *
145: * 22 23 24 20 21 22
146: * --
147: * 00 33 34 30 31 32
148: * -- --
149: * 01 11 44 40 41 42
150: *
1.6 ! bertrand 151: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
1.1 bertrand 152: * transpose of RFP A above. One therefore gets:
153: *
154: *
155: * RFP A RFP A
156: *
157: * -- -- -- -- -- -- -- -- --
158: * 02 12 22 00 01 00 10 20 30 40 50
159: * -- -- -- -- -- -- -- -- --
160: * 03 13 23 33 11 33 11 21 31 41 51
161: * -- -- -- -- -- -- -- -- --
162: * 04 14 24 34 44 43 44 22 32 42 52
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: * ..
168: * .. Local Scalars ..
169: LOGICAL LOWER, NISODD, NORMALTRANSR
170: INTEGER I, IJ, J, K, L, N1, N2, NT, NX2, NP1X2
171: * ..
172: * .. External Functions ..
173: LOGICAL LSAME
174: EXTERNAL LSAME
175: * ..
176: * .. External Subroutines ..
177: EXTERNAL XERBLA
178: * ..
179: * .. Intrinsic Functions ..
180: INTRINSIC DCONJG, MAX, MOD
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: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
196: INFO = -5
197: END IF
198: IF( INFO.NE.0 ) THEN
199: CALL XERBLA( 'ZTRTTF', -INFO )
200: RETURN
201: END IF
202: *
203: * Quick return if possible
204: *
205: IF( N.LE.1 ) THEN
206: IF( N.EQ.1 ) THEN
207: IF( NORMALTRANSR ) THEN
208: ARF( 0 ) = A( 0, 0 )
209: ELSE
210: ARF( 0 ) = DCONJG( A( 0, 0 ) )
211: END IF
212: END IF
213: RETURN
214: END IF
215: *
216: * Size of array ARF(1:2,0:nt-1)
217: *
218: NT = N*( N+1 ) / 2
219: *
220: * set N1 and N2 depending on LOWER: for N even N1=N2=K
221: *
222: IF( LOWER ) THEN
223: N2 = N / 2
224: N1 = N - N2
225: ELSE
226: N1 = N / 2
227: N2 = N - N1
228: END IF
229: *
230: * If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2.
231: * If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is
232: * N--by--(N+1)/2.
233: *
234: IF( MOD( N, 2 ).EQ.0 ) THEN
235: K = N / 2
236: NISODD = .FALSE.
237: IF( .NOT.LOWER )
1.6 ! bertrand 238: $ NP1X2 = N + N + 2
1.1 bertrand 239: ELSE
240: NISODD = .TRUE.
241: IF( .NOT.LOWER )
1.6 ! bertrand 242: $ NX2 = N + N
1.1 bertrand 243: END IF
244: *
245: IF( NISODD ) THEN
246: *
247: * N is odd
248: *
249: IF( NORMALTRANSR ) THEN
250: *
251: * N is odd and TRANSR = 'N'
252: *
253: IF( LOWER ) THEN
254: *
255: * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
256: * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
257: * T1 -> a(0), T2 -> a(n), S -> a(n1); lda=n
258: *
259: IJ = 0
260: DO J = 0, N2
261: DO I = N1, N2 + J
262: ARF( IJ ) = DCONJG( A( N2+J, I ) )
263: IJ = IJ + 1
264: END DO
265: DO I = J, N - 1
266: ARF( IJ ) = A( I, J )
267: IJ = IJ + 1
268: END DO
269: END DO
270: *
271: ELSE
272: *
273: * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
274: * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
275: * T1 -> a(n2), T2 -> a(n1), S -> a(0); lda=n
276: *
277: IJ = NT - N
278: DO J = N - 1, N1, -1
279: DO I = 0, J
280: ARF( IJ ) = A( I, J )
281: IJ = IJ + 1
282: END DO
283: DO L = J - N1, N1 - 1
284: ARF( IJ ) = DCONJG( A( J-N1, L ) )
285: IJ = IJ + 1
286: END DO
287: IJ = IJ - NX2
288: END DO
289: *
290: END IF
291: *
292: ELSE
293: *
294: * N is odd and TRANSR = 'C'
295: *
296: IF( LOWER ) THEN
297: *
298: * SRPA for LOWER, TRANSPOSE and N is odd
299: * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
300: * T1 -> A(0+0) , T2 -> A(1+0) , S -> A(0+n1*n1); lda=n1
301: *
302: IJ = 0
303: DO J = 0, N2 - 1
304: DO I = 0, J
305: ARF( IJ ) = DCONJG( A( J, I ) )
306: IJ = IJ + 1
307: END DO
308: DO I = N1 + J, N - 1
309: ARF( IJ ) = A( I, N1+J )
310: IJ = IJ + 1
311: END DO
312: END DO
313: DO J = N2, N - 1
314: DO I = 0, N1 - 1
315: ARF( IJ ) = DCONJG( A( J, I ) )
316: IJ = IJ + 1
317: END DO
318: END DO
319: *
320: ELSE
321: *
322: * SRPA for UPPER, TRANSPOSE and N is odd
323: * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
324: * T1 -> A(n2*n2), T2 -> A(n1*n2), S -> A(0); lda=n2
325: *
326: IJ = 0
327: DO J = 0, N1
328: DO I = N1, N - 1
329: ARF( IJ ) = DCONJG( A( J, I ) )
330: IJ = IJ + 1
331: END DO
332: END DO
333: DO J = 0, N1 - 1
334: DO I = 0, J
335: ARF( IJ ) = A( I, J )
336: IJ = IJ + 1
337: END DO
338: DO L = N2 + J, N - 1
339: ARF( IJ ) = DCONJG( A( N2+J, L ) )
340: IJ = IJ + 1
341: END DO
342: END DO
343: *
344: END IF
345: *
346: END IF
347: *
348: ELSE
349: *
350: * N is even
351: *
352: IF( NORMALTRANSR ) THEN
353: *
354: * N is even and TRANSR = 'N'
355: *
356: IF( LOWER ) THEN
357: *
358: * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
359: * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
360: * T1 -> a(1), T2 -> a(0), S -> a(k+1); lda=n+1
361: *
362: IJ = 0
363: DO J = 0, K - 1
364: DO I = K, K + J
365: ARF( IJ ) = DCONJG( A( K+J, I ) )
366: IJ = IJ + 1
367: END DO
368: DO I = J, N - 1
369: ARF( IJ ) = A( I, J )
370: IJ = IJ + 1
371: END DO
372: END DO
373: *
374: ELSE
375: *
376: * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
377: * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
378: * T1 -> a(k+1), T2 -> a(k), S -> a(0); lda=n+1
379: *
380: IJ = NT - N - 1
381: DO J = N - 1, K, -1
382: DO I = 0, J
383: ARF( IJ ) = A( I, J )
384: IJ = IJ + 1
385: END DO
386: DO L = J - K, K - 1
387: ARF( IJ ) = DCONJG( A( J-K, L ) )
388: IJ = IJ + 1
389: END DO
390: IJ = IJ - NP1X2
391: END DO
392: *
393: END IF
394: *
395: ELSE
396: *
397: * N is even and TRANSR = 'C'
398: *
399: IF( LOWER ) THEN
400: *
401: * SRPA for LOWER, TRANSPOSE and N is even (see paper, A=B)
402: * T1 -> A(0,1) , T2 -> A(0,0) , S -> A(0,k+1) :
403: * T1 -> A(0+k) , T2 -> A(0+0) , S -> A(0+k*(k+1)); lda=k
404: *
405: IJ = 0
406: J = K
407: DO I = K, N - 1
408: ARF( IJ ) = A( I, J )
409: IJ = IJ + 1
410: END DO
411: DO J = 0, K - 2
412: DO I = 0, J
413: ARF( IJ ) = DCONJG( A( J, I ) )
414: IJ = IJ + 1
415: END DO
416: DO I = K + 1 + J, N - 1
417: ARF( IJ ) = A( I, K+1+J )
418: IJ = IJ + 1
419: END DO
420: END DO
421: DO J = K - 1, N - 1
422: DO I = 0, K - 1
423: ARF( IJ ) = DCONJG( A( J, I ) )
424: IJ = IJ + 1
425: END DO
426: END DO
427: *
428: ELSE
429: *
430: * SRPA for UPPER, TRANSPOSE and N is even (see paper, A=B)
431: * T1 -> A(0,k+1) , T2 -> A(0,k) , S -> A(0,0)
432: * T1 -> A(0+k*(k+1)) , T2 -> A(0+k*k) , S -> A(0+0)); lda=k
433: *
434: IJ = 0
435: DO J = 0, K
436: DO I = K, N - 1
437: ARF( IJ ) = DCONJG( A( J, I ) )
438: IJ = IJ + 1
439: END DO
440: END DO
441: DO J = 0, K - 2
442: DO I = 0, J
443: ARF( IJ ) = A( I, J )
444: IJ = IJ + 1
445: END DO
446: DO L = K + 1 + J, N - 1
447: ARF( IJ ) = DCONJG( A( K+1+J, L ) )
448: IJ = IJ + 1
449: END DO
450: END DO
451: *
452: * Note that here J = K-1
453: *
454: DO I = 0, J
455: ARF( IJ ) = A( I, J )
456: IJ = IJ + 1
457: END DO
458: *
459: END IF
460: *
461: END IF
462: *
463: END IF
464: *
465: RETURN
466: *
467: * End of ZTRTTF
468: *
469: END