![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztfttp.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 ztfttp.f CVS log | Up to [local] / rpl / lapack / lapack |
1.10 bertrand 1: *> \brief \b ZTFTTP copies a triangular matrix from the rectangular full packed format (TF) to the standard packed format (TP).
1.7 bertrand 2: *
3: * =========== DOCUMENTATION ===========
1.1 bertrand 4: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
1.7 bertrand 8: *> \htmlonly
1.14 bertrand 9: *> Download ZTFTTP + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztfttp.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztfttp.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztfttp.f">
1.7 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.7 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
1.14 bertrand 22: *
1.7 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER TRANSR, UPLO
25: * INTEGER INFO, N
26: * ..
27: * .. Array Arguments ..
28: * COMPLEX*16 AP( 0: * ), ARF( 0: * )
29: * ..
1.14 bertrand 30: *
1.7 bertrand 31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> ZTFTTP copies a triangular matrix A from rectangular full packed
38: *> format (TF) to standard packed format (TP).
39: *> \endverbatim
40: *
41: * Arguments:
42: * ==========
43: *
44: *> \param[in] TRANSR
45: *> \verbatim
46: *> TRANSR is CHARACTER*1
47: *> = 'N': ARF is in Normal format;
48: *> = 'C': ARF is in Conjugate-transpose format;
49: *> \endverbatim
50: *>
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': A is upper triangular;
55: *> = 'L': A is lower triangular.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in] ARF
65: *> \verbatim
66: *> ARF is COMPLEX*16 array, dimension ( N*(N+1)/2 ),
67: *> On entry, the upper or lower triangular matrix A stored in
68: *> RFP format. For a further discussion see Notes below.
69: *> \endverbatim
70: *>
71: *> \param[out] AP
72: *> \verbatim
73: *> AP is COMPLEX*16 array, dimension ( N*(N+1)/2 ),
74: *> On exit, the upper or lower triangular matrix A, packed
75: *> columnwise in a linear array. The j-th column of A is stored
76: *> in the array AP as follows:
77: *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
78: *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
79: *> \endverbatim
80: *>
81: *> \param[out] INFO
82: *> \verbatim
83: *> INFO is INTEGER
84: *> = 0: successful exit
85: *> < 0: if INFO = -i, the i-th argument had an illegal value
86: *> \endverbatim
87: *
88: * Authors:
89: * ========
90: *
1.14 bertrand 91: *> \author Univ. of Tennessee
92: *> \author Univ. of California Berkeley
93: *> \author Univ. of Colorado Denver
94: *> \author NAG Ltd.
1.7 bertrand 95: *
96: *> \ingroup complex16OTHERcomputational
97: *
98: *> \par Further Details:
99: * =====================
100: *>
101: *> \verbatim
102: *>
103: *> We first consider Standard Packed Format when N is even.
104: *> We give an example where N = 6.
105: *>
106: *> AP is Upper AP is Lower
107: *>
108: *> 00 01 02 03 04 05 00
109: *> 11 12 13 14 15 10 11
110: *> 22 23 24 25 20 21 22
111: *> 33 34 35 30 31 32 33
112: *> 44 45 40 41 42 43 44
113: *> 55 50 51 52 53 54 55
114: *>
115: *>
116: *> Let TRANSR = 'N'. RFP holds AP as follows:
117: *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
118: *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
119: *> conjugate-transpose of the first three columns of AP upper.
120: *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
121: *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
122: *> conjugate-transpose of the last three columns of AP lower.
123: *> To denote conjugate we place -- above the element. This covers the
124: *> case N even and TRANSR = 'N'.
125: *>
126: *> RFP A RFP A
127: *>
128: *> -- -- --
129: *> 03 04 05 33 43 53
130: *> -- --
131: *> 13 14 15 00 44 54
132: *> --
133: *> 23 24 25 10 11 55
134: *>
135: *> 33 34 35 20 21 22
136: *> --
137: *> 00 44 45 30 31 32
138: *> -- --
139: *> 01 11 55 40 41 42
140: *> -- -- --
141: *> 02 12 22 50 51 52
142: *>
143: *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
144: *> transpose of RFP A above. One therefore gets:
145: *>
146: *>
147: *> RFP A RFP A
148: *>
149: *> -- -- -- -- -- -- -- -- -- --
150: *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
151: *> -- -- -- -- -- -- -- -- -- --
152: *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
153: *> -- -- -- -- -- -- -- -- -- --
154: *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
155: *>
156: *>
157: *> We next consider Standard Packed Format when N is odd.
158: *> We give an example where N = 5.
159: *>
160: *> AP is Upper AP is Lower
161: *>
162: *> 00 01 02 03 04 00
163: *> 11 12 13 14 10 11
164: *> 22 23 24 20 21 22
165: *> 33 34 30 31 32 33
166: *> 44 40 41 42 43 44
167: *>
168: *>
169: *> Let TRANSR = 'N'. RFP holds AP as follows:
170: *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
171: *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
172: *> conjugate-transpose of the first two columns of AP upper.
173: *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
174: *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
175: *> conjugate-transpose of the last two columns of AP lower.
176: *> To denote conjugate we place -- above the element. This covers the
177: *> case N odd and TRANSR = 'N'.
178: *>
179: *> RFP A RFP A
180: *>
181: *> -- --
182: *> 02 03 04 00 33 43
183: *> --
184: *> 12 13 14 10 11 44
185: *>
186: *> 22 23 24 20 21 22
187: *> --
188: *> 00 33 34 30 31 32
189: *> -- --
190: *> 01 11 44 40 41 42
191: *>
192: *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
193: *> transpose of RFP A above. One therefore gets:
194: *>
195: *>
196: *> RFP A RFP A
197: *>
198: *> -- -- -- -- -- -- -- -- --
199: *> 02 12 22 00 01 00 10 20 30 40 50
200: *> -- -- -- -- -- -- -- -- --
201: *> 03 13 23 33 11 33 11 21 31 41 51
202: *> -- -- -- -- -- -- -- -- --
203: *> 04 14 24 34 44 43 44 22 32 42 52
204: *> \endverbatim
205: *>
206: * =====================================================================
207: SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
1.1 bertrand 208: *
1.17 ! bertrand 209: * -- LAPACK computational routine --
1.1 bertrand 210: * -- LAPACK is a software package provided by Univ. of Tennessee, --
211: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212: *
213: * .. Scalar Arguments ..
214: CHARACTER TRANSR, UPLO
215: INTEGER INFO, N
216: * ..
217: * .. Array Arguments ..
218: COMPLEX*16 AP( 0: * ), ARF( 0: * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: * ..
225: * .. Local Scalars ..
226: LOGICAL LOWER, NISODD, NORMALTRANSR
227: INTEGER N1, N2, K, NT
228: INTEGER I, J, IJ
229: INTEGER IJP, JP, LDA, JS
230: * ..
231: * .. External Functions ..
232: LOGICAL LSAME
233: EXTERNAL LSAME
234: * ..
235: * .. External Subroutines ..
236: EXTERNAL XERBLA
237: * ..
238: * .. Intrinsic Functions ..
239: INTRINSIC DCONJG
240: * ..
241: * .. Intrinsic Functions ..
242: * ..
243: * .. Executable Statements ..
244: *
245: * Test the input parameters.
246: *
247: INFO = 0
248: NORMALTRANSR = LSAME( TRANSR, 'N' )
249: LOWER = LSAME( UPLO, 'L' )
250: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
251: INFO = -1
252: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
253: INFO = -2
254: ELSE IF( N.LT.0 ) THEN
255: INFO = -3
256: END IF
257: IF( INFO.NE.0 ) THEN
258: CALL XERBLA( 'ZTFTTP', -INFO )
259: RETURN
260: END IF
261: *
262: * Quick return if possible
263: *
264: IF( N.EQ.0 )
1.6 bertrand 265: $ RETURN
1.1 bertrand 266: *
267: IF( N.EQ.1 ) THEN
268: IF( NORMALTRANSR ) THEN
269: AP( 0 ) = ARF( 0 )
270: ELSE
271: AP( 0 ) = DCONJG( ARF( 0 ) )
272: END IF
273: RETURN
274: END IF
275: *
276: * Size of array ARF(0:NT-1)
277: *
278: NT = N*( N+1 ) / 2
279: *
280: * Set N1 and N2 depending on LOWER
281: *
282: IF( LOWER ) THEN
283: N2 = N / 2
284: N1 = N - N2
285: ELSE
286: N1 = N / 2
287: N2 = N - N1
288: END IF
289: *
290: * If N is odd, set NISODD = .TRUE.
291: * If N is even, set K = N/2 and NISODD = .FALSE.
292: *
293: * set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
294: * where noe = 0 if n is even, noe = 1 if n is odd
295: *
296: IF( MOD( N, 2 ).EQ.0 ) THEN
297: K = N / 2
298: NISODD = .FALSE.
299: LDA = N + 1
300: ELSE
301: NISODD = .TRUE.
302: LDA = N
303: END IF
304: *
305: * ARF^C has lda rows and n+1-noe cols
306: *
307: IF( .NOT.NORMALTRANSR )
1.6 bertrand 308: $ LDA = ( N+1 ) / 2
1.1 bertrand 309: *
310: * start execution: there are eight cases
311: *
312: IF( NISODD ) THEN
313: *
314: * N is odd
315: *
316: IF( NORMALTRANSR ) THEN
317: *
318: * N is odd and TRANSR = 'N'
319: *
320: IF( LOWER ) THEN
321: *
322: * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
323: * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
324: * T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
325: *
326: IJP = 0
327: JP = 0
328: DO J = 0, N2
329: DO I = J, N - 1
330: IJ = I + JP
331: AP( IJP ) = ARF( IJ )
332: IJP = IJP + 1
333: END DO
334: JP = JP + LDA
335: END DO
336: DO I = 0, N2 - 1
337: DO J = 1 + I, N2
338: IJ = I + J*LDA
339: AP( IJP ) = DCONJG( ARF( IJ ) )
340: IJP = IJP + 1
341: END DO
342: END DO
343: *
344: ELSE
345: *
346: * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
347: * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
348: * T1 -> a(n2), T2 -> a(n1), S -> a(0)
349: *
350: IJP = 0
351: DO J = 0, N1 - 1
352: IJ = N2 + J
353: DO I = 0, J
354: AP( IJP ) = DCONJG( ARF( IJ ) )
355: IJP = IJP + 1
356: IJ = IJ + LDA
357: END DO
358: END DO
359: JS = 0
360: DO J = N1, N - 1
361: IJ = JS
362: DO IJ = JS, JS + J
363: AP( IJP ) = ARF( IJ )
364: IJP = IJP + 1
365: END DO
366: JS = JS + LDA
367: END DO
368: *
369: END IF
370: *
371: ELSE
372: *
373: * N is odd and TRANSR = 'C'
374: *
375: IF( LOWER ) THEN
376: *
377: * SRPA for LOWER, TRANSPOSE and N is odd
378: * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
379: * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
380: *
381: IJP = 0
382: DO I = 0, N2
383: DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
384: AP( IJP ) = DCONJG( ARF( IJ ) )
385: IJP = IJP + 1
386: END DO
387: END DO
388: JS = 1
389: DO J = 0, N2 - 1
390: DO IJ = JS, JS + N2 - J - 1
391: AP( IJP ) = ARF( IJ )
392: IJP = IJP + 1
393: END DO
394: JS = JS + LDA + 1
395: END DO
396: *
397: ELSE
398: *
399: * SRPA for UPPER, TRANSPOSE and N is odd
400: * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
401: * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
402: *
403: IJP = 0
404: JS = N2*LDA
405: DO J = 0, N1 - 1
406: DO IJ = JS, JS + J
407: AP( IJP ) = ARF( IJ )
408: IJP = IJP + 1
409: END DO
410: JS = JS + LDA
411: END DO
412: DO I = 0, N1
413: DO IJ = I, I + ( N1+I )*LDA, LDA
414: AP( IJP ) = DCONJG( ARF( IJ ) )
415: IJP = IJP + 1
416: END DO
417: END DO
418: *
419: END IF
420: *
421: END IF
422: *
423: ELSE
424: *
425: * N is even
426: *
427: IF( NORMALTRANSR ) THEN
428: *
429: * N is even and TRANSR = 'N'
430: *
431: IF( LOWER ) THEN
432: *
433: * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
434: * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
435: * T1 -> a(1), T2 -> a(0), S -> a(k+1)
436: *
437: IJP = 0
438: JP = 0
439: DO J = 0, K - 1
440: DO I = J, N - 1
441: IJ = 1 + I + JP
442: AP( IJP ) = ARF( IJ )
443: IJP = IJP + 1
444: END DO
445: JP = JP + LDA
446: END DO
447: DO I = 0, K - 1
448: DO J = I, K - 1
449: IJ = I + J*LDA
450: AP( IJP ) = DCONJG( ARF( IJ ) )
451: IJP = IJP + 1
452: END DO
453: END DO
454: *
455: ELSE
456: *
457: * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
458: * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
459: * T1 -> a(k+1), T2 -> a(k), S -> a(0)
460: *
461: IJP = 0
462: DO J = 0, K - 1
463: IJ = K + 1 + J
464: DO I = 0, J
465: AP( IJP ) = DCONJG( ARF( IJ ) )
466: IJP = IJP + 1
467: IJ = IJ + LDA
468: END DO
469: END DO
470: JS = 0
471: DO J = K, N - 1
472: IJ = JS
473: DO IJ = JS, JS + J
474: AP( IJP ) = ARF( IJ )
475: IJP = IJP + 1
476: END DO
477: JS = JS + LDA
478: END DO
479: *
480: END IF
481: *
482: ELSE
483: *
484: * N is even and TRANSR = 'C'
485: *
486: IF( LOWER ) THEN
487: *
488: * SRPA for LOWER, TRANSPOSE and N is even (see paper)
489: * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
490: * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
491: *
492: IJP = 0
493: DO I = 0, K - 1
494: DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
495: AP( IJP ) = DCONJG( ARF( IJ ) )
496: IJP = IJP + 1
497: END DO
498: END DO
499: JS = 0
500: DO J = 0, K - 1
501: DO IJ = JS, JS + K - J - 1
502: AP( IJP ) = ARF( IJ )
503: IJP = IJP + 1
504: END DO
505: JS = JS + LDA + 1
506: END DO
507: *
508: ELSE
509: *
510: * SRPA for UPPER, TRANSPOSE and N is even (see paper)
511: * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
512: * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
513: *
514: IJP = 0
515: JS = ( K+1 )*LDA
516: DO J = 0, K - 1
517: DO IJ = JS, JS + J
518: AP( IJP ) = ARF( IJ )
519: IJP = IJP + 1
520: END DO
521: JS = JS + LDA
522: END DO
523: DO I = 0, K - 1
524: DO IJ = I, I + ( K+I )*LDA, LDA
525: AP( IJP ) = DCONJG( ARF( IJ ) )
526: IJP = IJP + 1
527: END DO
528: END DO
529: *
530: END IF
531: *
532: END IF
533: *
534: END IF
535: *
536: RETURN
537: *
538: * End of ZTFTTP
539: *
540: END