Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_lapack_dpstrf (f07kd)

## Purpose

nag_lapack_dpstrf (f07kd) computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.

## Syntax

[a, piv, rank, info] = f07kd(uplo, a, 'n', n, 'tol', tol)
[a, piv, rank, info] = nag_lapack_dpstrf(uplo, a, 'n', n, 'tol', tol)

## Description

nag_lapack_dpstrf (f07kd) forms the Cholesky factorization of a real symmetric positive semidefinite matrix $A$ either as ${P}^{\mathrm{T}}AP={U}^{\mathrm{T}}U$ if ${\mathbf{uplo}}=\text{'U'}$ or ${P}^{\mathrm{T}}AP=L{L}^{\mathrm{T}}$ if ${\mathbf{uplo}}=\text{'L'}$, where $P$ is a permutation matrix, $U$ is an upper triangular matrix and $L$ is lower triangular.
This algorithm does not attempt to check that $A$ 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 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:     $\mathrm{uplo}$ – string (length ≥ 1)
Specifies whether the upper or lower triangular part of $A$ is stored and how $A$ is to be factorized.
${\mathbf{uplo}}=\text{'U'}$
The upper triangular part of $A$ is stored and $A$ is factorized as ${U}^{\mathrm{T}}U$, where $U$ is upper triangular.
${\mathbf{uplo}}=\text{'L'}$
The lower triangular part of $A$ is stored and $A$ is factorized as $L{L}^{\mathrm{T}}$, where $L$ is lower triangular.
Constraint: ${\mathbf{uplo}}=\text{'U'}$ or $\text{'L'}$.
2:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The $n$ by $n$ symmetric positive semidefinite matrix $A$.
• If ${\mathbf{uplo}}=\text{'U'}$, the upper triangular part of $a$ must be stored and the elements of the array below the diagonal are not referenced.
• If ${\mathbf{uplo}}=\text{'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:     $\mathrm{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: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{tol}$ – double scalar
Default: $-1$
User defined tolerance. If ${\mathbf{tol}}<0$, then  will be used. The algorithm terminates at the $r$th step if the $\left(r+1\right)$th step pivot $\text{}<{\mathbf{tol}}$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{uplo}}=\text{'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 ${\mathbf{uplo}}=\text{'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:     $\mathrm{piv}\left({\mathbf{n}}\right)$int64int32nag_int array
piv is such that the nonzero entries of $P$ are $P\left({\mathbf{piv}}\left(\mathit{k}\right),\mathit{k}\right)=1$, for $\mathit{k}=1,2,\dots ,n$.
3:     $\mathrm{rank}$int64int32nag_int scalar
The computed rank of $A$ given by the number of steps the algorithm completed.
4:     $\mathrm{info}$int64int32nag_int scalar
${\mathbf{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.

${\mathbf{info}}<0$
If ${\mathbf{info}}=-i$, argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.
W  ${\mathbf{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

If ${\mathbf{uplo}}=\text{'L'}$ and ${\mathbf{rank}}=r$, the computed Cholesky factor $L$ and permutation matrix $P$ satisfy the following upper bound
 $A - PLLTPT 2 A2 ≤ 2r cr ε W 2 + 1 2 + Oε2 ,$
where
 $W = L 11 -1 L12 , L = L11 0 L12 0 , L11 ∈ ℝr×r ,$
$c\left(r\right)$ is a modest linear function of $r$, $\epsilon$ is machine precision, and
 $W2 ≤ 13 n-r 4r-1 .$
So there is no guarantee of stability of the algorithm for large $n$ and $r$, although ${‖W‖}_{2}$ is generally small in practice.

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= 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 .$
```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
```