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

## ◆ cpftrf()

 subroutine cpftrf ( character transr, character uplo, integer n, complex, dimension( 0: * ) a, integer info )

CPFTRF

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```

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: