This chapter is concerned with basic linear algebra methods which perform elementary algebraic operations involving scalars, vectors and matrices. It includes methods which conform to the specifications of the BLAS (Basic Linear Algebra Subprograms).

# Syntax

C# |
---|

public static class F06 |

Visual Basic |
---|

Public NotInheritable Class F06 |

Visual C++ |
---|

public ref class F06 abstract sealed |

F# |
---|

[<AbstractClassAttribute>] [<SealedAttribute>] type F06 = class end |

# Background to the Problems

A number of the methods in this chapter meet the specification of the
Basic Linear Algebra Subprograms (BLAS) as described in
Lawson et al. (1979), Dodson et al. (1991), Dongarra et al. (1988) and Dongarra et al. (1990). The first reference describes a set of methods concerned with operations on scalars and vectors: these will be referred to here as the Level-0 and the Level-1
BLAS; the second reference describes a set of methods concerned with operations on sparse vectors: these will be referred to here as the Level-1
Sparse BLAS; the third reference describes a set of methods concerned with matrix-vector operations: these will be referred to here as the Level-2 BLAS;
and the fourth reference describes a set of methods concerned with matrix-matrix operations: these will be referred to here as the Level-3 BLAS.

More generally we refer to the scalar methods in the chapter as
Level-0 methods, to the vector methods as Level-1 methods, to the matrix-vector and matrix methods as Level-2 methods, and to the matrix-matrix methods as Level-3 methods. The terminology reflects the number of operations involved. For example, a Level-2 method
involves $\mathit{O}\left({n}^{2}\right)$ operations for an $n\times n$ matrix.

# The Use of BLAS Names

Many of the methods in other chapters of the Library call the methods in this chapter, and in particular a number of the BLAS are called. These methods are usually called by the BLAS name and so, for correct operation of the Library, it is essential that you do not attempt to link your own versions of these methods. If you are in any doubt about how to avoid this, please consult your computer centre or the NAG
Response Centre.

The BLAS names are used in order to make use of efficient implementations of the methods when these exist. Such implementations are stringently tested before being used, to ensure that they correctly meet the specification of the BLAS, and that they return the desired accuracy (see, for example,
Dodson et al. (1991), Dongarra et al. (1988) and Dongarra et al. (1990)).

# Background Information

Most of the methods in this chapter implement straightforward scalar,
vector and matrix operations that need no further explanation beyond a statement of the purpose of the method. In this section we give some additional background information to those few cases where additional explanation may be necessary. A sub-section is devoted to each topic.

# Real plane rotations

There are a number of methods in the chapter concerned with setting up and applying plane rotations. This section discusses the real case and the next section looks at the complex case. For further background information see
Golub and Van Loan (1996).

A plane rotation matrix for the $\left(i,j\right)$ plane,
${R}_{ij}$,
is an orthogonal matrix that is different from the unit matrix only in the elements ${r}_{ii}$,
${r}_{jj}$,
${r}_{ij}$ and ${r}_{ji}$. If we put

then, in the real case, it is usual to choose ${R}_{ij}$ so that

An exception is method (F06FPF not in this release) which applies the so-called symmetric rotation for which

The application of plane rotations is straightforward and needs no further elaboration, so further comment is made only on the construction of plane rotations.

$$R=\left(\begin{array}{cc}{r}_{ii}& {r}_{ij}\\ {r}_{ji}& {r}_{jj}\end{array}\right)\text{,}$$ | (1) |

$$R=\left(\begin{array}{cc}\phantom{-}c& s\\ -s& c\end{array}\right)\text{, \hspace{1em}}c=\mathrm{cos}\u200a\theta \text{, \hspace{1em}}s=\mathrm{sin}\u200a\theta \text{.}$$ |

$$R=\left(\begin{array}{cc}c& \phantom{-}s\\ s& -c\end{array}\right)\text{.}$$ | (2) |

The most common use of plane rotations is to choose $c$ and $s$ so that for given $a$ and $b$,

In such an application the matrix $R$ is often termed a

$$\left(\begin{array}{cc}\phantom{-}c& s\\ -s& c\end{array}\right)\text{\hspace{1em}}\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}d\\ 0\end{array}\right)\text{.}$$ | (3) |

**Givens rotation**matrix. There are two approaches to the construction of real Givens rotations in**F06**class.The BLAS method (F06AAF not in this release),
see
Lawson et al. (1979) and Dodson and Grimes (1982),
computes $c$, $s$ and $d$ as

where $\sigma =\left\{\begin{array}{cc}\mathrm{sign}\u200aa\text{,}& \left|a\right|>\left|b\right|\\ \mathrm{sign}\u200ab\text{,}& \left|a\right|\le \left|b\right|\end{array}\right.\text{.}$

$$d=\sigma {\left({a}^{2}+{b}^{2}\right)}^{1/2}\text{,}$$ |

$$c=\left\{\begin{array}{ll}a/d\text{,}& d\ne 0\text{,}\\ 1\text{,}& d=0\text{,}\end{array}\right.\text{\hspace{1em} \hspace{1em}}s=\left\{\begin{array}{ll}b/d\text{,}& d\ne 0\text{,}\\ 0\text{,}& d=0\text{,}\end{array}\right.$$ | (4) |

The value $z$ defined as

is also computed and this enables $c$ and $s$ to be reconstructed from the single value $z$ as

The other

where $\mathit{flmax}=1/\mathit{flmin}$ and $\mathit{flmin}$ is the small positive value returned by x02am. The values of $c$ and $s$ are then computed or reconstructed via $t$ as

where $\mathrm{eps}$ is the machine precision. Note that $c$ is always non-negative in this scheme and that the same expressions are used in the initial computation of $c$ and $s$ from $a$ and $b$ as in any subsequent recovery of $c$ and $s$ via $t$. This is the approach used by many of the NAG Library methods that require plane rotations. $d$ is computed simply as

You need not be too concerned with the above detail, since methods are provided for setting up, recovering and applying such rotations.

$$z=\left\{\begin{array}{ll}s\text{,}& \left|s\right|<c\text{\hspace{1em} or \hspace{1em}}c=0\\ 1/c\text{,}& 0<\left|c\right|\le s\end{array}\right.$$ | (5) |

$$c=\left\{\begin{array}{ll}0\text{,}& z=1\\ {\left(1-{z}^{2}\right)}^{1/2}\text{,}& \left|z\right|<1\\ 1/z\text{,}& \left|z\right|>1\end{array}\right.\text{\hspace{1em} \hspace{1em}}s=\left\{\begin{array}{ll}1\text{,}& z=1\\ z\text{,}& \left|z\right|<1\\ {\left(1-{c}^{2}\right)}^{1/2}\text{,}& \left|z\right|>1\end{array}\right.\text{.}$$ |

**F06**class methods for constructing Givens rotations are based on the computation of the tangent, $t=\mathrm{tan}\u200a\theta $. $t$ is computed as$$t=\left\{\begin{array}{ll}0\text{,}& b=0\\ b/a\text{,}& \left|b\right|\le \left|a\right|.\mathit{flmax},b\ne 0\\ \mathrm{sign}\left(b/a\right).\mathit{flmax}\text{,}& \left|b\right|>\left|a\right|.\mathit{flmax}\\ \mathrm{sign}\left(b\right).\mathit{flmax}\text{,}& b\ne 0,a=0\end{array}\right.$$ | (6) |

$$c=\left\{\begin{array}{ll}1/{\left(1+{t}^{2}\right)}^{1/2}\text{,}& \sqrt{\mathrm{eps}}\le \left|t\right|\le 1/\sqrt{\mathrm{eps}}\\ 1\text{,}& \left|t\right|<\sqrt{\mathrm{eps}}\\ 1/\left|t\right|\text{,}& \left|t\right|>1/\sqrt{\mathrm{eps}}\end{array}\right.\text{\hspace{1em} \hspace{1em}}s=\left\{\begin{array}{ll}c.t\text{,}& \sqrt{\mathrm{eps}}\le \left|t\right|\le 1/\sqrt{\mathrm{eps}}\\ t\text{,}& \left|t\right|<\sqrt{\mathrm{eps}}\\ \mathrm{sign}\u200at\text{,}& \left|t\right|>1/\sqrt{\mathrm{eps}}\end{array}\right.$$ | (7) |

$$d=c.a+s.b\text{.}$$ |

Another use of plane rotations is to choose $c$ and $s$ so that for given $x$, $y$ and $z$

In such an application the matrix $R$ is often termed a

$$\left(\begin{array}{cc}\phantom{-}c& s\\ -s& c\end{array}\right)\text{\hspace{1em}}\left(\begin{array}{cc}x& y\\ y& z\end{array}\right)\text{\hspace{1em}}\left(\begin{array}{cc}c& -s\\ s& \phantom{-}c\end{array}\right)=\left(\begin{array}{cc}a& 0\\ 0& b\end{array}\right)\text{.}$$ | (8) |

**Jacobi rotation**matrix. The method that generates a Jacobi rotation ( (F06BEF not in this release)) first computes the tangent $t$ and then computes $c$ and $s$ via $t$ as described above for the Givens rotation.# Complex plane rotations

In the complex case a plane rotation matrix for the $\left(i,j\right)$ plane,
${R}_{ij}$ is a unitary matrix and, analogously to the real case,
it is usual to choose ${R}_{ij}$ so that

where $\stackrel{-}{a}$ denotes the complex conjugate of $a$.

$$R=\left(\begin{array}{cc}\phantom{-}\stackrel{-}{c}& \stackrel{-}{s}\\ -s& c\end{array}\right)\text{, \hspace{1em}}{\left|c\right|}^{2}+{\left|s\right|}^{2}=1\text{,}$$ | (9) |

The BLAS (see Lawson et al. (1979)) do not contain a method for the generation of complex rotations, and so the methods in

**F06**class are all based upon computing $c$ and $s$ via $t=b/a$ in an analogous manner to the real case. $R$ can be chosen to have either $c$ real, or $s$ real and there are methods for both cases.When $c$ is real then it is non-negative and the transformation

is such that if $a$ is real then $d$ is also real.

$$\left(\begin{array}{cc}\phantom{-}c& \stackrel{-}{s}\\ -s& c\end{array}\right)\text{\hspace{1em}}\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}d\\ 0\end{array}\right)$$ | (10) |

When $s$ is real then the transformation

is such that if $b$ is real then $d$ is also real.

$$\left(\begin{array}{cc}\phantom{-}\stackrel{-}{c}& s\\ -s& c\end{array}\right)\text{\hspace{1em}}\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}d\\ 0\end{array}\right)$$ | (11) |

# Elementary real (Householder) reflections

There are a number of methods in the chapter concerned with setting up and applying Householder transformations. This section discusses the real case and the next section looks at the complex case. For further background information see
Golub and Van Loan (1996).

A real elementary reflector, $P$, is a matrix of the form

where $\mu $ is a scalar and $u$ is a vector,
and $P$ is both symmetric and orthogonal. In the methods in

because in many applications $\zeta $ and $z$ are not contiguous elements. The usual use of elementary reflectors is to choose $\mu $ and $u$ so that for given $\alpha $ and $x$

Such a transformation is often termed a

$$P=I-\mu u{u}^{\mathrm{T}}\text{, \hspace{1em}}\mu {u}^{\mathrm{T}}u=2\text{,}$$ | (12) |

**F06**class, $u$ is expressed in the form$$u=\left(\begin{array}{c}\zeta \\ z\end{array}\right)\text{, \hspace{1em}}\zeta \text{ a scalar}$$ | (13) |

$$P\left(\begin{array}{c}\alpha \\ x\end{array}\right)=\left(\begin{array}{c}\beta \\ 0\end{array}\right)\text{, \hspace{1em}}\alpha \text{ and}\beta \text{ scalars.}$$ | (14) |

**Householder transformation**. There are two choices of $\mu $ and $u$ available in**F06**class.The first form of the Householder transformation is compatible with that used by LINPACK (see Dongarra et al. (1979)) and has

$$\mu =1/\zeta \text{.}$$ | (15) |

This choice makes $\zeta $ satisfy

The second form, and the form used by many of the NAG Library
methods, has

which makes

In both cases the special setting

is used by the methods to flag the case where $P=I$.

$$1\le \zeta \le 2\text{.}$$ |

$$\mu =1$$ | (16) |

$$1\le \zeta \le \sqrt{2}\text{.}$$ |

$$\zeta =0$$ | (17) |

Note that while there are methods to apply an elementary reflector to a vector, there are no methods available in

and so we can call a matrix-vector product method to form $b={A}^{\mathrm{T}}u$ and then call a rank-one update method to form $\left(A-\mu u{b}^{\mathrm{T}}\right)$. Of course, we must skip the transformation when $\zeta $ has been set to zero.

**F06**class to apply an elementary reflector to a matrix. This is because such transformations can readily and efficiently be achieved by calls to the matrix-vector Level 2 BLAS methods. For example, to form $PA$ for a given matrix$$\begin{array}{lcl}PA& =& \left(I-\mu u{u}^{\mathrm{T}}\right)A=A-\mu u{u}^{\mathrm{T}}A\\ & =& A-\mu u{b}^{\mathrm{T}}\text{, \hspace{1em}}b={A}^{\mathrm{T}}u\text{,}\end{array}$$ | (18) |

# Elementary complex (Householder) reflections

A complex elementary reflector, $P$, is a matrix of the form

where ${u}^{\mathrm{H}}$ denotes the complex conjugate of ${u}^{\mathrm{T}}$, and $P$ is both Hermitian and unitary. For convenience in a number of applications this definition can be generalized slightly by allowing $\mu $ to be complex and so defining the generalized elementary reflector as

for which $P$ is still unitary, but is no longer Hermitian.

$$P=I-\mu u{u}^{\mathrm{H}}\text{, \hspace{1em}}\mu {u}^{\mathrm{H}}u=2\text{, \hspace{1em}}\mu \text{ real,}$$ |

$$P=I-\mu u{u}^{\mathrm{H}}\text{, \hspace{1em}}{\left|\mu \right|}^{2}{u}^{\mathrm{H}}u=\mu +\stackrel{-}{\mu}$$ | (19) |

The

and this reduces to (12) with the choice
(16) when $\mu $ and $u$ are real. This choice is used because $\mu $ and $u$ can now be chosen so that in the Householder transformation
(14) we can make

and, as in the real case,

Rather than returning $\mu $ and $\zeta $ as separate parameters the

Obviously $\zeta $ and $\mu $ can be recovered as

The special setting

is used to flag the case where $P=I$, and

is used to flag the case where

and in this case $\theta $ actually contains the value of $\gamma $. Notice that with both (18) and (21) we merely have to supply $\stackrel{-}{\theta}$ rather than $\theta $ in order to represent ${P}^{\mathrm{H}}$.

**F06**class methods choose $\mu $ and $\zeta $ so that$$\mathrm{Re}\left(\mu \right)=1\text{, \hspace{1em}}\mathrm{Im}\left(\zeta \right)=0$$ | (20) |

$$\mathrm{Im}\left(\beta \right)=0$$ |

$$1\le \zeta \le \sqrt{2}\text{.}$$ |

**F06**class methods return the single complex value $\theta $ defined as$$\theta =\zeta +i.\mathrm{Im}\left(\mu \right)\text{, \hspace{1em}}i=\sqrt{-1}\text{.}$$ |

$$\zeta =\mathrm{Re}\left(\theta \right)\text{, \hspace{1em}}\mu =1+i.\mathrm{Im}\left(\theta \right)\text{.}$$ |

$$\theta =0$$ |

$$\mathrm{Re}\left(\theta \right)\le 0\text{, \hspace{1em}}\mathrm{Im}\left(\theta \right)\ne 0$$ |

$$P=\left(\begin{array}{cc}\gamma & 0\\ 0& I\end{array}\right)\text{, \hspace{1em}}\gamma \text{ a scalar}$$ | (21) |

# References

Dodson D S and Grimes R G (1982) Remark on Algorithm 539

*ACM Trans. Math. Software***8**403–404Dodson D S, Grimes R G and Lewis J G (1991) Sparse extensions to the Fortran basic linear algebra subprograms

*ACM Trans. Math. Software***17**253–263Dongarra J J, Du Croz J J, Duff I S and Hammarling S (1990) A set of Level 3 basic linear algebra subprograms

*ACM Trans. Math. Software***16**1–28Dongarra J J, Du Croz J J, Hammarling S and Hanson R J (1988) An extended set of FORTRAN basic linear algebra subprograms

*ACM Trans. Math. Software***14**1–32Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979)

*LINPACK Users' Guide*SIAM, PhiladelphiaGolub G H and Van Loan C F (1996)

*Matrix Computations*(3rd Edition) Johns Hopkins University Press, BaltimoreLawson C L, Hanson R J, Kincaid D R and Krogh F T (1979) Basic linear algebra supbrograms for Fortran usage

*ACM Trans. Math. Software***5**308–325