F07 Chapter Contents
F07 Chapter Introduction
NAG Library Manual

# NAG Library Routine DocumentF07PPF (ZHPSVX)

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

## 1  Purpose

F07PPF (ZHPSVX) uses the diagonal pivoting factorization
 $A=UDUH or A=LDLH$
to compute the solution to a complex system of linear equations
 $AX=B ,$
where $A$ is an $n$ by $n$ Hermitian matrix stored in packed format and $X$ and $B$ are $n$ by $r$ matrices. Error bounds on the solution and a condition estimate are also provided.

## 2  Specification

 SUBROUTINE F07PPF ( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
 INTEGER N, NRHS, IPIV(N), LDB, LDX, INFO REAL (KIND=nag_wp) RCOND, FERR(NRHS), BERR(NRHS), RWORK(N) COMPLEX (KIND=nag_wp) AP(*), AFP(*), B(LDB,*), X(LDX,*), WORK(2*N) CHARACTER(1) FACT, UPLO
The routine may be called by its LAPACK name zhpsvx.

## 3  Description

F07PPF (ZHPSVX) performs the following steps:
1. If ${\mathbf{FACT}}=\text{'N'}$, the diagonal pivoting method is used to factor $A$ as $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 routine 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}}=\mathbf{N}+{\mathbf{1}}$ is returned as a warning, but the routine 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.

## 4  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

## 5  Parameters

1:     FACT – CHARACTER(1)Input
On entry: specifies whether or not the factorized form of the matrix $A$ has been supplied.
${\mathbf{FACT}}=\text{'F'}$
AFP and IPIV contain the factorized form of the matrix $A$. AFP and IPIV will not be modified.
${\mathbf{FACT}}=\text{'N'}$
The matrix $A$ will be copied to AFP and factorized.
Constraint: ${\mathbf{FACT}}=\text{'F'}$ or $\text{'N'}$.
2:     UPLO – CHARACTER(1)Input
On entry: 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:     N – INTEGERInput
On entry: $n$, the number of linear equations, i.e., the order of the matrix $A$.
Constraint: ${\mathbf{N}}\ge 0$.
4:     NRHS – INTEGERInput
On entry: $r$, the number of right-hand sides, i.e., the number of columns of the matrix $B$.
Constraint: ${\mathbf{NRHS}}\ge 0$.
5:     AP($*$) – COMPLEX (KIND=nag_wp) arrayInput
Note: the dimension of the array AP must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}×\left({\mathbf{N}}+1\right)/2\right)$.
On entry: the $n$ by $n$ Hermitian matrix $A$, packed by columns.
More precisely,
• if ${\mathbf{UPLO}}=\text{'U'}$, the upper triangle of $A$ must be stored with element ${A}_{ij}$ in ${\mathbf{AP}}\left(i+j\left(j-1\right)/2\right)$ for $i\le j$;
• if ${\mathbf{UPLO}}=\text{'L'}$, the lower triangle of $A$ must be stored with element ${A}_{ij}$ in ${\mathbf{AP}}\left(i+\left(2n-j\right)\left(j-1\right)/2\right)$ for $i\ge j$.
6:     AFP($*$) – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the dimension of the array AFP must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}×\left({\mathbf{N}}+1\right)/2\right)$.
On entry: if ${\mathbf{FACT}}=\text{'F'}$, AFP contains the block diagonal matrix $D$ and the multipliers used to obtain the factor $U$ or $L$ from the factorization $A=UD{U}^{\mathrm{H}}$ or $A=LD{L}^{\mathrm{H}}$ as computed by F07PRF (ZHPTRF), stored as a packed triangular matrix in the same storage format as $A$.
On exit: if ${\mathbf{FACT}}=\text{'N'}$, AFP contains the block diagonal matrix $D$ and the multipliers used to obtain the factor $U$ or $L$ from the factorization $A=UD{U}^{\mathrm{H}}$ or $A=LD{L}^{\mathrm{H}}$ as computed by F07PRF (ZHPTRF), stored as a packed triangular matrix in the same storage format as $A$.
7:     IPIV(N) – INTEGER arrayInput/Output
On entry: if ${\mathbf{FACT}}=\text{'F'}$, IPIV contains details of the interchanges and the block structure of $D$, as determined by F07PRF (ZHPTRF).
• 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.
On exit: if ${\mathbf{FACT}}=\text{'N'}$, IPIV contains details of the interchanges and the block structure of $D$, as determined by F07PRF (ZHPTRF), as described above.
8:     B(LDB,$*$) – COMPLEX (KIND=nag_wp) arrayInput
Note: the second dimension of the array B must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On entry: the $n$ by $r$ right-hand side matrix $B$.
9:     LDB – INTEGERInput
On entry: the first dimension of the array B as declared in the (sub)program from which F07PPF (ZHPSVX) is called.
Constraint: ${\mathbf{LDB}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
10:   X(LDX,$*$) – COMPLEX (KIND=nag_wp) arrayOutput
Note: the second dimension of the array X must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On exit: if ${\mathbf{INFO}}={\mathbf{0}}$ or $\mathbf{N}+{\mathbf{1}}$, the $n$ by $r$ solution matrix $X$.
11:   LDX – INTEGERInput
On entry: the first dimension of the array X as declared in the (sub)program from which F07PPF (ZHPSVX) is called.
Constraint: ${\mathbf{LDX}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
12:   RCOND – REAL (KIND=nag_wp)Output
On exit: 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 . Otherwise, if RCOND is less than the machine precision, the matrix is singular to working precision. This condition is indicated by ${\mathbf{INFO}}=\mathbf{N}+{\mathbf{1}}$.
13:   FERR(NRHS) – REAL (KIND=nag_wp) arrayOutput
On exit: 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.
14:   BERR(NRHS) – REAL (KIND=nag_wp) arrayOutput
On exit: 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).
15:   WORK($2×{\mathbf{N}}$) – COMPLEX (KIND=nag_wp) arrayWorkspace
16:   RWORK(N) – REAL (KIND=nag_wp) arrayWorkspace
17:   INFO – INTEGEROutput
On exit: ${\mathbf{INFO}}=0$ unless the routine detects an error (see Section 6).

## 6  Error Indicators and Warnings

Errors or warnings detected by the routine:
${\mathbf{INFO}}<0$
If ${\mathbf{INFO}}=-i$, the $i$th argument had an illegal value. An explanatory message is output, and execution of the program is terminated.
If ${\mathbf{INFO}}\le {\mathbf{N}}$, $d\left(i,i\right)$ 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.
${\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.

## 7  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 routine is F07PBF (DSPSVX). The complex symmetric analogue of this routine is F07QPF (ZSPSVX).

## 9  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.

### 9.1  Program Text

Program Text (f07ppfe.f90)

### 9.2  Program Data

Program Data (f07ppfe.d)

### 9.3  Program Results

Program Results (f07ppfe.r)