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