LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cpftrf ( character  TRANSR,
character  UPLO,
integer  N,
complex, dimension( 0: * )  A,
integer  INFO 
)

CPFTRF

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

Purpose: 
 CPFTRF 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 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
November 2011

Definition at line 213 of file cpftrf.f.

213 *
214 * -- LAPACK computational routine (version 3.4.0) --
215 * -- LAPACK is a software package provided by Univ. of Tennessee, --
216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * November 2011
218 *
219 * .. Scalar Arguments ..
220  CHARACTER transr, uplo
221  INTEGER n, info
222 * ..
223 * .. Array Arguments ..
224  COMPLEX a( 0: * )
225 *
226 * =====================================================================
227 *
228 * .. Parameters ..
229  REAL one
230  COMPLEX cone
231  parameter ( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+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, cherk, cpotrf, ctrsm
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( 'CPFTRF', -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 cpotrf( 'L', n1, a( 0 ), n, info )
308  IF( info.GT.0 )
309  $ RETURN
310  CALL ctrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ), n,
311  $ a( n1 ), n )
312  CALL cherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
313  $ a( n ), n )
314  CALL cpotrf( '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 cpotrf( 'L', n1, a( n2 ), n, info )
325  IF( info.GT.0 )
326  $ RETURN
327  CALL ctrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ), n,
328  $ a( 0 ), n )
329  CALL cherk( 'U', 'C', n2, n1, -one, a( 0 ), n, one,
330  $ a( n1 ), n )
331  CALL cpotrf( '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 cpotrf( 'U', n1, a( 0 ), n1, info )
348  IF( info.GT.0 )
349  $ RETURN
350  CALL ctrsm( 'L', 'U', 'C', 'N', n1, n2, cone, a( 0 ), n1,
351  $ a( n1*n1 ), n1 )
352  CALL cherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1, one,
353  $ a( 1 ), n1 )
354  CALL cpotrf( '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 cpotrf( 'U', n1, a( n2*n2 ), n2, info )
365  IF( info.GT.0 )
366  $ RETURN
367  CALL ctrsm( 'R', 'U', 'N', 'N', n2, n1, cone, a( n2*n2 ),
368  $ n2, a( 0 ), n2 )
369  CALL cherk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
370  $ a( n1*n2 ), n2 )
371  CALL cpotrf( '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 cpotrf( 'L', k, a( 1 ), n+1, info )
394  IF( info.GT.0 )
395  $ RETURN
396  CALL ctrsm( 'R', 'L', 'C', 'N', k, k, cone, a( 1 ), n+1,
397  $ a( k+1 ), n+1 )
398  CALL cherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
399  $ a( 0 ), n+1 )
400  CALL cpotrf( '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 cpotrf( 'L', k, a( k+1 ), n+1, info )
411  IF( info.GT.0 )
412  $ RETURN
413  CALL ctrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
414  $ n+1, a( 0 ), n+1 )
415  CALL cherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
416  $ a( k ), n+1 )
417  CALL cpotrf( '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 cpotrf( 'U', k, a( 0+k ), k, info )
434  IF( info.GT.0 )
435  $ RETURN
436  CALL ctrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ), n1,
437  $ a( k*( k+1 ) ), k )
438  CALL cherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k, one,
439  $ a( 0 ), k )
440  CALL cpotrf( '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 cpotrf( 'U', k, a( k*( k+1 ) ), k, info )
451  IF( info.GT.0 )
452  $ RETURN
453  CALL ctrsm( 'R', 'U', 'N', 'N', k, k, cone,
454  $ a( k*( k+1 ) ), k, a( 0 ), k )
455  CALL cherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
456  $ a( k*k ), k )
457  CALL cpotrf( '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 CPFTRF
470 *
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ctrsm
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
cpotrf
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:109
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: