1 SUBROUTINE ZPFTRF( TRANSR, UPLO, N, A, 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 * ..
12 * .. Scalar Arguments ..
13 CHARACTER TRANSR, UPLO
14 INTEGER N, INFO
15 * ..
16 * .. Array Arguments ..
17 COMPLEX*16 A( 0: * )
18 *
19 * Purpose
20 * =======
21 *
22 * ZPFTRF computes the Cholesky factorization of a complex Hermitian
23 * positive definite matrix A.
24 *
25 * The factorization has the form
26 * A = U**H * U, if UPLO = 'U', or
27 * A = L * L**H, if UPLO = 'L',
28 * where U is an upper triangular matrix and L is lower triangular.
29 *
30 * This is the block version of the algorithm, calling Level 3 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * TRANSR (input) CHARACTER*1
36 * = 'N': The Normal TRANSR of RFP A is stored;
37 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of RFP A is stored;
41 * = 'L': Lower triangle of RFP A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * A (input/output) COMPLEX array, dimension ( N*(N+1)/2 );
47 * On entry, the Hermitian matrix A in RFP format. RFP format is
48 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
49 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
50 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
51 * the Conjugate-transpose of RFP A as defined when
52 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
53 * follows: If UPLO = 'U' the RFP A contains the nt elements of
54 * upper packed A. If UPLO = 'L' the RFP A contains the elements
55 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
56 * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
57 * is odd. See the Note below for more details.
58 *
59 * On exit, if INFO = 0, the factor U or L from the Cholesky
60 * factorization RFP A = U**H*U or RFP A = L*L**H.
61 *
62 * INFO (output) INTEGER
63 * = 0: successful exit
64 * < 0: if INFO = -i, the i-th argument had an illegal value
65 * > 0: if INFO = i, the leading minor of order i is not
66 * positive definite, and the factorization could not be
67 * completed.
68 *
69 * Further Notes on RFP Format:
70 * ============================
71 *
72 * We first consider Standard Packed Format when N is even.
73 * We give an example where N = 6.
74 *
75 * AP is Upper AP is Lower
76 *
77 * 00 01 02 03 04 05 00
78 * 11 12 13 14 15 10 11
79 * 22 23 24 25 20 21 22
80 * 33 34 35 30 31 32 33
81 * 44 45 40 41 42 43 44
82 * 55 50 51 52 53 54 55
83 *
84 *
85 * Let TRANSR = 'N'. RFP holds AP as follows:
86 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
87 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
88 * conjugate-transpose of the first three columns of AP upper.
89 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
90 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
91 * conjugate-transpose of the last three columns of AP lower.
92 * To denote conjugate we place -- above the element. This covers the
93 * case N even and TRANSR = 'N'.
94 *
95 * RFP A RFP A
96 *
97 * -- -- --
98 * 03 04 05 33 43 53
99 * -- --
100 * 13 14 15 00 44 54
101 * --
102 * 23 24 25 10 11 55
103 *
104 * 33 34 35 20 21 22
105 * --
106 * 00 44 45 30 31 32
107 * -- --
108 * 01 11 55 40 41 42
109 * -- -- --
110 * 02 12 22 50 51 52
111 *
112 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
113 * transpose of RFP A above. One therefore gets:
114 *
115 *
116 * RFP A RFP A
117 *
118 * -- -- -- -- -- -- -- -- -- --
119 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
120 * -- -- -- -- -- -- -- -- -- --
121 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
122 * -- -- -- -- -- -- -- -- -- --
123 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
124 *
125 *
126 * We next consider Standard Packed Format when N is odd.
127 * We give an example where N = 5.
128 *
129 * AP is Upper AP is Lower
130 *
131 * 00 01 02 03 04 00
132 * 11 12 13 14 10 11
133 * 22 23 24 20 21 22
134 * 33 34 30 31 32 33
135 * 44 40 41 42 43 44
136 *
137 *
138 * Let TRANSR = 'N'. RFP holds AP as follows:
139 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
140 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
141 * conjugate-transpose of the first two columns of AP upper.
142 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
143 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
144 * conjugate-transpose of the last two columns of AP lower.
145 * To denote conjugate we place -- above the element. This covers the
146 * case N odd and TRANSR = 'N'.
147 *
148 * RFP A RFP A
149 *
150 * -- --
151 * 02 03 04 00 33 43
152 * --
153 * 12 13 14 10 11 44
154 *
155 * 22 23 24 20 21 22
156 * --
157 * 00 33 34 30 31 32
158 * -- --
159 * 01 11 44 40 41 42
160 *
161 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
162 * transpose of RFP A above. One therefore gets:
163 *
164 *
165 * RFP A RFP A
166 *
167 * -- -- -- -- -- -- -- -- --
168 * 02 12 22 00 01 00 10 20 30 40 50
169 * -- -- -- -- -- -- -- -- --
170 * 03 13 23 33 11 33 11 21 31 41 51
171 * -- -- -- -- -- -- -- -- --
172 * 04 14 24 34 44 43 44 22 32 42 52
173 *
174 * =====================================================================
175 *
176 * .. Parameters ..
177 DOUBLE PRECISION ONE
178 COMPLEX*16 CONE
179 PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ) )
180 * ..
181 * .. Local Scalars ..
182 LOGICAL LOWER, NISODD, NORMALTRANSR
183 INTEGER N1, N2, K
184 * ..
185 * .. External Functions ..
186 LOGICAL LSAME
187 EXTERNAL LSAME
188 * ..
189 * .. External Subroutines ..
190 EXTERNAL XERBLA, ZHERK, ZPOTRF, ZTRSM
191 * ..
192 * .. Intrinsic Functions ..
193 INTRINSIC MOD
194 * ..
195 * .. Executable Statements ..
196 *
197 * Test the input parameters.
198 *
199 INFO = 0
200 NORMALTRANSR = LSAME( TRANSR, 'N' )
201 LOWER = LSAME( UPLO, 'L' )
202 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
203 INFO = -1
204 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
205 INFO = -2
206 ELSE IF( N.LT.0 ) THEN
207 INFO = -3
208 END IF
209 IF( INFO.NE.0 ) THEN
210 CALL XERBLA( 'ZPFTRF', -INFO )
211 RETURN
212 END IF
213 *
214 * Quick return if possible
215 *
216 IF( N.EQ.0 )
217 $ RETURN
218 *
219 * If N is odd, set NISODD = .TRUE.
220 * If N is even, set K = N/2 and NISODD = .FALSE.
221 *
222 IF( MOD( N, 2 ).EQ.0 ) THEN
223 K = N / 2
224 NISODD = .FALSE.
225 ELSE
226 NISODD = .TRUE.
227 END IF
228 *
229 * Set N1 and N2 depending on LOWER
230 *
231 IF( LOWER ) THEN
232 N2 = N / 2
233 N1 = N - N2
234 ELSE
235 N1 = N / 2
236 N2 = N - N1
237 END IF
238 *
239 * start execution: there are eight cases
240 *
241 IF( NISODD ) THEN
242 *
243 * N is odd
244 *
245 IF( NORMALTRANSR ) THEN
246 *
247 * N is odd and TRANSR = 'N'
248 *
249 IF( LOWER ) THEN
250 *
251 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
252 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
253 * T1 -> a(0), T2 -> a(n), S -> a(n1)
254 *
255 CALL ZPOTRF( 'L', N1, A( 0 ), N, INFO )
256 IF( INFO.GT.0 )
257 $ RETURN
258 CALL ZTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, A( 0 ), N,
259 $ A( N1 ), N )
260 CALL ZHERK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE,
261 $ A( N ), N )
262 CALL ZPOTRF( 'U', N2, A( N ), N, INFO )
263 IF( INFO.GT.0 )
264 $ INFO = INFO + N1
265 *
266 ELSE
267 *
268 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
269 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
270 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
271 *
272 CALL ZPOTRF( 'L', N1, A( N2 ), N, INFO )
273 IF( INFO.GT.0 )
274 $ RETURN
275 CALL ZTRSM( 'L', 'L', 'N', 'N', N1, N2, CONE, A( N2 ), N,
276 $ A( 0 ), N )
277 CALL ZHERK( 'U', 'C', N2, N1, -ONE, A( 0 ), N, ONE,
278 $ A( N1 ), N )
279 CALL ZPOTRF( 'U', N2, A( N1 ), N, INFO )
280 IF( INFO.GT.0 )
281 $ INFO = INFO + N1
282 *
283 END IF
284 *
285 ELSE
286 *
287 * N is odd and TRANSR = 'C'
288 *
289 IF( LOWER ) THEN
290 *
291 * SRPA for LOWER, TRANSPOSE and N is odd
292 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
293 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
294 *
295 CALL ZPOTRF( 'U', N1, A( 0 ), N1, INFO )
296 IF( INFO.GT.0 )
297 $ RETURN
298 CALL ZTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, A( 0 ), N1,
299 $ A( N1*N1 ), N1 )
300 CALL ZHERK( 'L', 'C', N2, N1, -ONE, A( N1*N1 ), N1, ONE,
301 $ A( 1 ), N1 )
302 CALL ZPOTRF( 'L', N2, A( 1 ), N1, INFO )
303 IF( INFO.GT.0 )
304 $ INFO = INFO + N1
305 *
306 ELSE
307 *
308 * SRPA for UPPER, TRANSPOSE and N is odd
309 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
310 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
311 *
312 CALL ZPOTRF( 'U', N1, A( N2*N2 ), N2, INFO )
313 IF( INFO.GT.0 )
314 $ RETURN
315 CALL ZTRSM( 'R', 'U', 'N', 'N', N2, N1, CONE, A( N2*N2 ),
316 $ N2, A( 0 ), N2 )
317 CALL ZHERK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE,
318 $ A( N1*N2 ), N2 )
319 CALL ZPOTRF( 'L', N2, A( N1*N2 ), N2, INFO )
320 IF( INFO.GT.0 )
321 $ INFO = INFO + N1
322 *
323 END IF
324 *
325 END IF
326 *
327 ELSE
328 *
329 * N is even
330 *
331 IF( NORMALTRANSR ) THEN
332 *
333 * N is even and TRANSR = 'N'
334 *
335 IF( LOWER ) THEN
336 *
337 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
338 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
339 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
340 *
341 CALL ZPOTRF( 'L', K, A( 1 ), N+1, INFO )
342 IF( INFO.GT.0 )
343 $ RETURN
344 CALL ZTRSM( 'R', 'L', 'C', 'N', K, K, CONE, A( 1 ), N+1,
345 $ A( K+1 ), N+1 )
346 CALL ZHERK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE,
347 $ A( 0 ), N+1 )
348 CALL ZPOTRF( 'U', K, A( 0 ), N+1, INFO )
349 IF( INFO.GT.0 )
350 $ INFO = INFO + K
351 *
352 ELSE
353 *
354 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
355 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
356 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
357 *
358 CALL ZPOTRF( 'L', K, A( K+1 ), N+1, INFO )
359 IF( INFO.GT.0 )
360 $ RETURN
361 CALL ZTRSM( 'L', 'L', 'N', 'N', K, K, CONE, A( K+1 ),
362 $ N+1, A( 0 ), N+1 )
363 CALL ZHERK( 'U', 'C', K, K, -ONE, A( 0 ), N+1, ONE,
364 $ A( K ), N+1 )
365 CALL ZPOTRF( 'U', K, A( K ), N+1, INFO )
366 IF( INFO.GT.0 )
367 $ INFO = INFO + K
368 *
369 END IF
370 *
371 ELSE
372 *
373 * N is even and TRANSR = 'C'
374 *
375 IF( LOWER ) THEN
376 *
377 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
378 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
379 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
380 *
381 CALL ZPOTRF( 'U', K, A( 0+K ), K, INFO )
382 IF( INFO.GT.0 )
383 $ RETURN
384 CALL ZTRSM( 'L', 'U', 'C', 'N', K, K, CONE, A( K ), N1,
385 $ A( K*( K+1 ) ), K )
386 CALL ZHERK( 'L', 'C', K, K, -ONE, A( K*( K+1 ) ), K, ONE,
387 $ A( 0 ), K )
388 CALL ZPOTRF( 'L', K, A( 0 ), K, INFO )
389 IF( INFO.GT.0 )
390 $ INFO = INFO + K
391 *
392 ELSE
393 *
394 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
395 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
396 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
397 *
398 CALL ZPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO )
399 IF( INFO.GT.0 )
400 $ RETURN
401 CALL ZTRSM( 'R', 'U', 'N', 'N', K, K, CONE,
402 $ A( K*( K+1 ) ), K, A( 0 ), K )
403 CALL ZHERK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE,
404 $ A( K*K ), K )
405 CALL ZPOTRF( 'L', K, A( K*K ), K, INFO )
406 IF( INFO.GT.0 )
407 $ INFO = INFO + K
408 *
409 END IF
410 *
411 END IF
412 *
413 END IF
414 *
415 RETURN
416 *
417 * End of ZPFTRF
418 *
419 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 * ..
12 * .. Scalar Arguments ..
13 CHARACTER TRANSR, UPLO
14 INTEGER N, INFO
15 * ..
16 * .. Array Arguments ..
17 COMPLEX*16 A( 0: * )
18 *
19 * Purpose
20 * =======
21 *
22 * ZPFTRF computes the Cholesky factorization of a complex Hermitian
23 * positive definite matrix A.
24 *
25 * The factorization has the form
26 * A = U**H * U, if UPLO = 'U', or
27 * A = L * L**H, if UPLO = 'L',
28 * where U is an upper triangular matrix and L is lower triangular.
29 *
30 * This is the block version of the algorithm, calling Level 3 BLAS.
31 *
32 * Arguments
33 * =========
34 *
35 * TRANSR (input) CHARACTER*1
36 * = 'N': The Normal TRANSR of RFP A is stored;
37 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of RFP A is stored;
41 * = 'L': Lower triangle of RFP A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * A (input/output) COMPLEX array, dimension ( N*(N+1)/2 );
47 * On entry, the Hermitian matrix A in RFP format. RFP format is
48 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
49 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
50 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is
51 * the Conjugate-transpose of RFP A as defined when
52 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as
53 * follows: If UPLO = 'U' the RFP A contains the nt elements of
54 * upper packed A. If UPLO = 'L' the RFP A contains the elements
55 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
56 * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N
57 * is odd. See the Note below for more details.
58 *
59 * On exit, if INFO = 0, the factor U or L from the Cholesky
60 * factorization RFP A = U**H*U or RFP A = L*L**H.
61 *
62 * INFO (output) INTEGER
63 * = 0: successful exit
64 * < 0: if INFO = -i, the i-th argument had an illegal value
65 * > 0: if INFO = i, the leading minor of order i is not
66 * positive definite, and the factorization could not be
67 * completed.
68 *
69 * Further Notes on RFP Format:
70 * ============================
71 *
72 * We first consider Standard Packed Format when N is even.
73 * We give an example where N = 6.
74 *
75 * AP is Upper AP is Lower
76 *
77 * 00 01 02 03 04 05 00
78 * 11 12 13 14 15 10 11
79 * 22 23 24 25 20 21 22
80 * 33 34 35 30 31 32 33
81 * 44 45 40 41 42 43 44
82 * 55 50 51 52 53 54 55
83 *
84 *
85 * Let TRANSR = 'N'. RFP holds AP as follows:
86 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
87 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
88 * conjugate-transpose of the first three columns of AP upper.
89 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
90 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
91 * conjugate-transpose of the last three columns of AP lower.
92 * To denote conjugate we place -- above the element. This covers the
93 * case N even and TRANSR = 'N'.
94 *
95 * RFP A RFP A
96 *
97 * -- -- --
98 * 03 04 05 33 43 53
99 * -- --
100 * 13 14 15 00 44 54
101 * --
102 * 23 24 25 10 11 55
103 *
104 * 33 34 35 20 21 22
105 * --
106 * 00 44 45 30 31 32
107 * -- --
108 * 01 11 55 40 41 42
109 * -- -- --
110 * 02 12 22 50 51 52
111 *
112 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
113 * transpose of RFP A above. One therefore gets:
114 *
115 *
116 * RFP A RFP A
117 *
118 * -- -- -- -- -- -- -- -- -- --
119 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
120 * -- -- -- -- -- -- -- -- -- --
121 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
122 * -- -- -- -- -- -- -- -- -- --
123 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
124 *
125 *
126 * We next consider Standard Packed Format when N is odd.
127 * We give an example where N = 5.
128 *
129 * AP is Upper AP is Lower
130 *
131 * 00 01 02 03 04 00
132 * 11 12 13 14 10 11
133 * 22 23 24 20 21 22
134 * 33 34 30 31 32 33
135 * 44 40 41 42 43 44
136 *
137 *
138 * Let TRANSR = 'N'. RFP holds AP as follows:
139 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
140 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
141 * conjugate-transpose of the first two columns of AP upper.
142 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
143 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
144 * conjugate-transpose of the last two columns of AP lower.
145 * To denote conjugate we place -- above the element. This covers the
146 * case N odd and TRANSR = 'N'.
147 *
148 * RFP A RFP A
149 *
150 * -- --
151 * 02 03 04 00 33 43
152 * --
153 * 12 13 14 10 11 44
154 *
155 * 22 23 24 20 21 22
156 * --
157 * 00 33 34 30 31 32
158 * -- --
159 * 01 11 44 40 41 42
160 *
161 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
162 * transpose of RFP A above. One therefore gets:
163 *
164 *
165 * RFP A RFP A
166 *
167 * -- -- -- -- -- -- -- -- --
168 * 02 12 22 00 01 00 10 20 30 40 50
169 * -- -- -- -- -- -- -- -- --
170 * 03 13 23 33 11 33 11 21 31 41 51
171 * -- -- -- -- -- -- -- -- --
172 * 04 14 24 34 44 43 44 22 32 42 52
173 *
174 * =====================================================================
175 *
176 * .. Parameters ..
177 DOUBLE PRECISION ONE
178 COMPLEX*16 CONE
179 PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ) )
180 * ..
181 * .. Local Scalars ..
182 LOGICAL LOWER, NISODD, NORMALTRANSR
183 INTEGER N1, N2, K
184 * ..
185 * .. External Functions ..
186 LOGICAL LSAME
187 EXTERNAL LSAME
188 * ..
189 * .. External Subroutines ..
190 EXTERNAL XERBLA, ZHERK, ZPOTRF, ZTRSM
191 * ..
192 * .. Intrinsic Functions ..
193 INTRINSIC MOD
194 * ..
195 * .. Executable Statements ..
196 *
197 * Test the input parameters.
198 *
199 INFO = 0
200 NORMALTRANSR = LSAME( TRANSR, 'N' )
201 LOWER = LSAME( UPLO, 'L' )
202 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
203 INFO = -1
204 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
205 INFO = -2
206 ELSE IF( N.LT.0 ) THEN
207 INFO = -3
208 END IF
209 IF( INFO.NE.0 ) THEN
210 CALL XERBLA( 'ZPFTRF', -INFO )
211 RETURN
212 END IF
213 *
214 * Quick return if possible
215 *
216 IF( N.EQ.0 )
217 $ RETURN
218 *
219 * If N is odd, set NISODD = .TRUE.
220 * If N is even, set K = N/2 and NISODD = .FALSE.
221 *
222 IF( MOD( N, 2 ).EQ.0 ) THEN
223 K = N / 2
224 NISODD = .FALSE.
225 ELSE
226 NISODD = .TRUE.
227 END IF
228 *
229 * Set N1 and N2 depending on LOWER
230 *
231 IF( LOWER ) THEN
232 N2 = N / 2
233 N1 = N - N2
234 ELSE
235 N1 = N / 2
236 N2 = N - N1
237 END IF
238 *
239 * start execution: there are eight cases
240 *
241 IF( NISODD ) THEN
242 *
243 * N is odd
244 *
245 IF( NORMALTRANSR ) THEN
246 *
247 * N is odd and TRANSR = 'N'
248 *
249 IF( LOWER ) THEN
250 *
251 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
252 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
253 * T1 -> a(0), T2 -> a(n), S -> a(n1)
254 *
255 CALL ZPOTRF( 'L', N1, A( 0 ), N, INFO )
256 IF( INFO.GT.0 )
257 $ RETURN
258 CALL ZTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, A( 0 ), N,
259 $ A( N1 ), N )
260 CALL ZHERK( 'U', 'N', N2, N1, -ONE, A( N1 ), N, ONE,
261 $ A( N ), N )
262 CALL ZPOTRF( 'U', N2, A( N ), N, INFO )
263 IF( INFO.GT.0 )
264 $ INFO = INFO + N1
265 *
266 ELSE
267 *
268 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
269 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
270 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
271 *
272 CALL ZPOTRF( 'L', N1, A( N2 ), N, INFO )
273 IF( INFO.GT.0 )
274 $ RETURN
275 CALL ZTRSM( 'L', 'L', 'N', 'N', N1, N2, CONE, A( N2 ), N,
276 $ A( 0 ), N )
277 CALL ZHERK( 'U', 'C', N2, N1, -ONE, A( 0 ), N, ONE,
278 $ A( N1 ), N )
279 CALL ZPOTRF( 'U', N2, A( N1 ), N, INFO )
280 IF( INFO.GT.0 )
281 $ INFO = INFO + N1
282 *
283 END IF
284 *
285 ELSE
286 *
287 * N is odd and TRANSR = 'C'
288 *
289 IF( LOWER ) THEN
290 *
291 * SRPA for LOWER, TRANSPOSE and N is odd
292 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
293 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
294 *
295 CALL ZPOTRF( 'U', N1, A( 0 ), N1, INFO )
296 IF( INFO.GT.0 )
297 $ RETURN
298 CALL ZTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, A( 0 ), N1,
299 $ A( N1*N1 ), N1 )
300 CALL ZHERK( 'L', 'C', N2, N1, -ONE, A( N1*N1 ), N1, ONE,
301 $ A( 1 ), N1 )
302 CALL ZPOTRF( 'L', N2, A( 1 ), N1, INFO )
303 IF( INFO.GT.0 )
304 $ INFO = INFO + N1
305 *
306 ELSE
307 *
308 * SRPA for UPPER, TRANSPOSE and N is odd
309 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
310 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
311 *
312 CALL ZPOTRF( 'U', N1, A( N2*N2 ), N2, INFO )
313 IF( INFO.GT.0 )
314 $ RETURN
315 CALL ZTRSM( 'R', 'U', 'N', 'N', N2, N1, CONE, A( N2*N2 ),
316 $ N2, A( 0 ), N2 )
317 CALL ZHERK( 'L', 'N', N2, N1, -ONE, A( 0 ), N2, ONE,
318 $ A( N1*N2 ), N2 )
319 CALL ZPOTRF( 'L', N2, A( N1*N2 ), N2, INFO )
320 IF( INFO.GT.0 )
321 $ INFO = INFO + N1
322 *
323 END IF
324 *
325 END IF
326 *
327 ELSE
328 *
329 * N is even
330 *
331 IF( NORMALTRANSR ) THEN
332 *
333 * N is even and TRANSR = 'N'
334 *
335 IF( LOWER ) THEN
336 *
337 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
338 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
339 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
340 *
341 CALL ZPOTRF( 'L', K, A( 1 ), N+1, INFO )
342 IF( INFO.GT.0 )
343 $ RETURN
344 CALL ZTRSM( 'R', 'L', 'C', 'N', K, K, CONE, A( 1 ), N+1,
345 $ A( K+1 ), N+1 )
346 CALL ZHERK( 'U', 'N', K, K, -ONE, A( K+1 ), N+1, ONE,
347 $ A( 0 ), N+1 )
348 CALL ZPOTRF( 'U', K, A( 0 ), N+1, INFO )
349 IF( INFO.GT.0 )
350 $ INFO = INFO + K
351 *
352 ELSE
353 *
354 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
355 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
356 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
357 *
358 CALL ZPOTRF( 'L', K, A( K+1 ), N+1, INFO )
359 IF( INFO.GT.0 )
360 $ RETURN
361 CALL ZTRSM( 'L', 'L', 'N', 'N', K, K, CONE, A( K+1 ),
362 $ N+1, A( 0 ), N+1 )
363 CALL ZHERK( 'U', 'C', K, K, -ONE, A( 0 ), N+1, ONE,
364 $ A( K ), N+1 )
365 CALL ZPOTRF( 'U', K, A( K ), N+1, INFO )
366 IF( INFO.GT.0 )
367 $ INFO = INFO + K
368 *
369 END IF
370 *
371 ELSE
372 *
373 * N is even and TRANSR = 'C'
374 *
375 IF( LOWER ) THEN
376 *
377 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
378 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
379 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
380 *
381 CALL ZPOTRF( 'U', K, A( 0+K ), K, INFO )
382 IF( INFO.GT.0 )
383 $ RETURN
384 CALL ZTRSM( 'L', 'U', 'C', 'N', K, K, CONE, A( K ), N1,
385 $ A( K*( K+1 ) ), K )
386 CALL ZHERK( 'L', 'C', K, K, -ONE, A( K*( K+1 ) ), K, ONE,
387 $ A( 0 ), K )
388 CALL ZPOTRF( 'L', K, A( 0 ), K, INFO )
389 IF( INFO.GT.0 )
390 $ INFO = INFO + K
391 *
392 ELSE
393 *
394 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
395 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
396 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
397 *
398 CALL ZPOTRF( 'U', K, A( K*( K+1 ) ), K, INFO )
399 IF( INFO.GT.0 )
400 $ RETURN
401 CALL ZTRSM( 'R', 'U', 'N', 'N', K, K, CONE,
402 $ A( K*( K+1 ) ), K, A( 0 ), K )
403 CALL ZHERK( 'L', 'N', K, K, -ONE, A( 0 ), K, ONE,
404 $ A( K*K ), K )
405 CALL ZPOTRF( 'L', K, A( K*K ), K, INFO )
406 IF( INFO.GT.0 )
407 $ INFO = INFO + K
408 *
409 END IF
410 *
411 END IF
412 *
413 END IF
414 *
415 RETURN
416 *
417 * End of ZPFTRF
418 *
419 END