LAPACK 3.12.0
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 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)
Definition cblat2.f:3285
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
Here is the call graph for this function:
Here is the caller graph for this function: