LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cpftrf()

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.

Definition at line 210 of file cpftrf.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 A( 0: * )
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  REAL ONE
227  COMPLEX CONE
228  parameter( one = 1.0e+0, cone = ( 1.0e+0, 0.0e+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, cherk, cpotrf, ctrsm
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( 'CPFTRF', -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 cpotrf( 'L', n1, a( 0 ), n, info )
305  IF( info.GT.0 )
306  $ RETURN
307  CALL ctrsm( 'R', 'L', 'C', 'N', n2, n1, cone, a( 0 ), n,
308  $ a( n1 ), n )
309  CALL cherk( 'U', 'N', n2, n1, -one, a( n1 ), n, one,
310  $ a( n ), n )
311  CALL cpotrf( '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 cpotrf( 'L', n1, a( n2 ), n, info )
322  IF( info.GT.0 )
323  $ RETURN
324  CALL ctrsm( 'L', 'L', 'N', 'N', n1, n2, cone, a( n2 ), n,
325  $ a( 0 ), n )
326  CALL cherk( 'U', 'C', n2, n1, -one, a( 0 ), n, one,
327  $ a( n1 ), n )
328  CALL cpotrf( '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 cpotrf( 'U', n1, a( 0 ), n1, info )
345  IF( info.GT.0 )
346  $ RETURN
347  CALL ctrsm( 'L', 'U', 'C', 'N', n1, n2, cone, a( 0 ), n1,
348  $ a( n1*n1 ), n1 )
349  CALL cherk( 'L', 'C', n2, n1, -one, a( n1*n1 ), n1, one,
350  $ a( 1 ), n1 )
351  CALL cpotrf( '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 cpotrf( 'U', n1, a( n2*n2 ), n2, info )
362  IF( info.GT.0 )
363  $ RETURN
364  CALL ctrsm( 'R', 'U', 'N', 'N', n2, n1, cone, a( n2*n2 ),
365  $ n2, a( 0 ), n2 )
366  CALL cherk( 'L', 'N', n2, n1, -one, a( 0 ), n2, one,
367  $ a( n1*n2 ), n2 )
368  CALL cpotrf( '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 cpotrf( 'L', k, a( 1 ), n+1, info )
391  IF( info.GT.0 )
392  $ RETURN
393  CALL ctrsm( 'R', 'L', 'C', 'N', k, k, cone, a( 1 ), n+1,
394  $ a( k+1 ), n+1 )
395  CALL cherk( 'U', 'N', k, k, -one, a( k+1 ), n+1, one,
396  $ a( 0 ), n+1 )
397  CALL cpotrf( '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 cpotrf( 'L', k, a( k+1 ), n+1, info )
408  IF( info.GT.0 )
409  $ RETURN
410  CALL ctrsm( 'L', 'L', 'N', 'N', k, k, cone, a( k+1 ),
411  $ n+1, a( 0 ), n+1 )
412  CALL cherk( 'U', 'C', k, k, -one, a( 0 ), n+1, one,
413  $ a( k ), n+1 )
414  CALL cpotrf( '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 cpotrf( 'U', k, a( 0+k ), k, info )
431  IF( info.GT.0 )
432  $ RETURN
433  CALL ctrsm( 'L', 'U', 'C', 'N', k, k, cone, a( k ), n1,
434  $ a( k*( k+1 ) ), k )
435  CALL cherk( 'L', 'C', k, k, -one, a( k*( k+1 ) ), k, one,
436  $ a( 0 ), k )
437  CALL cpotrf( '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 cpotrf( 'U', k, a( k*( k+1 ) ), k, info )
448  IF( info.GT.0 )
449  $ RETURN
450  CALL ctrsm( 'R', 'U', 'N', 'N', k, k, cone,
451  $ a( k*( k+1 ) ), k, a( 0 ), k )
452  CALL cherk( 'L', 'N', k, k, -one, a( 0 ), k, one,
453  $ a( k*k ), k )
454  CALL cpotrf( '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 CPFTRF
467 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:173
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
subroutine cpotrf(UPLO, N, A, LDA, INFO)
CPOTRF
Definition: cpotrf.f:107
Here is the call graph for this function:
Here is the caller graph for this function: