hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_lapack_zpstrf (f07kr)

Purpose

nag_lapack_zpstrf (f07kr) computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.

Syntax

[a, piv, rank, info] = f07kr(uplo, a, 'n', n, 'tol', tol)
[a, piv, rank, info] = nag_lapack_zpstrf(uplo, a, 'n', n, 'tol', tol)

Description

nag_lapack_zpstrf (f07kr) forms the Cholesky factorization of a complex Hermitian positive semidefinite matrix AA either as PTAP = UHUPTAP=UHU if uplo = 'U'uplo='U' or PTAP = LLHPTAP=LLH if uplo = 'L'uplo='L', where PP is a permutation matrix, UU is an upper triangular matrix and LL is lower triangular.
This algorithm does not attempt to check that AA is positive semidefinite.

References

Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia
Lucas C (2004) LAPACK-style codes for Level 2 and 3 pivoted Cholesky factorizations LAPACK Working Note 161. Technical Report CS-04-522 Department of Computer Science, University of Tennessee, 107 Ayres Hall, Knoxville, TN 37996-1301, USA

Parameters

Compulsory Input Parameters

1:     uplo – string (length ≥ 1)
Specifies whether the upper or lower triangular part of AA is stored and how AA is to be factorized.
uplo = 'U'uplo='U'
The upper triangular part of AA is stored and AA is factorized as UHUUHU, where UU is upper triangular.
uplo = 'L'uplo='L'
The lower triangular part of AA is stored and AA is factorized as LLHLLH, where LL is lower triangular.
Constraint: uplo = 'U'uplo='U' or 'L''L'.
2:     a(lda, : :) – complex array
The first dimension of the array a must be at least max (1,n)max(1,n)
The second dimension of the array must be at least max (1,n)max(1,n)
The nn by nn Hermitian positive semidefinite matrix AA.
  • If uplo = 'U'uplo='U', the upper triangular part of aa must be stored and the elements of the array below the diagonal are not referenced.
  • If uplo = 'L'uplo='L', the lower triangular part of aa must be stored and the elements of the array above the diagonal are not referenced.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array a The second dimension of the array a.
nn, the order of the matrix AA.
Constraint: n0n0.
2:     tol – double scalar
User defined tolerance. If tol < 0tol<0, then n × maxk = 1,n |Akk| × machine precisionn×maxk=1,n|Akk|×machine precision will be used. The algorithm terminates at the rrth step if the (r + 1)(r+1)th step pivot < tol<tol.
Default: -1-1

Input Parameters Omitted from the MATLAB Interface

lda work

Output Parameters

1:     a(lda, : :) – complex array
The first dimension of the array a will be max (1,n)max(1,n)
The second dimension of the array will be max (1,n)max(1,n)
ldamax (1,n)ldamax(1,n).
If uplo = 'U'uplo='U', the first rank rows of the upper triangle of AA are overwritten with the nonzero elements of the Cholesky factor UU, and the remaining rows of the triangle are destroyed.
If uplo = 'L'uplo='L', the first rank columns of the lower triangle of AA are overwritten with the nonzero elements of the Cholesky factor LL, and the remaining columns of the triangle are destroyed.
2:     piv(n) – int64int32nag_int array
piv is such that the nonzero entries of PP are P(piv(k),k) = 1P(piv(k),k)=1, for k = 1,2,,nk=1,2,,n.
3:     rank – int64int32nag_int scalar
The computed rank of AA given by the number of steps the algorithm completed.
4:     info – int64int32nag_int scalar
info = 0info=0 unless the function detects an error (see Section [Error Indicators and Warnings]).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  INFO > 0INFO>0
The leading minor of order __ is not positive definite and the factorization could not be completed. Hence AA itself is not positive definite. This may indicate an error in forming the matrix AA. To factorize a Hermitian matrix which is not positive definite, call nag_lapack_zhetrf (f07mr) instead.
W INFO = 1INFO=1
The matrix AA is either rank deficient with computed rank as returned in rank, or is indefinite, see Section [Further Comments].

Accuracy

If uplo = 'L'uplo='L' and rank = rrank=r, the computed Cholesky factor LL and permutation matrix PP satisfy the following upper bound
(APLLHPT2)/(A2) 2r c(r) ε (W2 + 1)2 + O(ε2) ,
A - PLLHPT 2 A2 2r c(r) ε ( W 2 + 1 ) 2 + O(ε2) ,
where
W = L111 L12 ,   L =
(L110)
L12 0
,   L11 r × r ,
W = L 11 -1 L12 ,   L = L11 0 L12 0 ,   L11 r×r ,
c(r)c(r) is a modest linear function of rr, εε is machine precision, and
W2 sqrt( (1/3) (nr) (4r1) ) .
W2 13 (n-r) (4r-1) .
So there is no guarantee of stability of the algorithm for large nn and rr, although W2W2 is generally small in practice.

Further Comments

The total number of real floating point operations is approximately 4nr28 / 3r34nr2-8/3r3, where rr is the computed rank of AA.
This algorithm does not attempt to check that AA is positive semidefinite, and in particular the rank detection criterion in the algorithm is based on AA being positive semidefinite. If there is doubt over semidefiniteness then you should use the indefinite factorization nag_lapack_zhetrf (f07mr). See Lucas (2004) for further information.
The real analogue of this function is nag_lapack_dpstrf (f07kd).

Example

function nag_lapack_zpstrf_example
a = [12.40,       2.39,       5.50+0.05i, 4.47,       11.89;
      2.39,       1.63,       1.04+0.10i, 1.14,        1.81;
      5.50+0.05i, 1.04+0.10i, 2.45,       1.98-0.03i,  5.28-0.02i;
      4.47,       1.14,       1.98-0.03i, 1.71,        4.14;
     11.89,       1.81,       5.28-0.02i, 4.14,       11.63];
uplo = 'l';
% Factorize a
[aOut, piv, rank, info] = nag_lapack_zpstrf(uplo, a);

fprintf('\nComputed rank: %d\n\n', rank);
% Zero out columns rank+1 onwards
aOut(:, double(rank)+1:5) = 0;
[ifail] = ...
    nag_file_print_matrix_complex_gen_comp(uplo, 'n', aOut, 'b', 'f5.2', 'Factor', 'i', 'i', int64(80), int64(0));
fprintf('\n piv:\n  ');
fprintf('             %d', piv);
fprintf('\n');
 
Warning: nag_lapack_zpstrf (f07kr) returned a warning indicator (1) 

Computed rank: 3

 Factor
                1             2             3             4             5
 1  ( 3.52, 0.00)
 2  ( 0.68, 0.00) ( 1.08, 0.00)
 3  ( 1.27, 0.00) ( 0.26, 0.00) ( 0.18, 0.00)
 4  ( 1.56, 0.01) (-0.02, 0.08) ( 0.01,-0.05) ( 0.00, 0.00)
 5  ( 3.38, 0.00) (-0.45, 0.00) (-0.17, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)

 piv:
               1             2             4             3             5

function f07kr_example
a = [12.40,       2.39,       5.50+0.05i, 4.47,       11.89;
      2.39,       1.63,       1.04+0.10i, 1.14,        1.81;
      5.50+0.05i, 1.04+0.10i, 2.45,       1.98-0.03i,  5.28-0.02i;
      4.47,       1.14,       1.98-0.03i, 1.71,        4.14;
     11.89,       1.81,       5.28-0.02i, 4.14,       11.63];
uplo = 'l';
% Factorize a
[aOut, piv, rank, info] = f07kr(uplo, a);

fprintf('\nComputed rank: %d\n\n', rank);
% Zero out columns rank+1 onwards
aOut(:, double(rank)+1:5) = 0;
[ifail] = x04db(uplo, 'n', aOut, 'b', 'f5.2', 'Factor', 'i', 'i', int64(80), int64(0));
fprintf('\n piv:\n  ');
fprintf('             %d', piv);
fprintf('\n');
 
Warning: nag_lapack_zpstrf (f07kr) returned a warning indicator (1) 

Computed rank: 3

 Factor
                1             2             3             4             5
 1  ( 3.52, 0.00)
 2  ( 0.68, 0.00) ( 1.08, 0.00)
 3  ( 1.27, 0.00) ( 0.26, 0.00) ( 0.18, 0.00)
 4  ( 1.56, 0.01) (-0.02, 0.08) ( 0.01,-0.05) ( 0.00, 0.00)
 5  ( 3.38, 0.00) (-0.45, 0.00) (-0.17, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)

 piv:
               1             2             4             3             5


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013