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