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