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_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 A$A$ either as PTAP = UHU${P}^{\mathrm{T}}AP={U}^{\mathrm{H}}U$ if uplo = 'U'${\mathbf{uplo}}=\text{'U'}$ or PTAP = LLH${P}^{\mathrm{T}}AP=L{L}^{\mathrm{H}}$ 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 UHU${U}^{\mathrm{H}}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 LLH$L{L}^{\mathrm{H}}$, where L$L$ is lower triangular.
Constraint: uplo = 'U'${\mathbf{uplo}}=\text{'U'}$ or 'L'$\text{'L'}$.
2:     a(lda, : $:$) – complex 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$ Hermitian 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, : $:$) – complex 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 Hermitian matrix which is not positive definite, call nag_lapack_zhetrf (f07mr) 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 − PLLHPT‖2)/(‖A‖2) ≤ 2r c(r) ε (‖W‖2 + 1)2 + O(ε2) , $‖ A - PLLHPT ‖ 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 precision, 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 real floating point operations is approximately 4nr28 / 3r3$4n{r}^{2}-8/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_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

```