\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts}

\newcommand{\x}{\boldsymbol{x}}
\newcommand{\y}{\boldsymbol{y}}
\newcommand{\z}{\boldsymbol{z}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\r}{\boldsymbol{r}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\newcommand{\e}{\boldsymbol{e}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}


\begin{document}

\title{Numerical Methods I -- Lab 4}
\author{Fall Semester 2005}
\date{October 11, 2005}
\maketitle

\noindent

\emph{Iterative refinement} is a technique for cheaply improving the
accuracy of the solution of $A\x=\b$ initially obtained through
$LU$-decomposition.  It was widely used in times when computations
where usually single precision, and the accuracy of the obtained
solutions was often insufficient.  It is sometimes useful even today.

The Algorithm is as follows.

\begin{enumerate}
\item (Usual $LU$-solve.) Given $A\x=\b$ with nonsingular $A\in
M(n\times n)$, find the $LU$-factorization of $A$.  Solve $L\y=\b$
(forward substitution), then $U\x=\y$ (backward substitution).

\item (Compute residual.) Compute $\r = \b-A\x$, called the residual.
If possible, this step should be done in higher precision.

\item (Update.) Compute $\x_{\text{new}} = \x + \z$, where $\z$ solves
$A\z=\r$; use the $LU$-factorization from part (1) to solve this system.

\end{enumerate}
Steps (2) and (3) can be repeated if necessary.

\paragraph{Tasks:}

\begin{enumerate}
\item Write a program that implements the above algorithm. For
obtaining the $LU$-factorization, you may use the \textsl{Octave}
function \texttt{lu}.  In order to mimic single precision, perturb the
entries of U and the off-diagonal entries of L by random small
relative error $\leq 10^{-7}$.  To mimic additional stability
problems, larger perturbations may be tried as well.

\item Test the effect of iterative refinement on the Vandermonde test
problem from the homework.  Give an estimate of the number of
operations in Step 2 and 3, and compare it with the operation count
for Step 1.

\item You may interpret iterative refinement as a linear iterative
method for solving $A\x = \b$ with the starting vector $\x_0$ obtained
in Step 1, or even with $\x_0 = 0$. What is its iteration matrix $B$?
Can you explain (based on your knowledge on the convergence theory of
linear iterative methods) why in practice one step of iterative
refinement is often sufficient to produce a high-accuracy result? This
question relates to the testing experience in part (2).

\end{enumerate}
\end{document}
