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_zhesvx (f07mp)

## Purpose

nag_lapack_zhesvx (f07mp) uses the diagonal pivoting factorization to compute the solution to a complex system of linear equations
 $AX=B ,$
where $A$ is an $n$ by $n$ Hermitian matrix and $X$ and $B$ are $n$ by $r$ matrices. Error bounds on the solution and a condition estimate are also provided.

## Syntax

[af, ipiv, x, rcond, ferr, berr, info] = f07mp(fact, uplo, a, af, ipiv, b, 'n', n, 'nrhs_p', nrhs_p)
[af, ipiv, x, rcond, ferr, berr, info] = nag_lapack_zhesvx(fact, uplo, a, af, ipiv, b, 'n', n, 'nrhs_p', nrhs_p)

## Description

nag_lapack_zhesvx (f07mp) performs the following steps:
 1 If ${\mathbf{fact}}=\text{'N'}$, the diagonal pivoting method is used to factor $A$. The form of the factorization is $A=UD{U}^{\mathrm{H}}$ if ${\mathbf{uplo}}=\text{'U'}$ or $A=LD{L}^{\mathrm{H}}$ if ${\mathbf{uplo}}=\text{'L'}$, where $U$ (or $L$) is a product of permutation and unit upper (lower) triangular matrices, and $D$ is Hermitian and block diagonal with $1$ by $1$ and $2$ by $2$ diagonal blocks. 2 If some ${d}_{ii}=0$, so that $D$ is exactly singular, then the function returns with ${\mathbf{info}}=i$. Otherwise, the factored form of $A$ is used to estimate the condition number of the matrix $A$. If the reciprocal of the condition number is less than machine precision, ${\mathbf{info}}\ge {\mathbf{n}}+1$ is returned as a warning, but the function still goes on to solve for $X$ and compute error bounds as described below. 3 The system of equations is solved for $X$ using the factored form of $A$. 4 Iterative refinement is applied to improve the computed solution matrix and to calculate error bounds and backward error estimates for it.

## References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{fact}$ – string (length ≥ 1)
Specifies whether or not the factorized form of the matrix $A$ has been supplied.
${\mathbf{fact}}=\text{'F'}$
af and ipiv contain the factorized form of the matrix $A$. af and ipiv will not be modified.
${\mathbf{fact}}=\text{'N'}$
The matrix $A$ will be copied to af and factorized.
Constraint: ${\mathbf{fact}}=\text{'F'}$ or $\text{'N'}$.
2:     $\mathrm{uplo}$ – string (length ≥ 1)
If ${\mathbf{uplo}}=\text{'U'}$, the upper triangle of $A$ is stored.
If ${\mathbf{uplo}}=\text{'L'}$, the lower triangle of $A$ is stored.
Constraint: ${\mathbf{uplo}}=\text{'U'}$ or $\text{'L'}$.
3:     $\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 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.
4:     $\mathrm{af}\left(\mathit{ldaf},:\right)$ – complex array
The first dimension of the array af must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array af must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{fact}}=\text{'F'}$, af contains the block diagonal matrix $D$ and the multipliers used to obtain the factor $U$ or $L$ from the factorization ${\mathbf{a}}=UD{U}^{\mathrm{H}}$ or ${\mathbf{a}}=LD{L}^{\mathrm{H}}$ as computed by nag_lapack_zhetrf (f07mr).
5:     $\mathrm{ipiv}\left(:\right)$int64int32nag_int array
The dimension of the array ipiv must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
If ${\mathbf{fact}}=\text{'F'}$, ipiv contains details of the interchanges and the block structure of $D$, as determined by nag_lapack_zhetrf (f07mr).
• if ${\mathbf{ipiv}}\left(i\right)=k>0$, ${d}_{ii}$ is a $1$ by $1$ pivot block and the $i$th row and column of $A$ were interchanged with the $k$th row and column;
• if ${\mathbf{uplo}}=\text{'U'}$ and ${\mathbf{ipiv}}\left(i-1\right)={\mathbf{ipiv}}\left(i\right)=-l<0$, $\left(\begin{array}{cc}{d}_{i-1,i-1}& {\stackrel{-}{d}}_{i,i-1}\\ {\stackrel{-}{d}}_{i,i-1}& {d}_{ii}\end{array}\right)$ is a $2$ by $2$ pivot block and the $\left(i-1\right)$th row and column of $A$ were interchanged with the $l$th row and column;
• if ${\mathbf{uplo}}=\text{'L'}$ and ${\mathbf{ipiv}}\left(i\right)={\mathbf{ipiv}}\left(i+1\right)=-m<0$, $\left(\begin{array}{cc}{d}_{ii}& {d}_{i+1,i}\\ {d}_{i+1,i}& {d}_{i+1,i+1}\end{array}\right)$ is a $2$ by $2$ pivot block and the $\left(i+1\right)$th row and column of $A$ were interchanged with the $m$th row and column.
6:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – complex array
The first dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$.
The $n$ by $r$ right-hand side matrix $B$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the arrays a, af, b and the second dimension of the arrays a, af, ipiv.
$n$, the number of linear equations, i.e., the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{nrhs_p}$int64int32nag_int scalar
Default: the second dimension of the array b.
$r$, the number of right-hand sides, i.e., the number of columns of the matrix $B$.
Constraint: ${\mathbf{nrhs_p}}\ge 0$.

### Output Parameters

1:     $\mathrm{af}\left(\mathit{ldaf},:\right)$ – complex array
The first dimension of the array af will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array af will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{fact}}=\text{'N'}$, af returns the block diagonal matrix $D$ and the multipliers used to obtain the factor $U$ or $L$ from the factorization ${\mathbf{a}}=UD{U}^{\mathrm{H}}$ or ${\mathbf{a}}=LD{L}^{\mathrm{H}}$.
2:     $\mathrm{ipiv}\left(:\right)$int64int32nag_int array
The dimension of the array ipiv will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
If ${\mathbf{fact}}=\text{'N'}$, ipiv contains details of the interchanges and the block structure of $D$, as determined by nag_lapack_zhetrf (f07mr), as described above.
3:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – complex array
The first dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$.
If ${\mathbf{info}}={\mathbf{0}}$ or $\mathbf{n}+{\mathbf{1}}$, the $n$ by $r$ solution matrix $X$.
4:     $\mathrm{rcond}$ – double scalar
The estimate of the reciprocal condition number of the matrix $A$. If ${\mathbf{rcond}}=0.0$, the matrix may be exactly singular. This condition is indicated by ${\mathbf{info}}>{\mathbf{0}} \text{and} {\mathbf{info}}\le \mathbf{n}$. Otherwise, if rcond is less than the machine precision, the matrix is singular to working precision. This condition is indicated by ${\mathbf{info}}\ge {\mathbf{n}}+1$.
5:     $\mathrm{ferr}\left(:\right)$ – double array
The dimension of the array ferr will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$
If ${\mathbf{info}}={\mathbf{0}}$ or $\mathbf{n}+{\mathbf{1}}$, an estimate of the forward error bound for each computed solution vector, such that ${‖{\stackrel{^}{x}}_{j}-{x}_{j}‖}_{\infty }/{‖{x}_{j}‖}_{\infty }\le {\mathbf{ferr}}\left(j\right)$ where ${\stackrel{^}{x}}_{j}$ is the $j$th column of the computed solution returned in the array x and ${x}_{j}$ is the corresponding column of the exact solution $X$. The estimate is as reliable as the estimate for rcond, and is almost always a slight overestimate of the true error.
6:     $\mathrm{berr}\left(:\right)$ – double array
The dimension of the array berr will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$
If ${\mathbf{info}}={\mathbf{0}}$ or $\mathbf{n}+{\mathbf{1}}$, an estimate of the component-wise relative backward error of each computed solution vector ${\stackrel{^}{x}}_{j}$ (i.e., the smallest relative change in any element of $A$ or $B$ that makes ${\stackrel{^}{x}}_{j}$ an exact solution).
7:     $\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}}>0 \text{and} {\mathbf{info}}\le {\mathbf{n}}$
Element $_$ of the diagonal is exactly zero. The factorization has been completed, but the factor $D$ is exactly singular, so the solution and error bounds could not be computed. ${\mathbf{rcond}}=0.0$ is returned.
W  ${\mathbf{info}}={\mathbf{n}}+1$
$D$ is nonsingular, but rcond is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of rcond would suggest.

## Accuracy

For each right-hand side vector $b$, the computed solution $\stackrel{^}{x}$ is the exact solution of a perturbed system of equations $\left(A+E\right)\stackrel{^}{x}=b$, where
 $E1 = Oε A1 ,$
where $\epsilon$ is the machine precision. See Chapter 11 of Higham (2002) for further details.
If $\stackrel{^}{x}$ is the true solution, then the computed solution $x$ satisfies a forward error bound of the form
 $x-x^∞ x^∞ ≤ wc condA,x^,b$
where $\mathrm{cond}\left(A,\stackrel{^}{x},b\right)={‖\left|{A}^{-1}\right|\left(\left|A\right|\left|\stackrel{^}{x}\right|+\left|b\right|\right)‖}_{\infty }/{‖\stackrel{^}{x}‖}_{\infty }\le \mathrm{cond}\left(A\right)={‖\left|{A}^{-1}\right|\left|A\right|‖}_{\infty }\le {\kappa }_{\infty }\left(A\right)$. If $\stackrel{^}{x}$ is the $j$th column of $X$, then ${w}_{c}$ is returned in ${\mathbf{berr}}\left(j\right)$ and a bound on ${‖x-\stackrel{^}{x}‖}_{\infty }/{‖\stackrel{^}{x}‖}_{\infty }$ is returned in ${\mathbf{ferr}}\left(j\right)$. See Section 4.4 of Anderson et al. (1999) for further details.

The factorization of $A$ requires approximately $\frac{4}{3}{n}^{3}$ floating-point operations.
For each right-hand side, computation of the backward error involves a minimum of $16{n}^{2}$ floating-point operations. Each step of iterative refinement involves an additional $24{n}^{2}$ operations. At most five steps of iterative refinement are performed, but usually only one or two steps are required. Estimating the forward error involves solving a number of systems of equations of the form $Ax=b$; the number is usually $4$ or $5$ and never more than $11$. Each solution involves approximately $8{n}^{2}$ operations.
The real analogue of this function is nag_lapack_dsysvx (f07mb). The complex symmetric analogue of this function is nag_lapack_zsysvx (f07np).

## Example

This example solves the equations
 $AX=B ,$
where $A$ is the Hermitian matrix
 $A = -1.84i+0.00 0.11-0.11i -1.78-1.18i 3.91-1.50i 0.11+0.11i -4.63i+0.00 -1.84+0.03i 2.21+0.21i -1.78+1.18i -1.84-0.03i -8.87i+0.00 1.58-0.90i 3.91+1.50i 2.21-0.21i 1.58+0.90i -1.36i+0.00$
and
 $B = 2.98-10.18i 28.68-39.89i -9.58+03.88i -24.79-08.40i -0.77-16.05i 4.23-70.02i 7.79+05.48i -35.39+18.01i .$
Error estimates for the solutions, and an estimate of the reciprocal of the condition number of the matrix $A$ are also output.
```function f07mp_example

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

% Hermitian indefinite matrix A (Upper triangular part stored)
uplo = 'Upper';
a = [-1.84 + 0i,  0.11 - 0.11i, -1.78 - 1.18i,  3.91 - 1.5i;
0    + 0i, -4.63 + 0i,    -1.84 + 0.03i,  2.21 + 0.21i;
0    + 0i,  0    + 0i,    -8.87 + 0i,     1.58 - 0.9i;
0    + 0i,  0    + 0i,     0    + 0i,    -1.36 + 0i];

% RHS
b = [ 2.98 - 10.18i,   28.68 - 39.89i;
-9.58 +  3.88i,  -24.79 -  8.40i;
-0.77 - 16.05i,    4.23 - 70.02i;
7.79 +  5.48i,  -35.39 + 18.01i];

% Input parameters
fact = 'Not factored';
af = a;
ipiv = zeros(size(b,1), 1, 'int64');

% Solve
[af, ipiv, x, rcond, ferr, berr, info] = ...
f07mp( ...
fact, uplo, a, af, ipiv, b);

disp('Solution(s)');
disp(x);
disp('Backward errors (machine-dependent)');
fprintf('%10.1e',berr);
fprintf('\n');
disp('Estimated forward error bounds (machine-dependent)');
fprintf('%10.1e',ferr);
fprintf('\n\n');
disp('Estimate of reciprocal condition number');
fprintf('%10.1e\n\n',rcond);

```
```f07mp example results

Solution(s)
2.0000 + 1.0000i  -8.0000 + 6.0000i
3.0000 - 2.0000i   7.0000 - 2.0000i
-1.0000 + 2.0000i  -1.0000 + 5.0000i
1.0000 - 1.0000i   3.0000 - 4.0000i

Backward errors (machine-dependent)
7.3e-17   8.1e-17
Estimated forward error bounds (machine-dependent)
2.6e-15   3.1e-15

Estimate of reciprocal condition number
1.5e-01

```