Lucas C (2004) LAPACK-style codes for Level 2 and 3 pivoted Cholesky factorizations LAPACK Working Note No. 161. Technical Report CS-04-522 Department of Computer Science, University of Tennessee, 107 Ayres Hall, Knoxville, TN 37996-1301, USA http://www.netlib.org/lapack/lawnspdf/lawn161.pdf

Parameters

Compulsory Input Parameters

1: $uplo$ – string (length ≥ 1)

Specifies whether the upper or lower triangular part of

A

is stored and how

A

is to be factorized.

$uplo ='U'$: The upper triangular part of $A$ is stored and $A$ is factorized as $U^{T} U$ , where $U$ is upper triangular.
$uplo ='L'$: The lower triangular part of $A$ is stored and $A$ is factorized as $L L^{T}$ , where $L$ is lower triangular.

Constraint:

uplo ='U'

'L'

2: $a (lda :)$ – double array

The first dimension of the array a must be at least

\max (1, n)

The second dimension of the array a must be at least

\max (1, n)

The

n

n

symmetric positive semidefinite matrix

A

If $uplo ='U'$ , the upper triangular part of $a$ must be stored and the elements of the array below the diagonal are not referenced.
If $uplo ='L'$ , the lower triangular part of $a$ 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 and the second dimension of the array a.
$n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $tol$ – double scalar: Default: $- 1$
User defined tolerance. If $tol < 0$ , then $n \times max_{k = 1, n} |A_{k k}| \times machine precision$ will be used. The algorithm terminates at the $r$ th step if the $(r + 1)$ th step pivot $< tol$ .

Output Parameters

1: $a (lda :)$ – double array: The first dimension of the array a will be $\max (1, n)$ .
The second dimension of the array a will be $\max (1, n)$ .

If $uplo ='U'$ , the first rank rows of the upper triangle of $A$ are overwritten with the nonzero elements of the Cholesky factor $U$ , and the remaining rows of the triangle are destroyed.
If $uplo ='L'$ , the first rank columns of the lower triangle of $A$ are overwritten with the nonzero elements of the Cholesky factor $L$ , and the remaining columns of the triangle are destroyed.
2: $piv (n)$ – int64int32nag_int array: piv is such that the nonzero entries of $P$ are $P (piv (k), k) = 1$ , for $k = 1, 2, \dots, n$ .
3: $rank$ – int64int32nag_int scalar: The computed rank of $A$ given by the number of steps the algorithm completed.
4: $info$ – int64int32nag_int scalar: $info = 0$ unless the function detects an error (see 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 < 0$: If $info = - i$ , argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.

W $info = 1$: The matrix $A$ is not positive definite. It is either positive semidefinite with computed rank as returned in rank and less than $n$ , or it may be indefinite, see Further Comments.

Accuracy

uplo ='L'

and

rank = r

, the computed Cholesky factor

L

and permutation matrix

P

satisfy the following upper bound

\frac{{‖A - P L L^{T} P^{T}‖}_{2}}{{‖A‖}_{2}} \leq 2 r c (r) ε {({‖W‖}_{2} + 1)}^{2} + O (ε^{2}),

where

W = L_{11}^{- 1} L_{12}, L = (\begin{matrix} L_{11} & 0 \\ L_{12} & 0 \end{matrix}), L_{11} \in ℝ^{r \times r},

c (r)

is a modest linear function of

r

ε

is machine precision, and

{‖W‖}_{2} \leq \sqrt{\frac{1}{3} (n - r) (4^{r} - 1)} .

So there is no guarantee of stability of the algorithm for large

n

and

r

, although

{‖W‖}_{2}

is generally small in practice.

Further Comments

The total number of floating-point operations is approximately

n r^{2} - 2 / 3 r^{3}

, where

r

is the computed rank of

A

This algorithm does not attempt to check that

A

is positive semidefinite, and in particular the rank detection criterion in the algorithm is based on

A

being positive semidefinite. If there is doubt over semidefiniteness then you should use the indefinite factorization nag_lapack_dsytrf (f07md). See Lucas (2004) for further information.

The complex analogue of this function is nag_lapack_zpstrf (f07kr).

Example

This example computes the Cholesky factorization of the matrix

A

, where

A = (\begin{matrix} 2.51 & 4.04 & 3.34 & 1.34 & 1.29 \\ 4.04 & 8.22 & 7.38 & 2.68 & 2.44 \\ 3.34 & 7.38 & 7.06 & 2.24 & 2.14 \\ 1.34 & 2.68 & 2.24 & 0.96 & 0.80 \\ 1.29 & 2.44 & 2.14 & 0.80 & 0.74 \end{matrix}) .

Open in the MATLAB editor: f07kd_example

function f07kd_example


fprintf('f07kd example results\n\n');

% Semidefinite matrix A
a = [2.51, 4.04, 3.34, 1.34, 1.29;
     4.04, 8.22, 7.38, 2.68, 2.44;
     3.34, 7.38, 7.06, 2.24, 2.14;
     1.34, 2.68, 2.24, 0.96, 0.80;
     1.29, 2.44, 2.14, 0.80, 0.74];

% Catch warnings about rank defficient matrix ifail=1
wstat = warning();
warning('OFF');

% Factorize a
uplo = 'l';
[afac, piv, rank, info] = f07kd( ...
                                 uplo, a);

if (info==0 || info==1) 
  fprintf('Computed rank: %d\n\n', rank);
  % Zero out columns rank+1 onwards
  afac(:, rank+1:5) = 0;

  [ifail] = x04ca( ...
                   uplo, 'n', afac, 'Factor');

  fprintf('\n piv:\n   ');
  fprintf('%11d', piv);
  fprintf('\n');
else
  fprintf('\nf07kd returned with error, info = %d.\n',info);
end

warning(wstat);

f07kd example results

Computed rank: 3

 Factor
             1          2          3          4          5
 1      2.8671
 2      1.4091     0.7242
 3      2.5741    -0.3965     0.5262
 4      0.9348     0.0315    -0.2920     0.0000
 5      0.8510     0.1254    -0.0018     0.0000     0.0000

 piv:
             2          1          3          4          5

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox