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$ either as ${P}^{\mathrm{T}}AP={U}^{\mathrm{H}}U$ if ${\mathbf{uplo}}=\text{'U'}$ or ${P}^{\mathrm{T}}AP=L{L}^{\mathrm{H}}$ 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{H}}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{H}}$, where $L$ is lower triangular.
Constraint: ${\mathbf{uplo}}=\text{'U'}$ or $\text{'L'}$.
2:     $\mathrm{a}\left(\mathit{lda},:\right)$ – complex 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$ Hermitian 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)$ – complex 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 - PLLHPT 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 real floating-point operations is approximately $4n{r}^{2}-8/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_zhetrf (f07mr). See Lucas (2004) for further information.
The real analogue of this function is nag_lapack_dpstrf (f07kd).

## Example

This example computes the Cholesky factorization of the matrix $A$, where
 $A= 12.40+0.00i 2.39+0.00i 5.50+0.05i 4.47+0.00i 11.89+0.00i 2.39+0.00i 1.63+0.00i 1.04+0.10i 1.14+0.00i 1.81+0.00i 5.50+0.05i 1.04+0.10i 2.45+0.00i 1.98-0.03i 5.28-0.02i 4.47+0.00i 1.14+0.00i 1.98-0.03i 1.71+0.00i 4.14+0.00i 11.89+0.00i 1.81+0.00i 5.28-0.02i 4.14+0.00i 11.63+0.00i .$
function f07kr_example

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

% Semidefinite matrix A
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];

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

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

fprintf('\nComputed rank: %d\n\n', rank);
% Zero out columns rank+1 onwards and display
afac(:, rank+1:5) = 0;

ncol   = int64(80);
indent = int64(0);
[ifail] = x04db( ...
uplo, 'non-unit', afac, 'bracketed', 'f5.2', 'Factor', ...
'Integer', 'Integer', ncol, indent);

fprintf('\n piv:\n   ');
fprintf('%14d', piv);
fprintf('\n');

warning(wstat);

f07kr example results

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