LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zpftrf()

subroutine zpftrf ( character transr,
character uplo,
integer n,
complex*16, dimension( 0: * ) a,
integer info )

ZPFTRF

Download ZPFTRF + dependencies [TGZ] [ZIP] [TXT]

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 principal minor of order i
!>                is not positive, 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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 208 of file zpftrf.f.

209*
210* -- LAPACK computational routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 CHARACTER TRANSR, UPLO
216 INTEGER N, INFO
217* ..
218* .. Array Arguments ..
219 COMPLEX*16 A( 0: * )
220*
221* =====================================================================
222*
223* .. Parameters ..
224 DOUBLE PRECISION ONE
225 COMPLEX*16 CONE
226 parameter( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ) )
227* ..
228* .. Local Scalars ..
229 LOGICAL LOWER, NISODD, NORMALTRANSR
230 INTEGER N1, N2, K
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 EXTERNAL lsame
235* ..
236* .. External Subroutines ..
237 EXTERNAL xerbla, zherk, zpotrf, ztrsm
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC mod
241* ..
242* .. Executable Statements ..
243*
244* Test the input parameters.
245*
246 info = 0
247 normaltransr = lsame( transr, 'N' )
248 lower = lsame( uplo, 'L' )
249 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
250 info = -1
251 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
252 info = -2
253 ELSE IF( n.LT.0 ) THEN
254 info = -3
255 END IF
256 IF( info.NE.0 ) THEN
257 CALL xerbla( 'ZPFTRF', -info )
258 RETURN
259 END IF
260*
261* Quick return if possible
262*
263 IF( n.EQ.0 )
264 $ RETURN
265*
266* If N is odd, set NISODD = .TRUE.
267* If N is even, set K = N/2 and NISODD = .FALSE.
268*
269 IF( mod( n, 2 ).EQ.0 ) THEN
270 k = n / 2
271 nisodd = .false.
272 ELSE
273 nisodd = .true.
274 END IF
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* start execution: there are eight cases
287*
288 IF( nisodd ) THEN
289*
290* N is odd
291*
292 IF( normaltransr ) THEN
293*
294* N is odd and TRANSR = 'N'
295*
296 IF( lower ) THEN
297*
298* SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
299* T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
300* T1 -> a(0), T2 -> a(n), S -> a(n1)
301*
302 CALL zpotrf( 'L', n1, a( 0 ), n, info )
303 IF( info.GT.0 )
304 $ RETURN
305 CALL ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ),
306 $ n,
307 $ a( n1 ), n )
308 CALL zherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
309 $ a( n ), n )
310 CALL zpotrf( 'U', n2, a( n ), n, info )
311 IF( info.GT.0 )
312 $ info = info + n1
313*
314 ELSE
315*
316* SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
317* T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
318* T1 -> a(n2), T2 -> a(n1), S -> a(0)
319*
320 CALL zpotrf( 'L', n1, a( n2 ), n, info )
321 IF( info.GT.0 )
322 $ RETURN
323 CALL ztrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ),
324 $ 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 ),
348 $ n1,
349 $ a( n1*n1 ), n1 )
350 CALL zherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1,
351 $ one,
352 $ a( 1 ), n1 )
353 CALL zpotrf( 'L', n2, a( 1 ), n1, info )
354 IF( info.GT.0 )
355 $ info = info + n1
356*
357 ELSE
358*
359* SRPA for UPPER, TRANSPOSE and N is odd
360* T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
361* T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
362*
363 CALL zpotrf( 'U', n1, a( n2*n2 ), n2, info )
364 IF( info.GT.0 )
365 $ RETURN
366 CALL ztrsm( 'R', 'U', 'N', 'N', n2, n1, cone,
367 $ a( n2*n2 ),
368 $ n2, a( 0 ), n2 )
369 CALL zherk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
370 $ a( n1*n2 ), n2 )
371 CALL zpotrf( 'L', n2, a( n1*n2 ), n2, info )
372 IF( info.GT.0 )
373 $ info = info + n1
374*
375 END IF
376*
377 END IF
378*
379 ELSE
380*
381* N is even
382*
383 IF( normaltransr ) THEN
384*
385* N is even and TRANSR = 'N'
386*
387 IF( lower ) THEN
388*
389* SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
390* T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
391* T1 -> a(1), T2 -> a(0), S -> a(k+1)
392*
393 CALL zpotrf( 'L', k, a( 1 ), n+1, info )
394 IF( info.GT.0 )
395 $ RETURN
396 CALL ztrsm( 'R', 'L', 'C', 'N', k, k, cone, a( 1 ),
397 $ n+1,
398 $ a( k+1 ), n+1 )
399 CALL zherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
400 $ a( 0 ), n+1 )
401 CALL zpotrf( 'U', k, a( 0 ), n+1, info )
402 IF( info.GT.0 )
403 $ info = info + k
404*
405 ELSE
406*
407* SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
408* T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
409* T1 -> a(k+1), T2 -> a(k), S -> a(0)
410*
411 CALL zpotrf( 'L', k, a( k+1 ), n+1, info )
412 IF( info.GT.0 )
413 $ RETURN
414 CALL ztrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
415 $ n+1, a( 0 ), n+1 )
416 CALL zherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
417 $ a( k ), n+1 )
418 CALL zpotrf( 'U', k, a( k ), n+1, info )
419 IF( info.GT.0 )
420 $ info = info + k
421*
422 END IF
423*
424 ELSE
425*
426* N is even and TRANSR = 'C'
427*
428 IF( lower ) THEN
429*
430* SRPA for LOWER, TRANSPOSE and N is even (see paper)
431* T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
432* T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
433*
434 CALL zpotrf( 'U', k, a( 0+k ), k, info )
435 IF( info.GT.0 )
436 $ RETURN
437 CALL ztrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ),
438 $ n1,
439 $ a( k*( k+1 ) ), k )
440 CALL zherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k,
441 $ one,
442 $ a( 0 ), k )
443 CALL zpotrf( 'L', k, a( 0 ), k, info )
444 IF( info.GT.0 )
445 $ info = info + k
446*
447 ELSE
448*
449* SRPA for UPPER, TRANSPOSE and N is even (see paper)
450* T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
451* T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
452*
453 CALL zpotrf( 'U', k, a( k*( k+1 ) ), k, info )
454 IF( info.GT.0 )
455 $ RETURN
456 CALL ztrsm( 'R', 'U', 'N', 'N', k, k, cone,
457 $ a( k*( k+1 ) ), k, a( 0 ), k )
458 CALL zherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
459 $ a( k*k ), k )
460 CALL zpotrf( 'L', k, a( k*k ), k, info )
461 IF( info.GT.0 )
462 $ info = info + k
463*
464 END IF
465*
466 END IF
467*
468 END IF
469*
470 RETURN
471*
472* End of ZPFTRF
473*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zpotrf
subroutine zpotrf(uplo, n, a, lda, info)
ZPOTRF
Definition zpotrf.f:105
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180
Here is the call graph for this function:
Here is the caller graph for this function: