# NAG CL InterfaceF07 (Lapacklin)Linear Equations (LAPACK)

Settings help

CL Name Style:

## 1Scope of the Chapter

This chapter provides functions for the solution of systems of simultaneous linear equations, and associated computations. It provides functions for
• matrix factorizations;
• solution of linear equations;
• estimating matrix condition numbers;
• computing error bounds for the solution of linear equations;
• matrix inversion.
Functions are provided for both real and complex data.
For a general introduction to the solution of systems of linear equations, you should turn first to the F04 Chapter Introduction. The decision trees, in Section 4 in the F04 Chapter Introduction, direct you to the most appropriate functions in Chapters F04 or F07 for solving your particular problem. In particular, Chapters F04 and F07 contain Black Box (or driver) functions which enable some standard types of problem to be solved by a call to a single function. Where possible, functions in Chapter F04 call Chapter F07 functions to perform the necessary computational tasks.
There are two types of driver functions in this chapter: simple drivers which just return the solution to the linear equations; and expert drivers which also return condition and error estimates and, in many cases, also allow equilibration. The simple drivers for real matrices have names of the form nag_d..sv (f07.ac) and for complex matrices have names of the form nag_z..sv (f07.nc). The expert drivers for real matrices have names of the form nag_d..svx (f07.bc) and for complex matrices have names of the form nag_z..svx (f07.pc).
The functions in this chapter (Chapter F07) handle only dense and band matrices (not matrices with more specialised structures, or general sparse matrices).
The functions in this chapter have all been derived from the LAPACK project (see Anderson et al. (1999)). They have been designed to be efficient on a wide range of high-performance computers, without compromising efficiency on conventional serial machines.

## 2Background to the Problems

This section is only a brief introduction to the numerical solution of systems of linear equations. Consult a standard textbook, for example Golub and Van Loan (1996) for a more thorough discussion.

### 2.1Notation

We use the standard notation for a system of simultaneous linear equations:
 $Ax=b$ (1)
where $A$ is the coefficient matrix, $b$ is the right-hand side, and $x$ is the solution. $A$ is assumed to be a square matrix of order $n$.
If there are several right-hand sides, we write
 $AX=B$ (2)
where the columns of $B$ are the individual right-hand sides, and the columns of $X$ are the corresponding solutions.
We also use the following notation, both here and in the function documents:
 $\stackrel{^}{x}$ a computed solution to $Ax=b$, (which usually differs from the exact solution $x$ because of round-off error) $r=b-A\stackrel{^}{x}$ the residual corresponding to the computed solution $\stackrel{^}{x}$ ${‖x‖}_{\infty }=\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}|{x}_{i}|$ the $\infty$-norm of the vector $x$ ${‖x‖}_{1}=\sum _{j=1}^{n}|{x}_{j}|$ the $1$-norm of the vector $x$ ${‖A‖}_{\infty }=\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\sum _{j}\phantom{\rule{0.25em}{0ex}}|{a}_{ij}|$ the $\infty$-norm of the matrix $A$ ${‖A‖}_{1}=\underset{j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\sum _{i=1}^{n}|{a}_{ij}|$ the $1$-norm of the matrix $A$ $|x|$ the vector with elements $|{x}_{i}|$ $|A|$ the matrix with elements $|{a}_{ij}|$
Inequalities of the form $|A|\le |B|$ are interpreted component-wise, that is $|{a}_{ij}|\le |{b}_{ij}|$ for all $i,j$.

### 2.2Matrix Factorizations

If $A$ is upper or lower triangular, $Ax=b$ can be solved by a straightforward process of backward or forward substitution.
Otherwise, the solution is obtained after first factorizing $A$, as follows.
General matrices (LU factorization with partial pivoting)
 $A=PLU$
where $P$ is a permutation matrix, $L$ is lower-triangular with diagonal elements equal to $1$, and $U$ is upper-triangular; the permutation matrix $P$ (which represents row interchanges) is needed to ensure numerical stability.
Symmetric positive definite matrices (Cholesky factorization)
 $A=UTU or A=LLT$
where $U$ is upper triangular and $L$ is lower triangular.
Symmetric positive semidefinite matrices (pivoted Cholesky factorization)
 $A=PUTUPT or A=PLLTPT$
where $P$ is a permutation matrix, $U$ is upper triangular and $L$ is lower triangular. The permutation matrix $P$ (which represents row-and-column interchanges) is needed to ensure numerical stability and to reveal the numerical rank of $A$.
Symmetric indefinite matrices (Bunch–Kaufman factorization)
 $A = PUD UT PT or A = PLD LT PT$
where $P$ is a permutation matrix, $U$ is upper triangular, $L$ is lower triangular, and $D$ is a block diagonal matrix with diagonal blocks of order $1$ or $2$; $U$ and $L$ have diagonal elements equal to $1$, and have $2×2$ unit matrices on the diagonal corresponding to the $2×2$ blocks of $D$. The permutation matrix $P$ (which represents symmetric row-and-column interchanges) and the $2×2$ blocks in $D$ are needed to ensure numerical stability. If $A$ is in fact positive definite, no interchanges are needed and the factorization reduces to $A=UD{U}^{\mathrm{T}}$ or $A=LD{L}^{\mathrm{T}}$ with diagonal $D$, which is simply a variant form of the Cholesky factorization.

### 2.3Solution of Systems of Equations

Given one of the above matrix factorizations, it is straightforward to compute a solution to $Ax=b$ by solving two subproblems, as shown below, first for $y$ and then for $x$. Each subproblem consists essentially of solving a triangular system of equations by forward or backward substitution; the permutation matrix $P$ and the block diagonal matrix $D$ introduce only a little extra complication:
General matrices ( LU factorization)
 $Ly=PTb Ux=y$
Symmetric positive definite matrices (Cholesky factorization)
 $UTy=b Ux=y or Ly=b LTx=y$
Symmetric indefinite matrices (Bunch–Kaufman factorization)
 $PUDy=b UTPTx=y or PLDy=b LTPTx=y$

### 2.4Sensitivity and Error Analysis

#### 2.4.1Normwise error bounds

Frequently, in practical problems the data $A$ and $b$ are not known exactly, and it is then important to understand how uncertainties or perturbations in the data can affect the solution.
If $x$ is the exact solution to $Ax=b$, and $x+\delta x$ is the exact solution to a perturbed problem $\left(A+\delta A\right)\left(x+\delta x\right)=\left(b+\delta b\right)$, then
 $‖δx‖ ‖x‖ ≤κ(A) (‖δA‖ ‖A‖ +‖δb‖ ‖b‖ )+⋯(second-order terms)$
where $\kappa \left(A\right)$ is the condition number of $A$ defined by
 $κ(A) = ‖A‖.‖A-1‖ .$ (3)
In other words, relative errors in $A$ or $b$ may be amplified in $x$ by a factor $\kappa \left(A\right)$. Section 2.4.2 discusses how to compute or estimate $\kappa \left(A\right)$.
Similar considerations apply when we study the effects of rounding errors introduced by computation in finite precision. The effects of rounding errors can be shown to be equivalent to perturbations in the original data, such that $\frac{‖\delta A‖}{‖A‖}$ and $\frac{‖\delta b‖}{‖b‖}$ are usually at most $p\left(n\right)\epsilon$, where $\epsilon$ is the machine precision and $p\left(n\right)$ is an increasing function of $n$ which is seldom larger than $10n$ (although in theory it can be as large as ${2}^{n-1}$).
In other words, the computed solution $\stackrel{^}{x}$ is the exact solution of a linear system $\left(A+\delta A\right)\stackrel{^}{x}=b+\delta b$ which is close to the original system in a normwise sense.

#### 2.4.2Estimating condition numbers

The previous section has emphasized the usefulness of the quantity $\kappa \left(A\right)$ in understanding the sensitivity of the solution of $Ax=b$. To compute the value of $\kappa \left(A\right)$ from equation (3) is more expensive than solving $Ax=b$ in the first place. Hence it is standard practice to estimate $\kappa \left(A\right)$, in either the $1$-norm or the $\infty$-norm, by a method which only requires $\mathit{O}\left({n}^{2}\right)$ additional operations, assuming that a suitable factorization of $A$ is available.
The method used in this chapter is Higham's modification of Hager's method (see Higham (1988)). It yields an estimate which is never larger than the true value, but which seldom falls short by more than a factor of $3$ (although artificial examples can be constructed where it is much smaller). This is acceptable since it is the order of magnitude of $\kappa \left(A\right)$ which is important rather than its precise value.
Because $\kappa \left(A\right)$ is infinite if $A$ is singular, the functions in this chapter actually return the reciprocal of $\kappa \left(A\right)$.

#### 2.4.3Scaling and Equilibration

The condition of a matrix and hence the accuracy of the computed solution, may be improved by scaling; thus if ${D}_{1}$ and ${D}_{2}$ are diagonal matrices with positive diagonal elements, then
 $B = D1 A D2$
is the scaled matrix. A general matrix is said to be equilibrated if it is scaled so that the lengths of its rows and columns have approximately equal magnitude. Similarly a general matrix is said to be row-equilibrated (column-equilibrated) if it is scaled so that the lengths of its rows (columns) have approximately equal magnitude. Note that row scaling can affect the choice of pivot when partial pivoting is used in the factorization of $A$.
A symmetric or Hermitian positive definite matrix is said to be equilibrated if the diagonal elements are all approximately equal to unity.
For further information on scaling and equilibration see Section 3.5.2 of Golub and Van Loan (1996), Section 7.2, 7.3 and 9.8 of Higham (1988) and Section 5 of Chapter 4 of Wilkinson (1965).
Functions are provided to return the scaling factors that equilibrate a matrix for general, general band, symmetric and Hermitian positive definite and symmetric and Hermitian positive definite band matrices.

#### 2.4.4Componentwise error bounds

A disadvantage of normwise error bounds is that they do not reflect any special structure in the data $A$ and $b$ – that is, a pattern of elements which are known to be zero – and the bounds are dominated by the largest elements in the data.
Componentwise error bounds overcome these limitations. Instead of the normwise relative error, we can bound the relative error in each component of $A$ and $b$:
 $maxijk (|δaij| |aij| ,|δbk| |bk| ) ≤ω$
where the component-wise backward error bound $\omega$ is given by
 $ω= maxi |ri| (|A|.|x^|+|b|)i .$
Functions are provided in this chapter which compute $\omega$, and also compute a forward error bound which is sometimes much sharper than the normwise bound given earlier:
 $‖x-x^‖∞ ‖x‖∞ ≤ ‖|A-1|.|r|‖∞ ‖x‖∞ .$
Care is taken when computing this bound to allow for rounding errors in computing $r$. The norm ${‖|{A}^{-1}|.|r|‖}_{\infty }$ is estimated cheaply (without computing ${A}^{-1}$) by a modification of the method used to estimate $\kappa \left(A\right)$.

#### 2.4.5Iterative refinement of the solution

If $\stackrel{^}{x}$ is an approximate computed solution to $Ax=b$, and $r$ is the corresponding residual, then a procedure for iterative refinement of $\stackrel{^}{x}$ can be defined as follows, starting with ${x}_{0}=\stackrel{^}{x}$:
• for $i=0,1,\dots \text{}$, until convergence
 compute ${r}_{i}=b-A{x}_{i}$ solve $A{d}_{i}={r}_{i}$ compute ${x}_{i+1}={x}_{i}+{d}_{i}$
In Chapter F04, functions are provided which perform this procedure using additional precision to compute $r$, and are thus able to reduce the forward error to the level of machine precision.
The functions in this chapter do not use additional precision to compute $r$, and cannot guarantee a small forward error, but can guarantee a small backward error (except in rare cases when $A$ is very ill-conditioned, or when $A$ and $x$ are sparse in such a way that $|A|.|x|$ has a zero or very small component). The iterations continue until the backward error has been reduced as much as possible; usually only one iteration is needed.

### 2.5Matrix Inversion

It is seldom necessary to compute an explicit inverse of a matrix. In particular, do not attempt to solve $Ax=b$ by first computing ${A}^{-1}$ and then forming the matrix-vector product $x={A}^{-1}b$; the procedure described in Section 2.3 is more efficient and more accurate.
However, functions are provided for the rare occasions when an inverse is needed, using one of the factorizations described in Section 2.2.

### 2.6Packed Storage Formats

Functions which handle symmetric matrices are usually designed so that they use either the upper or lower triangle of the matrix; it is not necessary to store the whole matrix. If the upper or lower triangle is stored conventionally in the upper or lower triangle of a two-dimensional array, the remaining elements of the array can be used to store other useful data.
However, that is not always convenient, and if it is important to economize on storage, the upper or lower triangle can be stored in a one-dimensional array of length $n\left(n+1\right)/2$ or a two-dimensional array with $n\left(n+1\right)/2$ elements; in other words, the storage is almost halved.
The one-dimensional array storage format is referred to as packed storage; it is described in Section 3.4.2. The two-dimensional array storage format is referred to as Rectangular Full Packed (RFP) format; it is described in Section 3.4.3. They may also be used for triangular matrices.
Functions designed for these packed storage formats perform the same number of arithmetic operations as functions which use conventional storage. Those using a packed one-dimensional array are usually less efficient, especially on high-performance computers, so there is then a trade-off between storage and efficiency. The RFP functions are as efficient as for conventional storage, although only a small subset of functions use this format.

### 2.7Band and Tridiagonal Matrices

A band matrix is one whose nonzero elements are confined to a relatively small number of subdiagonals or superdiagonals on either side of the main diagonal. A tridiagonal matrix is a special case of a band matrix with just one subdiagonal and one superdiagonal. Algorithms can take advantage of bandedness to reduce the amount of work and storage required. The storage scheme used for band matrices is described in Section 3.4.4.
The $LU$ factorization for general matrices, and the Cholesky factorization for symmetric and Hermitian positive definite matrices both preserve bandedness. Hence functions are provided which take advantage of the band structure when solving systems of linear equations.
The Cholesky factorization preserves bandedness in a very precise sense: the factor $U$ or $L$ has the same number of superdiagonals or subdiagonals as the original matrix. In the $LU$ factorization, the row-interchanges modify the band structure: if $A$ has ${k}_{l}$ subdiagonals and ${k}_{u}$ superdiagonals, then $L$ is not a band matrix but still has at most ${k}_{l}$ nonzero elements below the diagonal in each column; and $U$ has at most ${k}_{l}+{k}_{u}$ superdiagonals.
The Bunch–Kaufman factorization does not preserve bandedness, because of the need for symmetric row-and-column permutations; hence no functions are provided for symmetric indefinite band matrices.
The inverse of a band matrix does not in general have a band structure, so no functions are provided for computing inverses of band matrices.

### 2.8Block Partitioned Algorithms

Many of the functions in this chapter use what is termed a block partitioned algorithm. This means that at each major step of the algorithm a block of rows or columns is updated, and most of the computation is performed by matrix-matrix operations on these blocks. The matrix-matrix operations are performed by calls to the Level 3 BLAS (see Chapter F16), which are the key to achieving high performance on many modern computers. See Golub and Van Loan (1996) or Anderson et al. (1999) for more about block partitioned algorithms.
The performance of a block partitioned algorithm varies to some extent with the block size – that is, the number of rows or columns per block. This is a machine-dependent argument, which is set to a suitable value when the Library is implemented on each range of machines. You do not normally need to be aware of what value is being used. Different block sizes may be used for different functions. Values in the range $16$ to $64$ are typical.
On some machines there may be no advantage from using a block partitioned algorithm, and then the functions use an unblocked algorithm (effectively a block size of $1$), relying solely on calls to the Level 2 BLAS (see Chapter F16 again).

### 2.9Mixed Precision LAPACK Routines

Some LAPACK routines use mixed precision arithmetic in an effort to solve problems more efficiently on modern hardware. They work by converting a double precision problem into an equivalent single precision problem, solving it and then using iterative refinement in double precision to find a full precision solution to the original problem. The method may fail if the problem is too ill-conditioned to allow the initial single precision solution, in which case the functions fall back to solve the original problem entirely in double precision. The vast majority of problems are not so ill-conditioned, and in those cases the technique can lead to significant gains in speed without loss of accuracy. This is particularly true on machines where double precision arithmetic is significantly slower than single precision.

## 3Recommendations on Choice and Use of Available Functions

### 3.1Available Functions

Tables 1 and 8 in Section 3.6 show the functions which are provided for performing different computations on different types of matrices. Table 1 shows functions for real matrices; Table 8 shows functions for complex matrices. Each entry in the table gives the NAG function name and the LAPACK double precision name (see Section 3.3).
Functions are provided for the following types of matrix:
• general
• general band
• symmetric or Hermitian positive definite
• symmetric or Hermitian positive definite (packed storage)
• symmetric or Hermitian positive definite (RFP storage)
• symmetric or Hermitian positive definite band
• symmetric or Hermitian positive definite tridiagonal
• symmetric or Hermitian indefinite
• symmetric or Hermitian indefinite (packed storage)
• triangular
• triangular (packed storage)
• triangular (RFP storage)
• triangular band
For each of the above types of matrix (except where indicated), functions are provided to perform the following computations:
1. (a)(except for RFP matrices) solve a system of linear equations (driver functions);
2. (b)(except for RFP matrices) solve a system of linear equations with condition and error estimation (expert drivers);
3. (c)(except for triangular matrices) factorize the matrix (see Section 2.2);
4. (d)solve a system of linear equations, using the factorization (see Section 2.3);
5. (e)(except for RFP matrices) estimate the condition number of the matrix, using the factorization (see Section 2.4.2); these functions also require the norm of the original matrix (except when the matrix is triangular) which may be computed by a function in Chapter F16;
6. (f)(except for RFP matrices) refine the solution and compute forward and backward error bounds (see Sections 2.4.4 and 2.4.5); these functions require the original matrix and right-hand side, as well as the factorization returned from (a) and the solution returned from (b);
7. (g)(except for band and tridiagonal matrices) invert the matrix, using the factorization (see Section 2.5).
Thus, to solve a particular problem, it is usually only necessary to call a single driver function, but alternatively two or more functions may be called in succession. This is illustrated in the example programs in the function documents.

### 3.2NAG Names and LAPACK Names

As well as the NAG function name (beginning f07), Tables 1 and 8 show the LAPACK routine names in double precision.
The functions may be called either by their NAG short names or by their NAG long names. The NAG long names for a function is simply the LAPACK name (in lower case) prepended by nag_, for example, nag_lapacklin_dpotrf is the long name for f07fdc.
References to Chapter F07 functions in the manual normally include the LAPACK double precision names, for example, f07adc.
The LAPACK routine names follow a simple scheme. Most names have the structure xyyzzz, where the components have the following meanings:
• the initial letter x indicates the data type (real or complex) and precision:  s real, single precision d real, double precision c complex, single precision z complex, double precision
• exceptionally, the mixed precision LAPACK routines described in Section 2.9 replace the initial first letter by a pair of letters, as:  ds double precision function using single precision internally zc double complex function using single precision complex internally
• the letters yy indicate the type of the matrix $A$ (and in some cases its storage scheme):  ge general gb general band po symmetric or Hermitian positive definite pf symmetric or Hermitian positive definite (rectangular full packed (RFP) storage) pp symmetric or Hermitian positive definite (packed storage) pb symmetric or Hermitian positive definite band sy symmetric indefinite sf symmetric indefinite (rectangular full packed (RFP) storage) sp symmetric indefinite (packed storage) he (complex) Hermitian indefinite hf (complex) Hermitian indefinite (rectangular full packed (RFP) storage) hp (complex) Hermitian indefinite (packed storage) gt general tridiagonal pt symmetric or Hermitian positive definite tridiagonal tr triangular tf triangular (rectangular full packed (RFP) storage) tp triangular (packed storage) tb triangular band
• the last two or three letters zz or zzz indicate the computation performed. Examples are:  trf triangular factorization trs solution of linear equations, using the factorization con estimate condition number rfs refine solution and compute error bounds tri compute inverse, using the factorization
Thus the function nag_lapacklin_dgetrf performs a triangular factorization of a real general matrix in double precision; the corresponding function for a complex general matrix is nag_lapacklin_zgetrf.

### 3.3Matrix Storage Schemes

In this chapter the following different storage schemes are used for matrices:
• conventional storage;
• packed storage for symmetric, Hermitian or triangular matrices;
• rectangular full packed (RFP) storage for symmetric, Hermitian or triangular matrices;
• band storage for band matrices.
These storage schemes are compatible with those used in Chapter F16 (especially in the BLAS) and Chapter F08, but different schemes for packed or band storage are used in a few older functions in Chapters F01, F02, F03 and F04.
In the examples below, $*$ indicates an array element which need not be set and is not referenced by the functions. The examples illustrate only the relevant part of the arrays; array arguments may of course have additional rows or columns, according to the usual rules for passing array arguments.

#### 3.3.1Conventional storage

Matrices may be stored column-wise or row-wise as described in Section 3.1.4 in the Introduction to the NAG Library CL Interface: a matrix $A$ is stored in a one-dimensional array a, with matrix element ${a}_{i,j}$ stored column-wise in array element $\mathbf{a}\left[\left(j-1\right)×\mathbf{pda}+i-1\right]$ or row-wise in array element $\mathbf{a}\left[\left(i-1\right)×\mathbf{pda}+j-1\right]$ where pda is the principle dimension of the array (i.e., the stride separating row or column elements of the matrix respectively). Most functions in this chapter contain the order argument which can be set to Nag_ColMajor for column-wise storage or Nag_RowMajor for row-wise storage of matrices. Where groups of functions are intended to be used together, the value of the order argument passed must be consistent throughout.
If a matrix is triangular (upper or lower, as specified by the argument uplo), only the elements of the relevant triangle are stored; the remaining elements of the array need not be set. Such elements are indicated by * or $⌴$ in the examples below.
For example, when $n=3$:
order uplo Triangular matrix $\mathbit{A}$ Storage in array a
Nag_ColMajor Nag_Upper $\left(\begin{array}{lll}{a}_{11}& {a}_{12}& {a}_{13}\\ & {a}_{22}& {a}_{23}\\ & & {a}_{33}\end{array}\right)$ ${a}_{11}⌴⌴{a}_{12}{a}_{22}⌴{a}_{13}{a}_{23}{a}_{33}$
Nag_RowMajor Nag_Upper $\left(\begin{array}{lll}{a}_{11}& {a}_{12}& {a}_{13}\\ & {a}_{22}& {a}_{23}\\ & & {a}_{33}\end{array}\right)$ ${a}_{11}{a}_{12}{a}_{13}⌴{a}_{22}{a}_{23}⌴⌴{a}_{33}$
Nag_ColMajor Nag_Lower $\left(\begin{array}{lll}{a}_{11}& & \\ {a}_{21}& {a}_{22}& \\ {a}_{31}& {a}_{32}& {a}_{33}\end{array}\right)$ ${a}_{11}{a}_{21}{a}_{31}⌴{a}_{22}{a}_{32}⌴⌴{a}_{33}$
Nag_RowMajor Nag_Lower $\left(\begin{array}{lll}{a}_{11}& & \\ {a}_{21}& {a}_{22}& \\ {a}_{31}& {a}_{32}& {a}_{33}\end{array}\right)$ ${a}_{11}⌴⌴{a}_{21}{a}_{22}⌴{a}_{31}{a}_{32}{a}_{33}$
Functions which handle symmetric or Hermitian matrices allow for either the upper or lower triangle of the matrix (as specified by uplo) to be stored in the corresponding elements of the array; the remaining elements of the array need not be set.
For example, when $n=3$:
order uplo Hermitian matrix $\mathbit{A}$ Storage in array a
Nag_ColMajor Nag_Upper $\left(\begin{array}{lll}{a}_{11}& {a}_{12}& {a}_{13}\\ {\overline{a}}_{12}& {a}_{22}& {a}_{23}\\ {\overline{a}}_{13}& {\overline{a}}_{23}& {a}_{33}\end{array}\right)$ ${a}_{11}⌴⌴{a}_{12}{a}_{22}⌴{a}_{13}{a}_{23}{a}_{33}$
Nag_RowMajor Nag_Upper $\left(\begin{array}{lll}{a}_{11}& {a}_{12}& {a}_{13}\\ {\overline{a}}_{12}& {a}_{22}& {a}_{23}\\ {\overline{a}}_{13}& {\overline{a}}_{23}& {a}_{33}\end{array}\right)$ ${a}_{11}{a}_{12}{a}_{13}⌴{a}_{22}{a}_{23}⌴⌴{a}_{33}$
Nag_ColMajor Nag_Lower $\left(\begin{array}{lll}{a}_{11}& {\overline{a}}_{21}& {\overline{a}}_{31}\\ {a}_{21}& {a}_{22}& {\overline{a}}_{32}\\ {a}_{31}& {a}_{32}& {a}_{33}\end{array}\right)$ ${a}_{11}{a}_{21}{a}_{31}⌴{a}_{22}{a}_{32}⌴⌴{a}_{33}$
Nag_RowMajor Nag_Lower $\left(\begin{array}{lll}{a}_{11}& {\overline{a}}_{21}& {\overline{a}}_{31}\\ {a}_{21}& {a}_{22}& {\overline{a}}_{32}\\ {a}_{31}& {a}_{32}& {a}_{33}\end{array}\right)$ ${a}_{11}⌴⌴{a}_{21}{a}_{22}⌴{a}_{31}{a}_{32}{a}_{33}$

#### 3.3.2Packed storage

Symmetric, Hermitian or triangular matrices may be stored more compactly, if the relevant triangle (again as specified by uplo) is packed by columns in a one-dimensional array. In this chapter, as in Chapters F08 and F16, arrays which hold matrices in packed storage, have names ending in P. For a matrix of order $n$, the array must have at least $n\left(n+1\right)/2$ elements. So:
• if $\mathbf{uplo}=\mathrm{Nag_Upper}$, ${a}_{ij}$ is stored in $\mathbf{ap}\left[i-1+j\left(j-1\right)/2\right]$ for $i\le j$;
• if $\mathbf{uplo}=\mathrm{Nag_Lower}$, ${a}_{ij}$ is stored in $\mathbf{ap}\left[i-1+\left(2n-j\right)\left(j-1\right)/2\right]$ for $j\le i$.
For example:
uplo Triangle of matrix $\mathbit{A}$ Packed storage in array ap
Nag_Upper $\left(\begin{array}{llll}{a}_{11}& {a}_{12}& {a}_{13}& {a}_{14}\\ & {a}_{22}& {a}_{23}& {a}_{24}\\ & & {a}_{33}& {a}_{34}\\ & & & {a}_{44}\end{array}\right)$ ${a}_{11}\underbrace{{a}_{12}{a}_{22}}\underbrace{{a}_{13}{a}_{23}{a}_{33}}\underbrace{{a}_{14}{a}_{24}{a}_{34}{a}_{44}}$
Nag_Lower $\left(\begin{array}{llll}{a}_{11}& & & \\ {a}_{21}& {a}_{22}& & \\ {a}_{31}& {a}_{32}& {a}_{33}& \\ {a}_{41}& {a}_{42}& {a}_{43}& {a}_{44}\end{array}\right)$ $\underbrace{{a}_{11}{a}_{21}{a}_{31}{a}_{41}}\underbrace{{a}_{22}{a}_{32}{a}_{42}}\underbrace{{a}_{33}{a}_{43}}{a}_{44}$
Note that for real symmetric matrices, packing the upper triangle by columns is equivalent to packing the lower triangle by rows; packing the lower triangle by columns is equivalent to packing the upper triangle by rows. (For complex Hermitian matrices, the only difference is that the off-diagonal elements are conjugated.)

#### 3.3.3Rectangular Full Packed (RFP) Storage

The rectangular full packed (RFP) storage format offers the same savings in storage as the packed storage format (described in Section 3.4.2), but is likely to be much more efficient in general since the block structure of the matrix is maintained. This structure can be exploited using block partition algorithms (see Section 2.8) in a similar way to matrices that use conventional storage.
Figure 1 gives a graphical representation of the key idea of RFP for the particular case of a lower triangular matrix of even dimensions. In all cases the original triangular matrix of stored elements is separated into a trapezoidal part and a triangular part. The number of columns in these two parts is equal when the dimension of the matrix is even, $n=2k$, while the trapezoidal part has $k+1$ columns when $n=2k+1$. The smaller part is then transposed and fitted onto the trapezoidal part forming a rectangle. The rectangle has dimensions $2k+1$ and $q$, where $q=k$ when $n$ is even and $q=k+1$ when $n$ is odd.
For functions using RFP there is the option of storing the rectangle as described above ($\mathbf{transr}=\mathrm{Nag_RFP_Normal}$) or its transpose ($\mathbf{transr}=\mathrm{Nag_RFP_Trans}$, for real a) or its conjugate transpose ($\mathbf{transr}=\mathrm{Nag_RFP_ConjTrans}$, for complex a). It is useful to note that the storage ordering for Nag_RowMajor is the same as that for Nag_ColMajor with the value of $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ switched to $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$ or vice versa.
As an example, we first consider RFP for the case $n=2k$ with $k=3$.
If $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$, then ar holds a as follows:
• For $\mathbf{uplo}=\mathrm{Nag_Upper}$ the upper trapezoid $\mathbf{AR}\left(1:6,1:3\right)$ consists of the last three columns of a upper. The lower triangle $\mathbf{AR}\left(5:7,1:3\right)$ consists of the transpose of the first three columns of a upper.
• For $\mathbf{uplo}=\mathrm{Nag_Lower}$ the lower trapezoid $\mathbf{AR}\left(2:7,1:3\right)$ consists of the first three columns of a lower. The upper triangle $\mathbf{AR}\left(1:3,1:3\right)$ consists of the transpose of the last three columns of a lower.
If $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$, then ar in both uplo cases is just the transpose of ar as defined when $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$.
uplo Triangle of matrix $\mathbit{A}$ Rectangular Full Packed matrix $\mathbit{AR}$
$\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$
Nag_Upper $\left(\begin{array}{llllll}\mathbf{00}& \mathbf{01}& \mathbf{02}& 03& 04& 05\\ & \mathbf{11}& \mathbf{12}& 13& 14& 15\\ & & \mathbf{22}& 23& 24& 25\\ & & & 33& 34& 35\\ & & & & 44& 45\\ & & & & & 55\end{array}\right)$ $\begin{array}{ccc}03& 04& 05\\ 13& 14& 15\\ 23& 24& 25\\ 33& 34& 35\\ \mathbf{00}& 44& 45\\ \mathbf{01}& \mathbf{11}& 55\\ \mathbf{02}& \mathbf{12}& \mathbf{22}\end{array}$ $\begin{array}{ccccccc}03& 13& 23& 33& \mathbf{00}& \mathbf{01}& \mathbf{02}\\ 04& 14& 24& 34& 44& \mathbf{11}& \mathbf{12}\\ 05& 15& 25& 35& 45& 55& \mathbf{22}\end{array}$
Nag_Lower $\left(\begin{array}{l}00\\ 10& 11\\ 20& 21& 22\\ 30& 31& 32& \mathbf{33}\\ 40& 41& 42& \mathbf{43}& \mathbf{44}\\ 50& 51& 52& \mathbf{53}& \mathbf{54}& \mathbf{55}\end{array}\right)$ $\begin{array}{ccc}\mathbf{33}& \mathbf{43}& \mathbf{53}\\ 00& \mathbf{44}& \mathbf{54}\\ 10& 11& \mathbf{55}\\ 20& 21& 22\\ 30& 31& 32\\ 40& 41& 42\\ 50& 51& 52\end{array}$ $\begin{array}{ccccccc}\mathbf{33}& 00& 10& 20& 30& 40& 50\\ \mathbf{43}& \mathbf{44}& 11& 21& 31& 41& 51\\ \mathbf{53}& \mathbf{54}& \mathbf{55}& 22& 32& 42& 52\end{array}$
Now we consider RFP for the case $n=2k+1$ and $k=2$.
If $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$. ar holds a as follows:
• if $\mathbf{uplo}=\mathrm{Nag_Upper}$ the upper trapezoid $\mathbf{AR}\left(1:5,1:3\right)$ consists of the last three columns of a upper. The lower triangle $\mathbf{AR}\left(4:5,1:2\right)$ consists of the transpose of the first two columns of a upper;
• if $\mathbf{uplo}=\mathrm{Nag_Lower}$ the lower trapezoid $\mathbf{AR}\left(1:5,1:3\right)$ consists of the first three columns of a lower. The upper triangle $\mathbf{AR}\left(1:2,2:3\right)$ consists of the transpose of the last two columns of a lower.
If $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$. ar in both uplo cases is just the transpose of ar as defined when $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$.
uplo Triangle of matrix $\mathbit{A}$ Rectangular Full Packed matrix $\mathbit{AR}$
$\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$
Nag_Upper $\left(\begin{array}{lllll}\mathbf{00}& \mathbf{01}& 02& 03& 04\\ & \mathbf{11}& 12& 13& 14\\ & & 22& 23& 24\\ & & & 33& 34\\ & & & & 44\end{array}\right)$ $\begin{array}{ccc}02& 03& 04\\ 12& 13& 14\\ 22& 23& 24\\ \mathbf{00}& 33& 34\\ \mathbf{01}& \mathbf{11}& 44\end{array}$ $\begin{array}{ccccc}02& 12& 22& \mathbf{00}& \mathbf{01}\\ 03& 13& 23& 33& \mathbf{11}\\ 04& 14& 24& 34& 44\end{array}$
Nag_Lower $\left(\begin{array}{l}00\\ 10& 11\\ 20& 21& 22\\ 30& 31& 32& \mathbf{33}\\ 40& 41& 42& \mathbf{43}& \mathbf{44}\end{array}\right)$ $\begin{array}{ccc}00& \mathbf{33}& \mathbf{43}\\ 10& 11& \mathbf{44}\\ 20& 21& 22\\ 30& 31& 32\\ 40& 41& 42\end{array}$ $\begin{array}{cccccc}00& 10& 20& 30& 40& 50\\ \mathbf{33}& 11& 21& 31& 41& 51\\ \mathbf{43}& \mathbf{44}& 22& 32& 42& 52\end{array}$
Explicitly, in the real matrix case, ar is a one-dimensional array of length $n\left(n+1\right)/2$ and, when Nag_ColMajor, contains the elements of a as follows:
for $\mathbf{uplo}=\mathrm{Nag_Upper}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(\mathit{i}-1\right)+\mathit{j}+k\right]$, for $1\le \mathit{j}\le k$ and $1\le i\le \mathit{j}$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(\mathit{j}-k-1\right)+i-1\right]$, for $k and $1\le i\le j$;
for $\mathbf{uplo}=\mathrm{Nag_Upper}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(j+k\right)+i-1\right]$, for $1\le j\le k$ and $1\le i\le j$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(i-1\right)+\mathit{j}-k-1\right]$, for $k and $1\le i\le j$;
for $\mathbf{uplo}=\mathrm{Nag_Lower}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(j-1\right)+i+k-q\right]$, for $1\le j\le q$ and $j\le i\le n$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(i-k-1\right)+j-q-1\right]$, for $q and $j\le i\le n$;
for $\mathbf{uplo}=\mathrm{Nag_Lower}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(i+k-q\right)+j-1\right]$, for $1\le j\le q$ and $1\le i\le n$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(j-1-q\right)+i-k-1\right]$, for $q and $1\le i\le n$.
When Nag_RowMajor the above storage formulae can be used by looking up the opposite case for transr, i.e., when $\mathbf{transr}=\mathrm{Nag_RFP_Trans}$ look up the storage order above for the cases when $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ and vice versa.
In the case of complex matrices, the assumption is that the full matrix, if it existed, would be Hermitian. Thus, when $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$, the triangular portion of a that is, in the real case, transposed into the notional $\left(2k+1\right)×q$ RFP matrix is also conjugated. When $\mathbf{transr}=\mathrm{Nag_RFP_ConjTrans}$ the notional $q×\left(2k+1\right)$ RFP matrix is the conjugate transpose of the corresponding $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ RFP matrix. Explicitly, for complex a, the array ar contains the elements (or conjugated elements) of a as follows:
for $\mathbf{uplo}=\mathrm{Nag_Upper}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$,
${\overline{a}}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(\mathit{i}-1\right)+\mathit{j}+k\right]$, for $1\le \mathit{j}\le k$ and $1\le i\le \mathit{j}$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(\mathit{j}-k-1\right)+i-1\right]$, for $k and $1\le i\le j$;
for $\mathbf{uplo}=\mathrm{Nag_Upper}$ and $\mathbf{transr}=\mathrm{Nag_RFP_ConjTrans}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(j+k\right)+i-1\right]$, for $1\le j\le k$ and $1\le i\le j$, and
${\overline{a}}_{ij}$ is stored in $\mathbf{ar}\left[q\left(i-1\right)+\mathit{j}-k-1\right]$, for $k and $1\le i\le j$;
for $\mathbf{uplo}=\mathrm{Nag_Lower}$ and $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$,
${a}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(j-1\right)+i+k-q\right]$, for $1\le j\le q$ and $j\le i\le n$, and
${\overline{a}}_{ij}$ is stored in $\mathbf{ar}\left[\left(2k+1\right)\left(i-k-1\right)+j-q-1\right]$, for $q and $j\le i\le n$;
for $\mathbf{uplo}=\mathrm{Nag_Lower}$ and $\mathbf{transr}=\mathrm{Nag_RFP_ConjTrans}$,
${\overline{a}}_{ij}$ is stored in $\mathbf{ar}\left[q\left(i+k-q\right)+j-1\right]$, for $1\le j\le q$ and $1\le i\le n$, and
${a}_{ij}$ is stored in $\mathbf{ar}\left[q\left(j-1-q\right)+i-k-1\right]$, for $q and $1\le i\le n$.
When Nag_RowMajor the above storage formulae can be used by looking up the opposite case for transr, i.e., when $\mathbf{transr}=\mathrm{Nag_RFP_ConjTrans}$ look up the storage order above for the cases when $\mathbf{transr}=\mathrm{Nag_RFP_Normal}$ and vice versa.

#### 3.3.4Band storage

A band matrix with ${k}_{l}$ subdiagonals and ${k}_{u}$ superdiagonals may be stored compactly in a notional two-dimensional array with ${k}_{l}+{k}_{u}+1$ rows and $n$ columns if stored column-wise or $n$ rows and ${k}_{l}+{k}_{u}+1$ columns if stored row-wise. In column-major order, elements of a column of the matrix are stored contiguously in the array, and elements of the diagonals of the matrix are stored with constant stride (i.e., in a row of the two-dimensional array). In row-major order, elements of a row of the matrix are stored contiguously in the array, and elements of a diagonal of the matrix are stored with constant stride (i.e., in a column of the two-dimensional array). These storage schemes should only be used in practice if ${k}_{l}$, ${k}_{u}\ll n$, although the functions in Chapters F07 and F08 work correctly for all values of ${k}_{l}$ and ${k}_{u}$. In Chapters F07 and F08 arrays which hold matrices in band storage have names ending in $\mathrm{B}$.
To be precise, elements of matrix elements ${a}_{ij}$ are stored as follows:
• if $\mathbf{order}=\mathrm{Nag_ColMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(j-1\right)×\mathbf{pdab}+{k}_{u}+i-j\right]$;
• if $\mathbf{order}=\mathrm{Nag_RowMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(i-1\right)×\mathbf{pdab}+{k}_{l}+j-i\right]$,
where $\mathbf{pdab}\ge {k}_{l}+{k}_{u}+1$ is the stride between diagonal elements and where $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,i-{k}_{l}\right)\le j\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,i+{k}_{u}\right)$.
For example, when $n=5$, ${k}_{l}=2$ and ${k}_{u}=1$:
Band matrix $\mathbit{A}$ Band storage in array ab
$\mathbf{order}=\mathrm{Nag_ColMajor}$ $\mathbf{order}=\mathrm{Nag_RowMajor}$
$\begin{array}{lllll}{a}_{11}& {a}_{12}& & & \\ {a}_{21}& {a}_{22}& {a}_{23}& & \\ {a}_{31}& {a}_{32}& {a}_{33}& {a}_{34}& \\ & {a}_{42}& {a}_{43}& {a}_{44}& {a}_{45}\\ & & {a}_{53}& {a}_{54}& {a}_{55}\end{array}$ $\begin{array}{lllll}\text{*}& {a}_{12}& {a}_{23}& {a}_{34}& {a}_{45}\\ {a}_{11}& {a}_{22}& {a}_{33}& {a}_{44}& {a}_{55}\\ {a}_{21}& {a}_{32}& {a}_{43}& {a}_{54}& \text{*}\\ {a}_{31}& {a}_{42}& {a}_{53}& \text{*}& \text{*}\end{array}$ $\begin{array}{llll}\text{*}& \text{*}& {a}_{11}& {a}_{12}\\ \text{*}& {a}_{21}& {a}_{22}& {a}_{23}\\ {a}_{31}& {a}_{32}& {a}_{33}& {a}_{34}\\ {a}_{42}& {a}_{43}& {a}_{44}& {a}_{45}\\ {a}_{53}& {a}_{54}& {a}_{55}& \text{*}\end{array}$
The elements marked $*$ in the upper left and lower right corners of the array ab need not be set, and are not referenced by the functions. In this example, if $\mathbf{order}=\mathrm{Nag_ColMajor}$ and ldab takes the minimum value of $4$, then $\mathbf{ab}\left[0\right]$ need not be set, $\mathbf{ab}\left[1\right]={a}_{11},\mathbf{ab}\left[2\right]={a}_{21},\dots ,\mathbf{ab}\left[17\right]={a}_{55}$. On the other hand, if $\mathbf{order}=\mathrm{Nag_RowMajor}$ ($\mathbf{ldab}=4$), then $\mathbf{ab}\left[0\right]$ and $\mathbf{ab}\left[1\right]$ need not be set, $\mathbf{ab}\left[2\right]={a}_{11},\mathbf{ab}\left[3\right]={a}_{12},\dots ,\mathbf{ab}\left[18\right]={a}_{55}$.
Note:  when a general band matrix is supplied for $LU$ factorization, space must be allowed to store an additional ${k}_{l}$ superdiagonals, generated by fill-in as a result of row interchanges. This means that the matrix is stored according to the above scheme, but with ${k}_{l}+{k}_{u}$ superdiagonals; it also means that the principal dimension has the constraint $\mathbf{ldab}\ge 2{k}_{l}+{k}_{u}+1$.
Triangular band matrices are stored in the same format, with either ${k}_{l}=0$ if upper triangular, or ${k}_{u}=0$ if lower triangular.
For symmetric or Hermitian band matrices with $k$ subdiagonals or superdiagonals, only the upper or lower triangle (as specified by uplo) need be stored:
• if $\mathbf{uplo}=\mathrm{Nag_Upper}$ then
• if $\mathbf{order}=\mathrm{Nag_ColMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(j-1\right)×\mathbf{pdab}+k+i-j\right]$;
• if $\mathbf{order}=\mathrm{Nag_RowMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(i-1\right)×\mathbf{pdab}+j-i\right]$.
for $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j-k\right)\le i\le j$;
• if $\mathbf{uplo}=\mathrm{Nag_Lower}$ then
• if $\mathbf{order}=\mathrm{Nag_ColMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(j-1\right)×\mathbf{pdab}+i-j\right]$;
• if $\mathbf{order}=\mathrm{Nag_RowMajor}$, ${a}_{ij}$ is stored in $\mathbf{ab}\left[\left(i-1\right)×\mathbf{pdab}+k+j-i\right]$.
for $j\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+k\right)$,
where $\mathbf{pdab}\ge k+1$ is the stride separating diagonal matrix elements in the array ab.
For example, when $n=5$ and $k=2$:
uplo Hermitian band matrix $\mathbit{A}$ Band storage in array ab
$\mathbf{order}=\mathrm{Nag_ColMajor}$ $\mathbf{order}=\mathrm{Nag_RowMajor}$
Nag_Upper $\left(\begin{array}{lllll}{a}_{11}& {a}_{12}& {a}_{13}& & \\ {\overline{a}}_{12}& {a}_{22}& {a}_{23}& {a}_{24}& \\ {\overline{a}}_{13}& {\overline{a}}_{23}& {a}_{33}& {a}_{34}& {a}_{35}\\ & {\overline{a}}_{24}& {\overline{a}}_{34}& {a}_{44}& {a}_{45}\\ & & {\overline{a}}_{35}& {\overline{a}}_{45}& {a}_{55}\end{array}\right)$ $\begin{array}{lllll}\text{*}& \text{*}& {a}_{13}& {a}_{24}& {a}_{35}\\ \text{*}& {a}_{12}& {a}_{23}& {a}_{34}& {a}_{45}\\ {a}_{11}& {a}_{22}& {a}_{33}& {a}_{44}& {a}_{55}\end{array}$ $\begin{array}{lll}{a}_{11}& {a}_{12}& {a}_{13}\\ {a}_{22}& {a}_{23}& {a}_{24}\\ {a}_{33}& {a}_{34}& {a}_{35}\\ {a}_{44}& {a}_{45}& \text{*}\\ {a}_{55}& \text{*}& \text{*}\end{array}$
Nag_Lower $\left(\begin{array}{lllll}{a}_{11}& {\overline{a}}_{21}& {\overline{a}}_{31}& & \\ {a}_{21}& {a}_{22}& {\overline{a}}_{32}& {\overline{a}}_{42}& \\ {a}_{31}& {a}_{32}& {a}_{33}& {\overline{a}}_{43}& {\overline{a}}_{53}\\ & {a}_{42}& {a}_{43}& {a}_{44}& {\overline{a}}_{54}\\ & & {a}_{53}& {a}_{54}& {a}_{55}\end{array}\right)$ $\begin{array}{lllll}{a}_{11}& {a}_{22}& {a}_{33}& {a}_{44}& {a}_{55}\\ {a}_{21}& {a}_{32}& {a}_{43}& {a}_{54}& \text{*}\\ {a}_{31}& {a}_{42}& {a}_{53}& \text{*}& \text{*}\end{array}$ $\begin{array}{lll}\text{*}& \text{*}& {a}_{11}\\ \text{*}& {a}_{21}& {a}_{22}\\ {a}_{31}& {a}_{32}& {a}_{33}\\ {a}_{42}& {a}_{43}& {a}_{44}\\ {a}_{53}& {a}_{54}& {a}_{55}\end{array}$
Note that different storage schemes for band matrices are used by some functions in Chapters F01, F02, F03 and F04. In the above example, if $\mathbf{order}=\mathrm{Nag_ColMajor}$ and $\mathbf{pdab}=3$, then for $\mathbf{uplo}=\mathrm{Nag_Upper}$, $\mathbf{ab}\left[2\right]={a}_{11},\mathbf{ab}\left[4\right]={a}_{12},\dots ,\mathbf{ab}\left[14\right]={a}_{55}$; while for $\mathbf{uplo}=\mathrm{Nag_Lower}$, $\mathbf{ab}\left[0\right]={a}_{11},\mathbf{ab}\left[1\right]={a}_{21},\dots ,\mathbf{ab}\left[12\right]={a}_{55}$. If $\mathbf{order}=\mathrm{Nag_RowMajor}$ ($\mathbf{pdab}=3$), then for $\mathbf{uplo}=\mathrm{Nag_Upper}$, $\mathbf{ab}\left[0\right]={a}_{11},\mathbf{ab}\left[1\right]={a}_{12},\dots ,\mathbf{ab}\left[12\right]={a}_{55}$; while for $\mathbf{uplo}=\mathrm{Nag_Lower}$, $\mathbf{ab}\left[2\right]={a}_{11},\mathbf{ab}\left[4\right]={a}_{21},\dots ,\mathbf{ab}\left[14\right]={a}_{55}$.

#### 3.3.5Unit triangular matrices

Some functions in this chapter have an option to handle unit triangular matrices (that is, triangular matrices with diagonal elements $\text{}=1$). This option is specified by an argument diag. If $\mathbf{diag}=\mathrm{Nag_UnitDiag}$ (Unit triangular), the diagonal elements of the matrix need not be stored, and the corresponding array elements are not referenced by the functions. The storage scheme for the rest of the matrix (whether conventional, packed or band) remains unchanged.

#### 3.3.6Real diagonal elements of complex matrices

Complex Hermitian matrices have diagonal elements that are by definition purely real. In addition, complex triangular matrices which arise in Cholesky factorization are defined by the algorithm to have real diagonal elements.
If such matrices are supplied as input to functions in Chapters F07 and F08, the imaginary parts of the diagonal elements are not referenced, but are assumed to be zero. If such matrices are returned as output by the functions, the computed imaginary parts are explicitly set to zero.

### 3.4Parameter Conventions

#### 3.4.1Option arguments

In addition to the order argument of type Nag_OrderType, most functions in this Chapter have one or more option arguments of various types; only options of the correct type may be supplied.
For example,
`nag_lapacklin_dpotrf(Nag_RowMajor,Nag_Upper,...)`

#### 3.4.2Problem dimensions

It is permissible for the problem dimensions (for example, m in f07adc, n or nrhs in f07aec) to be passed as zero, in which case the computation (or part of it) is skipped. Negative dimensions are regarded as an error.

### 3.5Tables of Driver and Computational Functions

#### 3.5.1Real matrices

Each entry gives:
• the NAG function short name
• the LAPACK routine name from which the NAG function long name is derived by prepending nag_.
Table 1
Functions for real general matrices
Type of matrix and storage scheme
Operation general general band general tridiagonal
driver f07aac f07bac f07cac
expert driver f07abc f07bbc f07cbc
mixed precision driver f07acc
solve f07aec f07bec f07cec
scaling factors f07afc f07bfc
condition number f07agc f07bgc f07cgc
error estimate f07ahc f07bhc f07chc
invert f07ajc
Table 2
Functions for real symmetric positive definite matrices
Type of matrix and storage scheme
Operation symmetric positive definite symmetric positive definite (packed storage) symmetric positive definite band symmetric positive definite tridiagonal
driver f07fac f07gac f07hac f07jac
expert driver f07fbc f07gbc f07hbc f07jbc
factorize f07fdc f07gdc f07hdc f07jdc
solve f07fec f07gec f07hec f07jec
scaling factors f07ffc f07gfc f07hfc
condition number f07fgc f07ggc f07hgc f07jgc
error estimate f07fhc f07ghc f07hhc f07jhc
invert f07fjc f07gjc
Table 3
Functions for real symmetric indefinite matrices
Type of matrix and storage scheme
Operation symmetric indefinite symmetric indefinite (packed storage)
driver f07mac f07pac
expert driver f07mbc f07pbc
factorize f07mdc f07pdc
solve f07mec f07pec
condition number f07mgc f07pgc
error estimate f07mhc f07phc
invert f07mjc f07pjc
Table 4
Functions for real triangular matrices
Type of matrix and storage scheme
Operation triangular triangular (packed storage) triangular band
solve f07tec f07uec f07vec
condition number f07tgc f07ugc f07vgc
error estimate f07thc f07uhc f07vhc
invert f07tjc f07ujc

#### 3.5.2Complex matrices

Each entry gives:
• the NAG function short name
• the LAPACK routine name from which the NAG function long name is derived by prepending nag_.
Table 5
Functions for complex general matrices
Type of matrix and storage scheme
Operation general general band general tridiagonal
driver f07anc f07bnc f07cnc
expert driver f07apc f07bpc f07cpc
mixed precision driver f07aqc
factorize f07arc f07brc f07crc
solve f07asc f07bsc f07csc
scaling factors f07atc f07btc
condition number f07auc f07buc f07cuc
error estimate f07avc f07bvc f07cvc
invert f07awc
Table 6
Functions for complex Hermitian positive definite matrices
Type of matrix and storage scheme
Operation Hermitian positive definite Hermitian positive definite (packed storage) Hermitian positive definite band Hermitian positive definite tridiagonal
driver f07fnc f07gnc f07hnc f07jnc
expert driver f07fpc f07gpc f07hpc f07jpc
factorize f07frc f07grc f07hrc f07jrc
solve f07fsc f07gsc f07hsc f07jsc
scaling factors f07ftc f07gtc
condition number f07fuc f07guc f07huc f07juc
error estimate f07fvc f07gvc f07hvc f07jvc
invert f07fwc f07gwc
Table 7
Functions for complex Hermitian and symmetric indefinite matrices
Type of matrix and storage scheme
Operation Hermitian indefinite symmetric indefinite (packed storage) Hermitian indefinite band symmetric indefinite tridiagonal
driver f07mnc f07nnc f07pnc f07qnc
expert driver f07mpc f07npc f07ppc f07qpc
factorize f07mrc f07nrc f07prc f07qrc
solve f07msc f07nsc f07psc f07qsc
condition number f07muc f07nuc f07puc f07quc
error estimate f07mvc f07nvc f07pvc f07qvc
invert f07mwc f07nwc f07pwc f07qwc
Table 8
Functions for complex triangular matrices
Type of matrix and storage scheme
Operation triangular triangular (packed storage) triangular band
solve f07tsc f07usc f07vsc
condition number f07tuc f07uuc f07vuc
error estimate f07tvc f07uvc f07vvc
invert f07twc f07uwc

## 4Functionality Index

 Apply iterative refinement to the solution and compute error estimates,
 after factorizing the matrix of coefficients,
 complex band matrix f07bvc
 complex Hermitian indefinite matrix f07mvc
 complex Hermitian indefinite matrix, packed storage f07pvc
 complex Hermitian positive definite band matrix f07hvc
 complex Hermitian positive definite matrix f07fvc
 complex Hermitian positive definite matrix, packed storage f07gvc
 complex Hermitian positive definite tridiagonal matrix f07jvc
 complex matrix f07avc
 complex symmetric indefinite matrix f07nvc
 complex symmetric indefinite matrix, packed storage f07qvc
 complex tridiagonal matrix f07cvc
 real band matrix f07bhc
 real matrix f07ahc
 real symmetric indefinite matrix f07mhc
 real symmetric indefinite matrix, packed storage f07phc
 real symmetric positive definite band matrix f07hhc
 real symmetric positive definite matrix f07fhc
 real symmetric positive definite matrix, packed storage f07ghc
 real symmetric positive definite tridiagonal matrix f07jhc
 real tridiagonal matrix f07chc
 Compute error estimates,
 complex triangular band matrix f07vvc
 complex triangular matrix f07tvc
 complex triangular matrix, packed storage f07uvc
 real triangular band matrix f07vhc
 real triangular matrix f07thc
 real triangular matrix, packed storage f07uhc
 Compute row and column scalings,
 complex band matrix f07btc
 complex Hermitian positive definite band matrix f07htc
 complex Hermitian positive definite matrix f07ftc
 complex Hermitian positive definite matrix, packed storage f07gtc
 complex matrix f07atc
 real band matrix f07bfc
 real matrix f07afc
 real symmetric positive definite band matrix f07hfc
 real symmetric positive definite matrix f07ffc
 real symmetric positive definite matrix, packed storage f07gfc
 Condition number estimation,
 after factorizing the matrix of coefficients,
 complex band matrix f07buc
 complex Hermitian indefinite matrix f07muc
 complex Hermitian indefinite matrix, packed storage f07puc
 complex Hermitian positive definite band matrix f07huc
 complex Hermitian positive definite matrix f07fuc
 complex Hermitian positive definite matrix, packed storage f07guc
 complex Hermitian positive definite tridiagonal matrix f07juc
 complex matrix f07auc
 complex symmetric indefinite matrix f07nuc
 complex symmetric indefinite matrix, packed storage f07quc
 complex tridiagonal matrix f07cuc
 real band matrix f07bgc
 real matrix f07agc
 real symmetric indefinite matrix f07mgc
 real symmetric indefinite matrix, packed storage f07pgc
 real symmetric positive definite band matrix f07hgc
 real symmetric positive definite matrix f07fgc
 real symmetric positive definite matrix, packed storage f07ggc
 real symmetric positive definite tridiagonal matrix f07jgc
 real tridiagonal matrix f07cgc
 complex triangular band matrix f07vuc
 complex triangular matrix f07tuc
 complex triangular matrix, packed storage f07uuc
 real triangular band matrix f07vgc
 real triangular matrix f07tgc
 real triangular matrix, packed storage f07ugc
 $LD{L}^{\mathrm{T}}$ factorization,
 complex Hermitian positive definite tridiagonal matrix f07jrc
 real symmetric positive definite tridiagonal matrix f07jdc
 $L{L}^{\mathrm{T}}$ or ${U}^{\mathrm{T}}U$ factorization,
 complex Hermitian positive definite band matrix f07hrc
 complex Hermitian positive definite matrix f07frc
 complex Hermitian positive definite matrix, packed storage f07grc
 complex Hermitian positive definite matrix, RFP storage f07wrc
 complex Hermitian positive semidefinite matrix f07krc
 real symmetric positive definite band matrix f07hdc
 real symmetric positive definite matrix f07fdc
 real symmetric positive definite matrix, packed storage f07gdc
 real symmetric positive definite matrix, RFP storage f07wdc
 real symmetric positive semidefinite matrix f07kdc
 $LU$ factorization,
 complex band matrix f07brc
 complex matrix f07arc
 complex tridiagonal matrix f07crc
 real band matrix f07bdc
 real tridiagonal matrix f07cdc
 Matrix inversion,
 after factorizing the matrix of coefficients,
 complex Hermitian indefinite matrix f07mwc
 complex Hermitian indefinite matrix, packed storage f07pwc
 complex Hermitian positive definite matrix f07fwc
 complex Hermitian positive definite matrix, packed storage f07gwc
 complex Hermitian positive definite matrix, RFP storage f07wwc
 complex matrix f07awc
 complex symmetric indefinite matrix f07nwc
 complex symmetric indefinite matrix, packed storage f07qwc
 real matrix f07ajc
 real symmetric indefinite matrix f07mjc
 real symmetric indefinite matrix, packed storage f07pjc
 real symmetric positive definite matrix f07fjc
 real symmetric positive definite matrix, packed storage f07gjc
 real symmetric positive definite matrix, RFP storage f07wjc
 complex triangular matrix f07twc
 complex triangular matrix, packed storage f07uwc
 complex triangular matrix, RFP storage,
 expert driver f07wxc
 real triangular matrix f07tjc
 real triangular matrix, packed storage f07ujc
 real triangular matrix, RFP storage,
 expert driver f07wkc
 $PLD{L}^{\mathrm{T}}{P}^{\mathrm{T}}$ or $PUD{U}^{\mathrm{T}}{P}^{\mathrm{T}}$ factorization,
 complex Hermitian indefinite matrix f07mrc
 complex Hermitian indefinite matrix, packed storage f07prc
 complex symmetric indefinite matrix f07nrc
 complex symmetric indefinite matrix, packed storage f07qrc
 real symmetric indefinite matrix f07mdc
 real symmetric indefinite matrix, packed storage f07pdc
 Solution of simultaneous linear equations,
 after factorizing the matrix of coefficients,
 complex band matrix f07bsc
 complex Hermitian indefinite matrix f07msc
 complex Hermitian indefinite matrix, packed storage f07psc
 complex Hermitian positive definite band matrix f07hsc
 complex Hermitian positive definite matrix f07fsc
 complex Hermitian positive definite matrix, packed storage f07gsc
 complex Hermitian positive definite matrix, RFP storage f07wsc
 complex Hermitian positive definite tridiagonal  matrix f07jsc
 complex matrix f07asc
 complex symmetric indefinite matrix f07nsc
 complex symmetric indefinite matrix, packed storage f07qsc
 complex tridiagonal matrix f07csc
 real band matrix f07bec
 real matrix f07aec
 real symmetric indefinite matrix f07mec
 real symmetric indefinite matrix, packed storage f07pec
 real symmetric positive definite band matrix f07hec
 real symmetric positive definite matrix f07fec
 real symmetric positive definite matrix, packed storage f07gec
 real symmetric positive definite matrix, RFP storage f07wec
 real symmetric positive definite tridiagonal matrix f07jec
 real tridiagonal matrix f07cec
 expert drivers (with condition and error estimation):
 complex band matrix f07bpc
 complex Hermitian indefinite matrix f07mpc
 complex Hermitian indefinite matrix, packed storage f07ppc
 complex Hermitian positive definite band matrix f07hpc
 complex Hermitian positive definite matrix f07fpc
 complex Hermitian positive definite matrix, packed storage f07gpc
 complex Hermitian positive definite tridiagonal matrix f07jpc
 complex matrix f07apc
 complex symmetric indefinite matrix f07npc
 complex symmetric indefinite matrix, packed storage f07qpc
 complex tridiagonal matrix f07cpc
 real band matrix f07bbc
 real matrix f07abc
 real symmetric indefinite matrix f07mbc
 real symmetric indefinite matrix, packed storage f07pbc
 real symmetric positive definite band matrix f07hbc
 real symmetric positive definite matrix f07fbc
 real symmetric positive definite matrix, packed storage f07gbc
 real symmetric positive definite tridiagonal matrix f07jbc
 real tridiagonal matrix f07cbc
 simple drivers,
 complex band matrix f07bnc
 complex Hermitian indefinite matrix f07mnc
 complex Hermitian indefinite matrix, packed storage f07pnc
 complex Hermitian positive definite band matrix f07hnc
 complex Hermitian positive definite matrix f07fnc
 complex Hermitian positive definite matrix, packed storage f07gnc
 complex Hermitian positive definite tridiagonal matrix f07jnc
 complex matrix f07anc
 complex matrix, using mixed precision f07aqc
 complex symmetric indefinite matrix f07nnc
 complex symmetric indefinite matrix, packed storage f07qnc
 complex triangular band matrix f07vsc
 complex triangular matrix f07tsc
 complex triangular matrix, packed storage f07usc
 complex tridiagonal matrix f07cnc
 real band matrix f07bac
 real matrix f07aac
 real matrix, using mixed precision f07acc
 real symmetric indefinite matrix f07mac
 real symmetric indefinite matrix, packed storage f07pac
 real symmetric positive definite band matrix f07hac
 real symmetric positive definite matrix f07fac
 real symmetric positive definite matrix, packed storage f07gac
 real symmetric positive definite tridiagonal matrix f07jac
 real triangular band matrix f07vec
 real triangular matrix f07tec
 real triangular matrix, packed storage f07uec
 real tridiagonal matrix f07cac

None.

None.

## 7References

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
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Higham N J (1988) Algorithm 674: Fortran codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation ACM Trans. Math. Software 14 381–396
Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford