LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zpftrf()

 subroutine zpftrf ( character TRANSR, character UPLO, integer N, complex*16, dimension( 0: * ) A, integer INFO )

ZPFTRF

Purpose:
``` ZPFTRF computes the Cholesky factorization of a complex Hermitian
positive definite matrix A.

The factorization has the form
A = U**H * U,  if UPLO = 'U', or
A = L  * L**H,  if UPLO = 'L',
where U is an upper triangular matrix and L is lower triangular.

This is the block version of the algorithm, calling Level 3 BLAS.```
Parameters
 [in] TRANSR ``` TRANSR is CHARACTER*1 = 'N': The Normal TRANSR of RFP A is stored; = 'C': The Conjugate-transpose TRANSR of RFP A is stored.``` [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of RFP A is stored; = 'L': Lower triangle of RFP A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is COMPLEX*16 array, dimension ( N*(N+1)/2 ); On entry, the Hermitian matrix A in RFP format. RFP format is described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is the Conjugate-transpose of RFP A as defined when TRANSR = 'N'. The contents of RFP A are defined by UPLO as follows: If UPLO = 'U' the RFP A contains the nt elements of upper packed A. If UPLO = 'L' the RFP A contains the elements of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N is odd. See the Note below for more details. On exit, if INFO = 0, the factor U or L from the Cholesky factorization RFP A = U**H*U or RFP A = L*L**H.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i is not positive definite, and the factorization could not be completed. Further Notes on RFP Format: ============================ We first consider Standard Packed Format when N is even. We give an example where N = 6. AP is Upper AP is Lower 00 01 02 03 04 05 00 11 12 13 14 15 10 11 22 23 24 25 20 21 22 33 34 35 30 31 32 33 44 45 40 41 42 43 44 55 50 51 52 53 54 55 Let TRANSR = 'N'. RFP holds AP as follows: For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last three columns of AP upper. The lower triangle A(4:6,0:2) consists of conjugate-transpose of the first three columns of AP upper. For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first three columns of AP lower. The upper triangle A(0:2,0:2) consists of conjugate-transpose of the last three columns of AP lower. To denote conjugate we place -- above the element. This covers the case N even and TRANSR = 'N'. RFP A RFP A -- -- -- 03 04 05 33 43 53 -- -- 13 14 15 00 44 54 -- 23 24 25 10 11 55 33 34 35 20 21 22 -- 00 44 45 30 31 32 -- -- 01 11 55 40 41 42 -- -- -- 02 12 22 50 51 52 Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- transpose of RFP A above. One therefore gets: RFP A RFP A -- -- -- -- -- -- -- -- -- -- 03 13 23 33 00 01 02 33 00 10 20 30 40 50 -- -- -- -- -- -- -- -- -- -- 04 14 24 34 44 11 12 43 44 11 21 31 41 51 -- -- -- -- -- -- -- -- -- -- 05 15 25 35 45 55 22 53 54 55 22 32 42 52 We next consider Standard Packed Format when N is odd. We give an example where N = 5. AP is Upper AP is Lower 00 01 02 03 04 00 11 12 13 14 10 11 22 23 24 20 21 22 33 34 30 31 32 33 44 40 41 42 43 44 Let TRANSR = 'N'. RFP holds AP as follows: For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last three columns of AP upper. The lower triangle A(3:4,0:1) consists of conjugate-transpose of the first two columns of AP upper. For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first three columns of AP lower. The upper triangle A(0:1,1:2) consists of conjugate-transpose of the last two columns of AP lower. To denote conjugate we place -- above the element. This covers the case N odd and TRANSR = 'N'. RFP A RFP A -- -- 02 03 04 00 33 43 -- 12 13 14 10 11 44 22 23 24 20 21 22 -- 00 33 34 30 31 32 -- -- 01 11 44 40 41 42 Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- transpose of RFP A above. One therefore gets: RFP A RFP A -- -- -- -- -- -- -- -- -- 02 12 22 00 01 00 10 20 30 40 50 -- -- -- -- -- -- -- -- -- 03 13 23 33 11 33 11 21 31 41 51 -- -- -- -- -- -- -- -- -- 04 14 24 34 44 43 44 22 32 42 52```

Definition at line 210 of file zpftrf.f.

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*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:173
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
Here is the call graph for this function:
Here is the caller graph for this function: