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

lda work

### Output Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldamax (1,n)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If uplo = 'U'${\mathbf{uplo}}=\text{'U'}$, the first rank rows of the upper triangle of A$A$ are overwritten with the nonzero elements of the Cholesky factor U$U$, and the remaining rows of the triangle are destroyed.
If uplo = 'L'${\mathbf{uplo}}=\text{'L'}$, the first rank columns of the lower triangle of A$A$ are overwritten with the nonzero elements of the Cholesky factor L$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$P$ are P(piv(k),k) = 1$P\left({\mathbf{piv}}\left(\mathit{k}\right),\mathit{k}\right)=1$, for k = 1,2,,n$\mathit{k}=1,2,\dots ,n$.
3:     rank – int64int32nag_int scalar
The computed rank of A$A$ given by the number of steps the algorithm completed.
4:     info – int64int32nag_int scalar
info = 0${\mathbf{info}}=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 > 0${\mathbf{INFO}}>0$
The leading minor of order _$_$ is not positive definite and the factorization could not be completed. Hence A$A$ itself is not positive definite. This may indicate an error in forming the matrix A$A$. To factorize a symmetric matrix which is not positive definite, call nag_lapack_dsytrf (f07md) instead.
W INFO = 1${\mathbf{INFO}}=1$
The matrix A$A$ is either rank deficient with computed rank as returned in rank, or is indefinite, see Section [Further Comments].

## Accuracy

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

The total number of floating point operations is approximately nr22 / 3r3$n{r}^{2}-2/3{r}^{3}$, where r$r$ is the computed rank of A$A$.
This algorithm does not attempt to check that A$A$ is positive semidefinite, and in particular the rank detection criterion in the algorithm is based on A$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

```function nag_lapack_dpstrf_example
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];
uplo = 'l';
% Factorize a
[aOut, piv, rank, info] = nag_lapack_dpstrf(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_real_gen(uplo, 'n', aOut, 'Factor');
fprintf('\n piv:\n   ');
fprintf('          %d', piv);
fprintf('\n');
```
```
Warning: nag_lapack_dpstrf (f07kd) returned a warning indicator (1)

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

```
```function f07kd_example
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];
uplo = 'l';
% Factorize a
[aOut, piv, rank, info] = f07kd(uplo, a);

fprintf('\nComputed rank: %d\n\n', rank);
% Zero out columns rank+1 onwards
aOut(:, double(rank)+1:5) = 0;
[ifail] = x04ca(uplo, 'n', aOut, 'Factor');
fprintf('\n piv:\n   ');
fprintf('          %d', piv);
fprintf('\n');
```
```
Warning: nag_lapack_dpstrf (f07kd) returned a warning indicator (1)

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

```