\chapter[ Scientific Computations]{Scientific Computing}
This title refers to a fast developing interdisciplinary area between
mathematics, computers and applications. The subject is also often called as
Computational Science and Engineering. Its aim is the efficient use of
computer algorithms to solve engineering and scientific problems. One can say
with a certain simplification that our subject is related to numerical
mathematics, software engineering, computer graphics and applications. Here we
can deal only with some basic elements of the subject such as the fundamentals
of the \emph{floating point computer arithmetic,} \emph{error analysis,} the basic numerical
methods of \emph{linear algebra} and related mathematical software. \\


\section{Floating point arithmetic and error analysis}

\subsection{Classical error analysis}
Let $x$ be the exact value and let $a$ be an approximation of $x$ ($a\approx
x$). The error of the approximation $a$ is defined by the formula $\Delta
a=x-a$ (sometimes with opposite sign). The quantity $\delta a\geq0$ is called
an \ki{(absolute) error (bound)}\inddef{error bound} of approximation $a$, if $\left|  x-a\right|
=\left|  \Delta a\right|  \leq$ $\delta a$. For example, the error of the
approximation $\sqrt{2}\approx1.41$ is at most $0.01$. In other words, the
error bound of the approximation is $0.01$. The quantities $x$ and $a$ (and
\ accordingly $\Delta a$ and $\delta a$) may be vectors or matrices. In such
cases the absolute value and relation operators must be understood
componentwise. We also measure the error by using matrix and vector norms. In
such cases, the quantity $\delta a\in\mathbb{R}$ is an error bound, if the
inequality $\left\|  \Delta a\right\|  \leq$ $\delta a$ holds.

The absolute error bound can be irrelevant in many cases. For example, an
approximation with error bound $0.05$ has no value in estimating a quantity of
order $0.001$. The goodness of an approximation is measured by the \ki{relative
error}\inddef{relative
error} $\delta a/\left|  x\right|  $ ($\delta a/\left\|  x\right\|  $ for
vectors and matrices), which compares the error bound to the approximated
quantity. Since the exact value is generally unknown, we use the approximate
relative error $\delta a/\left|  a\right|  $ ($\delta a/\left\|  a\right\|
$). The committed error is proportional to the quantity $\left(  \delta
a\right)  ^{2}$, which can be neglected, if the absolute value (norm) of $x$
and $a$ is much greater than $\left(  \delta a\right)  ^{2}$. The relative
error is often expressed in percentages.


In practice, the (absolute) error bound is used as a substitute for the
generally unknown true error.

In the \emph{classical error analysis} \index{classical error analysis} we assume input data with given error bounds,
exact computations (operations) and seek for the error bound of the final
result. Let $x$ and $y$ be exact values with approximations $a$ and $b$,
respectively. Assume that the absolute error bounds of approximations $a$ and
$b$ are $\delta a$ and $\delta b$, respectively. Using the classical error
analysis approach we obtain the following error bounds for the four basic
arithmetic operations:%

\begin{align*}
\delta\left(  a+b\right)   &  =\delta a+\delta b,\qquad\frac{\delta
(a+b)}{\left|  a+b\right|  }=\max\left\{  \frac{\delta a}{\left|  a\right|
},\frac{\delta b}{\left|  b\right|  }\right\}  \quad\left(  ab>0\right) \koz ,\\
\delta\left(  a-b\right)   &  =\delta a+\delta b,\qquad\frac{\delta
(a-b)}{\left|  a-b\right|  }=\frac{\delta a+\delta b}{\left|  a-b\right|
}\quad\left(  ab>0\right)  \koz ,\\
\delta\left(  ab\right)   &  \approx\left|  a\right|  \delta b+\left|
b\right|  \delta a\qquad\frac{\delta(ab)}{\left|  ab\right|  }\approx
\frac{\delta a}{\left|  a\right|  }+\frac{\delta b}{\left|  b\right|  }%
\quad\left(  ab\neq0\right)  \koz ,\\
\delta(a/b)  &  \approx\frac{\left|  a\right|  \delta b+\left|  b\right|
\delta a}{\left|  b\right|  ^{2}}\qquad\frac{\delta(a/b)}{\left|  a/b\right|
}\approx\frac{\delta a}{\left|  a\right|  }+\frac{\delta b}{\left|  b\right|
}\quad\left(  ab\neq0\right)  \koz .
\end{align*}


We can see that the division with a number near to $0$ can make the absolute
error arbitrarily big. Similarly, if the result of subtraction is near to $0$,
then its relative error can become arbitrarily big. One has to avoid these
cases. Especially the subtraction operation can be quite dangerous.

\begin{pelda}
Calculate the quantity $\sqrt{1996}-\sqrt{1995}$ with approximations
$\sqrt{1996}\approx44.67$ and $\sqrt{1995}\approx44.66$ whose common absolute
and relative error bounds are $0.01$ and $0.022\%$, respectively. One obtains
the approximate value $\sqrt{1996}-\sqrt{1995}\approx0.01$, whose relative
error bound is%
\[
\frac{0.01+0.01}{0.01}=2 \koz ,
\]
that is $200\%$. The true relative error is about $10.66\%$. Yet it is too big,
since it is approximately $5\times10^{2}$ times bigger than the relative error
of the initial data. \ We can avoid the subtraction operation by using the
following trick%
\[
\sqrt{1996}-\sqrt{1995}=\frac{1996-1995}{\sqrt{1996}+\sqrt{1995}}=\frac
{1}{\sqrt{1996}+\sqrt{1995}}\approx\frac{1}{89.33}\approx0.01119 \koz .
\]
Here the nominator is exact, while the absolute error of the denominator is
$0.02$. Hence the relative error (bound) of the quotient is about
$0.02/89.33\approx0.00022=0.022\%$. The latter result is in agreement with the
relative error of the initial data and it is substantially smaller than the
one obtained with direct subtraction operation.
\end{pelda}

The first order error terms of twice differentiable functions can be obtained
by their first order \emph{Taylor polynomial}\nevindex{Taylor, Brook}:%
\begin{align*}
\delta\left(  f\left(  a\right)  \right)   &  \approx\left|  f^{\prime
}(a)\right|  \delta a,\quad f:\mathbb{R}\rightarrow\mathbb{R} \koz ,\\
\delta\left(  f\left(  a\right)  \right)   &  \approx\sum_{i=1}^{n}\left|
\frac{\partial f(a)}{\partial x_{i}}\right|  \delta a_{i},\quad f:\mathbb{R}%
^{n}\rightarrow\mathbb{R} \koz .
\end{align*}
The numerical sensitivity of functions at a given point is characterised by
the \ki{condition number,}\inddef{condition number} which is the ratio of the relative errors of approximate
function value and the input data (the \emph{Jacobian matrix} of functions
$F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is denoted by $F^{\prime}(a)$ at
the point $a\in\mathbb{R}^{n}$):%
\begin{align*}
c(f,a)  &  =\frac{\left|  f^{\prime}(a)\right|  \left|  a\right|  }{\left|
f(a)\right|  },\quad f:\mathbb{R}\rightarrow\mathbb{R} \koz ,\\
c(F,a)  &  =\frac{\left\|  a\right\|  \left\|  F^{\prime}(a)\right\|
}{\left\|  F(a)\right\|  },\quad F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m} \koz .
\end{align*}


We can consider the condition number as the magnification number of the input
relative error. Therefore the functions is considered \ki{numerically stable}\inddef{numerically stable} (or
\ki{well-conditioned}\inddef{well-conditioned}) at the point $a$, if $c\left(  f,a\right)  $ is
\idez{small}.
Otherwise $f$ is considered as \ki{numerically unstable}\inddef{numerically unstable}
(\ki{ill-conditioned.}\inddef{ill-conditioned}) The
condition number depends on the point $a$. A function can be well-conditioned
at point $a$, while it is ill-conditioned at point $b$. The term
\idez{small} is
relative. It depends on the problem, the computer and the required precision.

The condition number of matrices can be defined as the upper bound of a
function condition number. Let us define the mapping $F:\mathbb{R}%
^{n}\rightarrow\mathbb{R}^{n}$ by the solution of the equation $Ay=x$
($A\in\mathbb{R}^{n\times n}$, $\det(A)\neq0$), that is, let $F(x)=A^{-1}x$ .
Then $F^{\prime}\equiv A^{-1}$ and
\[
c(F,a)=\frac{\left\|  a\right\|  \left\|  A^{-1}\right\|  }{\left\|
A^{-1}a\right\|  }=\frac{\left\|  Ay\right\|  \left\|  A^{-1}\right\|
}{\left\|  y\right\|  }\leq\left\|  A\right\|  \left\|  A^{-1}\right\|
\quad\left(  Ay=a\right)  \koz .
\]


The upper bound of the right side is called \ the \ki{condition number of the
matrix}\inddef{condition number of the
matrix} $A$. This bound is sharp, since there exists a vector $a\in
\mathbb{R}^{n\text{ }}$ such that $c(F,a)=\left\|  A\right\|  \left\|
A^{-1}\right\|.$

\subsection{Forward and backward errors}
Let us investigate the calculation of the function value $f\left(  x\right)
$. If we calculate the approximation $\hat{y}$ instead of the exact value
$y=f\left(  x\right)  $, then the forward error $\Delta y=\hat{y}-y$. If for a
value $x+\Delta x$ the equality $\hat{y}=f\left(  x+\Delta x\right)  $ holds,
that is, $\hat{y}$ is the exact function value of the perturbed input data
$\hat{x}=x+\Delta x$, then $\Delta x$ is called the \ki{backward error.}\inddef{backward error} The
connection of the two concepts is shown on the Figure \ref{fig:numerr}.%

%TCIMACRO{\FRAME{dtbpF}{10.8601cm}{2.5788cm}{0pt}{}{}{numerr.wmf}%
%{\special{ language "Scientific Word";  type "GRAPHIC";
%maintain-aspect-ratio TRUE;  display "USEDEF";  valid_file "F";
%width 10.8601cm;  height 2.5788cm;  depth 0pt;  original-width 3.8035in;
%original-height 0.8579in;  cropleft "-0.0046";  croptop "1.0205";
%cropright "1.0046";  cropbottom "-0.0205";
%filename 'C:/Documents and Settings/matjy/Local Settings/Temporary Internet Files/Content.IE5/SPAJG5QN/numerr.wmf';file-properties "NPEU";}%
%}}%
%BeginExpansion
%\begin{center}
%\includegraphics[
%trim=-0.017496in -0.017587in -0.017496in -0.017587in,
%natheight=0.857900in,
%natwidth=3.803500in,
%height=2.5788cm,
%width=10.8601cm
%]%
%{C:/Documents and Settings/matjy/Local Settings/Temporary Internet Files/Content.IE5/SPAJG5QN/numerr.wmf}%
%\end{center}
%EndExpansion

\epskep{figs/Numerr}{Forward and backward error.\label{fig:numerr}}

The continuous line shows exact value, while the dashed one indicates computed
value. The analysis of the backward error is called the \ki{backward error analysis.} 
\inddef{backward error analysis} If there exist more than one backward error, then the estimation of
the smallest one is the most important.

An algorithm for computing the value $y=f\left(  x\right)  $ is called
\ki{backward stable,}\inddef{backward stable} if for any $x$ it gives a computed value $\hat{y}$ with small
 backward error $\Delta x$. Again, the term \idez{small} is relative to the
problem environment.

The connection of the forward and backward errors is described by the
approximate\ thumb rule%

\begin{equation}
\frac{\delta\hat{y}}{\left|  y\right|  }\lessapprox c\left(  f,x\right)
\frac{\delta\hat{x}}{\left|  x\right|  } \koz ,
\end{equation}
which means that
\[
\text{relative forward error}\leq\text{condition number}\times
\text{relative backward error.}%
\]


This inequality indicates that the computed solution of an ill-conditioned
problem may have a big relative forward error. An algorithm is said to be
\ki{forward stable}\inddef{forward stable} if the forward error is small. A forward stable method is not
necessarily backward stable. If the forward error and the condition number are
small, then the algorithm is forward stable.

\begin{pelda}
Consider the function $f\left(  x\right)  =\log x$ the condition number of
which is $c\left(  f,x\right)  =c\left(  x\right)  =1/\left|  \log x\right|
$. For $x\approx1$ the condition number $c\left(  f,x\right)  $ is big.
Therefore the relative forward error is big for $x\approx1$.
\end{pelda}

\subsection{Rounding errors and floating point arithmetic}
The classical error analysis investigates only the effects of the input data
errors and assumes exact arithmetic operations. The digital computers however
are representing the numbers with a finite number of digits, the arithmetic
computations are carried out on the elements of a finite set $F$ of such
numbers and the results of operations belong to $F$. Hence the computer
representation of the numbers may add further errors to the input data and the
results of arithmetic operations may also be subject to further rounding. If
the result of operation belongs to $F$, then we have the exact result.
Otherwise we have three cases:

(i) rounding to representable (nonzero) number;

(ii) underflow (rounding to $0$);

(iii) overflow (in case of results whose moduli too large).

The most of the scientific-engineering calculations are done in \emph{floating point
arithmetic}\index{floating point arithmetic} whose generally accepted model is the following:

\begin{defi}
The set of floating point numbers is given by%
\[%
\begin{array}
[c]{c}%
F(\beta,t,L,U)=\\
=\left\{  \pm m\times\beta^{e}\mid\frac{1}{\beta}\leq m<1,\text{ }%
m=0.d_{1}d_{2}\ldots d_{t},\text{ }L\leq e\leq U\right\}  \cup\left\{
0\right\}  \koz ,
\end{array}
\]
where\newline\quad{- }$\beta$\ is the base (or radix) of the number
system,\newline\quad- $m$\ is the mantissa in the number system with base
$\beta$,\newline\quad- $e$\ is the exponent,\newline{\quad- }$t$\ is the
length of mantissa (the precision of arithmetic),\newline{\quad- }$L$\ is the
smallest exponent (underflow exponent),\newline{\quad- }$U$\ is the biggest
exponent (overflow exponent).
\end{defi}

The parameters of the three most often used number systems are indicated in
the following table
\[%
\begin{tabular}
[c]{|c|c|c|}\hline
Name & $\beta$ & Machines\\\hline
binary & 2 & most computer\\\hline
decimal & 10 & most calculators\\\hline
hexadecimal & 16 & IBM mainframe computers\\\hline
\end{tabular}
\]
The mantissa can be written in the form%
\begin{equation}
m=0.d_{1}d_{2}\ldots d_{t}=\frac{d_{1}}{\beta}+\frac{d_{2}}{\beta^{2}}%
+\cdots+\frac{d_{t}}{\beta^{t}}. \label{mant1}%
\end{equation}
We can observe that condition $1/\beta\leq m<1$ implies the inequality $1\leq
d_{1}\leq\beta-1$ for the first digit $d_{1}$. The remaining digits must
satisfy $0\leq d_{i}\leq\beta-1$ ($i=2,\ldots,t$). \ Such arithmetic systems
are called \ki{normalized.}\inddef{normalized} \ The zero digit and the dot is not represented. If
$\beta=2$, then the first digit is $1$, which is also unrepresented. Using the
representation (\ref{mant1}) we can give the set $F=F(\beta,t,L,U)$ in the
form%
\begin{equation}
F=\left\{  \pm\left(  \frac{d_{1}}{\beta}+\frac{d_{2}}{\beta^{2}}+\cdots
+\frac{d_{t}}{\beta^{t}}\right)  \beta^{e}\mid\text{ }L\leq e\leq U\right\}
\cup\left\{  0\right\}  \koz ,
\end{equation}
where $0\leq d_{i}\leq\beta-1\ (i=1,\ldots,t)$ and $1\leq d_{1}$.

\begin{pelda}
The set $F\left(  2,3,-1,2\right)  $ contains $33$ elements and its positive
elements are given by
\[
\left\{  \frac{1}{4},\frac{5}{16},\frac{6}{16},\frac{7}{16},\frac{1}{2}%
,\frac{5}{8},\frac{6}{8},\frac{7}{8},1,\frac{10}{8},\frac{12}{8},\frac{14}%
{8},2,\frac{20}{8},3,\frac{28}{8}\right\} \koz .
\]
\end{pelda}

The elements of $F$ are not equally distributed on the real line. The distance
of two consecutive numbers in $\left[  1/\beta,1\right]  \cap F$ is
$\beta^{-t}$. Since the elements of $F$ are of the form $\pm m\times\beta^{e}%
$, the distance of two consecutive numbers in $F$ is changing with the
exponent. The maximum distance of two consecutive floating point numbers is
$\beta^{U-t}$, while the minimum distance is $\beta^{L-t}$.

For the mantissa we have $m\in\left[  1/\beta,1-1/\beta^{t}\right]  ,$ since
\[
\frac{1}{\beta}\leq m=\frac{d_{1}}{\beta}+\frac{d_{2}}{\beta^{2}}+\cdots
+\frac{d_{t}}{\beta^{t}}\leq\frac{\beta-1}{\beta}+\frac{\beta-1}{\beta^{2}%
}+\cdots+\frac{\beta-1}{\beta^{t}}=1-\frac{1}{\beta^{t}}\text{.}%
\]


Using this observation we can easily prove the following result on the range
of floating point numbers.

\begin{tetel}
If $a\in F$, $a\neq0$, then $M_{L}\leq\left|  a\right|  \leq M_{U}$, where
\[
M_{L}=\beta^{L-1}\text{,\quad}M_{U}=\beta^{U}(1-\beta^{-t}).
\]

\end{tetel}

Let $a,b\in F$ and denote $\square$ any of the four arithmetic operations
$(+,-,\ast,/).$ The following cases are possible:

(1) $a\square b\in F$ (exact result),

(2) $\left|  a\square b\right|  >M_{U}$ (arithmetic overflow),

(3) $0<\left|  a\square b\right|  <M_{L}$ (arithmetic underflow),

(4) $a\square b\notin F$, $M_{L}<\left|  a\square b\right|  <M_{U}$ (not
representable result).

In the last two cases the \emph{floating point arithmetic}\index{floating point
arithmetic} is rounding the result
$a\square b$ to the nearest floating point number in $F$. If two consecutive
floating point numbers are equally distant from $a\square b$, then we
generally round to the greater number. For example, in a five digit decimal
arithmetic, the number $2.6457513$ is rounded to the number $2.6458$.

Let $G=\left[  -M_{U},M_{U}\right]  $. It is clear that $F\subset G$. Let
$x\in G$. The $fl\left(  x\right)  $ denotes an element of $F$ nearest to $x$.
The mapping $x\rightarrow fl\left(  x\right)  $ is called rounding. The
quantity $\left|  x-fl\left(  x\right)  \right|  $ is called the rounding
error. If $fl\left(  x\right)  =1$, then the rounding error is at most
$\beta^{1-t}/2$. The quantity $u=\beta^{1-t}/2$ is called the \ki{unit roundoff.}\inddef{unit roundoff} 
The quantity $u$ is the relative error bound\index{relative error bound} of $fl\left(  x\right)  $.

\begin{tetel}
If $x\in G$, then
\[
fl(x)=x(1+\varepsilon),\quad\left|  \varepsilon\right|  \leq u \koz .
\]

\end{tetel}

\begin{biz}
Without loss of generality we can assume that $x>0$. Let $m_{1}\beta^{e}$,
$m_{2}\beta^{e}\in F$ be two consecutive numbers such that%
\[
m_{1}\beta^{e}\leq x\leq m_{2}\beta^{e} \koz .
\]
Either $1/\beta\leq m_{1}<m_{2}\leq1-\beta^{-t}$ or $1-\beta^{-t}=m_{1}%
<m_{2}=1$ holds. Since $m_{2}-m_{1}=\beta^{-t}$ holds in both cases, we have
\[
\left|  fl\left(  x\right)  -x\right|  \leq\frac{\left|  m_{2}-m_{1}\right|
}{2}\beta^{e}=\frac{\beta^{e-t}}{2}%
\]
either $fl\left(  x\right)  =m_{1}\beta^{e}$ or $fl\left(  x\right)
=m_{2}\beta^{e}$. It follows that%
\[
\frac{\left|  fl\left(  x\right)  -x\right|  }{\left|  x\right|  }\leq
\frac{\left|  fl\left(  x\right)  -x\right|  }{m_{1}\beta^{e}}\leq\frac
{\beta^{e-t}}{2m_{1}\beta^{e}}=\frac{\beta^{-t}}{2m_{1}}\leq\frac{1}{2}%
\beta^{1-t}=u \koz .
\]
Hence $fl\left(  x\right)  -x=\lambda xu$, where $\left|  \lambda\right|
\leq1$. A simple arrangement yields
\[
fl\left(  x\right)  =x(1+\varepsilon)\quad\left(  \varepsilon=\lambda
u\right)
\]
Since $\left|  \varepsilon\right|  \leq u,$ we proved the claim.
\end{biz}

Thus we proved that the relative error of the rounding is bounded in floating
point arithmetic\index{floating
point arithmetic} and the bound is the unit roundoff $u$.

Another quantity used to measure the rounding errors is the so called the
\ki{machine epsilon}\inddef{machine epsilon} $\epsilon_{M}=2u=\beta^{1-t}$ ($\epsilon_{M}=2u$). The number
$\epsilon_{M}$ is the distance of $1$ and its nearest neighbour greater than
$1$. The following algorithm determines $\epsilon_{M}$ in the case of binary base.

\begin{alg}{Machine-Epsilon}\inddef{\textsc{Machine-Epsilon}}
1 \' $ x\leftarrow 1$\\
2 \'  \key{while} \= $1+x>1$\\
3 \'  \> \key{do} \=$ x\leftarrow x/2$\\
4 \'  $\epsilon _{M}\leftarrow 2x$ \\
5 \' \key{return} $\epsilon _{M}$
\end{alg}

\noindent\ In the MATLAB system $\epsilon_{M}\approx2.2204\times10^{-16}$.

For the results of floating point arithmetic operations we assume the
following (standard model):%
\begin{equation}
fl(a\square b)=\left(  a\square b\right)  \left(  1+\varepsilon\right)
,\quad\left|  \varepsilon\right|  \leq u\quad\left(  a,b\in F\right)  \koz .
\label{standmod}%
\end{equation}


The IEEE arithmetic standard satisfies this assumption. It is an important
consequence of the assumption that for $a\square b\neq0$ the relative error of
arithmetic operations satisfies
\[
\frac{\left|  fl(a\square b)-(a\square b)\right|  }{\left|  a\square b\right|
}\leq u.
\]
Hence the relative error of the floating point arithmetic operations is small.

There exist computer floating point arithmetics that do not comply with the
standard model (\ref{standmod}). The usual reason for this is that the
arithmetic lacks a guard digit in subtraction. For simplicity we investigate
the subtraction $1-0.111$ in a three digit binary arithmetic. In the first
step we equate the exponents:
\[%
\begin{array}
[c]{ccccccccc}
& 2 & \times & 0 & . & 1 & 0 & 0 & \\
- & 2 & \times & 0 & . & 0 & 1 & 1 & 1
\end{array}
 \koz .
\]
If the computation is done with four digits, the result is the following
\[%
\begin{array}
[c]{ccccccccc}
& 2^{1} & \times & 0 & . & 1 & 0 & 0 & \\
- & 2^{1} & \times & 0 & . & 0 & 1 & 1 & 1\\\hline
& 2^{1} & \times & 0 & . & 0 & 0 & 0 & 1
\end{array}
\koz ,
\]
from which the normalized result is $2^{-2}\times0.100$. Observe that the
subtracted number is unnormalised. The temporary fourth digit of the mantissa
is called a guard digit. Without a guard digit the computations are the
following:
\[%
\begin{array}
[c]{cccccccc}
& 2^{1} & \times & 0 & . & 1 & 0 & 0\\
- & 2^{1} & \times & 0 & . & 0 & 1 & 1\\\hline
& 2^{1} & \times & 0 & . & 0 & 0 & 1
\end{array}
\koz .
\]
Hence the normalized result is $2^{-1}\times0.100$ with a relative error of
$100\%$. Several CRAY computers and pocket calculators lack guard digits.

Without the guard digit the floating point arithmetic operations satisfy only
the weaker conditions
\begin{align}
fl\left(  x\pm y\right)   &  =x\left(  1+\alpha\right)  \pm y\left(
1+\beta\right)  ,\quad\left|  \alpha\right|  ,\left|  \beta\right|  \leq u,\\
fl\left(  x\square y\right)   &  =\left(  x\square y\right)  \left(
1+\delta\right)  ,\quad\left|  \delta\right|  \leq u,\quad\square=\ast,/\, \koz .
\end{align}


Assume that we have a guard digit and the arithmetic complies with standard
model (\ref{standmod}). Introduce the following notations:
\begin{equation}
\left|  z\right|  =\left[  \left|  z_{1}\right|  ,\ldots,\left|  z_{n}\right|
\right]  ^{T}\quad\left(  z\in\mathbb{R}^{n}\right)  \koz ,
\end{equation}%
\begin{equation}
\left|  A\right|  =\left[  \left|  a_{ij}\right|  \right]  _{i,j=1}^{m,n}%
\quad\left(  A\in\mathbb{R}^{m\times n}\right)  ,
\end{equation}%
\begin{equation}
A\leq B\Leftrightarrow a_{ij}\leq b_{ij}\quad\left(  A,B\in\mathbb{R}^{m\times
n}\right)  \koz .
\end{equation}
The following results hold:
\begin{equation}
\left|  fl\left(  x^{T}y\right)  -x^{T}y\right|  \leq1.01nu\left|  x\right|
^{T}\left|  y\right|  \quad\left(  nu\leq0.01\right)  \koz ,
\end{equation}%
\begin{equation}
fl\left(  \alpha A\right)  =\alpha A+E\quad\left(  \left|  E\right|  \leq
u\left|  \alpha A\right|  \right)  \koz ,
\end{equation}%
\begin{equation}
fl\left(  A+B\right)  =\left(  A+B\right)  +E\quad\left(  \left|  E\right|
\leq u\left|  A+B\right|  \right)  \koz ,
\end{equation}%
\begin{equation}
fl\left(  AB\right)  =AB+E\quad\left(  \left|  E\right|  \leq nu\left|
A\right|  \left|  B\right|  +O\left(  u^{2}\right)  \right) \koz   ,
\end{equation}
where $E$ denotes the error (matrix) of the actual operation.

The standard \emph{floating point arithmetics}\index{floating point
arithmetic} have many special properties. It is an
important property that the addition is not associative because of the rounding.

\begin{pelda}
If $a=1$, $b=c=3\times10^{-16}$, then using MATLAB and AT386 type PC we
obtain
\[
1.000000000000000e+000=(a+b)+c\neq a+(b+c)=1.000000000000001e+000 \koz .
\]
We can have a similar result on Pentium1 machine with the choice
$b=c=1.15\times10^{-16}$.
\end{pelda}

The example also indicates that for different (numerical) processors may
produce different computational results for the same calculations. The
commutativity can also be lost in addition. Consider the computation of the
sum $\sum_{i=1}^{n}x_{i}$. The usual algorithm is the recursive summation.

\begin{alg}{Recursive-Summation($n,x$)}\inddef{\textsc{Recursive-Summation}}
1 \' $s\leftarrow 0$ \\
2 \'  \key{for} $i\leftarrow 1$ to $n$ \\
3 \'  \> \key{do} \= $s\leftarrow s+x_{i}$ \\
4 \'  \key{return} $s$
\end{alg}


\begin{pelda}
Compute the sum
\[
s_{n}=1+\sum_{i=1}^{n}\frac{1}{i^{2}+i}%
\]
for $n=4999$. The recursive summation algorithm (and MATLAB) gives the result%
\[
1.999800000000002e+000 \koz .
\]
If the summation is done in the reverse (increasing) order, then the result
is
\[
1.999800000000000e+000 \koz .
\]
If the two values are compared with the exact result $s_{n}=2-1/(n+1)$, then
we can see that the second summation gives better result. In this case the sum
of smaller numbers gives significant digits to the final result unlike in the
first case.
\end{pelda}

The last example indicates that the summation of a large number of data
varying in modulus and sign is a complicated task. The following algorithm of
W. Kahan\nevindex{Kahan, W.} is one of the most interesting procedures to solve the problem.

\begin{alg}{Compensated-Summation($n,x$)}\inddef{\textsc{Compensated-Summation}}
1 \' $s\leftarrow 0$ \\
2 \' $e\leftarrow 0$ \\
3 \' \key{for } \= $i\leftarrow 1$ \key{to} $n$ \\
4 \'  \> \key{do} \= $t \leftarrow s$ \\
5 \'  \> \> $y\leftarrow x_{i}  +e$ \\
6 \'  \> \> $s\leftarrow t+y$ \\
7 \'  \> \> $e\leftarrow \left(  t-s\right)  +y$ \\
8 \' \key{return} $s$
\end{alg}


\subsection{The floating point arithmetic standard}
The ANSI/IEEE Standard 754-1985 of a binary ($\beta=2$) \emph{floating point
arithmetic system}\index{floating point
arithmetic system} was published in 1985. The standard specifies the basic
arithmetic operations, comparisons, rounding modes, the arithmetic exceptions
and their handling, and conversion between the different arithmetic formats.
The square root is included as a basic operation. The standard does not deal
with the exponential and transcendent functions. The standard defines two main
floating point formats:
\[%
\begin{tabular}
[c]{|l|l|l|l|l|l|}\hline
Type & Size & Mantissa & $e$ & $u$ & $\left[  M_{L},M_{U}\right]  \approx
$\\\hline
Single & $32$ bits & $23+1$ bits & $8$ bits & $2^{-24}\approx5.96\times
10^{-8}$ & $10^{\pm38}$\\\hline
Double & $64$ bits & $52+1$ bits & $11$ bits & $2^{-53}\approx1.11\times
10^{-16}$ & $10^{\pm308}$\\\hline
\end{tabular}
\]
In both formats one bit is reserved as a sign bit. Since the floating point
numbers are normalized and the first digit is always $1$, this bit is not
stored. This hidden bit is denoted by the \idez{$+1$} in the table.

The arithmetic standard contains the handling of arithmetic exceptions.%

\[%
\begin{tabular}
[c]{|l|l|l|}\hline
Exception type & Example & Default result\\\hline
Invalid operation & $0/0$, $0\times\infty$, $\sqrt{-1}$ & NaN (Not a
Number)\\\hline
Overflow & $\left|  x\square y\right|  >M_{U}$ & $\pm\infty$\\\hline
Divide by zero & Finite nonzero$/0$ & $\pm\infty$\\\hline
Underflow & $0<\left|  x\square y\right|  <M_{L}$ & Subnormal numbers\\\hline
Inexact & $fl\left(  x\square y\right)  \neq x\square y$ & Correctly rounded
result\\\hline
\end{tabular}
\]
(The numbers of the form $\pm m\times\beta^{L-t},$ $0<m<\beta^{t-1}$ are
called \ki{subnormal numbers.}\inddef{subnormal number}) The IEEE arithmetic is a closed system. Every
arithmetic operations has a result, whether it is expected mathematically or
not. The exceptional operations raise a signal and continue. The arithmetic
standard conforms with the standard model (\ref{standmod}).

The first hardware implementation of the IEEE standard was the Intel 8087
mathematical coprocessor. Since then it is generally accepted and used.

\textit{Remark.}
In the single precision we have about 7 significant digit precision in the
decimal system. For double precision we have approximately 16 digit precision
in decimals. There also exists an extended precision format of 80 bits, where
$t=63$ and the exponential has $15$ bits.

\begin{gyak}
\ujgyak The measured values of two resistors are $R_{1}=110.2\pm0.3\Omega$ and
$R_{2}=65.6\pm0.2\Omega$. We connect the two resistors parallel and obtain the
circuit resistance $R_{e}=R_{1}R_{2}/(R_{1}+R_{2})$. Calculate the relative
error bounds of the initial data and the approximate value of the resistance
$R_{e}$. Evaluate the absolute and relative error bounds $\delta R_{e}$ and
$\delta R_{e}/R_{e}$, respectively in the following three ways: \newline(i)
Estimate first $\delta R_{e}$ using only the absolute error bounds of the
input data, then estimate the relative error bound $\delta R_{e}/R_{e}$.
\newline(ii) Estimate first the relative error bound $\delta R_{e}/R_{e}$
using only the relative error bounds of the input data, then estimate the
absolute error bound $\delta R_{e}$.\newline(iii) Consider the circuit
resistance as a two variable function $R_{e}=F\left(  R_{1},R_{2}\right)  $.
\ujgyak Assume that $\sqrt{2}$ is calculated with the absolute error bound $10^{-8}$.
The following two expressions are theoretically equal:\newline(i) $1/\left(
1+\sqrt{2}\right)  ^{6};$\newline(ii) $99-70\sqrt{2}.$\newline Which
expression can be calculated with less relative error and why?
\ujgyak Consider the arithmetic operations as two variable functions of the form
$f\left(  x,y\right)  =x\square y$, where $\square\in\left\{  +,-,\ast
,/\right\}  $.\newline(i) Derive the \emph{error bounds}\gyakindex{error bound} of the arithmetic operations
from the error formula of two variable functions.\newline(ii) Derive the
\emph{condition numbers}\gyakindex{condition number} of these functions. When are they ill-conditioned?\newline%
(iii) Derive \emph{error bounds}\gyakindex{error bound} for the power function assuming that both the base
and the exponent have errors. What is the result if the exponent is
exact?\newline(iv) Let $y=16x^{2}$, $x\approx a$ and $y\approx b=16a^{2}$.
Determine the smallest and the greatest value of $a$ as a function of $x$ such
that the \emph{relative error bound}\gyakindex{relative error bound} of $b$ should be at most $0.01$.
\ujgyak Assume that the number $C=\exp(4\pi^{2}/\sqrt{83})$ ($=76.1967868\ldots$) is
calculated in a 24 bit long mantissa and the exponential function is also
calculated with 24 significant bits. Estimate the \emph{absolute error}\gyakindex{absolute error} of the
result. Estimate the \emph{relative error}\gyakindex{relative error} without using the actual value of $C$.
\ujgyak Consider the emph{floating point number set}\gyakindex{floating point number set} $F(\beta,t,L,U)$ and show
that\newline(i) Every arithmetic operation can result \emph{arithmetic
overflow}\gyakindex{arithmetic overflow};\newline(ii) Every arithmetic operation can result \emph{arithmetic underflow}.
\gyakindex{arithmetic underflow}
\ujgyak Show that the following expressions are \emph{numerically unstable}\gyakindex{numerically unstable} for $x\approx
0$:\newline(i) $(1-\cos x)/\sin^{2}x;$\newline(ii) $\sin(100\pi+x)-\sin
(x)$;\newline(iii) $2-\sin x-\cos x-e^{-x}.$\newline Calculate the values of
the above expressions for $x=10^{-3},10^{-5},10^{-7}$ and estimate the error.
Manipulate the expressions into \emph{numerically stable}\gyakindex{numerically stable} ones and estimate the error
as well.
\ujgyak How many elements does the set $F=F\left(  \beta,t,L,U\right)  $ have? How
many \emph{subnormal numbers}\gyakindex{subnormal number} can we find?
\ujgyak If $x,y\geq0$, then $\left(  x+y\right)  /2\geq\sqrt{xy}$ and equality holds
if and only if $x=y$. Is it true numerically? Check the inequality
experimentally for various data (small and large numbers, numbers close to
each other or different in magnitude).
\end{gyak}

\section{Linear systems of equations}
The general form of linear algebraic systems with $n$ unknowns and $m$ equations is given by
\begin{equation}%
\begin{array}
[c]{ccc}%
a_{11}x_{1}+\cdots+a_{1j}x_{j}+\cdots+a_{1n}x_{n} & = & b_{1}\\
& \vdots & \\
a_{i1}x_{1}+\cdots+a_{ij}x_{j}+\cdots+a_{in}x_{n} & = & b_{i}\\
& \vdots & \\
a_{m1}x_{1}+\cdots+a_{mj}x_{j}+\cdots+a_{mn}x_{n} & = & b_{m}%
\end{array}
\end{equation}
This system can be written in the more compact form
\begin{equation}
Ax=b \koz ,
\end{equation}
where
\[
A=\left[  a_{ij}\right]  _{i,j=1}^{m,n}\in\mathbb{R}^{m\times n}%
,\ x\in\mathbb{R}^{n},\ b\in\mathbb{R}^{m} \koz .
\]
The systems is called \emph{underdetermined} if $m<n$. For $m>n$, the systems is
called \emph{overdetermined}. Here we investigate only the case $m=n$, when the
coefficient matrix $A$ is square. We also assume that the inverse matrix
$A^{-1}$ exists (or equivalently $\det\left(  A\right)  \neq0$). Under this
assumption the linear system $Ax=b$ has exactly one solution: $x=A^{-1}b$.

\subsection{Direct methods for solving linear systems}

\subsubsection{Triangular linear systems}

\begin{defi}
The matrix $A=[a_{ij}]_{i,j=1}^{n}$ is \ki{upper triangular} if \ $a_{ij}=0$ for
all $i>j$. The matrix $A$ is \ki{lower triangular} if $a_{ij}=0$ for all $i<j$.
\end{defi}

For example the general form of the upper triangular matrices is the
following:
\[
\left[
\begin{array}
[c]{ccccc}%
\ast & \ast & \cdots & \cdots & \ast\\
0 & \ast &  &  & \vdots\\
\vdots & \ddots & \ddots &  & \vdots\\
\vdots &  & \ddots & \ast & \ast\\
0 & \cdots & \cdots & 0 & \ast
\end{array}
\right]  \koz .
\]
We note that the diagonal matrices are both lower and upper triangular. It is
easy to show that $\det(A)=a_{11}a_{22}\ldots a_{nn}$ holds for the upper or
lower triangular matrices. It is easy to solve linear systems with triangular
coefficient matrices. Consider the following upper triangular linear system:%
\[%
\begin{array}
[c]{ccccccc}%
a_{11}x_{1}+ & \cdots & +a_{1i}x_{i}+ & \cdots & +a_{1n}x_{n} & = & b_{1}\\
& \ddots & \vdots &  & \vdots &  & \vdots\\
&  & a_{ii}x_{i}+ & \cdots & +a_{in}x_{n} & = & b_{i}\\
&  &  & \ddots & \vdots &  & \vdots\\
&  &  &  & a_{nn}x_{n} & = & b_{n}%
\end{array}
\]
This can be solved by the so called \emph{back substitution} algorithm.

\begin{alg}{Back-Substitution($A,b,n$)}\inddef{\textsc{Back-Substitution}}
1 \' $x_{n}\leftarrow b_{n}/a_{nn}$ \\
2 \'  \key{for} \= $i\leftarrow n-1$ \= \key{downto} $1$ \\
3 \' \> \key{do} \=$x_{i}\leftarrow (b_{i}-\sum_{j=i+1}^{n}a_{ij}x_{j})/a_{ii}$ \\
4 \'  \key{return} $x$
\end{alg}


\noindent The solution of lower triangular systems is similar.

\textbf{The Gauss method.} The Gauss method \nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
or Gaussian elimination (GE) consists of two phases:

\noindent I. The linear system $Ax=b$ is transformed to an equivalent upper
triangular system using elementary operations (see Figure \ref{fig:Gauss}). %



\epskep{figs/Gauss}{Gaussian elimination.\label{fig:Gauss}}

\noindent II. The obtained upper triangular system is then solved by the back
substitution algorithm.

The first phase is often called the elimination or forward phase. The second
phase of GE is called the backward phase. The elementary operations are of the
following three types:

1. Add a multiple of one equation to another equation.

2. Interchange two equations.

3. Multiply an equation by a nonzero constant.

The elimination phase of GE is based on the following observation. Multiply
equation $k$ by $\gamma\neq0$ and subtract it from equation $i$:%
\[
\left(  a_{i1}-\gamma a_{k1}\right)  x_{1}+\cdots+\left(  a_{ij}-\gamma
a_{kj}\right)  x_{j}+\cdots+\left(  a_{in}-\gamma a_{kn}\right)  x_{n}%
=b_{i}-\gamma b_{k} \koz.
\]
If $a_{kj}\neq0$, then by choosing $\gamma=a_{ij}/a_{kj}$, the coefficient of
$x_{j}$ becomes $0$ in the new equivalent equation, which replaces equation
$i$. Thus we can eliminate variable $x_{j}$ (or coefficient $a_{ij}$) from
equation $i$.

The Gauss method eliminates the coefficients (variables) under the main
diagonal of $A$ in a systematic way. First variable $x_{1}$ is eliminated from
equations $i=2,\ldots,n$ using equation $1$, then $x_{2}$ is eliminated from
equations $i=3,\ldots,n$ using equation $2$, and so on.

Assume that the unknowns are eliminated in the first $\left(  k-1\right)  $
columns under the main diagonal and the resulting linear system has the form%

\[%
\begin{array}
[c]{cccccccccc}%
a_{11}x_{1}+ & \cdots & \cdots & +a_{1k}x_{k} & + & \cdots & + & a_{1n}x_{n} &
= & b_{1}\\
& \ddots &  & \vdots &  &  &  & \vdots &  & \vdots\\
&  & \ddots & \vdots &  &  &  & \vdots &  & \vdots\\
&  &  & a_{kk}x_{k} & + & \cdots & + & a_{kn}x_{n} & = & b_{k}\\
&  &  & \vdots &  &  &  & \vdots &  & \vdots\\
&  &  & a_{ik}x_{k} & + & \cdots & + & a_{in}x_{n} & = & b_{i}\\
&  &  & \vdots &  &  &  & \vdots &  & \vdots\\
&  &  & a_{nk}x_{k} & + & \cdots & + & a_{nn}x_{n} & = & b_{n}%
\end{array} \koz .
\]
If $a_{kk}\neq0$, then multiplying row $k$ by $\gamma$ and subtracting it from
equation $i$ we obtain%

\[
(a_{ik}-\gamma a_{kk})x_{k}+(a_{i,k+1}-\gamma a_{k,k+1})x_{k+1}+\cdots
+(a_{in}-\gamma a_{kn})x_{n}=b_{i}-\gamma b_{k} \koz .
\]
Since $a_{ik}-\gamma a_{kk}=0$ for $\gamma=a_{ik}/a_{kk}$, we eliminated the
coefficient $a_{ik}$ (variable $x_{k}$) from equation $i>k$. Repeating this
process for $i=k+1,\ldots,n$ we can eliminate the coefficients under the main
diagonal entry $a_{kk}$. Next we denote by $A\left[  i,j\right]  $ the element
$a_{ij}$ of matrix $A$ and by $A\left[  i,j:n\right]  $ the vector $\left[
a_{ij},a_{i,j+1},\ldots,a_{in}\right]  $. The Gauss method has the following
form (where the \emph{pivoting} discussed later is also included):

\begin{algN}{Gauss-Method($A,b$)}\inddef{\textsc{Gauss-Method}}
1 \'  \` $\rhd$ Forward phase: \\
2 \'  $n\leftarrow \textit{rows}[A]$ \\
3 \'  \key{for} \= $k\leftarrow 1$ \= \key{to} \= $n-1$ \\
4 \'   \> \key{do} \= \'  \{pivoting and interchange of rows and columns\} \\
5 \'  \> \>  \key{for} \= $i\leftarrow k+1$ \= \key{to} \= $n $ \\
6 \'  \>  \> \> \key{do} \= 
$\gamma _{ik}\leftarrow A\left[ i,k\right] /A\left[ k,k\right] $ \\
7 \'  \>  \> \> \>   $A\left[ i,k+1:n\right] \leftarrow
A\left[ i,k+1:n\right] -\gamma _{ik}\ast A\left[ k,k+1:n\right]$  \\
8 \'  \>  \>  \>  \>  $b_{i}\leftarrow b_{i}-\gamma _{ik} b_{k}$ \\
9 \'  \` $\rhd$ Backward phase: see the back substitution algorithm. \\
10 \' \key{return} $x$
\end{algN}

The algorithm overwrites the original matrix $A$ and vector $b$. It does not
write however the zero entries under the main diagonal since these elements
are not necessary for the second phase of the algorithm. Hence the lower
triangular part of matrix $A$ can be used to store information for the $LU$
decomposition of matrix $A$.

The above version of the Gauss method \nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
can be performed only if the elements
$a_{kk}$ occurring in the computation are not zero. For this and numerical
stability reasons we use the Gaussian elimination with \emph{pivoting}.

\subsubsection{The Gauss method with pivoting.}
If $a_{kk}=0$, then we can interchange row $k$ with another row, say $i$, so
that the new entry ($a_{ki}$) at position $\left(  k,k\right)  $ should be
nonzero. If this is not possible, then all the coefficients $a_{kk}%
,a_{k+1,k},\ldots,a_{nk}$ are zero and $\det\left(  A\right)  =0$. In the
latter case $Ax=b$ has no unique solution. The element $a_{kk}$ is called the
$k^{\text{th}}$ \ki{pivot element.}\inddef{pivot element} We can always select new pivot elements by
interchanging the rows. The selection of the pivot element has a great
influence on the reliability of the computed results. The simple fact that we
divide by the pivot element indicates this influence. We recall that
$\delta(a/b)$ is proportional to $1/\left|  b\right|  ^{2}$. It is considered
advantageous if the pivot element is selected so that it has the greatest
possible modulus. The process of selecting the pivot element is called
\emph{pivoting}. We mention the following two pivoting processes.

\textbf{Partial pivoting:} At the $k^{\text{th}}$ step, interchange the rows
of the matrix so the largest remaining element, say $a_{ik}$, in the
$k^{\text{th}}$ column is used as pivot. After the pivoting we have%
\[
\left|  a_{kk}\right|  =\max_{k\leq i\leq n}\left|  a_{ik}\right|  \koz .
\]


\textbf{Complete pivoting:} At the $k^{\text{th}}$ step, interchange both the
rows and columns of the matrix so that the largest element, say $a_{ij}$, in
the remaining matrix is used as pivot After the pivoting we have%
\[
\left|  a_{kk}\right|  =\max_{k\leq i,j\leq n}\left|  a_{ij}\right|  \koz .
\]


Note that the interchange of two columns implies the interchange of the
corresponding unknowns. The significance of pivoting is well illustrated by
the following

\begin{pelda}
The exact solution of the linear system%
\[%
\begin{array}
[c]{ccccc}%
10^{-17}x & + & y & = & 1\\
x & + & y & = & 2
\end{array}
\]
is $x=1/(1-10^{-17})$ and $y=1-10^{-17}/(1-10^{-17})$. The MATLAB program
gives the result $x=1$, $y=1$ and this is the best available result in
standard double precision arithmetic. Solving this system with the Gaussian
elimination without pivoting (also in double precision) we obtain the
catastrophic result $x=0$ and $y=1$. Using partial pivoting with the Gaussian
elimination we obtain the best available numerical result $x=y=1$.
\end{pelda}

\begin{megj}
Theoretically we do not need pivoting in the following cases: 1. If $A$ is
symmetric and positive definite ($A\in\mathbb{R}^{n\times n}$ is positive
definite $\Leftrightarrow$ $x^{T}Ax>0$, $\forall x\in\mathbb{R}^{n}$, $x\neq
0$). 2. If $A$ is diagonally dominant in the following sense:%
\[
\left|  a_{ii}\right|  >\sum_{j\neq i}\left|  a_{ij}\right|  \quad\left(
1\leq i\leq n\right)  \koz .
\]
In case of symmetric and positive definite matrices we use the Cholesky method
\nevindex{Cholesky, André-Louis (1875--1918)} which is a special version of the Gauss-type methods.
\end{megj}

During the Gaussian elimination we obtain a sequence of equivalent linear
systems%
\[
A^{\left(  0\right)  }x=b^{\left(  0\right)  }\rightarrow A^{\left(  1\right)
}x=b^{\left(  1\right)  }\rightarrow\cdots\rightarrow A^{\left(  n-1\right)
}x=b^{\left(  n-1\right)  } \koz ,
\]
where%
\[
A^{\left(  0\right)  }=A,\quad A^{\left(  k\right)  }=\left[  a_{ij}^{\left(
k\right)  }\right]  _{i,j=1}^{n}.
\]
Note that matrices $A^{\left(  k\right)  }$ are stored in the place of
$A=A^{\left(  0\right)  }$. The last coefficient matrix of phase I has the
form
\[
A^{\left(  n-1\right)  }=\left[
\begin{array}
[c]{cccc}%
a_{11}^{\left(  0\right)  } & a_{12}^{\left(  0\right)  } & \cdots &
a_{1n}^{\left(  0\right)  }\\
0 & a_{22}^{\left(  1\right)  } & \cdots & a_{2n}^{\left(  1\right)  }\\
\vdots &  & \ddots & \vdots\\
0 & \cdots & \cdots & a_{nn}^{\left(  n-1\right)  }%
\end{array}
\right]  \koz ,
\]
where $a_{kk}^{\left(  k-1\right)  }$ is the $k^{\text{th}}$ pivot element.
The \ki{growth factor of pivot elements}\inddef{growth factor of pivot elements} is given by
\[
\rho=\rho_{n}=\max_{1\leq k\leq n}\left|  a_{kk}^{\left(  k-1\right)  }%
/a_{11}^{\left(  0\right)  }\right|  \koz .
\]
Wilkinson\nevindex{Wilkinson, James H.} proved that the error of the computed solution is proportional to
the growth factor $\rho$ and the bounds%

\[
\rho\leq\sqrt{n}\left(  2\cdot3^{\frac{1}{2}}\cdots n^{\frac{1}{n-1}}\right)
^{\frac{1}{2}}\sim cn^{\frac{1}{2}}n^{\frac{1}{4}log\left(  n\right)  }%
\]
and
\[
\rho\leq2^{n-1}%
\]
hold for complete and partial pivoting, respectively. Wilkinson\nevindex{Wilkinson, James H.} conjectured
that $\rho\leq n$ for complete pivoting. This has been proved by researchers
for small values of $n$. Statistical investigations on random matrices
($n\leq1024$) indicate that the average of $\rho$ is $\Theta\left(
n^{2/3}\right)  $ for the partial pivoting and $\Theta\left(  n^{1/2}\right)
$ for the complete pivoting. Hence the case $\rho>n$ hardly occurs in the
statistical sense.

We remark that Wilkinson\nevindex{Wilkinson, James H.} constructed a linear system on which $\rho=2^{n-1}$
for the partial pivoting. Hence Wilkinson's\nevindex{Wilkinson, James H.} bound for $\rho$ is sharp in the
case of partial pivoting. There also exist examples of linear systems
concerning discretisations of differential and integral equations, where
$\rho$ is increasing exponentially if Gaussian elimination is used with
partial pivoting.

The growth factor $\rho$ can be very large, if the 
Gaussian\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} elimination is
used without pivoting. For example, $\rho=\rho_{4}\left(  A\right)
=1.23\times10^{5}$, if
\[
A=\left[
\begin{array}
[c]{rrrr}%
1.7846 & -0.2760 & -0.2760 & -0.2760\\
-3.3848 & 0.7240 & -0.3492 & -0.2760\\
-0.2760 & -0.2760 & 1.4311 & -0.2760\\
-0.2760 & -0.2760 & -0.2760 & 0.7240
\end{array}
\right]  \koz .
\]


\textbf{Operations counts.} The Gauss method\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
gives the solution of the linear system $Ax=b$ ($A\in
\mathbb{R}^{n\times n}$) in a finite number of steps and arithmetic operations
$(\,+,-,$ $\ast,/\,)$. The amount of necessary arithmetic operations is an
important characteristic of the direct linear system solvers, since the CPU
time is largely proportional to the number of arithmetic operations. It was
also observed that the number of additive and multiplicative operations are
nearly the same in the numerical algorithms of linear algebra. For measuring
the cost of such algorithms C. B. Moler\nevindex{Moler, Cleve B.} introduced the concept of flop.

\begin{defi}
\ki{One (old) flop} \inddef{One (old) flop} is the computational work necessary for the operation
$s=s+x\ast y$ (1 addition + 1 multiplication). \ki{One (new)}\inddef{One (new) flop} flop is the
computational work necessary for any of the arithmetic operations $+,-,\ast,/$.
\end{defi}

The new flop can be used if the computational time of additive and
multiplicative operations are approximately the same. Two new flops equals to
one old flop. Here we use the notion of old flop.

For the Gauss method\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
a simple counting gives the number of additive and
multiplicative operations.

\begin{tetel}
\label{muvig} The computational cost of the Gauss 
method\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} is $n^{3}%
/3+\Theta(n^{2})$ flops.
\end{tetel}

V. V. Klyuyev\nevindex{Klyuyev, v. v.} and N. Kokovkin-Shcherbak\nevindex{Kokovkin-Shcherbak, N.}
proved that if only elementary row and
column operations (multiplication of row or column by a number, interchange of
rows or columns, addition of a multiple of row or column to another row or
column) are allowed, then the linear system $Ax=b$ cannot be solved in less
than $n^{3}/3+\Omega(n^{2})$ flops.

Using fast matrix inversion procedures we can solve the $n\times n$ linear
system $Ax=b$ in $O(n^{2.808})$ flops. These theoretically interesting
algorithms are not used in practice since they are considered as numerically unstable.

\textbf{The $LU$-decomposition.}
In many cases it is easier to solve a linear system if the coefficient matrix
can be decomposed into the product of two triangular matrices.

\begin{defi}
The matrix $A\in\mathbb{R}^{n\times n}$ has an 
\ki{$LU$-decomposition,}\inddef{LU-decomposition@$LU$-decomposition} if $A=LU$,
where $L\in\mathbb{R}^{n\times n}$ is lower and $U\in\mathbb{R}^{n\times n}$
is upper triangular matrix.
\end{defi}

The $LU$-decomposition is not unique. If a nonsingular matrix has an
$LU$-decomposition, then it has a particular $LU$-decomposition, where the
main diagonal of a given component matrix consists of $1$'s. Such triangular
matrices are called \ki{unit (upper or lower) triangular}\inddef{unit (upper or lower) triangular} matrices. The $LU$
decomposition is unique, if $L$ is set to be lower unit triangular or $U$ is
set to be unit upper triangular.

The $LU$-decomposition of nonsingular matrices is closely related to the
Gaussian\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
elimination method. If $A=LU$, where $L$ is unit lower triangular,
then $l_{ik}=\gamma_{ik}$ ($i>k$), where $\gamma_{ik}$ is given by the Gauss 
algorithm. The matrix $U$ is the upper triangular part of the matrix we obtain
at the end of the forward phase. The matrix $L$ can also be derived from this
matrix, if the columns of the lower triangular part are divided by the
corresponding main diagonal elements. We remind that the first phase of the
Gaussian elimination does not annihilate the matrix elements under the main
diagonal. It is clear that a nonsingular matrix \ has $LU$-decomposition if
and only if $a_{kk}^{\left(  k-1\right)  }\neq0$ holds for each pivot element
for the Gauss method without pivoting.

\begin{defi}
A matrix $P\in\mathbb{R}^{n\times n}$ whose every row and column has one and
only one non-zero element, that element being $1$, is called a \ki{permutation matrix.}\inddef{permutation matrix} 
\end{defi}

In case of partial pivoting we permute the rows of the coefficient matrix
(multiply $A$ by a permutation matrix on the left) so that $a_{kk}^{\left(
k-1\right)  }\neq0$ $(k=1,\ldots,n)$ holds for a nonsingular matrix. Hence we have

\begin{tetel}
If $A\in\mathbb{R}^{n\times n}$ is nonsingular then there exists a permutation
matrix $P$ such that $PA$ has an $LU$-decomposition.
\end{tetel}

The the algorithm of $LU$-decomposition is essentially the 
Gaussian\nevindex{Gauss, Johann Carl Friedrich (1777--1855)}
elimination method. If pivoting is used then the interchange of rows must also
be executed on the elements under the main diagonal and the permutation matrix
$P$ must be recorded. A vector containing the actual order of the original
matrix rows is obviously sufficient for this purpose.

\textbf{The $LU$- and Cholesky-methods.}\nevindex{Cholesky, André-Louis (1875--1918)}
Let $A=LU$ and consider the equation $Ax=b$. Since $Ax=LUx=L(Ux)=b$, we can
decompose $Ax=b$ into the equivalent linear system $Ly=b$ and $Ux=b$, where
$L$ is lower triangular and $U$ is upper triangular.

\begin{alg}{\textit{LU}-Method$(A,b)$}\inddef{\textsc{LU-Method}}
1 \' Determine the $LU$-decomposition $A=LU$.\\
2 \' Solve $Ly=b$. \\
3 \' Solve $Ux=y$. \\
4 \' \key{return} $x$
\end{alg}


\textit{Remark.} In case of partial pivoting we obtain the decomposition
$\hat{A}=PA=LU$ and we set $\hat{b}=Pb$ instead of $b$.

In the first phase of the Gauss method\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} 
we produce decomposition $A=LU$ and the
equivalent linear system $Ux=L^{-1}b$ with upper triangular coefficient
matrix. The latter is solved in the second phase. In the $LU$-method we
decompose the first phase of the Gauss method into two steps. In the first
step we obtain only the decomposition $A=LU$. In the second step we produce
the vector $y=L^{-1}b$. The third step of the algorithm is identical with the
second phase of the original Gauss method.

The $LU$-method is especially advantageous if we have to solve several linear
systems with the same coefficient matrix:
\[
Ax=b_{1},\ Ax=b_{2},\ldots,\ Ax=b_{k}\text{.}%
\]
In such a case we determine the $LU$-decomposition of matrix $A$ only once,
and then we solve the linear systems $Ly_{i}=b_{i}$, $Ux_{i}=y_{i}$
($x_{i},y_{i,}b_{i}\in\mathbf{R}^{n},$ $i=1,\ldots,k$). The computational cost
of this process is $n^{3}/3+kn^{2}+\Theta\left(  kn\right)  $ flops.

The \emph{inversion of a matrix}\index{inversion of a matrix} $A\in\mathbb{R}^{n\times n}$ can be done as follows:

1. Determine the $LU$-decomposition $A=LU$. .

2. Solve $Ly_{i}=e_{i}$, $Ux_{i}=y_{i}$ ($e_{i}$ is the $i^{\text{th}}$ unit
vector $i=1,\ldots,n$).

\noindent The inverse of $A$ is given by $A^{-1}=\left[  x_{1},\ldots
,x_{n}\right]  $. The computational cost of the algorithm is $4n^{3}%
/3+\Theta\left(  n^{2}\right)  $ flops.

\subsubsection{The \emph{LU}-method with pointers.}
This implementation of the $LU$-method is known since the 60's. Vector $P$
contains the indices of the rows. At the start we set $P\left[  i\right]  =i$
($1\leq i\leq n$). When exchanging rows we exchange only those components of
vector $P$ that correspond to the rows.

\inddef{LU-method-with-Pointers@\textsc{LU-Method-with-Pointers}}
\begin{algN}{\textsc{LU-Method-with-Pointers}($A,b$)}
1 \' $n\leftarrow \textit{rows}[A]$ \\
2 \' $P\leftarrow [1,2,\ldots,n]$ \\
3 \'  \key{for} \= $k\leftarrow  1$ \=  \key{to} \= $n-1$  \\
4 \'  \>  \key{do} \= 
compute index $t$ such that $\left| A\left[ P\left[ t\right]
,k\right] \right| = \max_{k\leq i\leq n}\left| A\left[ P\left[
i\right] ,k\right] \right| $ . \\
5 \'  \> \>  \key{if} \= $k<t$ \\
6 \' \>  \>  \>   \key{then} \= 
exchange the components $P\left[ k\right] $ and $P\left[ t\right] $. \\
7 \'  \>  \>   \key{for} \= $i\leftarrow k+1$ \= \key{to} \= $n$ \\
8 \'  \> \>  \>   \key{do} \=  $A\left[ P\left[ i\right] ,k\right] 
  \leftarrow A\left[ P\left[ i\right], k\right] /A\left[ P\left[ k\right],k\right]$ \\
9 \'  \> \> \> \> $A\left[ P\left[ i\right] ,k+1:n\right]$ \\ 
  \'  \> \> \> \> $\leftarrow A\left[ P\left[ i\right],k+1:n\right] -A\left[ P\left[ i\right] 
,k\right] \ast A\left[ P\left[ k\right] ,k+1:n\right] $ 
\algnewpage
10 \'    \key{for} \=  $i\leftarrow 1$  \= \key{to} \= $n$ \\
11 \'  \>  \key{do} \= $s\leftarrow 0$ \\
12 \'  \> \>  \key{for} \= $j\leftarrow 1$ \=  \key{to} \= $i-1$ \\
13 \'  \> \> \>  \key{do} \=  $s\leftarrow
s+A\left[ P\left[ i\right] ,j\right] \ast x\left[ j\right] $ \\
14 \'  \> \> $x\left[ i\right] \leftarrow b[P[i]]-s$ \\
15 \'  \key{for} \= $i\leftarrow n$ \=  \key{downto} $1$ \\
16 \'  \>  \key{do} \= $s\leftarrow 0$ \\
17 \'  \> \>  \key{for} \= $j\leftarrow i+1$  \key{to} \= $n$ \\
18 \'  \> \> \> $s\leftarrow s+A\left[ P\left[ i\right] ,j\right] 
\ast x\left[ j\right] $ \\
19 \'  \> \> $x\left[ i\right] \leftarrow 
\left( x\left[ i\right] -s\right) /A\left[
P\left[ i\right] ,i\right] $ \\
20 \'  \key{return} $x$
\end{algN}

If $A\in\mathbb{R}^{n\times n}$ is \emph{symmetric and positive definite}, then it
can be decomposed in the form $A=LL^{T}$, where $L$ is lower triangular
matrix. The $LL^{T}$-decomposition is called the \ki{Cholesky-decomposition.}\inddef{Cholesky-decomposition}
\nevindex{Cholesky, André-Louis (1875--1918)} In
this case we can save approximately half of the storage place for $A$ and half
of the computational cost of the $LU$-decomposition ($LL^{T}$-decomposition).
Let
\[
A=\left[
\begin{array}
[c]{ccc}%
a_{11} & \cdots & a_{1n}\\
a_{21} & \cdots & a_{2n}\\
\vdots &  & \vdots\\
a_{n1} & \cdots & a_{nn}%
\end{array}
\right]  =\left[
\begin{array}
[c]{cccc}%
l_{11} & 0 & \cdots & 0\\
l_{21} & l_{22} & \ddots & \vdots\\
\vdots & \vdots & \ddots & 0\\
l_{n1} & l_{n2} & \cdots & l_{nn}%
\end{array}
\right]  \left[
\begin{array}
[c]{cccc}%
l_{11} & l_{21} & \cdots & l_{n1}\\
0 & l_{22} & \cdots & l_{n2}\\
\vdots & \ddots & \ddots & \vdots\\
0 & \cdots & 0 & l_{nn}%
\end{array}
\right]  .
\]
Observing that only the first $k$ elements may be nonzero in the
$k^{\text{th}}$ column of $L^{T}$ we obtain that%
\begin{align*}
a_{kk}  &  =l_{k1}^{2}+l_{k2}^{2}+\cdots+l_{k,k-1}^{2}+l_{kk}^{2},\\
a_{ik}  &  =l_{i1}l_{k1}+l_{i2}l_{k2}+\cdots+l_{i,k-1}l_{k,k-1}+l_{ik}%
l_{kk}\quad\left(  i=k+1,\ldots,n\right)  \koz .
\end{align*}
This gives the formulae%
\begin{align*}
l_{kk}  &  =(a_{kk}-\sum_{j=1}^{k-1}l_{kj}^{2})^{1/2},\\
l_{ik}  &  =(a_{ik}-\sum_{j=1}^{k-1}l_{ij}l_{kj})/l_{kk}\quad\left(
i=k+1,\ldots,n\right)  \koz .
\end{align*}

Using the notation $\sum_{j=i}^{k}s_{j}=0$ ($k<i$) we can formulate the
Cholesky-method \nevindex{Cholesky, André-Louis (1875--1918)}as follows.

\begin{alg}{Cholesky-Method($A$)}\inddef{\textsc{Cholesky-Method}}
1 \' $n\leftarrow \textit{rows}[A]$ \\
2 \'  \key{for} \= $k\leftarrow 1$ \= \key{ to }\= $n$ \\
3 \' \> \key{do } \= $a_{kk}\leftarrow 
(a_{kk}-\sum_{j=1}^{k-1}a_{kj}^{2})^{1/2}$  \\
4 \' \> \> \key{for} \= $i\leftarrow k+1$ \= \key{ to } \= $n$ \\
5 \' \> \> \> \key{do } \= $a_{ik}\leftarrow 
(a_{ik}-\sum_{j=1}^{k-1}a_{ij}a_{kj})/a_{kk}$ \\
6 \'  \key{return } $A$ 
\end{alg}

The lower triangular part of $A$ contains $L$. The computational cost of the
algorithm is $n^{3}/6+\Theta(n^{2})$ flops and $n$ square roots. The
algorithm, which can be considered as a special case of the 
Gauss-methods\nevindex{Gauss, Johann Carl Friedrich (1777--1855)},
does not require pivoting, at least in principle.

\subsubsection{The \emph{LU}- and Cholesky-methods on banded matrices.}\nevindex{Cholesky, André-Louis (1875--1918)}
It often happens that linear systems have \emph{banded coefficient matrices}\index{banded matrix}.

\begin{defi}
Matrix $A\in\mathbb{R}^{n\times n}$ is banded with lower bandwidth $p$ and
upper bandwidth $q$ if
\[
a_{ij}=0,\quad\text{if }i>j+p\text{ or }j>i+q \koz .
\]
\end{defi}

The possibly non-zero elements $a_{ij}$ ($i-p\leq j\leq i+q$) form a band like
structure. Schematically $A$ has the form%
\[
A=\left[
\begin{array}
[c]{ccccccccc}%
a_{11} & a_{12} & \cdots & \cdots & a_{1,1+q} & 0 & \cdots & \cdots & 0\\
a_{21} & a_{22} &  &  &  & \ddots &  &  & \vdots\\
\vdots &  & \ddots &  &  &  & \ddots &  & \vdots\\
a_{1+p,1} &  &  & \ddots &  &  &  & \ddots & 0\\
0 & \ddots &  &  & \ddots &  &  &  & a_{n-q,n}\\
\vdots &  & \ddots &  &  & \ddots &  &  & \vdots\\
\vdots &  &  & \ddots &  &  & \ddots &  & \vdots\\
\vdots &  &  &  & \ddots &  &  & \ddots & a_{n-1,n}\\
0 & \cdots & \cdots & \cdots & 0 & a_{n,n-p} & \cdots & a_{n,n-1} & a_{nn}%
\end{array}
\right]  \koz .
\]
The \emph{banded matrices}\index{banded matrix} yield very efficient algorithms if $p$ and $q$ are
significantly less than $n$. If a banded matrix $A$ with lower bandwidth $p$
and upper bandwidth $q$ has an $LU$-decomposition, then both $L$ and $U$ are
banded with lower bandwidth $p$ and upper bandwidth $q$, respectively.

Next we give the $LU$-method for banded matrices in three parts.

\begin{alg}{The-LU-Decomposition-of-Banded-Matrix($A,n,p,q$)}
\inddef{The-LU-Decomposition-of-Banded-Matrix@\textsc{The-LU-Decomposition-of-Banded-Matrix}}
1 \' \key{for } \= $k\leftarrow 1$ \= \key{ to } \= $n-1$ \\
2 \' \> \key{do} \= \key{for} \= $i\leftarrow k+1$ \= \key{ to } \=
$\min \left\{ k+p,n\right\}$  \\
3 \' \> \> \> \key{do} \= $a_{ik}\leftarrow a_{ik}/a_{kk}$ \\
4 \' \> \> \> \> \key{for} \= $j\leftarrow k+1$ \= \key{ to } \=
$\min \left\{ k+q,n\right\}$  \\
5 \'  \> \> \> \> \> \key{do} \= $a_{ij}\leftarrow a_{ij}-a_{ik}a_{kj}$ \\
6 \'  \key{return} $A$
\end{alg}


Entry $a_{ij}$ is overwritten by $l_{ij}$, if $i>j$ and by $u_{ij}$, if $i\leq
j$. The computational cost of \ is $c\left(  p,q\right)  $ flops, where
\[
c\left(  p,q\right)  =\left\{
\begin{array}
[c]{l}%
npq-\frac{1}{2}pq^{2}-\frac{1}{6}p^{3}+pn,\quad p\leq q\\
npq-\frac{1}{2}qp^{2}-\frac{1}{6}q^{3}+qn,\quad p>q
\end{array} \koz .
\right.
\]


The following algorithm overwrites $b$ by the solution of equation $Ly=b$.

\begin{alg}{Solution-of-Banded-Unit-Lower-Triangular-System$(L,b,n,p)$}
\inddef{\textsc{Solution-of-Banded-Unit-Lower-Triangular-System}}
1 \' \key{for} \= $i\leftarrow 1$ \= \key{ to } \= $n$ \\
2 \' \> \key{do } \= $b_{i}\leftarrow b_{i}-
\sum_{j=\max \left\{ 1,i-p\right\}}^{i-1}l_{ij}b_{j}$ \\
3 \'  \key{return} $b$
\end{alg}


The total cost of the algorithm is $np-p^{2}/2$ flops. The next algorithm
overwrites vector $b$ by the solution of $Ux=b$.

\begin{alg}{Solution-of-Banded-Upper-Triangular-System$(U,b,n,q)$}
\inddef{\textsc{Solution-of-Banded-Upper-Triangular-System}}
1 \' \key{for} \= $i\leftarrow n$ \= \key{downto} \= $1$ \\
2 \' \> \key{do } \= 
$b_{i}\leftarrow \left( b_{i}-\sum_{j=i+1}^{\min \left\{
i+q,n\right\} }u_{ij}b_{j}\right) /u_{ii}$ \\
3 \' \key{return} $b$
\end{alg}


The computational cost is $n\left(  q+1\right)  -q^{2}/2$ flops.

Assume that $A\in\mathbb{R}^{n\times n}$ is symmetric, positive definite and
banded with lower bandwidth $p$. The banded version of the 
Cholesky-methods\nevindex{Cholesky, André-Louis (1875--1918)} is
given by

\begin{alg}{Cholesky-decomposition-of-Banded-Matrices($A,n,p$)}
\inddef{\textsc{Cholesky-Decomposition-of-Banded-Matrices}}
1 \' \key{for} \= $i\leftarrow 1$ \= \key{ to } \= $n$ \\
2 \'  \> \key{do} \= 
\key{for } \= $j\leftarrow \max \left\{ 1,i-p\right\}$ \=
\key{ to } \= $i-1$ \\
3 \'  \> \> \> \key{do} \= $a_{ij}\leftarrow
\left( a_{ij}-\sum_{k=\max \left\{ 1,i-p\right\}
}^{j-1}a_{ik}a_{jk}\right) /a_{jj}$  \\
4 \'  \>  \> $a_{ii}\leftarrow 
\left( a_{ii}-\sum_{k=\max \left\{ 1,i-p\right\}
}^{i-1}a_{ik}^{2}\right) ^{1/2}$ \\
5 \'  \key{return} $A$
\end{alg}


The elements $a_{ij}$ are overwritten by $l_{ij}$ ($i\geq j$). The total
amount of work is given by $\left(  np^{2}/2\right)  -\left(  p^{3}/3\right)
+\left(  3/2\right)  \left(  np-p^{2}\right)  $ flops \'{e}s $n$ square roots.

\textit{Remark.} If $A\in\mathbb{R}^{n\times n}$ has lower bandwidth $p$ and
upper bandwidth $q$ and partial pivoting takes place, then the upper bandwidth
of $U$ increases up to $\hat{q}=p+q$.

\subsection{Iterative methods for linear systems}
There are several iterative methods for solving linear systems of algebraic
equations. The best known iterative algorithms are the classical 
\emph{Jacobi}-,\nevindex{Jacobi, Carl Gustav Jacob (1804--1851)} the
\emph{Gauss-Seidel}\nevindex{Gauss, Johann Carl 
Friedrich (1777--1855)}\nevindex{Seidel, Philipp Ludwig, von (1821--1896)}- and the 
\emph{relaxation methods}. The greatest advantage of these
iterative algorithms is their easy implementation to large systems. At the
same time they usually have slow convergence. However for parallel computers
the \emph{multisplitting} iterative algorithms seem to be efficient.

Consider the iteration
\[
x_{i}=Gx_{i-1}+b\quad\left(  i=1,2,\ldots\right)
\]
where $G\in\mathbb{R}^{n\times n}$ \'{e}s $x_{0},b\in\mathbb{R}^{n}$. It is
known that $\left\{  x_{i}\right\}  _{i=0}^{\infty}$ converges for all $x_{0}%
$, $b\in\mathbb{R}^{n}$ if and only if the spectral radius of $G$ satisfies
$\rho\left(  G\right)  <1$ ($\rho\left(  G\right)  =\max\left|  \lambda
\right|  \mid\lambda$ is an eigenvalue of $G$). In case of convergence
$x_{i}\rightarrow x^{\ast}=\left(  I-G\right)  ^{-1}b$, that is we obtain the
solution of the equation $\left(  I-G\right)  x=b$. The speed of convergence
depends on the spectral radius $\rho\left(  G\right)  $. Smaller the spectral
radius$\rho\left(  G\right)  $, faster the convergence.

Consider now the linear system
\[
Ax=b \koz ,
\]
where $A\in\mathbb{R}^{n\times n}$ is nonsingular. The matrices $M_{l}%
,N_{l},E_{l}\in\mathbf{R}^{n\times n}$ form a multisplitting of $A$ {if}

(i) $A=M_{i}-N_{i}$, $i=1,2,\ldots,L,$

(ii) $M_{i}$ is nonsingular, $i=1,2,\ldots,L,$

(iii) $E_{i}$ is non-negative diagonal matrix, $i=1,2,\ldots,L,$

(iv) $\sum_{i=1}^{L}E_{i}=I$.

\noindent Let $x_{0}\in\mathbb{R}^{n}$ be a given initial vector. The
multisplitting iterative method is the following.

\begin{alg}{Multisplitting-Iteration($x_{0},b,L,M_{l},N_{l},E_{l},
l=1,\ldots ,L$)}\inddef{\textsc{Multisplitting-Iteration}}
1 \'  $i\leftarrow 0$ \\
2 \' \key{while} \= \emph{exit condition} = \textsc{false} \\
3 \' \>\key{do} \= $i\leftarrow i+1$ \\
4 \' \> \> \key{for} \= $l\leftarrow 1$ \key{to} $L$ \\
5 \' \> \> \> \key{do} \= $M_{l}y_{l}\leftarrow N_{l}x_{i-1}+b$ \\
6 \' \> \> $x_{i}\leftarrow \sum_{l=1}^{L}E_{l}y_{l}$ \\
7 \' \key{return} $x_{i}$
\end{alg}

It is easy to show that $y_{l}=M_{l}^{-1}N_{l}x_{i-1}+M_{l}^{-1}b$ and
\begin{align*}
x_{i}  &  =\sum_{l=1}^{L}E_{l}y_{l}=\sum_{l=1}^{L}E_{l}M_{l}^{-1}N_{l}%
x_{i-1}+\sum_{l=1}^{L}E_{l}M_{l}^{-1}b\\
&  =Hx_{i-1}+c \koz .
\end{align*}
Thus the condition of convergence is $\rho\left(  H\right)  <1$. The
\emph{multisplitting iteration} is a true \emph{parallel algorithm} because we
can solve $L$ linear systems parallel in each iteration (synchronised
parallelism). The bottleneck of the algorithm is the computation of iterate
$x_{i}$.

The selection of matrices $M_{i}$ and $E_{i}$ is such that the solution of the
linear system $M_{i}y=c$ should be cheap. Let $S_{1},S_{2},\ldots,S_{L}$ be a
partition of $\left\{  1,\ldots,n\right\}  $, that is $S_{i}\neq\emptyset$,
$S_{i}\cap S_{j}=\emptyset$ ($i\neq j$) and $\cup_{i=1}^{L}S_{i}=\left\{
1,\ldots,n\right\}  $. Furthermore let $S_{i}\subseteq T_{i}\subseteq\left\{
1,\ldots,n\right\}  $ ($i=1,\ldots,L$) be such that $S_{l}\neq T_{l}$ for at
least one $l$.

The \ki{non-overlapping block Jacobi splitting}\inddef{non-overlapping block Jacobi splitting}
\nevindex{Jacobi, Carl Gustav Jacob (1804--1851)} of $A$ is given by
\begin{align*}
M_{l}  &  =\left[  M_{ij}^{\left(  l\right)  }\right]  _{i,j=1}^{n},\quad
M_{ij}^{\left(  l\right)  }=\left\{
\begin{array}
[c]{l}%
a_{ij},\quad\text{if }i,j\in S_{l} \koz,\\
a_{ii},\quad\text{if }i=j \koz,\\
0,\quad\text{otherwise \koz,}%
\end{array}
\right. \\
& \\
N_{l}  &  =M_{l}-A,\\
& \\
E_{l}  &  =\left[  E_{ij}^{\left(  l\right)  }\right]  _{i,j=1}^{n},\quad
E_{ij}^{\left(  l\right)  }=\left\{
\begin{array}
[c]{l}%
1,\quad\text{if }i=j\in S_{l}\,\\
0,\quad\text{otherwise}%
\end{array}
\right.
\end{align*}
for $l=1,\ldots,L$.

Define now the simple splitting%
\[
A=M-N \koz ,
\]
where $M$ is nonsingular,
\[
M=\left[  M_{ij}\right]  _{i,j=1}^{n},\quad M_{ij}=\left\{
\begin{array}
[c]{l}%
a_{ij},\quad\text{if }i,j\in S_{l}\text{ for some }l\in\left\{  1,\ldots
,n\right\}\koz, \\
0,\quad\text{otherwise.}%
\end{array}
\right.  
\]
It can be shown that
\[
H=\sum_{l=1}^{L}E_{l}M_{l}^{-1}N_{l}=M^{-1}N
\]
holds for the non-overlapping block Jacobi multisplitting.

The \ki{overlapping block Jacobi multisplitting}\inddef{overlapping block Jacobi multisplitting}
 \nevindex{Jacobi, Carl Gustav Jacob (1804--1851)} of $A$ \ is defined by
\begin{align*}
\widetilde{M}_{l}  &  =\left[  \widetilde{M}_{ij}^{\left(  l\right)  }\right]
_{i,j=1}^{n},\quad\widetilde{M}_{ij}^{\left(  l\right)  }=\left\{
\begin{array}
[c]{l}%
a_{ij},\quad\text{if }i,j\in T_{l}\\
a_{ii},\quad\text{if }i=j\\
0,\quad\text{otherwise}%
\end{array}
\right.  ,\\
& \\
\widetilde{N}_{l}  &  =\widetilde{M}_{l}-A,\\
& \\
\widetilde{e}_{l}  &  =\left[  \widetilde{e}_{ij}^{\left(  l\right)  }\right]
_{i,j=1}^{n},\quad E_{ii}^{\left(  l\right)  }=0\text{, if }i\notin T_{l}%
\end{align*}
for $l=1,\ldots,L$.

A nonsingular matrix $A\in\mathbb{R}^{n\times n}$ is called an $M$-matrix, if
$a_{ij}\leq0$ ($i\neq j$) and all the elements of $A^{-1}$ are nonnegative.

\begin{tetel}
Assume that $A\in\mathbb{R}^{n\times n}$ is nonsingular $M$-matrix, $\left\{
M_{i},N_{i},E_{i}\right\}  _{i=1}^{L}$\emph{\ }is a non-overlapping, $\left\{
\widetilde{M}_{i},\widetilde{N}_{i},E_{i}\right\}  _{i=1}^{L}$ is an
overlapping block Jacobi multisplitting of \ $A$, where the weighting matrices
$E_{i}$ are the same. The we have%
\[
\rho\left(  \widetilde{H}\right)  \leq\rho\left(  H\right)  <1 \koz ,
\]
where\emph{ }$H=\sum_{l=1}^{L}E_{l}M_{l}^{-1}N_{l}$\emph{\ }and\emph{
}$\widetilde{H}=\sum_{l=1}^{L}E_{l}\widetilde{M}_{l}^{-1}\widetilde{N}_{l}$
\emph{. }
\end{tetel}

We can observe that both iteration procedures are convergent and the
convergence of the overlapping multisplitting is not slower than that of the
non-overlapping procedure. The theorem remains true if we use block
Gauss-Seidel multisplittings instead of the block Jacobi multisplittings. In
this case we replace the above defined matrices $M_{i}$ and $\widetilde{M}%
_{i}$ with their lower triangular parts.

The multisplitting algorithm has multi-stage and asynchronous variants as well.

\subsection{Error analysis of linear algebraic systems}
We analyse the direct and inverse errors. We use the following notations and
concepts. The exact (theoretical) solution of $Ax=b$ is denoted by $x$, while
any approximate solution is denoted by $\hat{x}$. The direct error of the
approximate solution is given by $\Delta x=\hat{x}-x$. The quantity
$r=r\left(  y\right)  =Ay-b$ is called the \ki{residual error.}\inddef{residual error} For the exact
solution $r\left(  x\right)  =0$, while for the approximate solution
\[
r\left(  \hat{x}\right)  =A\hat{x}-b=A\left(  \hat{x}-x\right)  =A\Delta x \koz .
\]
We use various models to estimate the inverse error. In the most general case
we assume that the computed solution $\hat{x}$ satisfies the linear
system$\hat{A}\hat{x}=\hat{b}$, where $\hat{A}=A+\Delta A$ and $\hat
{b}=b+\Delta b$. The quantities $\Delta A$ and $\Delta b$ are called the
inverse errors.

One has to distinguish between the sensitivity of the problem and the
stability of the solution algorithm. By \ki{sensitivity of a problem}\inddef{sensitivity of a problem} we mean the
sensitivity of the solution to changes in the input parameters (data). By the
\ki{stability (or sensitivity) of an algorithm}\inddef{stability (or sensitivity) of an algorithm}
we mean the influence of
computational errors on the computed solution. We measure the sensitivity of a
problem or algorithm in various ways. One such characterization is the
\emph{condition number}\idez{condition number}, which compares the relative errors of the input and output values.

The following general principles are used when applying any algorithm:

\qquad- We use only stable or well-conditioned algorithms.

\qquad- We cannot solve an unstable (ill-posed or ill-conditioned) problem
with a general purpose algorithm, in general.

\subsubsection{Sensitivity analysis.}
Assume that we solve the perturbed equation
\begin{equation}
A\hat{x}=b+\Delta b
\end{equation}
instead of the original $Ax=b$. Let $\widehat{x}=x+\Delta x$ and investigate
the difference of the two solutions.

\begin{tetel}
\label{amodell} If $A$ is nonsingular and $b\neq0$, then
\begin{equation}
\frac{\left\|  \Delta x\right\|  }{\left\|  x\right\|  }\leq\text{\emph{cond}%
}(A)\frac{\left\|  \Delta b\right\|  }{\left\|  b\right\|  }=\text{\emph{cond}%
}(A)\frac{\left\|  r\left(  \hat{x}\right)  \right\|  }{\left\|  b\right\|  } \koz ,
\end{equation}
where \emph{cond}$(A)=\left\|  A\right\|  \left\|  A^{-1}\right\|  $\ is the
condition number of $A$.
\end{tetel}

Here we can see that the condition number of $A$ may strongly influence the
relative error of the perturbed solution $\widehat{x}$. A linear algebraic
system is said to be \emph{well-conditioned} if cond$(A)$ is small, and
\emph{ill-conditioned}, if cond$(A)$ is big. \ It is clear that the terms \idez{small}
and \idez{big} are relative and the condition number depends on the norm chosen.
We identify the applied norm if it is essential for some reason. For example
cond$_{\infty}\left(  A\right)  =\left\|  A\right\|  _{\infty}\left\|
A^{-1}\right\|  _{\infty}$. The next example gives possible geometric
characterization of the condition number.

\begin{pelda}
The linear system
\[%
\begin{array}
[c]{ccccc}%
1000x_{1} & + & 999x_{2} & = & b_{1}\\
999x_{1} & + & 998x_{2} & = & b_{2}%
\end{array}
\]
is ill-conditioned (cond$_{\infty}(A)=3.99\times10^{6}$). The two lines, whose
meshpoint defines the system, are almost parallel. Therefore if we perturb the
right hand side, the new meshpoint of the two lines will be far from the
previous meshpoint.
\end{pelda}

The inverse error is $\Delta b$ in the sensitivity model under investigation.
Theorem \ref{amodell} gives an estimate of the direct error which conforms
with the thumb rule. It follows that we can expect a small relative error of
\ the perturbed solution $\widehat{x}$, if the condition number of $A$ is small.

\begin{pelda}
Consider the linear system $Ax=b$ with
\[
A=\left[
\begin{array}
[c]{cc}%
1+\epsilon & 1\\
1 & 1
\end{array}
\right]  ,\quad b=\left[
\begin{array}
[c]{c}%
1\\
1
\end{array}
\right]  ,\quad x=\left[
\begin{array}
[c]{c}%
0\\
1
\end{array}
\right]  .
\]
Let $\hat{x}=\left[
\begin{array}
[c]{r}%
2\\
-1
\end{array}
\right]  $. Then $r=\left[
\begin{array}
[c]{c}%
2\epsilon\\
0
\end{array}
\right]  $ and $\left\|  r\right\|  _{\infty}/\left\|  b\right\|  _{\infty
}=2\epsilon$, but $\left\|  \hat{x}-x\right\|  _{\infty}/\left\|  x\right\|
_{\infty}=2$.
\end{pelda}

Consider now the perturbed linear system
\begin{equation}
\left(  A+\Delta A\right)  \hat{x}=b
\end{equation}
instead of $Ax=b$. It can be proved that for this perturbation model there
exist more than one \emph{inverse errors}\idez{inverse error} among which $\Delta A=-r\left(  \widehat
{x}\right)  \widehat{x}^{T}/\widehat{x}^{T}\widehat{x}$ is the inverse error
with minimal spectral norm, provided that $\widehat{x},r\left(  \widehat
{x}\right)  \neq0$.

The following theorem establish that for small relative residual error the
relative inverse error is also small.

\begin{tetel}
Assume that $\hat{x}\neq0$\ is the approximate solution of $Ax=b$,
$\det\left(  A\right)  \neq0$\ and $b\neq0$. If $\left\|  r\left(  \hat
{x}\right)  \right\|  _{2}/\left\|  b\right\|  _{2}=\alpha<1$, the the matrix
$\Delta A=-r\left(  \hat{x}\right)  \hat{x}^{T}/\hat{x}^{T}\hat{x}$\ satisfies
$\left(  A+\Delta A\right)  \hat{x}=b$\ and $\left\|  \Delta A\right\|
_{2}/\left\|  A\right\|  _{2}\leq\alpha/\left(  1-\alpha\right) \koz .$
\end{tetel}

If the relative inverse error and the condition number of $A$ are small, then
the relative residual error is small.

\begin{tetel}
\label{bmodell}If $r\left(  \hat{x}\right)  $=$A\hat{x}-b$, $\left(  A+\Delta
A\right)  \hat{x}$=$b$, $A\neq0$, $b\neq0$\ and $\text{\emph{cond}}\left(
A\right)  \frac{\left\|  \Delta A\right\|  }{\left\|  A\right\|  }$
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
$1$, then
\begin{equation}
\frac{\left\|  r\left(  \hat{x}\right)  \right\|  }{\left\|  b\right\|  }%
\leq\frac{\text{\emph{cond}}(A)\frac{\left\|  \Delta A\right\|  }{\left\|
A\right\|  }}{1-\text{\emph{cond}}(A)\frac{\left\|  \Delta A\right\|
}{\left\|  A\right\|  }} \koz .
\end{equation}

\end{tetel}

If $A$ is ill-conditioned, then Theorem \ref{bmodell} is not true.

\begin{pelda}
Let $A=\left[
\begin{array}
[c]{cc}%
1+\epsilon & 1\\
1 & 1-\epsilon
\end{array}
\right]  $, $\Delta A=\left[
\begin{array}
[c]{cc}%
0 & 0\\
0 & \epsilon^{2}%
\end{array}
\right]  $ and $b=\left[
\begin{array}
[c]{r}%
1\\
-1
\end{array}
\right]  $, ($0<\epsilon\ll1$). Then $\text{cond}_{\infty}\left(  A\right)
=\left(  2+\epsilon\right)  ^{2}/\epsilon^{2}\approx4/\epsilon^{2}$ and
$\left\|  \Delta A\right\|  _{\infty}/\left\|  A\right\|  _{\infty}%
=\epsilon^{2}/\left(  2+\epsilon\right)  \approx\epsilon^{2}/2$. Let
\[
\hat{x}=\left(  A+\Delta A\right)  ^{-1}b=\frac{1}{\epsilon^{3}}\left[
\begin{array}
[c]{c}%
2-\epsilon+\epsilon^{2}\\
-2-\epsilon
\end{array}
\right]  \approx\left[
\begin{array}
[c]{r}%
2/\epsilon^{3}\\
-2/\epsilon^{3}%
\end{array}
\right]  .
\]
Then $r\left(  \hat{x}\right)  =A\hat{x}-b=\left[
\begin{array}
[c]{c}%
0\\
2/{\epsilon}+1
\end{array}
\right]  $ and $\left\|  r\left(  \hat{x}\right)  \right\|  _{\infty}/\left\|
b\right\|  _{\infty}=2/{\epsilon}+1$, which is not small.
\end{pelda}

In the most general case we solve the perturbed equation
\begin{equation}
\left(  A+\Delta A\right)  \hat{x}=b+\Delta b
\end{equation}
instead of $Ax=b$. The following general result holds.

\begin{tetel}
\label{cmodell}If $A$ is nonsingular, $\text{\emph{cond}}\left(  A\right)
\frac{\left\|  \Delta A\right\|  }{\left\|  A\right\|  }<1$\ and $b\neq0$,
then
\begin{equation}
\frac{\left\|  \Delta x\right\|  }{\left\|  x\right\|  }\leq\frac
{\text{\emph{cond}}(A)\left(  \frac{\left\|  \Delta A\right\|  }{\left\|
A\right\|  }+\frac{\left\|  \Delta b\right\|  }{\left\|  b\right\|  }\right)
}{1-\text{\emph{cond}}(A)\frac{\left\|  \Delta A\right\|  }{\left\|
A\right\|  }} \koz .
\end{equation}
\end{tetel}

This theorem implies the following \idez{thumb rule}.

\noindent\textbf{Thumb rule. }\inddef{thumb rule} \emph{Assume that }$Ax=b$\emph{. If the entries
of }$A$\emph{\ and} $b$ \emph{are accurate to about }$s$\emph{\ decimal places
and }cond$(A)\sim10^{t}$\emph{, where }$t<s$\emph{, then the entries of the
computed solution are accurate to about }$s-t$\emph{\ decimal places.}

The assumption cond$\left(  A\right)  \left\|  \Delta A\right\|  /\left\|
A\right\|  <1$ of Theorem \ref{cmodell} guarantees that that matrix $A+\Delta
A$ is nonsingular. \ The inequality cond$\left(  A\right)  \left\|  \Delta
A\right\|  /\left\|  A\right\|  <1$ is equivalent with the inequality
$\left\|  \Delta A\right\|  <\frac{1}{\left\|  A^{-1}\right\|  }$ and the
distance of $A$ from the nearest singular matrix is just $1/\left\|
A^{-1}\right\|  $. Thus we can give a new characterization of the condition
number:
\begin{equation}
\frac{1}{\text{cond}\left(  A\right)  }=\min_{A+\Delta A\text{ is singular}%
}\frac{\left\|  \Delta A\right\|  }{\left\|  A\right\|  } \koz .
\end{equation}
Thus if a matrix is ill-conditioned, then it is close to a singular matrix.
Earlier we defined the condition numbers of matrices as the condition number
of the mapping $F\left(  x\right)  =A^{-1}x$.

Let us introduce the following definition.

\begin{defi}
A linear system solver is said to be \ki{weakly stable}\inddef{weakly stable} on a matrix class $H$, if
\ for all well-conditioned $A\in H$ and for all $b$, the computed solution
$\widehat{x}$ of the linear system $Ax=b$ has small relative error $\left\|
\hat{x}-x\right\|  /\left\|  x\right\|  $.
\end{defi}

Putting together Theorems \ref{amodell}--\ref{cmodell} we obtain the following.

\begin{tetel}[Bunch] \label{bunch}\nevindex{Bunch, James R.}
A linear system solver is weakly stable on a matrix
class $H$, if for all well-conditioned $A\in H$ and for all $b$, the computed
solution $\widehat{x}$ of the linear system $Ax=b$ satisfies any of the
following conditions: \newline(1) $\left\|  \hat{x}-x\right\|  /\left\|
x\right\|  $\ is small; \newline(2) $\left\|  r\left(  \hat{x}\right)
\right\|  /\left\|  b\right\|  $\ is small; \newline(3) There exists $\Delta
A$\ such that $\left(  A+\Delta A\right)  \hat{x}=b$\ and $\left\|  \Delta
A\right\|  /\left\|  A\right\|  $\ are small.
\end{tetel}

The estimate of Theorem \ref{cmodell} can be used in practice if \ we know
estimates of $\Delta b,\Delta A$ and cond$(A)$. If no estimates are available,
then we can only make a posteriori error estimates.

In the following we study the componentwise error estimates. We first give an
estimate for the absolute error of the approximate solution using the
components of the inverse error.

\begin{tetel}[Bauer, Skeel]
\label{bask}
\nevindex{Bauer, F. L.}\nevindex{Skeel, Robert D.}
Let $A\in\mathbb{R}^{n}$\ be nonsingular and
assume that the approximate solution $\widehat{x}$ of $Ax=b$ satisfies the
linear system $\left(  A+E\right)  \widehat{x}=b+e$. If $S\in\mathbb{R}%
^{n\times n}$ , $s\in\mathbb{R}^{n}$\ and $\varepsilon>0$\ are such that
$S\geq 0$, $s\geq 0$, $\left|  E\right|  \leq\varepsilon S$, $\left|  e\right|
\leq\varepsilon s$ and $\varepsilon\left\|  \left|  A^{-1}\right|  S\right\|
_{\infty}<1$, then
\begin{equation}
\left\|  \widehat{x}-x\right\|  _{\infty}\leq\frac{\varepsilon\left\|  \left|
A^{-1}\right|  \left(  S\left|  x\right|  +s\right)  \right\|  _{\infty}%
}{1-\varepsilon\left\|  \left|  A^{-1}\right|  S\right\|  _{\infty}}\koz.
\end{equation}

\end{tetel}

If $e=0$ ($s=0$), $S=\left|  A\right|  $ and%

\begin{equation}
k_{r}\left(  A\right)  =\left\|  \left|  A^{-1}\right|  \left|  A\right|
\right\|  _{\infty}<1 \koz,
\end{equation}
then we obtain the estimate%
\begin{equation}
\left\|  \widehat{x}-x\right\|  _{\infty}\leq\frac{\varepsilon k_{r}\left(
A\right)  }{1-\varepsilon k_{r}\left(  A\right)  }\koz.
\end{equation}
The quantity $k_{r}\left(  A\right)  $ is said to be \ki{Skeel-norm}\inddef{Skeel-norm}
\nevindex{Skeel, Robert D.}, although it
is not a norm in the earlier defined sense. The Skeel-norm satisfies the
inequality%
\begin{equation}
k_{r}\left(  A\right)  \leq\text{cond}_{\infty}\left(  A\right)  =\left\|
A\right\|  _{\infty}\left\|  A^{-1}\right\|  _{\infty} \koz.
\end{equation}
Therefore the above estimate is not worse than the traditional one that uses
the standard condition number.

The inverse error can be estimated componentwise by the following result of
Oettli and Prager. Let $A,\delta A\in\mathbb{R}^{n\times n}$ and $b,\delta
b\in\mathbb{R}^{n}$. Assume that $\delta A\geq0$ and $\delta b$ $\geq0$.
Furthermore let%
\[
D=\left\{  \Delta A\in\mathbb{R}^{n\times n}:\left|  \Delta A\right|
\leq\delta A\right\}  ,\quad G=\left\{  \Delta b\in\mathbb{R}^{n}:\left|
\Delta b\right|  \leq\delta b\right\}  \koz.
\]


\begin{tetel}[Oettli, Prager]\label{oettlp}\nevindex{Oettli, W.}\nevindex{Prager, W.}
The computed solution $\hat{x}$ satisfies
a perturbed equation $\left(  A+\Delta A\right)  \hat{x}=b+\Delta b$\ with
$\Delta A\in D$ and $\Delta b\in G$, if
\begin{equation}
\left|  r\left(  \hat{x}\right)  \right|  =\left|  A\hat{x}-b\right|
\leq\delta A\left|  \hat{x}\right|  +\delta b\koz. \label{ottli1}%
\end{equation}
\end{tetel}

We do not need the condition number to apply this theorem. In practice the
entries $\delta A$ and $\delta b$ are proportional to the machine epsilon.

\begin{tetel}[Wilkinson]
\label{wilkinson} \nevindex{Wilkinson, James H.}The approximate solution $\widehat{x}${ } of
$Ax=b$ obtained by the Gauss method in floating point arithmetic satisfies the
perturbed linear equation
\begin{equation}
\left(  A+\Delta A\right)  \widehat{x}=b
\end{equation}
with
\begin{equation}
\left\|  \Delta A\right\|  _{\infty}\leq8n^{3}\rho_{n}\left\|  A\right\|
_{\infty}u+O(u^{2})\koz,
\end{equation}
where $\rho_{n}$ denotes the groth factor of the pivot elements and $u$ is the
unit roundoff.
\end{tetel}

Since $\rho_{n}$ is small in practice, the relative error%
\[
\frac{\left\|  \Delta A\right\|  _{\infty}}{\left\|  A\right\|  _{\infty}}%
\leq8n^{3}\rho_{n}u+O(u^{2})
\]
is also small. Therefore Theorem\ref{bunch} implies that the Gauss method is
weakly stable both for full and partial pivoting.

Wilkinson's\nevindex{Wilkinson, James H.} theorem implies that%
\[
\text{cond}_{\infty}(A)\frac{\left\|  \Delta A\right\|  _{\infty}}{\left\|
A\right\|  _{\infty}}\leq8n^{3}\rho_{n}\text{cond}_{\infty}\left(  A\right)
u+O\left(  u^{2}\right)  \koz .
\]
For a small condition number we can assume that $1-$cond$_{\infty}(A)\left\|
\Delta A\right\|  _{\infty}/\left\|  A\right\|  _{\infty}\approx1$. Using
Theorems \ref{wilkinson} and \ref{cmodell} (case $\Delta b=0$) we obtain the
following estimate of the direct error:
\begin{equation}
\frac{\left\|  \Delta x\right\|  _{_{\infty}}}{\left\|  x\right\|  _{_{\infty
}}}\leq8n^{3}\rho_{n}\text{cond}_{\infty}\left(  A\right)  u\koz.
\end{equation}
The obtained result supports the thumb rule in the case of the Gauss method.

\begin{pelda}
Consider the following linear system whose coefficients can be represented
exactly:%
\[%
\begin{array}
[c]{c}%
888445x_{1}+887112x_{2}=1\koz ,\\
887112x_{1}+885781x_{2}=0\koz .
\end{array}
\]
Here $\text{cond}(A)_{\infty}$ is big, but $\text{cond}_{\infty}(A){\left\|
\Delta A\right\|  _{\infty}}/{\left\|  A\right\|  _{\infty}}$ is negligible.
The exact solution of the problem is $x_{1}=885781$, $x_{2}=-887112$. The
MATLAB gives the approximate solution $\hat{x}_{1}=885827.23$, $\hat{x}%
_{2}=-887158.30$ with the relative error
\[
\frac{\left\|  x-\hat{x}\right\|  _{_{\infty}}}{\left\|  x\right\|
_{_{\infty}}}=5.22\times10^{-5}\koz.
\]
Since $s\approx16$ and $\text{cond}\left(  A\right)  _{_{\infty}}%
\approx3.15\times10^{12}$, the result essentially corresponds to the Wilkinson\nevindex{Wilkinson, James H.}
theorem or the thumb rule. The Wilkinson\nevindex{Wilkinson, James H.} theorem gives the bound
\[
\left\|  \Delta A\right\|  _{\infty}\leq1.26\times10^{-8}%
\]
for the inverse error. If we use the Oettli-Prager theorem with the choice
$\delta A=\epsilon_{M}\left|  A\right|  $ and $\delta b=\epsilon_{M}\left|
b\right|  $, then we obtain the estimate $\left|  r\left(  \hat{x}\right)
\right|  \leq\delta A\left|  \hat{x}\right|  +\delta b$. Since $\left\|
\left|  \delta A\right|  \right\|  _{\infty}=3.94\times10^{-10}$, this
estimate is better than that of Wilkinson\nevindex{Wilkinson, James H.}.
\end{pelda}

\textbf{Scaling and preconditioning.}
Several matrices that occur in applications are ill-conditioned if their order
$n$ is large. For example the famous Hilbert-matrix\nevindex{Hilbert, David}
\begin{equation}
H_{n}=\left[  \frac{1}{i+j-1}\right]  _{i,j=1}^{n}%
\end{equation}
has cond$_{2}\left(  H_{n}\right)  \approx e^{3.5n}$, if $n\rightarrow\infty$.
There exist $2n\times2n$ matrices with integer entries that can be represented
exactly in standard IEEE754 floating point arithmetic while their condition
number is approximately $4\times10^{32n}$.

We have two main techniques to solve linear systems with large condition
numbers. Either we use multiple precision arithmetic or decrease the condition
number. There are two known forms of decreasing the condition number.

\textbf{1. Scaling.} We replace the linear system $Ax=b$ with the equation%
\begin{equation}
\left(  RAC\right)  y=\left(  Rb\right)  \koz,
\end{equation}
where $R$ and $C$ are diagonal matrices.

We apply the Gauss method to this scaled system and get the solution $y$. The
quantity $x=Cy$ defines the requested solution. If the condition number of the
matrix $RAC$ is smaller then we expect a smaller error in $y$ and consequently
in $x$. Various strategies are given to choose the scaling matrices $R$ and
$C$. One of the best known strategies is the \ki{balancing}\inddef{balancing} which forces every
column and row of $RAC$ to have approximately the same norm. For example, if
\[
\hat{D}=\text{diag}\left(  \frac{1}{\left\|  a_{1}^{T}\right\|  _{2}}%
,\ldots,\frac{1}{\left\|  a_{n}^{T}\right\|  _{2}}\right)
\]
where $a_{i}^{T}$ is the $i^{\text{th}}$ row vector of $A$, the Euclidean
norms of the rows of $\hat{D}A$ will be $1$ and the estimate
\[
\text{cond}_{2}\left(  \hat{D}A\right)  \leq\sqrt{n}\min_{D\in D_{+}}\text{
cond}_{2}\left(  DA\right)
\]
holds with $D_{+}=\left\{  \text{diag}\left(  d_{1},\ldots,d_{n}\right)
\ |\ d_{1},\ldots,d_{n}>0\right\}  $. This means that $\hat{D}$ optimally
scales the rows of $A$ in an approximate sense.

The next example shows that the scaling may lead to bad results.

\begin{pelda}
Consider the matrix%
\[
A=\left[
\begin{array}
[c]{ccc}%
\epsilon/2 & 1 & 1\\
1 & 1 & 1\\
1 & 1 & 2
\end{array}
\right]
\]
for $0<\epsilon\ll1$. It is easy to show that $\text{cond}_{\infty}\left(
A\right)  =12$. Let
\[
R=C=\left[
\begin{array}
[c]{lll}%
2/\sqrt{\epsilon} & 0 & 0\\
0 & \sqrt{\epsilon}/2 & 0\\
0 & 0 & \sqrt{\epsilon}/2
\end{array}
\right]  \koz .
\]
Then the scaled matrix%
\[
RAR=\left[
\begin{array}
[c]{lll}%
2 & 1 & 1\\
1 & \epsilon/4 & \epsilon/4\\
1 & \epsilon/4 & \epsilon/2
\end{array}
\right]  ,
\]
has the condition number $\text{cond}_{\infty}\left(  RAR\right)
=32/\epsilon$, which a very large value for small $\epsilon$.
\end{pelda}

\textbf{2. Preconditioning.} The preconditioning is very close to scaling. We
rewrite the linear system $Ax=b$ with the equivalent form
\begin{equation}
\widetilde{A}x=\left(  MA\right)  x=Mb=\tilde{b},
\end{equation}
where matrix $M$ is such that cond$\left(  M^{-1}A\right)  $ is smaller and
$Mz=y$ is easily solvable.

The preconditioning is often used with iterative methods on linear systems
with symmetric and positive definite matrices.

\textbf{A posteriori error estimates.}
The a posteriori estimate of the error of an approximate solution is necessary
to get some information on the reliability of the obtained result. There are
plenty of such estimates. Here we show three estimates whose computational
cost is $\Theta\left(  n^{2}\right)  $ flops. This cost is acceptable when
comparing to the cost of direct or iterative methods ($\Theta\left(
n^{3}\right)  $ or $\Theta\left(  n^{2}\right)  $ per iteration step).

\paragraph{The estimate of the direct error with the residual error.}

\begin{tetel}[Auchmuty]
\label{auchm}
\nevindex{Auchmuty, Giles}
Let $\widehat{x}$\ be the approximate solution of
$Ax=b$. Then
\[
\left\|  x-\hat{x}\right\|  _{2}=\frac{c\left\|  r(\hat{x})\right\|  _{2}^{2}%
}{\left\|  A^{T}r(\hat{x})\right\|  _{2}} \koz ,
\]
where $c\geq1$.
\end{tetel}

The error constant $c$ depends on $A$ and the direction of error vector
$\hat{x}-x$. Furthermore
\[
\frac{1}{2}\text{cond}_{2}\left(  A\right)  \approx C_{2}\left(  A\right)
=\frac{1}{2}\left(  \text{cond}_{2}\left(  A\right)  +\frac{1}{\text{cond}%
_{2}\left(  A\right)  }\right)  \leq\text{cond}_{2}\left(  A\right)  \koz .
\]


The error constant $c$ takes the upper value $C_{2}\left(  A\right)  $ only in
exceptional cases. The computational experiments indicate that the average
value of $c$ grows slowly with the order of $A$ and it depends more strongly
on $n$ than the condition number of $A$. The following experimental estimate
\begin{equation}
\left\|  x-\widehat{x}\right\|  _{2}\lessapprox0.5\dim\left(  A\right)
\left\|  r\left(  \widehat{x}\right)  \right\|  _{2}^{2}/\left\|
A^{T}r\left(  \widehat{x}\right)  \right\|  _{2} \label{pracest}%
\end{equation}
seems to hold with a high degree of probability.

\paragraph{The LINPACK estimate of $\left\|  A^{-1}\right\|.$}
The famous LINPACK program package uses the following process to estimate
$\left\|  A^{-1}\right\|  $. We solve the linear systems $A^{T}y=d$ and
$Aw=y$. Then the estimate of $\left\|  A^{-1}\right\|  $ is given by%
\begin{equation}
\left\|  A^{-1}\right\|  \approx\frac{\left\|  w\right\|  }{\left\|
y\right\|  }\quad\left(  \leq\left\|  A^{-1}\right\|  \right)  .
\end{equation}
Since
\[
\frac{\left\|  w\right\|  }{\left\|  y\right\|  }=\frac{\left\|  A^{-1}\left(
A^{-T}d\right)  \right\|  }{\left\|  A^{-T}d\right\|  } \koz ,
\]
we can interpret the process as an application of the power method of the
eigenvalue problem. The estimate can be used with the $1-,2-$ and $\infty
$-norms. The entries of vector $d$ are $\pm1$ possibly with random signs.

If the linear system $Ax=b$ is solved by the $LU$-method, then the solution of
further linear systems costs $\Theta(n^{2})$ flops per system. Thus the total
cost of the LINPACK estimate remains small. Having the estimate $\left\|
A^{-1}\right\|  $ we can easily estimate cond$(A)$ and the error of the
approximate solution (cf. Theorem \ref{cmodell} or the thumb rule). We remark
that several similar processes are known in the literature.

\paragraph{The Oettli-Prager estimate of the inverse error.}
We use the Oettli-Prager\nevindex{Oettli, W.}\nevindex{Prager, W.} theorem in the following form. Let $r\left(  \hat
{x}\right)  =A\hat{x}-b$ be the residual error, $E\in\mathbb{R}^{n\times n}$
and $f\in\mathbb{R}^{n}$ are given such that $E\geq0$ and $f\geq0$. Let
\[
\omega=\max_{i}\frac{\left|  r\left(  \hat{x}\right)  _{i}\right|  }{\left(
E\left|  \hat{x}\right|  +f\right)  _{i}} \koz ,
\]
where $0/0$ is set to $0$-nak, $\rho/0$ is set to $\infty$, if $\rho\neq0$.
Symbol $\left(  y\right)  _{i}$ denotes the $i^{\text{th}}$ component of the
vector $y$. If $\omega\neq\infty$, then there exist a matrix $\Delta A$ and a
vector $\Delta b$ for which
\[
\left|  \Delta A\right|  \leq\omega E,\quad\left|  \Delta b\right|  \leq\omega
f
\]
holds and
\[
\left(  A+\Delta A\right)  \hat{x}=b+\Delta b \koz .
\]
Moreover $\omega$ is the smallest number for which $\Delta A$ and $\Delta b$
exist with the above properties. The quantity $\omega$ measures the relative
inverse error in terms of $E$ and $f$. If for a given $E$, $f$ and $\hat{x}$,
the quantity $\omega$ is small, then the perturbed problem (and its solution)
are close to the original problem (and its solution). In practice, the choice
$E=\left|  A\right|  $ and $f=\left|  b\right|  $ is preferred

\subsubsection{Iterative refinement.}
Denote by $\hat{x}$ the approximate solution of $Ax=b$ and let $r(y)=Ay-b$ be
the residual error at the point $y$. The precision of the approximate solution
$\hat{x}$ can be improved with the following method.

\begin{alg}{Iterative-Refinement($A,b,\hat{x},tol$)}\inddef{\textsc{Iterative-Refinement}}
1 \' $k\leftarrow 1$ \\
2 \'$x_{1}\leftarrow \hat{x}$ \\
3 \' $\hat{d}\leftarrow \mathrm{inf}$ \\
4 \' \key{while} $\left\| \hat{d}\right\| /\left\| x_{k}\right\| >tol$ \\
5 \'  \> \key{do} \= $r\leftarrow Ax_{k}-b$ \\
6 \' \> \> Compute the approximate solution $\hat{d}$ of $Ad=r$ with
the $LU$-method.\\
7 \' \> \> $x_{k+1}\leftarrow x_{k}-\hat{d}$ \\
8 \'  \> \> $k\leftarrow k+1$ \\
9 \'  \key{return} $x_{k}$
\end{alg}


There are other variants of this process. We can use other linear solvers
instead of the $LU$-method.

Let $\eta$ be the smallest bound of relative inverse error with
\[
\left(  A+\Delta A\right)  \hat{x}=b+\Delta b,\quad\left|  \Delta A\right|
\leq\eta\left|  A\right|  ,\quad\left|  \Delta b\right|  \leq\eta\left|
b\right|  \koz .
\]
Furthermore let%
\[
\sigma\left(  A,x\right)  =\max_{k}\left(  \left|  A\right|  \left|  x\right|
\right)  _{k}/\min_{k}\left(  \left|  A\right|  \left|  x\right|  \right)
_{k},\quad\min_{k}\left(  \left|  A\right|  \left|  x\right|  \right)
_{k}>0 \koz .
\]


\begin{tetel}[Skeel]
\label{skeel}
\nevindex{Skeel, Robert D.}
If $k_{r}\left(  A^{-1}\right)  \sigma\left(
A,x\right)  \leq c_{1}<1/\epsilon_{M},$ then for sufficiently large $k$ we
have%
\begin{equation}
\left(  A+\Delta A\right)  x_{k}=b+\Delta b,\quad\left|  \Delta A\right|
\leq4\eta\epsilon_{M}\left|  A\right|  ,\quad\left|  \Delta b\right|
\leq4\eta\epsilon_{M}\left|  b\right|  \koz .
\end{equation}

\end{tetel}

This result often holds after the first iteration, i.e. for $k=2$. Jankowski\nevindex{Jankowski, T.}
and  Wozniakowski\nevindex{Wozniakowski, T.} investigated the iterative refinement for any method
$\phi$ which produces an approximate solution $\widehat{x}${ with relative
error less than }$1$. They showed that the iterative refinement improves the
precision of the approximate solution even in single precision arithmetic and
makes method $\phi$ to be weakly stable.

\begin{gyak}
\ujgyak Prove Theorem \ref{muvig}.
\ujgyak Consider the linear systems $Ax=b$ and $Bx=b$, where

\[
A=\left[
\begin{array}
[c]{rr}%
1 & 1/2\\
1/2 & 1/3
\end{array}
\right]  ,\quad B=\left[
\begin{array}
[c]{rr}%
1 & -1/2\\
1/2 & 1/3
\end{array}
\right]
\]
and $b\in\mathbb{R}^{2}$. Which equation is more sensitive to the perturbation
of $b$? What should be the \emph{relative error}\gyakindex{relative error} of $b$ in the more sensitive
equation in order to get the solutions of both equations with the same precision?
\ujgyak Let $\chi=3/2^{29}$, $\zeta=2^{14}$ and
\[
A=\left[
\begin{array}
[c]{ccc}%
\chi\zeta & -\zeta & \zeta\\
\zeta^{-1} & \zeta^{-1} & 0\\
\zeta^{-1} & -\chi\zeta^{-1} & \zeta^{-1}%
\end{array}
\right]  ,\quad b=\left[
\begin{array}
[c]{r}%
1\\
1+\varepsilon\\
1
\end{array}
\right]  \ .
\]
Solve the linear systems $Ax=b$ for $\varepsilon=10^{-1},10^{-3}%
,10^{-5},10^{-7},10^{-10}$. Explain the results.
\ujgyak Let $A$ be a $10\times10$ matrix and choose the band matrix consisting of the
main and the neighbouring two subdiagonals of $A$ as a preconditioning matrix.
How much does the \emph{condition number}\gyakindex{condition number} of $A$ improves if (i) $A$ is a random
matrix; (ii) $A$ is a Hilbert matrix?\nevindex{Hilbert, David}
\ujgyak Let
\[
A=\left[
\begin{array}
[c]{ccc}%
1/2 & 1/3 & 1/4\\
1/3 & 1/4 & 1/5\\
1/4 & 1/5 & 1/6
\end{array}
\right]  ,
\]
and assume that $\varepsilon$ is the common error bound of every component of
$b\in\mathbb{R}^{3}$. Give the sharpest possible error bounds for the solution
$\left[  x_{1},x_{2},x_{3}\right]  ^{T}$ of the equation $Ax=b$ and for the
sum $(x_{1}+x_{2}+x_{3})$.
\ujgyak Consider the linear system $Ax=b$ with the approximate solution $\hat{x}$.
\newline(i) Give an error bound for $\hat{x}$, if ($A+E)\hat{x}=b$ holds
exactly and both $A$ and $A+E$ is nonsingular.\newline(ii) Let
\[
A=\left[
\begin{array}
[c]{rrr}%
10 & 7 & 8\\
7 & 5 & 6\\
8 & 6 & 10
\end{array}
\right]  ,\,\quad b=\left[
\begin{array}
[c]{r}%
25\\
18\\
24
\end{array}
\right]
\]
and consider the solution of $Ax=b$. Give (if possible) a relative error bound
for the entries of $A$ such that the integer part of every solution component
remains constant within the range of this relative error bound.
\end{gyak}

\section{Eigenvalue problems}
The set of complex $n$-vectors will be denoted by $\mathbb{C}^{n}$.
Similarly,
$\mathbb{C}^{m\times n}$  denotes the set of 
\emph{complex $m\times n$ matrices}.

\begin{defi}
Let $A\in\mathbb{C}^{n\times n}$\ be an arbitrary matrix. The number
$\lambda\in\mathbb{C}$ is the \ki{eigenvalue}
\inddef{eigenvalue}
 of $A$ if there is vector
$x\in\mathbb{C}^{n}$\ ($x\neq0$) such that
\begin{equation}
Ax=\lambda x. \label{sdefegy}%
\end{equation}
Vector $x$ is called the (right) eigenvector of $A$ that belongs to the
eigenvalue $\lambda$.
\end{defi}

Equation $Ax=\lambda x$ can be written in the equivalent form $(A-\lambda
I)x=0$, where $I$ is the unit matrix of appropriate size. The latter
homogeneous linear system has a nonzero solution $x$ if and only if
\begin{equation}
\phi(\lambda)=\det(A-\lambda I)=\det\left(  \left[
\begin{array}
[c]{cccc}%
a_{11}-\lambda & a_{12} & \ldots & a_{1n}\\
a_{21} & a_{22}-\lambda & \ldots & a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n1} & a_{n2} & \ldots & a_{nn}-\lambda
\end{array}
\right]  \right)  =0\ . \label{karegy}%
\end{equation}
Equation (\ref{karegy}) is called the
\ki{characteristic equation}\inddef{characteristic equation} of matrix $A$.
The roots of this equation are the eigenvalues of matrix $A$. Expanding
$\det\left(  A-\lambda I\right)  $ we obtain a polynomial of degree $n$:
\[
\phi(\lambda)=(-1)^{n}(\lambda^{n}-p_{1}\lambda^{n-1}-\ldots-p_{n-1}%
\lambda-p_{n}) \koz .
\]
This polynomial called the 
\ki{characteristic polynomial} \inddef{characteristic polynomial} of $A$.
It follows from
the fundamental theorem of algebra that any matrix $A\in\mathbb{C}^{n\times
n}$ has exactly $n$ eigenvalues with multiplicities. The eigenvalues may be
complex or real. Therefore one needs to use complex arithmetic for eigenvalue
calculations. If the matrix is real and the computations are done in real
arithmetic, the complex eigenvalues and eigenvectors can be determined only
with special techniques.

If $x\neq 0$ is an eigenvector, $t\in\mathbb{C}$ ($t\neq0$), then $tx$ is also
eigenvector. The number of linearly independent eigenvectors that belong to an
eigenvalue $\lambda_{k}$ does not exceed the multiplicity of $\lambda_{k}$ in
the characteristic equation (\ref{karegy}). The eigenvectors that belong to
different eigenvalues are linearly independent.

The following results give estimates for the size and location of the eigenvalues.

\begin{tetel}
Let $\lambda$\ be any eigenvalue of matrix $A$. The upper estimate $\left|
\lambda\right|  \leq\left\|  A\right\|  $ holds in any induced matrix norm.
\end{tetel}

\begin{tetel}[Gersgorin]
\nevindex{Gersgorin, S. A.} 
Let $A\in\mathbb{C}^{n\times n}$,
\[
r_{i}=\sum_{j=1,j\neq i}^{n}\left|  a_{ij}\right|  \quad\left(  i=1,\ldots
,n\right)
\]
and
\[
D_{i}=\left\{  z\in\mathbf{C}|\left|  z-a_{ii}\right|  \leq r_{i}\right\}
\quad\left(  i=1,\ldots,n\right)  .
\]
Then for any eigenvalue $\lambda$ of $A$ we have $\lambda\in\cup_{i=1}%
^{n}D_{i}$.
\end{tetel}

For certain matrices the solution of the characteristic equation
(\ref{karegy}) is very easy. For example, if $A$ is a triangular matrix, then
its eigenvalues are entries of the main diagonal. In most cases however the
computation of all eigenvalues and eigenvectors is a very difficult task.
Those transformations of matrices that keeps the eigenvalues unchanged have
practical significance for this problem. Later we see that the eigenvalue
problem of transformed matrices is simpler.

\begin{defi}
The matrices $A,B\in\mathbb{C}^{n\times n}$ are similar if there is a matrix
$T$ such that $B=T^{-1}AT$. The mapping $A\rightarrow T^{-1}AT$ is said to be
\ki{similarity transformation} \inddef{similarity transformation} of $A$.
\end{defi}

\begin{tetel}
Assume that $\det(T)\neq0$. Then the eigenvalues of $A$\ and $B=T^{-1}AT$ are
the same. If $x$ is the eigenvector of $A$, then $y=T^{-1}x$\ is the
eigenvector of $B$.
\end{tetel}

Similar matrices have the same eigenvalues.

The difficulty of the eigenvalue problem also stems from the fact that the
eigenvalues and eigenvectors are very sensitive (unstable) to changes in the
matrix entries. The eigenvalues of $A$ and the perturbed matrix $A+\delta A$
may differ from each other significantly. Besides the multiplicity of the
eigenvalues may also change under perturbation. The following theorems and
examples show the very sensitivity of the eigenvalue problem.

\begin{tetel}[Ostrowski, Elsner]
\label{oseltetel}
\nevindex{Ostrowski, Alexander R.}
\nevindex{Elsner, Ludwig} For every eigenvalue $\lambda_{i}$ of
matrix $A\in\mathbb{C}^{n\times n}$ there exists an eigenvalue $\mu_{k}$ of
the perturbed matrix $A+\delta A$ such that%
\[
\left|  \lambda_{i}-\mu_{k}\right|  \leq\left(  2n-1\right)  \left(  \left\|
A\right\|  _{2}+\left\|  A+\delta A\right\|  _{2}\right)  ^{1-\frac{1}{n}%
}\left\|  \delta A\right\|  _{2}^{\frac{1}{n}} \koz .
\]

\end{tetel}

We can observe that the eigenvalues are changing continuously and the size of
change is proportional to the $n^{\text{th}}$ root of $\left\|  \delta
A\right\|  _{2}$.

\begin{pelda}
Consider the following perturbed 
\emph{Jordan matrix}\nevindex{Jordan, Camille} 
 of the size $r\times r$:
\[
\left[
\begin{array}
[c]{ccccc}%
\mu & 1 & 0 & \ldots & 0\\
0 & \mu & 1 & \ddots & \vdots\\
\vdots & \ddots & \ddots & \ddots & 0\\
0 &  & \ddots & \mu & 1\\
\epsilon & 0 & \ldots & 0 & \mu
\end{array}
\right]  \koz .
\]
The characteristic equation\index{characteristic equation}
is $\left(  \lambda-\mu\right)  ^{r}=\epsilon$,
which gives the $r$ different eigenvalues
\[
\lambda_{s}=\mu+\epsilon^{1/r}\left(  \cos\left(  2s\pi/r\right)
+i\sin\left(  2s\pi/r\right)  \right)  \quad(s=0,\ldots,r-1)
\]
instead of the original eigenvalue $\mu$ with multiplicity $r$. The size of
\ change is $\epsilon^{1/r}$, which corresponds to Theorem (\ref{oseltetel} ).
If $\left|  \mu\right|  \approx1$, $r=16$ and $\epsilon=\epsilon_{M}%
\approx2.2204\times10^{-16}$, then the perturbation size of the eigenvalues is
$\approx$ $0.1051$. This is a significant change relative to the input
perturbation $\epsilon$.
\end{pelda}

For special matrices and perturbations we may have much better
perturbation bounds.

\begin{tetel}[Bauer, Fike]
\nevindex{Bauer, F. L.}
Assume that $A\in\mathbb{C}^{n\times n}$ is 
\ki{diagonalisable,}\inddef{diagonalisable} 
that is a matrix $X$ exists such that $X^{-1}AX=\text{\emph{diag}}(\lambda
_{1},\ldots,\lambda_{n}).$ Denote $\mu$ an eigenvalue of $A+\delta A$. Then%
\begin{equation}
\min_{1\leq i\leq n}\left|  \lambda_{i}-\mu\right|  \leq\text{\emph{cond}}%
_{2}(X)\left\|  \delta A\right\|  _{2}. \label{bf1}%
\end{equation}

\end{tetel}

This result is better than that of Ostrowski and Elsner. 
\nevindex{Ostrowski, Alexander R.}
\nevindex{Elsner, Ludwig} Nevertheless
cond$_{2}\left(  X\right)  $, which is generally unknown, can be very big.

The eigenvalues are continuous functions of the matrix entries. This is also
true for the normalized eigenvectors if the eigenvalues are simple. The
following example shows that this property does not hold for multiple eigenvalues.

\begin{pelda}
Let
\[
A\left(  t\right)  =\left[
\begin{array}
[c]{rr}%
1+t\cos\left(  2/t\right)  & -t\sin\left(  2/t\right) \\
-t\sin\left(  2/t\right)  & 1-t\cos\left(  2/t\right)
\end{array}
\right]  \quad\left(  t\neq0\right)  .
\]
The eigenvalues of $A\left(  t\right)  $ are $\lambda_{1}=1+t$ and
$\lambda_{2}=1-t$. Vector $\left[  \sin\left(  1/t\right)  ,\cos\left(
1/t\right)  \right]  ^{T}$ is the eigenvector belonging to $\lambda_{1}$.
Vector $\left[  \cos\left(  1/t\right)  ,-\sin\left(  1/t\right)  \right]
^{T}$ is the eigenvector belonging to $\lambda _{2}$. If $t\rightarrow 0$, then
\[
A\left(  t\right)  \rightarrow I=\left[
\begin{array}
[c]{cc}%
1 & 0\\
0 & 1
\end{array}
\right]  ,\quad\lambda _{1},\lambda _{2}\rightarrow 1 \koz,
\]
while the eigenvectors do not have limit.
\end{pelda}

We study the numerical solution of the eigenvalue problem in the\ next
section. Unfortunately it is very difficult to estimate the goodness of
numerical approximations. From the fact that $Ax-\lambda x=0$ holds with a
certain error we cannot conclude anything in general.

\begin{pelda}
\label{respelda}Consider the matrix%
\[
A\left(  \epsilon\right)  =\left[
\begin{array}
[c]{cc}%
1 & 1\\
\epsilon & 1
\end{array}
\right]  \koz ,
\]
where $\epsilon\approx 0$ is small. The eigenvalues of $A\left(  \epsilon
\right)  $ are $1\pm\sqrt{\epsilon}$, while the corresponding eigenvectors are
$[1,\pm\sqrt{\epsilon}]^{T}$. Let $\mu=1$ be an approximation of the
eigenvalues and let $x=\left[  1,0\right]  ^{T}$ be the approximate
eigenvector. Then%
\[
\left\|  Ax-\mu x\right\|  _{2}=\left\|  \left[
\begin{array}
[c]{c}%
0\\
\epsilon
\end{array}
\right]  \right\|  _{2}=\epsilon \koz .
\]
If $\epsilon=10^{-10}$, then the residual error under estimate the true error
$10^{-5}$ by five order.
\end{pelda}

\begin{megj}
\label{kondmj}We can define the
\ki{condition number of eigenvalues}
\inddef{condition number of eigenvalues}  for simple
eigenvalues:
\[
\nu\left(  \lambda_{1}\right)  \approx\frac{\left\|  x\right\|  _{2}\left\|
y\right\|  _{2}}{\left|  x^{H}y\right|  } \koz ,
\]
where $x$ and $y$ are the right and left eigenvectors, respectively. For
multiple eigenvalues the condition number is not finite.
\end{megj}

\subsection{Iterative solutions of the eigenvalue problem}
We investigate only the real eigenvalues and eigenvectors of real matrices.
The methods under consideration can be extended to the complex case with
appropriate modifications.

\paragraph{The power method.}
This method is due to von Mieses.\nevindex{von Mieses, Richard (1853--1953)}
Assume that $A\in\mathbb{R}^{n\times n}$ has
exactly $n$ different real eigenvalues. Then the eigenvectors $x_{1}%
,\ldots,x_{n}$ belonging to the corresponding eigenvalues $\lambda_{1}%
,\ldots,\lambda_{n}$ are linearly independent. Assume that the eigenvalues
satisfy the condition%
\[
\left|  \lambda_{1}\right|  >\left|  \lambda_{2}\right|  \geq\cdots\geq\left|
\lambda_{n}\right|
\]
and let $v^{(0)}\in\mathbb{R}^{n}$ be a given vector. This vector is a unique
linear combination of the eigenvectors, that is $v^{(0)}=\alpha_{1}%
x_{1}+\alpha_{2}x_{2}+\cdots+\alpha_{n}x_{n}$. Assume that $\alpha_{1}\neq0$
and compute the sequence $v^{(k)}=Av^{(k-1)}=A^{k}v^{(0)}$ ($k=1,2,\ldots$).
The initial assumptions imply that%
\[%
\begin{array}
[c]{ccccl}%
v^{(k)} & = & Av^{(k-1)} & = & A(\alpha_{1}\lambda_{1}^{k-1}x_{1}+\alpha
_{2}\lambda_{2}^{k-1}x_{2}+\ldots+\alpha_{n}\lambda_{n}^{k-1}x_{n})\\
&  &  & = & \alpha_{1}\lambda_{1}^{k}x_{1}+\alpha_{2}\lambda_{2}^{k}%
x_{2}+\cdots+\alpha_{n}\lambda_{n}^{k}x_{n}\\
&  &  & = & \lambda_{1}^{k}\left(  \alpha_{1}x_{1}+\alpha_{2}\left(
\frac{\lambda_{2}}{\lambda_{1}}\right)  ^{k}x_{2}+\cdots+\alpha_{n}\left(
\frac{\lambda_{n}}{\lambda_{1}}\right)  ^{k}x_{n}\right)  \koz .
\end{array}
\]
Let $y\in\mathbb{R}^{n}$ be an arbitrary vector such that $y^{T}x_{1}\neq0$. Then%

\[
\frac{y^{T}Av^{\left(  k\right)  }}{y^{T}v^{\left(  k\right)  }}=\frac{
y^{T}v^{\left(  k+1\right)  }}{y^{T}v^{\left(  k\right)  }}=\frac{\lambda
_{1}^{k+1}\left(  \alpha_{1}y^{T}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(  \frac{
\lambda_{i}}{\lambda_{1}}\right)  ^{k+1}y^{T}x_{i}\right)  }{\lambda_{1}%
^{k}\left(  \alpha_{1}y^{T}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(  \frac{
\lambda_{i}}{\lambda_{1}}\right)  ^{k}y^{T}x_{i}\right)  }\rightarrow
\lambda_{1}\ \koz .
\]


Given the initial vector $v^{\left(  0\right)  }\in\mathbb{R}^{n}$, the power
method has the following form.

\begin{alg}{Power-Method($A,v^{\left(0\right)}$)}\inddef{\textsc{Power-Method}}
1 \' $k\leftarrow 0$ \\
2 \' \key{while} \= \emph{exit condition} = \textsc{false} \\
3 \' \> \key{do} \= $k\leftarrow k+1$ \\
4 \' \> \>  $z^{(k)}\leftarrow Av^{(k-1)}$ \\
5 \' \> \> Select vector $y$
such that $y^{T}v^{\left(  k-1\right)  }\neq0$ \\
6 \' \> \> $\gamma _{k}\leftarrow y^{T}z^{\left( k\right) }/y^{T}v^{
\left( k-1\right) }$ \\
7 \' \> \> $v^{(k)}\leftarrow z^{(k)}/\left\| z^{(k)}\right\| _{\infty }$ \\
8 \' \key{return} $\gamma _{k},v^{(k)}$

\end{alg}

It is clear that%
\[
v^{(k)}\rightarrow x_{1},\quad\gamma_{k}\rightarrow\lambda_{1} \koz .
\]
The convergence $v^{\left(  k\right)  }\rightarrow x_{1}$ here means that
$\left(  v^{\left(  k\right)  },x_{1}\right)  _{\measuredangle}\rightarrow0$,
that is the action line of $v^{\left(  k\right)  }$ tends to the action line
of $x_{1}$. There are various strategies to select $y$. We can select
$y=e_{i}$, where $i$ is defined by $\left|  v_{i}^{\left(  k\right)  }\right|
=\left\|  v^{\left(  k\right)  }\right\|  _{\infty}$. If we select
$y=v^{\left(  k-1\right)  }$, then $\gamma_{k}=v^{\left(  k-1\right)
T}Av^{\left(  k-1\right)  }/\left(  v^{\left(  k-1\right)  T}v^{\left(
k-1\right)  }\right)  $ will be identical with the 
\emph{Rayleigh quotient}\index{Rayleigh quotient}
\nevindex{Rayleigh, John William Strutt}
$\mathbf{R}\left(  v^{\left(  k-1\right)  }\right)  $. This choice gives an
approximation of $\lambda_{1}$ that have the minimal residual norm (Example
\ref{respelda} shows that this choice is not necessarily the best option).

The speed of convergence depends on the quotient $\left|  \lambda_{2}%
/\lambda_{1}\right|  $. The method is very sensitive to the choice of the
initial vector $v^{(0)}$. If $\alpha_{1}=0$, then the process does not
converge to the dominant eigenvalue $\lambda_{1}$. For certain matrix classes
the \ power method converges with probability $1$ if the initial vector
$v^{\left(  0\right)  }$ is randomly chosen. In case of complex eigenvalues or
multiple $\lambda_{1}$ we have to use modifications of the algorithm. The
speed of convergence can be accelerated if the method is applied to the
shifted matrix $A-\sigma I$, where $\sigma$ is an appropriately chosen number.
The shifted matrix $A-\sigma I$ has the eigenvalues $\lambda_{1}%
-\sigma,\lambda_{2}-\sigma,\ldots,\lambda_{n}-\sigma$ and the corresponding
convergence factor $\left|  \lambda_{2}-\sigma\right|  /\left|  \lambda
_{1}-\sigma\right|  $. The latter quotient can be made smaller than $\left|
\lambda_{2}/\lambda_{1}\right|  $ with the proper selection of $\sigma$.

The usual exit condition of the power method is%
\[
\left\|  E_{k}\right\|  _{2}=\frac{\left\|  r_{k}\right\|  _{2}}{\left\|
v^{\left(  k\right)  }\right\|  _{2}}=\frac{\left\|  Av^{\left(  k\right)
}-\gamma_{k}v^{\left(  k\right)  }\right\|  _{2}}{\left\|  v^{\left(
k\right)  }\right\|  _{2}}\leq\epsilon \koz .
\]
If we simultaneously apply the power method to the transposed matrix $A^{T}$
and $w_{k}=\left(  A^{T}\right)  ^{k}w_{0}$, then the quantity
\[
\nu\left(  \lambda_{1}\right)  \approx\frac{\left\|  w^{\left(  k\right)
}\right\|  _{2}\left\|  v^{\left(  k\right)  }\right\|  _{2}}{\left|
w_{k}^{T}v_{k}\right|  }%
\]
gives an estimate for the condition number of $\lambda_{1}$ (see Remark 
\ref{kondmj}). In such a case we use the exit condition
\[
\nu\left(  \lambda_{1}\right)  \left\|  E_{k}\right\|  _{2}\leq\epsilon \koz .
\]


The power method is very useful for large sparse matrices. It is often used to
determine the largest and the smallest eigenvalue. We can approximate the
smallest eigenvalue as follows. The eigenvalues of $A^{-1}$ are $1/\lambda
_{1},\ldots,1/\lambda_{n}$. The eigenvalue $1/\lambda_{n}$ will be the
eigenvalue with the largest modulus. We can approximate this value by applying
the power method to $A^{-1}$. This requires only a small modification of the
algorithm. We replace line 4. with the following:

\begin{center}
Solve equation{ $Az^{(k)}=v^{(k-1)}$ for $z^{(k)}$}
\end{center}

The modified algorithm is called the 
\ki{inverse power method.}\inddef{inverse power method} It is clear that
$\gamma_{k}\rightarrow1/\lambda_{n}$ and $v^{(k)}\rightharpoonup x_{n}$ hold
under appropriate conditions. If we use the $LU$-method to solve
$Az^{(k)}=v^{(k-1)}$, we can avoid the inversion of $A$.

If the inverse power method is applied to the shifted matrix $A-\mu I$, then
the eigenvalues of $\left(  A-\mu I\right)  ^{-1}$ are $\left(  \lambda
_{i}-\mu\right)  ^{-1}$. If $\mu$ approaches, say, to $\lambda_{t}$, then
$\lambda_{i}-\mu\rightarrow\lambda_{i}-\lambda_{t}$. Hence the inequality
\[
\left|  \lambda_{t}-\mu\right|  ^{-1}>\left|  \lambda_{i}-\mu\right|
^{-1}\quad\left(  i\neq t\right)
\]
holds for the eigenvalues of the shifted matrix. The speed of convergence is
determined by the quotient
\[
q=\left|  \lambda_{t}-\mu\right|  /\left\{  \max\left|  \lambda_{i}%
-\mu\right|  \right\}  \koz .
\]
If $\mu$ is close enough to $\lambda_{t}$, then $q$ is very small and the
inverse power iteration converges very fast. This property can be exploited in
the calculation of approximate eigenvectors if an approximate eigenvalue, say
$\mu$, is known. Assuming that $\det\left(  A-\mu I\right)  \neq0$, we apply
the inverse power method to the shifted matrix $A-\mu I$. In spite of the fact
that matrix $A-\mu I$ is nearly singular and the linear equation $\left(
A-\mu I\right)  z^{\left(  k\right)  }=v^{\left(  k\right)  }$ cannot be
solved with high precision, the algorithm gives very often good approximations
of the eigenvectors.

Finally we note that in principle the von Mieses method
\index{von Mieses method}
\nevindex{von Mieses, Richard (1853--1953)}
can be modified to
determine all eigenvalues and eigenvectors.

\subsubsection{Orthogonalisation processes.}
We need the following definition and theorem.

\begin{defi}
The matrix $Q\in\mathbb{R}^{n\times n}$  is said to be
\ki{orthogonal} \inddef{orthogonal} if
$Q^{T}Q=I$.
\end{defi}

\begin{tetel}[$QR$-decomposition]
Every matrix $A\in\mathbb{R}^{n\times m}$ having
linearly independent \ column vectors can be decomposed in the product form
$A=QR$, where $Q$ is orthogonal and $R$ is upper triangular matrix. \ 
\end{tetel}

We note that the $QR$-decomposition can be applied for solving linear systems
of equations, similarly to the $LU$-decomposition. If the $QR$-decomposition
of $A$ is known, then the equation $Ax=QRx=b$ can be written in the equivalent
form $Rx=Q^{T}b$. Thus we have to solve only an upper triangular linear system.

There are several methods to determine the $QR$-decomposition of a matrix. In
practice the Givens-, the Householder- 
\nevindex{Givens, Wallace J.}\nevindex{Householder, Alston Scott}
and the MGS-methods are used.

The MGS (Modified Gram-Schmidt)
\nevindex{Gram, Jorgen Pedersen (1850--1916)}\nevindex{Schmidt, Erhard (1876--1959)}
method is a stabilised, but algebraically
equivalent version of the classical Gram-Schmidt orthogonalisation algorithm.
The basic problem is the following: We seek for an 
\emph{orthonormal basis} $\left\{
q_{j}\right\}  _{j=1}^{m}$ of the subspace
\[
\mathcal{L}\left\{  a_{1},\ldots,a_{m}\right\}  =\left\{  \sum_{j=1}%
^{m}\lambda_{j}a_{j}\,\,|\ \lambda_{j}\in\mathbb{R},\ j=1,\ldots,m\right\}  \koz ,
\]
where $a_{1},\ldots,a_{m}\in\mathbb{R}^{n}$ ($m\leq n$) are linearly
independent vectors. That is we determine the linearly independent vectors
$q_{1},\ldots,q_{m}$ such that
\[
q_{i}^{T}q_{j}=0\ \left(  i\neq j\right)  ,\left\|  q_{i}\right\|
_{2}=1\ \left(  i=1,\ldots,m\right)
\]
and
\[
\mathcal{L}\left\{  a_{1},\ldots,a_{m}\right\}  =\mathcal{L}\left\{
q_{1},\ldots,q_{m}\right\}  \koz .
\]
The basic idea of the \emph{classical Gram-Schmidt-method}
\index{classical Gram-Schmidt-method} is the following:

Let $r_{11}=\left\|  a_{1}\right\|  _{2}$ and $q_{1}=a_{1}/r_{11}$. Assume
that vectors $q_{1},\ldots,q_{k-1}$ are already computed and orthonormal.
Assume that vector $\tilde{q}_{k}=a_{k}-\sum_{j=1}^{k-1}r_{jk}q_{j}$ is such
that $\tilde{q}_{k}\,\bot\,q_{i}$, that is $\tilde{q}_{k}^{T}q_{i}=a_{k}%
^{T}q_{i}-\sum_{j=1}^{k-1}r_{jk}q_{j}^{T}q_{i}=0$ holds for $i=1,\ldots,k-1$.
Since $q_{1},\ldots,q_{k-1}$ are orthonormal, $q_{j}^{T}q_{i}=0$ ($i\neq j$)
and $r_{ik}=a_{k}^{T}q_{i}\quad\left(  i=1,\ldots,k-1\right)  $. After
normalisation we obtain $q_{k}=\tilde{q}_{k}/\left\|  \tilde{q}_{k}\right\|
_{2}$.

The algorithm is formalised as follows.

\begin{alg}{CGS-Orthogonalization($m,a_{1},\ldots ,a_{m}$)}
\inddef{\textsc{CGS-Orthogonalization}}
1 \' \key{for} \= $k\leftarrow 1$ \= \key{to} \= $m$ \\
2 \' \> \key{do} \= \key{for} \= $i\leftarrow 1$ \= \key{to} \= $k-1$ \\
3 \'  \> \> \> \key{do} \= $r_{ik}\leftarrow a_{k}^{T}a_{i}$ \\
4 \' \> \> \> \> $a_{k}\leftarrow a_{k}-r_{ik}a_{i}$ \\
5 \' \> \> $r_{kk}\leftarrow \left\| a_{k}\right\| _{2}$ \\
6 \' \> \> $a_{k}\leftarrow a_{k}/r_{kk}$ \\
7 \' \key{return} $a_{1},\ldots ,a_{m}$
\end{alg}


The algorithm overwrites vectors $a_{i}$ by the orthonormal vectors $q_{i}$.
The connection with the $QR$-decomposition follows from the relation
$a_{k}=\sum_{j=1}^{k-1}r_{jk}q_{j}+r_{kk}q_{k}$. Since
\[%
\begin{array}
[c]{rcl}%
a_{1} & = & q_{1}r_{11},\\
a_{2} & = & q_{1}r_{12}+q_{2}r_{22},\\
a_{3} & = & q_{1}r_{13}+q_{2}r_{23}+q_{3}r_{33,}\\
& \vdots & \\
a_{m} & = & q_{1}r_{1m}+q_{2}r_{2m}+\ldots+q_{m}r_{mm} \koz ,
\end{array}
\]
we can write that%

\[
A=\left[  a_{1},\ldots,a_{m}\right]  =\underset{Q}{\underbrace{\left[
q_{1},\ldots,q_{m}\right]  }}\underset{R}{\underbrace{\left[
\begin{array}
[c]{ccccc}%
r_{11} & r_{12} & r_{13} & \ldots & r_{1m}\\
0 & r_{22} & r_{23} & \ldots & r_{2m}\\
0 & 0 & r_{33} & \ldots & r_{3m}\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \ldots & r_{mm}%
\end{array}
\right]  }}=QR\ .
\]


The \emph{numerically stable MGS method}
\index{numerically stable MGS method}
is given in the following form

\nevindex{Gram, Jorgen Pedersen (1850--1916)} \nevindex{Schmidt, Erhard (1876--1959)}
\begin{alg}{MGS-Orthogonalisation($m,a_{1},\ldots ,a_{m}$)}
\inddef{MGS-Orthogonalisation@\textsc{MGS-Orthogonalisation}}
1 \' \key{for} \= $k\leftarrow 1$ \= \key{to} \= $m$ \\
2 \' \> \key{do} \= $r_{kk}\leftarrow \left\| a_{k}\right\| _{2}$ \\
3 \' \> \> $a_{k}\leftarrow a_{k}/r_{kk}$ \\
4 \' \> \> \key{for} \= $j\leftarrow k+1$ \= \key{to} \= $m$ \\
5 \' \> \> \> \key{do} \= $r_{kj}\leftarrow a_{j}^{T}a_{k}$ \\
6 \' \> \> \> \> $a_{j}\leftarrow a_{j}-r_{kj}a_{k}$ \\
7 \' \key{return} $a_{1},\ldots ,a_{m}$
\end{alg}

The algorithm overwrites vectors $a_{i}$ by the orthonormal vectors $q_{i}$.
The MGS method is more stable than the CGS algorithm. Björck
\nevindex{Björck, Ake}proved that
for $m=n$ the computed matrix $\hat{Q}$ satisfies
\[
\hat{Q}^{T}\hat{Q}=I+E,\qquad\left\|  E\right\|  _{2}\cong\text{cond}\left(
A\right)  u \koz ,
\]
where $u$ is the unit roundoff.

\paragraph{The $QR$-method.}
Today the $QR$-method\index{QR-method@$QR$-method} is the most important numerical algorithm to compute all
eigenvalues of a general matrix. It can be shown that the $QR$-method is a
generalisation of the power method. The basic idea of the method is the
following: Starting from $A_{1}=A$ we compute the sequence $A_{k+1}=Q_{k}%
^{-1}A_{k}Q_{k}=Q_{k}^{T}A_{k}Q_{k}$, where $Q_{k}$ is orthogonal, $A_{k+1}$
is orthogonally similar to $A_{k}$ ($A$) and the lower triangular part of
$A_{k}$ tends to a diagonal matrix, whose entries will be the eigenvalues of
$A$. Here $Q_{k}$ is the orthogonal factor of the $QR$-decomposition
$A_{k}=Q_{k}R_{k}$. Therefore $A_{k+1}=Q_{k}^{T}(Q_{k}R_{k})Q_{k}=R_{k}Q_{k}$.
The basic algorithm is given in the following form.

\begin{alg}{\textsc{QR-Method}($A$)}\inddef{QR-method@\textsc{QR-Method}}
1 \' $k\leftarrow 1$ \\
2 \' $A_{1}\leftarrow A$ \\
3 \' \key{while} \= \emph{exit condition} = \textsc{false} \\
4 \'             \> \key{do} \= Compute the $QR$-decomposition $A_{k}=Q_{k}R_{k}$ \\
5 \'             \>          \> $A_{k+1}\leftarrow R_{k}Q_{k}$ \\
6 \'             \>          \> $k\leftarrow k+1$ \\
7 \' \key{return} $A_{k}$
\end{alg}


The following result holds.

\begin{tetel}[Parlett]
\label{parltetel}
\nevindex{Parlett, Beresford} If the matrix $A$ is diagonalisable, $X^{-1}%
AX=$diag$(\lambda_{1},\lambda_{2},\ldots,\lambda_{n})$, the eigenvalues
satisfy
\[
\left|  \lambda_{1}\right|  >\left|  \lambda_{2}\right|  >\ldots>\left|
\lambda_{n}\right|  >0
\]
and $X$ has an $LU$-decomposition, then the lower triangular part of $A_{k}$
converges to a diagonal matrix whose entries are the eigenvalues of $A$.
\end{tetel}

In general, matrices $A_{k}$ do not necessarily converge to a given matrix. If
$A$ has $p$ eigenvalues of the same modulus, the form of matrices $A_{k}$
converge to the form
\begin{equation}
\left[
\begin{array}
[c]{ccccccccc}%
\times &  &  &  &  &  &  &  & \times\\
0 & \ddots &  &  &  &  &  &  & \\
&  & \times &  &  &  &  &  & \\
0 &  & 0 & \ast & \cdots & \ast &  &  & \\
&  &  & \vdots &  & \vdots &  &  & \\
&  &  & \ast & \cdots & \ast &  &  & \\
0 &  &  &  &  & 0 & \times &  & \\
&  &  &  &  &  &  & \ddots & \\
0 &  &  &  &  &  &  & 0 & \times
\end{array}
\right]  , \label{qra}%
\end{equation}
where the entries of the submatrix denoted by $\ast$ do not converge. However
the eigenvalues of this submatrix will converge. This submatrix can be
identified and properly handled. A real matrix may have real and complex
eigenvalues. If there is a complex eigenvalues, than there is a corresponding
conjugate eigenvalue as well. For pairs of complex conjugated eigenvalues $p$
is at least $2$. Hence the sequence $A_{k}$ will show this phenomenon .

The $QR$-decomposition is very expensive. Its cost is $\Theta(n^{3})$ flops
for general $n\times n$ matrices. If $A$ has upper Hessenberg form, the cost
of $QR$-decomposition is $\Theta(n^{2})$ flops.

\begin{defi}
The matrix $A\in\mathbb{R}^{n\times n}$\ has 
\ki{upper Hessenberg form,}\inddef{upper Hessenberg form} 
\nevindex{Hessenberg, Gerhard}if
\[
A=\left[
\begin{array}
[c]{cccccc}%
a_{11} &  &  & \ldots &  & a_{1n}\\
a_{21} &  &  &  &  & \\
0 & a_{32} &  &  &  & \vdots\\
\vdots & 0 & \ddots &  &  & \\
&  & \ddots & a_{n-1,n-2} & a_{n-1,n-1} & a_{n-1,n}\\
0 &  & \ldots & 0 & a_{n,n-1} & a_{nn}%
\end{array}
\right]  .
\]

\end{defi}

The following theorem guarantees that if $A$ has upper Hessenberg form, then
every $A_{k}$ of the $QR$-method has also upper Hessenberg form.

\begin{tetel}
If $A$ has upper Hessenberg form and $A=QR$, then $RQ$ has also upper
Hessenberg form.
\end{tetel}

We can transform a matrix $A$ to a similar matrix of upper Hessenberg form in
many ways. One of the cheapest ways, that costs about $5/6n^{3}$ flops, is
based on the Gauss elimination method.
\nevindex{Gauss, Johann Carl Friedrich (1777--1855)} Considering the advantages of the upper
Hessenberg form
\nevindex{Hessenberg, Gerhard} the efficient implementation of the $QR$-method requires first
the similarity transformation of $A$ to upper Hessenberg form.

The convergence of the $QR$-method, similarly to the power method, depends on
the quotients $\left|  \lambda_{i+1}/\lambda_{i}\right|  $. The eigenvalues of
the shifted matrix $A-\sigma I$ are $\lambda_{1}-\sigma,\lambda_{2}%
-\sigma,\ldots,\lambda_{n}-\sigma$. The corresponding eigenvalue ratios are
$\left|  \left(  \lambda_{i+1}-\sigma\right)  /\left(  \lambda_{i}%
-\sigma\right)  \right|  $. A proper selection of $\sigma$ can fasten the convergence.

The usual form of the $QR$-method includes the transformation to upper
Hessenberg form and the shifting.

\begin{alg}{Shifted-$QR$-Method($A$)}\inddef{Shifter-QR-Method\textsc{Shifted-QR-Method}}
1 \' $H_{1}\leftarrow U^{-1}AU\quad$($H_{1}$
is of upper Hessenberg form) \\
2 \' $k\leftarrow 1$ \\
3 \' \key{while} \= \emph{exit condition} = \textsc{false} \\
4 \' \> \key{do} \= compute the $QR$-decomposition
$H_{k}-\sigma_{k}I=Q_{k}R_{k}$ \\
5 \' \> \> $H_{k+1}\leftarrow R_{k}Q_{k}+\sigma _{k}I$ \\
6 \' \> \> $k\leftarrow k+1$ \\
7 \' \key{return} $H_{k}$
\end{alg}

In practice the $QR$-method is used in shifted form. There are various
strategies to select $\sigma_{i}$. The most often used selection is given by
$\sigma_{k}=h_{nn}^{\left(  k\right)  }\quad\left(  H_{k}=\left[
h_{ij}^{\left(  k\right)  }\right]  _{i,j=1}^{n}\right)  $.

The eigenvectors of $A$ can also be determined by the $QR$-method. For this we
refer to the literature.

\begin{gyak}
\ujgyak Apply the \emph{power method} \gyakindex{power method}
 to the matrix $A=\left[
\begin{array}
[c]{cc}%
1 & 1\\
0 & 2
\end{array}
\right]  $ with the initial vector $v^{(0)}=\left[
\begin{array}
[c]{c}%
1\\
1
\end{array}
\right]  $. What is the result of the $20$th step?
\ujgyak Apply the \emph{power method},
\gyakindex{power method}
the \emph{inverse power method}\gyakindex{inverse power method}
and the $QR$-method \gyakindex{QR-method@$QR$-method}
to the
matrix%
\[
\left[
\begin{array}
[c]{rrr}%
-4 & -3 & -7\\
2 & 3 & 2\\
4 & 2 & 7
\end{array}
\right]  \koz .
\]
\ujgyak Apply the shifted $QR$-method\gyakindex{shifted $QR$-method}
to the matrix of the previous exercise with the
choice $\sigma_{i}=\sigma$ ($\sigma$ is fixed).
\end{gyak}

\section{Numerical program libraries and software tools}
We have plenty of devices and tools that support efficient coding and
implementation of numerical algorithms. One aim of such developments is to
free the programmers from writing the programs of frequently occurring
problems. This is usually done by writing safe, reliable and standardised
routines that can be downloaded from (public) program libraries. We just
mention the LINPACK, EISPACK, LAPACK, VISUAL\ NUMERICS (former IMSL) and NAG
libraries. Another way of developments is to produce software that work as a
programming language and makes the programming very easy. Such software
systems are the MATLAB and the SciLab.

\subsection{Standard linear algebra subroutines}
The main purpose of the BLAS (Basic Linear Algebra Subprograms) programs is
the standardisation and efficient implementation the most frequent
matrix-vector operations. Although the BLAS routines were published in FORTRAN
they can be accessed in optimised machine code form as well. The BLAS routines
have three levels:

- BLAS 1 (1979),

- BLAS 2 (1988),

- BLAS 3 (1989).

These levels corresponds to the computation cost of the implemented matrix
operations. The BLAS routines are considered as the best implementations of
the given matrix operations. The selection of the levels and individual BLAS
routines strongly influence the efficiency of the program. A sparse version of
BLAS also exists.

We note that the BLAS 3 routines were developed mainly for block parallel
algorithms. The standard linear algebra packages LINPACK, EISPACK and LAPACK
are built from BLAS routines. The parallel versions can be found in the
SCALAPACK package. These programs can be found in the public NETLIB library:

\begin{center}
\texttt{http:/www.netlib.org/index.html}
\end{center}

\paragraph{BLAS 1 routines.}
Let $\alpha\in\mathbb{R}$, $x,y\in\mathbb{R}^{n}$. The BLAS 1 routines are the
programs of the most important vector operations ($z=\alpha x$, $z=x+y$,
$dot=x^{T}y$), the computation of $\left\|  x\right\|  _{2}$, the swapping of
variables, rotations and the \emph{saxpy} operation which is defined by
\[
z=\alpha x+y \koz .
\]
The word \emph{saxpy} means that \idez{scalar alpha $x$ plus $y$}. The
\emph{saxpy} operation is implemented in the following way.

\begin{alg}{Saxpy$\left(\alpha,x,y\right)$}\inddef{\textsc{Saxpy}}
1 \'  $n\leftarrow \textit{elements}\left[ x\right] $ \\
2 \' \key{for} \= $i\leftarrow 1$ to $n$ \\
3 \' \> \key{do} \= $z\left[ i\right] =\alpha x\left[ i\right] +y\left[ i\right] $ \\
4 \' \key{return} $z$
\end{alg}


The \emph{saxpy} is a software driven operation. The cost of BLAS 1 routines
is $\Theta\left(  n\right)  $ flops.

\paragraph{BLAS 2 routines.}
The matrix-vector operations of BLAS 2 requires $\Theta\left(  n^{2}\right)  $
flops. These operations are $y=\alpha Ax+\beta y$, $y=Ax$, $y=A^{-1}x$,
$y=A^{T}x$, $A\leftarrow A+xy^{T}$ and their variants. Certain operations work
only with triangular matrices. We analyse two operations in detail. 
The \idez{outer or dyadic product} update
\[
A\leftarrow A+xy^{T}\quad\left(  A\in\mathbb{R}^{m\times n}\text{, }%
x\in\mathbb{R}^{m}\text{, }y\in\mathbb{R}^{n}\right)
\]
can be implemented in two ways.

The rowwise or \idez{$ij$} variant:

\begin{alg}{Outer-Product-Update-Version''ij`` $(A,x,y)$}
\inddef{\textsc{Outer-Product-Update-Version''ij``}}
1 \' $m\leftarrow \textit{rows}[A]$ \\
2 \'  \key{for} \= $i\leftarrow 1$ \= \key{to} \= $m$ \\
3 \'  \> \key{do} \=  $A\left[ i,:\right] \leftarrow 
A\left[ i,:\right] +x\left[ i \right] y^{T}$ \\
4 \'  \key{return} $A$
\end{alg}


The notation \idez{:} \ denotes all allowed indices. In our case this means the
indices $1\leq j\leq n$. Thus $A\left[  i,:\right]  $ denotes the
$i^{\text{th}}$ row of matrix $A$.

The columnwise or \idez{$ji$} variant\emph{:}

\begin{alg}{Outer-Product-Update-Version''ji``$(A,x,y)$}
\inddef{\textsc{Outer-Product-Update-Version''ji``}}
1 \' $n\leftarrow \textit{columns}[A]$ \\
2 \'  \key{for} \= $j\leftarrow 1$ \= \key{to} $n$ \\
3 \'  \> \key{do} \= $A\left[ :,j\right] \leftarrow A\left[ :,j\right] +
y\left[ j\right]x$ \\
4 \'  \key{return} $A$
\end{alg}


Here $A\left[  :,j\right]  $ denotes the $j^{\text{th}}$ column of matrix $A$.
Observe that both variants are based on the saxpy operation.

The \emph{gaxpy} operation is defined by
\[
z=y+Ax\quad\left(  x\in\mathbb{R}^{n}\text{, }y\in\mathbb{R}^{m}\text{, }%
A\in\mathbb{R}^{m\times n}\right)  \koz .
\]
The word \emph{gaxpy} means that \idez{general $A$ $x$ plus $y$}.
The \emph{gaxpy} operation
is also software driven and implemented in the following way:

\begin{alg}{Gaxpy$(A,x,y)$}\inddef{\textsc{Gaxpy}}
1 \' $n\leftarrow \textit{columns}[A]$ \\
2 \'  $z\leftarrow y$ \\
3 \'  \key{for} \= $j\leftarrow 1$ \= \key{to} \=  $n$ \\
4 \'  \> \key{do} \= $z\leftarrow z+x\left[ j\right] A\left[ :,j\right] $ \\
5 \'  \key{return} $z$
\end{alg}


Observe that the computation is done columnwise and the \emph{gaxpy}
operation is
essentially a generalised \emph{saxpy}.

\paragraph{BLAS 3 routines.}
These routines are the implementations of $\Theta\left(  n^{3}\right)  $
matrix-matrix and matrix-vector operations such as the operations
$C\leftarrow\alpha AB+\beta C$, $C\leftarrow\alpha AB^{T}+\beta C$,
$B\leftarrow\alpha T^{-1}B$ ($T$ \ is upper triangular) and their variants.
BLAS 3 operations can be implemented in several forms. For example, the matrix
product $C=AB$ can be implemented at least in three ways. Let $A\in
\mathbb{R}^{m\times r}$, $B\in\mathbb{R}^{r\times n}$.

\begin{alg}{Matrix-Product-Dot-Version($A,B$)}
\inddef{Matrix-Product-Dot-Version@\textsc{Matrix-Product-Dot-Version}}
1 \' $m\leftarrow \textit{rows}[A]$ \\
2 \' $r\leftarrow \textit{columns}[A]$ \\
3 \' $n\leftarrow \textit{columns}[B]$ \\
4 \'  $C\left[ 1:m,1:n\right] \leftarrow 0$ \\
5 \' \key{for} \= $i\leftarrow 1$ \= \key{to} \= $m$ \\
6 \' \> \key{do} \= \key{for} \= $j\leftarrow 1$ \= \key{to} \= $n$ \\
7 \' \> \> \> \key{do} \= \key{for} \= $k\leftarrow 1$ \=\key{to} \= $r$ \\
8 \' \> \> \> \> \> \key{do} \= $C\left[ i,j\right] \leftarrow
C\left[ i,j\right] +A\left[ i,k\right] B\left[k,j\right]$ \\
9 \' \key{return} $C$
\end{alg}

This algorithm computes $c_{ij}$ as the dot (inner) product of the
$i^{\text{th}}$ row of $A$ and the $j^{\text{th}}$ column of $B$. This
corresponds to the original definition of matrix products.

Now let $A,B$ and $C$ be partitioned columnwise as follows%
\[%
\begin{array}
[c]{c}%
A=\left[  a_{1},\ldots,a_{r}\right]  \quad\left(  a_{i}\in\mathbb{R}%
^{m}\right)  ,\\
B=\left[  b_{1},\ldots,b_{n}\right]  \quad\left(  b_{i}\in\mathbb{R}%
^{r}\right)  ,\\
C=\left[  c_{1},\ldots,c_{n}\right]  \quad\left(  c_{i}\in\mathbb{R}%
^{m}\right)  .
\end{array}
\]
Then we can write $c_{j}$ as the linear combination of the columns of $A$,
that is
\[
c_{j}=\sum_{k=1}^{r}b_{kj}a_{k}\quad\left(  j=1,\ldots,n\right)  .
\]
Hence the product can be implemented with \emph{saxpy} operations.

\begin{alg}{Matrix-Product-Gaxpy-Variant($A,B$)}
\inddef{\textsc{Matrix-Product-Gaxpy-Variant}}
1 \' $m\leftarrow \textit{rows}[A]$ \\
2 \' $r\leftarrow \textit{columns}[A]$ \\
3 \' $n\leftarrow \textit{columns}[B]$ \\
4 \' $C\left[ 1:m,1:n\right] \leftarrow 0$ \\
5 \' \key{for} \= $j\leftarrow 1$ \= \key{to} \= $n$ \\
6 \' \> \key{do} \= \key{for} \= $k\leftarrow 1$ \= \key{to} \= $r$ \\
7 \' \> \> \> \key{do} \= \key{for} \= $i\leftarrow 1$ \=\key{to} \= $m$ \\
8 \' \> \> \> \> \> \key{do} \= $C\left[ i,j\right] \leftarrow
C\left[ i,j\right] +A\left[ i,k\right] B\left[k,j\right]$ \\
9 \' \key{return} $C$
\end{alg}


The following equivalent form of the \idez{$jki$}-algorithm shows that it is indeed a
\emph{gaxpy} based process.

\begin{alg}{Matrix-Product-with-Gaxpy-Call($A,B$)}
\inddef{\textsc{Matrix-Product-with-Gaxpy-Call}}
1 \' $m\leftarrow \textit{rows}[A]$ \\
2 \' $n\leftarrow \textit{columns}[B]$ \\
3 \' $C\left[ 1:m,1:n\right]\leftarrow 0$ \\
4 \' \key{for} \= $j\leftarrow 1$ \= \key{to} \= $n$ \\
5 \' \> \key{do} \= $C\left[ :,j\right] =
gaxpy\left( A,B\left[ :,j\right] ,C\left[ :,j\right] \right)$  \\
6 \' \key{return} $C$
\end{alg}


Consider now the partitions $A=\left[  a_{1},\ldots,a_{r}\right]  $ ($a_{i}\in
R^{m}$) and
\[
B=\left[
\begin{array}
[c]{c}%
b_{1}^{T}\\
\vdots\\
b_{r}^{T}%
\end{array}
\right]  \quad\left(  b_{i}\in R^{n}\right)  .
\]
Then $C=AB=\sum_{k=1}^{r}a_{k}b_{k}^{T}.$

\begin{alg}{Matrix-Product-Outer-Product-Variant($A,B$)}\inddef{\textsc{Matrix-Product-Outer-Product-Variant}}
1 \' $m\leftarrow \textit{rows}[A]$ \\
2 \' $r\leftarrow \textit{columns}[A]$ \\
3 \' $n\leftarrow \textit{columns}[B]$ \\
4 \'  $C\left[ 1:m,1:n\right] =0$ \\
5 \' \key{for} \= $k\leftarrow 1$ \= \key{to} \= $r$ \\
6 \' \>\key{do} \= \key{for} \= $j\leftarrow 1$ \= \key{to} \= $n$ \\
7 \'  \> \> \> \key{do} \= \key{for} \= $i\leftarrow 1$ \=\key{to} \= $m$ \\
8 \'  \> \> \>          \>           \> \key{do} $C\left[ i,j\right] \leftarrow 
C\left[ i,j\right] +A\left[ i,k\right] B\left[k,j\right] $ \\
9 \' \key{return} $C$
\end{alg}

The inner loop realizes a \emph{saxpy} operation: it gives the multiple of
$a_{k}$ to the $j^{\text{th}}$ column of matrix $C$.

\subsection{Mathematical software}
These are those programming tools that help easy programming in concise
(possibly mathematical) form within an integrated program development system.
Such systems were developed primarily for solving mathematical problems. By
now they have been extended so that they can be applied in many other fields.
For example, Nokia uses MATLAB in the testing and quality control of mobile
phones. We give a short review on MATLAB in the next section. We also mention
the widely used MAPLE, DERIVE and MATEMATICA systems.

\paragraph{The MATLAB system.}
The MATLAB software was named after the expression MATrix LABoratory. The
name indicates that the matrix operations are very easy to make. The initial
versions of MATLAB had only one data type: the complex matrix. In the later
versions high dimension arrays, cells, records and objects also appeared. The
MATLAB can be learned quite easily and even a beginner can write programs for
relatively complicated problems.

The coding of matrix operations is similar to their standard mathematical
form. For example if $A$ and $B$ are two matrices of the same size, then their
sum is given by the command $C=A+B$. As a programming language the MATLAB
contains only four control structures known from other programming languages:

-- the simple statement $Z=$\emph{expression},
-- the if statement of the form

$\qquad\qquad\mathbf{if}$ \emph{expression, commands} \{$\mathbf{else/elseif}$
\emph{commands}\} $\mathbf{end}$,

-- the for loop of the form

\qquad\qquad$\mathbf{for}$ \emph{the values of the loop variable, commands}
$\mathbf{end}$

-- the while loop of the form

\qquad\qquad$\mathbf{while}$ \emph{expression, commands} $\mathbf{end}$.

\noindent The MATLAB has an extremely large number of built in functions that
help efficient programming. We mention the following ones as a sample.

-- $\max(A)$ selects the maximum element in every column of $A$,

-- $[v,s]=$eig$(A)$ returns the approximate eigenvalues and eigenvectors of
$A$,

-- The command $A\backslash b$ returns the numerical solution of the linear
system $Ax=b$.

The entrywise operations and partitioning of matrices can be done very
efficiently in MATLAB. For example, the statement%
\[
A([2,3],:)=1./A([3,2],:)
\]
exchange the second and third rows of $A$ while it takes the reciprocal of
each element.

The above examples only illustrate the possibilities and easy programming of
MATLAB. These examples require much more programming effort in other
languages, say e.g. in PASCAL. The built in functions of MATLAB can be easily
supplemented by other programs.

The higher number versions of MATLAB include more and more functions and
special libraries (tool boxes) to solve special problems such as optimisation,
statistics and so on.

There is a built in automatic technique to store and handle sparse matrices
that makes the MATLAB competitive in solving large computational problems. The
recent versions of MATLAB offer very rich graphic capabilities as well. There
is an extra interval arithmetic package that can be downloaded from the WEB site

\begin{center}
\texttt{http:/www.ti3.tu-harburg.de%
%TCIMACRO{\TEXTsymbol{\backslash}}%
%BeginExpansion
$\backslash$%
%EndExpansion
\%7Erump%
%TCIMACRO{\TEXTsymbol{\backslash}}%
%BeginExpansion
$\backslash$%
%EndExpansion
intlab}
\end{center}

\noindent There is a possibility to build certain C and FORTRAN programs into
MATLAB. Finally we mention that the system has an extremely well written help system.


\begin{fld}
\ujfld{Without overflow}
{Write a MATLAB program that computes the norm $\left\|
x\right\|  _{2}=\left(  \sum_{i=1}^{n}x_{i}^{2}\right)  ^{1/2}$ without
\emph{overflow} \felindex{overflow}
 in all cases when the result does not make overflow. It is also
required that the error of the final result can not be greater than that of
the original formula.}
\ujfld{Estimate}
{Equation $x^{3}-3.330000x^{2}+3.686300x-1.356531=0$ has the
solution $x_{1}=1.01.$ The perturbed equation $x^{3}-3.3300x^{2}%
+3.6863x-1.3565=0$ has the solutions $y_{1},\,\,y_{2},\,\,y_{3}$. Give an
estimate for the perturbation $\min_{i}\left|  x_{1}-y_{i}\right|  $.}
\ujfld{Double word length}
{Consider an arithmetic system that has double word
length such that every number represented with $2t$ digits are stored in two
$t$ digit word. Assume that the computer can only add numbers with $t$ digits.
Furthermore assume that the machine can recognise overflow. \newline(i) Find
an algorithm that add two positive numbers of $2t$ digit length. \newline(ii)
If the representation of numbers requires the sign digit for all numbers, then
modify algorithm (i) so that it can add negative and positive numbers both of
the same sign. We can assume that the sum does not overflow.}
\ujfld{Auchmuty theorem}
{Write a MATLAB program for the Auchmuty\nevindex{Auchmuty, Giles} error estimate
(see Theorem \ref{auchm}) and perform the following numerical testing.\newline%
(i) Solve the linear systems $Ax=b_{i}$, where $A\in\mathbb{R}^{n\times n}$ is
a given matrix, $b_{i}=Ay_{i}$, $y_{i}\in\mathbb{R}^{n}$ ($i=1,\ldots,N$) are
random vectors such that $\left\|  y_{i}\right\|  _{\infty}\leq\beta$. Compare
the true errors $\left\|  \widetilde{x}_{i}-y_{i}\right\|$ ($i=1,\ldots,N$),
and the estimated errors $EST_{i}=\left\|  r(\widetilde{x_{i}})\right\|
_{2}^{2}/\left\|  A^{T}r(\widetilde{x}_{i})\right\|  _{2}$, where
$\widetilde{x}_{i}$ is the approximate solution of $Ax=b_{i}$. What is the
minimum, maximum and average of numbers $c_{i}$? Use graphic for the
presentation of the results. Suggested values are $n\leq200$, $\beta=200$ and
$N=40$. \newline(ii) Analyse the effect of
condition number\felindex{condition number} and size.
\newline(iii) Repeat problems (i) and (ii) using LINPACK and BLAS.}
\ujfld{Hilbert matrix}
{Consider the linear system $Ax=b$, where $b=[1,1,1,1]^{T}$
and $A$ is the fourth order Hilbert matrix,
\nevindex{Hilbert, David}
that is $a_{i,j}=1/(i+j).$ $A$ is
ill-conditioned. The inverse of $A$ is approximated by
\newpage
\[
B=\left[
\begin{array}
[c]{rrrr}%
202 & -1212 & 2121 & -1131\\
-1212 & 8181 & -15271 & 8484\\
2121 & -15271 & 29694 & -16968\\
-1131 & 8484 & -16968 & 9898
\end{array}
\right]  \koz .
\]
Thus an $x_{0}$ approximation of the true solution $x$ is given by $x_{0}=Bb.$
Although the true solution is also integer $x_{0}$ is not an acceptable
approximation. Apply the \emph{iterative refinement}
\felindex{iterative refinement}
with $B$ instead of $A^{-1}$ to
find an acceptable integer solution.}
\ujfld{Consistent norm}
{Let $\left\|  A\right\|  $ be a \emph{consistent norm}
\felindex{consistent norm}
and consider the linear system $Ax=b$\newline(i) Prove that if $A+\Delta A$ is
singular, then $\text{cond}(A)\geq\left\|  A\right\|  /\left\|  \Delta
A\right\|  $. \newline(ii) Show that for the {2}-norm equality holds in (i),
if $\Delta A=-bx^{T}/(b^{t}x)$ and $\left\|  A^{-1}\right\|  _{2}\left\|
b\right\|  _{2}=\left\|  A^{-1}b\right\|  _{2}$. \newline(iii) Using the
result of (i) give a lower bound to $\text{cond}_{\infty}(A)$, if
\[
A=\left[
\begin{array}
[c]{rrr}%
1 & -1 & 1\\
-1 & \varepsilon & \varepsilon\\
1 & \varepsilon & \varepsilon
\end{array}
\right] \koz .
\]}
\ujfld{Cholesky-method}
{Use the Cholesky-method\nevindex{Cholesky, André-Louis (1875--1918)} to solve the linear system $Ax=b$,
where
\[
A=\left[
\begin{array}
[c]{cccccc}%
5.5 & 0 & 0 & 0 & 0 & 3.5\\
0 & 5.5 & 0 & 0 & 0 & 1.5\\
0 & 0 & 6.25 & 0 & 3.75 & 0\\
0 & 0 & 0 & 5.5 & 0 & 0.5\\
0 & 0 & 3.75 & 0 & 6.25 & 0\\
3.5 & 1.5 & 0 & 0.5 & 0 & 5.5
\end{array}
\right]  ,\,\,b=\left[
\begin{array}
[c]{c}%
1\\
1\\
1\\
1\\
1\\
1
\end{array}
\right]  \koz .
\]
Also give the exact \emph{Cholesky-decomposition}
\felindex{Cholesky-decomposition}
\nevindex{Cholesky, André-Louis (1875--1918)} $A=LL^{T}$  and the true solution
of $Ax=b$. The approximate Cholesky-factor\nevindex{Cholesky, André-Louis (1875--1918)} 
$\widetilde{L}$ satisfies the relation $\tilde{L}\tilde{L}^{T}=A+F.$ It can proved that in a
\emph{floating point arithmetic}\felindex{floating point arithmetic}
with $t$-digit mantissa and base $\beta$ the entries of $F$ satisfy
the inequality $\left|  f_{i,j}\right|  \leq e_{i,j}$, where
\[
E=\beta^{1-t}\left[
\begin{array}
[c]{cccccc}%
11 & 0 & 0 & 0 & 0 & 3.5\\
0 & 11 & 0 & 0 & 0 & 1.5\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 11 & 0 & 0.5\\
0 & 0 & 0 & 0 & 0 & 0\\
3.5 & 1.5 & 0 & 0.5 & 0 & 11
\end{array}
\right]  \koz .
\]
Give a bound for the relative error of the approximate solution $\tilde{x}$,
if $\beta=16$ and $t=14$ (IBM3033).}
\newpage
\ujfld{Bauer-Fike theorem}\nevindex{Bauer, F. L.}
{Let
\[
A=\left[
\begin{array}
[c]{cccccc}%
10 & 10 &  &  &  & \\
& 9 & 10 &  &  & \\
&  & 8 & 10 &  & \\
&  &  & \ddots & \ddots & \\
&  &  &  & 2 & 10\\
\varepsilon &  &  &  &  & 1
\end{array}
\right]
\]
(i) Analyze the perturbation of the eigenvalues for $\varepsilon
=10^{-5},10^{-6},10^{-7},0$. \newline(ii) Compare the estimate of 
Bauer-Fike\nevindex{Bauer, F. L.}
theorem to the matrix $A=A(0)$. }
\ujfld{Eigenvalues}
{Using the MATLAB eig routine compute the eigenvalues\felindex{eigenvalue}
of$B=AA^{T}$ for various (random) matrices $A\in\mathbb{R}^{n\times n}$ and
order $n$. Also compute the eigenvalues of the perturbed matrices $B+R_{i}$,
where $R_{i}$ are random matrices with entries from the interval $\left[
-10^{-5},10^{-5}\right]  $ ($i=1,\ldots,N$). What is the maximum perturbation
of the eigenvalues? How precise is the Bauer-Fike estimate?
\nevindex{Bauer, F. L.}
Suggested values are $N=10$ and $5\leq n\leq200$. How do the results depend on the condition
number and the order $n$? Display the maximum perturbations and the
Bauer-Fike estimates graphically.}
\end{fld}

\begin{fejmegj}
The a posteriori error estimates of linear algebraic systems are not
completely reliable. Demmel, Diament és 
Malajovich\nevindex{Ðemel, J.}\nevindex{Diament, B.}\nevindex{Malajovich, Diement G.} \cite{DemmelDM} showed
that for the $\Theta\left(n^{2}\right)$ number estimators there are always
cases when the estimate is unreliable (the error of the estimate exceeds a
given order). The first appearance of the iterative improvement is due to Fox,
Goodwin, Turing and Wilkinson (1946)\nevindex{Wilkinson, James H.}. 
\nevindex{Fox, L.}\nevindex{Goodwin, E. T.}
\nevindex{Turing, Alan (1912--1954)}
The experiences show that the decrease
of the residual error is not monotone.

Young \cite{Young}\nevindex{Young, L. A.}, Hageman and Young
\nevindex{Hageman, L. A.} \cite{HagemanY} give an excellent survey
of the theory and application of iterative methods. Barett, Berry et al.
\nevindex{Barett, R.}\nevindex{Berry, M.}
\cite{Barett} give a software oriented survey of the subject.
Frommer\nevindex{Frommer, Andreas}
\cite{Frommer} concentrates on the parallel computations.

The convergence of the $QR$-method is a delicate matter. It is analyzed in
great depth and much better results than Theorem \ref{parltetel} exist in the
literature. There are $QR$-like methods that involve double shifting.
Batterson \cite{Batterson}\nevindex{Batterson, Steve}
showed that there exists a $3\times3$
Hessenberg\nevindex{Hessenberg, Gerhard}
matrix with complex eigenvalues such that convergence cannot be achieved even
with multiple shifting.

Several other methods are known for solving the eigenvalue problems (see, e.g.
\cite{Watki1,Wilk1}). The $LR$-method is one of the best known ones.
It is very effective on positive definite Hermitian matrices. The $LR$-method
computes the Cholesky-decomposition\nevindex{Cholesky, André-Louis (1875--1918)} 
$A_{k}=LL^{\ast}$ and sets $A_{k+1}%
= L^{\ast}L .$
\end{fejmegj}

