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

◆ 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 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 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)
Definition cblat2.f:3285
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine cpotrf(uplo, n, a, lda, info)
CPOTRF
Definition cpotrf.f:107
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
Here is the call graph for this function:
Here is the caller graph for this function: