LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 213 of file zpftrf.f.

213 *
214 * -- LAPACK computational routine (version 3.6.1) --
215 * -- LAPACK is a software package provided by Univ. of Tennessee, --
216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * June 2016
218 *
219 * .. Scalar Arguments ..
220  CHARACTER transr, uplo
221  INTEGER n, info
222 * ..
223 * .. Array Arguments ..
224  COMPLEX*16 a( 0: * )
225 *
226 * =====================================================================
227 *
228 * .. Parameters ..
229  DOUBLE PRECISION one
230  COMPLEX*16 cone
231  parameter ( one = 1.0d+0, cone = ( 1.0d+0, 0.0d+0 ) )
232 * ..
233 * .. Local Scalars ..
234  LOGICAL lower, nisodd, normaltransr
235  INTEGER n1, n2, k
236 * ..
237 * .. External Functions ..
238  LOGICAL lsame
239  EXTERNAL lsame
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL xerbla, zherk, zpotrf, ztrsm
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC mod
246 * ..
247 * .. Executable Statements ..
248 *
249 * Test the input parameters.
250 *
251  info = 0
252  normaltransr = lsame( transr, 'N' )
253  lower = lsame( uplo, 'L' )
254  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'C' ) ) THEN
255  info = -1
256  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
257  info = -2
258  ELSE IF( n.LT.0 ) THEN
259  info = -3
260  END IF
261  IF( info.NE.0 ) THEN
262  CALL xerbla( 'ZPFTRF', -info )
263  RETURN
264  END IF
265 *
266 * Quick return if possible
267 *
268  IF( n.EQ.0 )
269  $ RETURN
270 *
271 * If N is odd, set NISODD = .TRUE.
272 * If N is even, set K = N/2 and NISODD = .FALSE.
273 *
274  IF( mod( n, 2 ).EQ.0 ) THEN
275  k = n / 2
276  nisodd = .false.
277  ELSE
278  nisodd = .true.
279  END IF
280 *
281 * Set N1 and N2 depending on LOWER
282 *
283  IF( lower ) THEN
284  n2 = n / 2
285  n1 = n - n2
286  ELSE
287  n1 = n / 2
288  n2 = n - n1
289  END IF
290 *
291 * start execution: there are eight cases
292 *
293  IF( nisodd ) THEN
294 *
295 * N is odd
296 *
297  IF( normaltransr ) THEN
298 *
299 * N is odd and TRANSR = 'N'
300 *
301  IF( lower ) THEN
302 *
303 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
304 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
305 * T1 -> a(0), T2 -> a(n), S -> a(n1)
306 *
307  CALL zpotrf( 'L', n1, a( 0 ), n, info )
308  IF( info.GT.0 )
309  $ RETURN
310  CALL ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ), n,
311  $ a( n1 ), n )
312  CALL zherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
313  $ a( n ), n )
314  CALL zpotrf( 'U', n2, a( n ), n, info )
315  IF( info.GT.0 )
316  $ info = info + n1
317 *
318  ELSE
319 *
320 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
321 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
322 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
323 *
324  CALL zpotrf( 'L', n1, a( n2 ), n, info )
325  IF( info.GT.0 )
326  $ RETURN
327  CALL ztrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ), n,
328  $ a( 0 ), n )
329  CALL zherk( 'U', 'C', n2, n1, -one, a( 0 ), n, one,
330  $ a( n1 ), n )
331  CALL zpotrf( 'U', n2, a( n1 ), n, info )
332  IF( info.GT.0 )
333  $ info = info + n1
334 *
335  END IF
336 *
337  ELSE
338 *
339 * N is odd and TRANSR = 'C'
340 *
341  IF( lower ) THEN
342 *
343 * SRPA for LOWER, TRANSPOSE and N is odd
344 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
345 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
346 *
347  CALL zpotrf( 'U', n1, a( 0 ), n1, info )
348  IF( info.GT.0 )
349  $ RETURN
350  CALL ztrsm( 'L', 'U', 'C', 'N', n1, n2, cone, a( 0 ), n1,
351  $ a( n1*n1 ), n1 )
352  CALL zherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1, one,
353  $ a( 1 ), n1 )
354  CALL zpotrf( 'L', n2, a( 1 ), n1, info )
355  IF( info.GT.0 )
356  $ info = info + n1
357 *
358  ELSE
359 *
360 * SRPA for UPPER, TRANSPOSE and N is odd
361 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
362 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
363 *
364  CALL zpotrf( 'U', n1, a( n2*n2 ), n2, info )
365  IF( info.GT.0 )
366  $ RETURN
367  CALL ztrsm( 'R', 'U', 'N', 'N', n2, n1, cone, 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 ), n+1,
397  $ a( k+1 ), n+1 )
398  CALL zherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
399  $ a( 0 ), n+1 )
400  CALL zpotrf( 'U', k, a( 0 ), n+1, info )
401  IF( info.GT.0 )
402  $ info = info + k
403 *
404  ELSE
405 *
406 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
407 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
408 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
409 *
410  CALL zpotrf( 'L', k, a( k+1 ), n+1, info )
411  IF( info.GT.0 )
412  $ RETURN
413  CALL ztrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
414  $ n+1, a( 0 ), n+1 )
415  CALL zherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
416  $ a( k ), n+1 )
417  CALL zpotrf( 'U', k, a( k ), n+1, info )
418  IF( info.GT.0 )
419  $ info = info + k
420 *
421  END IF
422 *
423  ELSE
424 *
425 * N is even and TRANSR = 'C'
426 *
427  IF( lower ) THEN
428 *
429 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
430 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
431 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
432 *
433  CALL zpotrf( 'U', k, a( 0+k ), k, info )
434  IF( info.GT.0 )
435  $ RETURN
436  CALL ztrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ), n1,
437  $ a( k*( k+1 ) ), k )
438  CALL zherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k, one,
439  $ a( 0 ), k )
440  CALL zpotrf( 'L', k, a( 0 ), k, info )
441  IF( info.GT.0 )
442  $ info = info + k
443 *
444  ELSE
445 *
446 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
447 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
448 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
449 *
450  CALL zpotrf( 'U', k, a( k*( k+1 ) ), k, info )
451  IF( info.GT.0 )
452  $ RETURN
453  CALL ztrsm( 'R', 'U', 'N', 'N', k, k, cone,
454  $ a( k*( k+1 ) ), k, a( 0 ), k )
455  CALL zherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
456  $ a( k*k ), k )
457  CALL zpotrf( 'L', k, a( k*k ), k, info )
458  IF( info.GT.0 )
459  $ info = info + k
460 *
461  END IF
462 *
463  END IF
464 *
465  END IF
466 *
467  RETURN
468 *
469 * End of ZPFTRF
470 *
zpotrf
subroutine zpotrf(UPLO, N, A, LDA, INFO)
ZPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
Definition: zpotrf.f:102
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175
ztrsm
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: