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