\chapter{Linear algebra functions}
%--------------------------------------------------------------------------
\label{sec:lin_alg}

\starthistory
  020502 & Created and written by Claudia Emde.\\
\stophistory

%
% Introduction
%

Solving the vector radiative transfer equation requires the
computation of linear equation systems and the matrix
exponential. This section describes the functions which are implemented
in ARTS and it gives instructions how these functions can be used, also
for other purposes than the radiative transfer calculations.

\section{Implementation files}
%-------------------------------------------------------------------------
\label{sec:lin_alg:files}

All the functions described below can be found in the files:
\begin{itemize}
\item \artsstyle{lin\_alg.h}
\item \artsstyle{lin\_alg.cc}
\end{itemize}
The template class \artsstyle{Array} and the classes \artsstyle{Matrix} and
\artsstyle{Vector} are used, therefore the linear algebra functions require
the files:
\begin{itemize}
\item \artsstyle{matpackI.h}
\item \artsstyle{make\_vector.h}
\item \artsstyle{array.h}
\item \artsstyle{matpackI.cc}
\item \artsstyle{make\_vector.cc}
\item \artsstyle{array.cc}
\end{itemize}
Furthermore logical functions contained in
\begin{itemize}
\item \artsstyle{logic.h}
\item \artsstyle{logic.cc}
\end{itemize}
are used to check the dimensions of input matrices for various functions.


\section{Linear Equation Systems}
%------------------------------------------------------------------------
\label{sec:lin_alg:lineqsys}


For solving a set of linear equations 
\begin{eqnarray}
\label{eq:lin_equ_sys}
{\bf A} {\bf x} = {\bf b}
\end{eqnarray}
the LU decomposition method is implemented.A slightly modified version
of the algorithm described in
[\cite{numerical_recipes_C:97}] is used here. An alternative method
is the Gauss-Jordan elimination, but this method is three times
slower than the LU decomposition method
[\cite{numerical_recipes_C:97}, p.36]. The LU decomposition method
reqires two functions, \artsstyle{ludcmp} and \artsstyle{lubacksub},
which will be decribed below.
\vspace{0.5cm}\\
The following example for a three dimensional equation sytem 
demonstrates how to solve a linear
equation sytem of the type
(\ref{eq:lin_equ_sys}):
\begin{itemize}
\item Create matrix A, vector b: \\
  \artsstyle{A = Matrix(3,3);} \\
  \artsstyle{A(1,1) = 4;}\\
  \artsstyle{A(2,1) = 3;}\\
  $\cdots$\\
  \artsstyle{b = Vector(3);}\\
  \artsstyle{b(1) = 7;}\\
  $\cdots$
\item Initialize solution vector x and two other variables needed for
  storing intermediate results:\\
  \artsstyle{x = Vector(3);}\\
  \artsstyle{LU = Matrix(3,3);}\\
  \artsstyle{indx = ArrayOfIndex(3);}
\item Call LU decomposition function (see Section \ref{sec:lin_alg:lu_decomp}): \\
  \artsstyle{ludcmp(LU, indx, A);}
\item Call LU backsubstitution function (see Section \ref{sec:lin_alg:backsub}): \\
  \artsstyle{lubacksub(x, LU, b, indx);}
\item Print the solution vector:\\
  \artsstyle{cout << x;}
\end{itemize}

\subsection{LU Decomposition}
%------------------------------------------------------------------------
\label{sec:lin_alg:lu_decomp}

A LU decomposition is a procedure for decomposing a square matrix {\bf
  A} with
dimension $n$ into a product of a lower triangular matrix {\bf L} (has
elements only on the diagonal elements and below) and
an upper triangular matrix {\bf U} (has elements only on the diagonal
and above):
\begin{equation}
  \label{eq:lu_decomp}
  {\bf L}\cdot{\bf U} ={\bf A} 
\end{equation}
For a 3 x 3 matrix equation \ref{eq:lu_decomp} would look like this: 
\[ \left(
  \begin{array}{ccc}
    l_{11} & 0 & 0 \\
    l_{21} & l_{22} & 0 \\
    l_{31} & l_{32} & l_{33}
    \end{array} \right)
\cdot
\left(
  \begin{array}{ccc}
    u{11} & u_{12} & u_{13} \\
    0 & u_{22} & u_{23}\\
    0 & 0 & u_{33}
    \end{array} \right)
=
\left(
  \begin{array}{ccc}
    a_{11} & a_{12} & a_{13} \\
    a_{21} & a_{22} & a_{23} \\
    a_{31} & a_{32} & a_{33}
    \end{array} \right)
\]
The decomposition can be used to rewrite the linear set of equations
(\ref{eq:lin_equ_sys}) in the following way:
\begin{eqnarray}
  {\bf A}\cdot{\bf x} = ({\bf L}\cdot{\bf U})\cdot{\bf x} = {\bf
    L}\cdot({\bf U}\cdot{\bf x}) = {\bf b}
\end{eqnarray}
First 
\begin{equation}
  {\bf L} \cdot{\bf y} = {\bf b}
\end{equation}
is solved for the vector ${\bf y}$ which can be done by
forward substitution (see section \ref{sec:lin_alg:backsub}). Then 
\begin{equation}
  {\bf U} \cdot{\bf x} = {\bf y}
\end{equation}
is solved again by backsubstitution. 
The advantage in breaking up one linear set into two successive ones
is that the solution of a triangular set of equations is quite trivial.

The function \artsstyle{ludcmp} requires a square matrix of arbitrary
dimension $n$ as input and performs the LU decomposition. It returns one
matrix which contains both matrices, ${\bf L}$ and ${\bf U}$. 
For the lower triangular matrix  ${\bf L}$ the diagonal elements 
are chosen to be 1, then the other elements of ${\bf L}$ and ${\bf U}$
are determined. This is possible, as the LU decomposition is an under
determined equation sytem with $n^2$ equations for $n^2+n$ unknowns. 
The output matrix does not include the diagonal of ${\bf L}$, in the
three-dimensional case it has the following elements:
\[ \left(
  \begin{array}{ccc}
    u_{11} & u_{12} & u_{13} \\
    l_{21} & u_{22} & u_{23} \\
    l_{31} & l_{32} & u_{33}
    \end{array} \right)
\]
This special arrangement of the LU decomposition is named {\sl
Crout's algorithm} and a matrix arranged in this form is named {\sl
Crout matrix} in this context.
  

Another output variable of the function \artsstyle{ludcmp} is an index
vector which contains information about pivoting which is absolutely
essential for the stability of
Crout's algorithm. Here partial pivoting,
i.e. interchange of rows is implemented. That means that not {\bf A} is
decomposed into $LU$-form but a rowwise permutation of {\bf A}. If the
index vector contains for example the elements $(2,1,0)$ the first and
the last row of a three dimensional matrix would be exchanged.


\subsection{Forward- and Backsubstitution}
%---------------------------------------------------------------
\label{sec:lin_alg:backsub}
An equation system of the form
\[ 
\left(
  \begin{array}{ccc}
    a{11} & a_{12} & a_{13} \\
    0 & a_{22} & a_{23}\\
    0 & 0 & a_{33}
    \end{array} \right)
\cdot
\left(
  \begin{array}{c}
    x_1\\x_2\\x_3
 \end{array} \right)
=
\left(
  \begin{array}{c}
    b_1\\b_2\\b_3
 \end{array} \right)
\]
can be solved very easy. The last element, here $x_3$, is already isolated,
namely
\begin{eqnarray}
  x_3 = b_3/a_{33}
\end{eqnarray}
As $x_3$ is known $x_2$ can be calculated using the second row of the
eqautions. Then, finally, $x_1$ can be calculated as well using the
first row. This procedure
is called backsubtitution. The same
method  applied for an equation system including a
lower triangular matrix is named forward substitution.   

The function \artsstyle{lubacksub} does forward and backward
substitution to solve the equation system described in
\ref{sec:lin_alg:lu_decomp}. As input it requires the output variables of
\artsstyle{ludcmp} which are the {\sl Crout matrix} and the index
vector. Output of the function is the solution vector ${\bf x}$ to the
equation system.


\subsection{More Applications of the LU Decomposition}
%-------------------------------------------------------------------
\label{lu_applications}

\begin{itemize}
\item Inverse of a matrix:\\
  To compute $(\bf K)^{-1}\cdot{\bf b}$, which is a part of the
solution to the vector radiative transfer equation (\ref{eq:scattering:RTE}) the LU
decomposition method can be used. The following equations show, that
the problem is equivalent to  solving a linear equation system of the type
\ref{eq:lin_equ_sys}.
\begin{eqnarray}
  {\bf K}^{-1}\cdot{\bf b} &=& {\bf x}\\
\Leftrightarrow \qquad  {\bf K}\cdot{\bf x} &=& {\bf b}
\end{eqnarray}

\item To solve the equation system
  \begin{eqnarray}
    {\bf A}\cdot{\bf X} &=& {\bf B}
  \end{eqnarray}
where {\bf A}, {\bf B} and  {\bf X} are matrices of dimension
$n$, the LU decomposition functions can be applied as well. Assume
that {\bf A} and {\bf B} are known and you want to solve for {\bf
 X}.
First you should do a LU decomposition of  {\bf A} and then
backsubstitute with the columns of B and you get the columns of {\bf
  X} as solution vectors.

\end{itemize}

\section{Matrix Exponential Function}
%----------------------------------------------------------------
\label{sec:lin_alg:mat_exp}

A very important function for solving differential equations is the
matrix exponential:
\begin{eqnarray}
  \label{eq:mat_exp}
  e^{{\bf A}s} = \sum_{k=0}^\infty{\frac{({\bf A} s)^k}{k!}}
\end{eqnarray}
In principle it could be computed using the Taylor power series but 
 this method is not efficient. {\sc Moler} and {\sc Van
  Loan} have shown for the simple example [\cite{Moler_Loan:79}]
\[ {\bf A} =
\left(
  \begin{array}{cc}
    -49 & 24\\
    -64 & 31
    \end{array} \right) \]
that convergence is obtained not until 59 terms. And if a relative
accuracy of only 10$^{-5}$ is taken, the method even leads to a wrong
result due to rounding errors.

\subsection{Pad\'e Approximation}
%----------------------------------------------------------------------
\label{sec:pade_approximation}

One of the better algorithms for computing the matrix exponential is
the Pad\'e approximation which is also shortly described in
[\cite{Moler_Loan:79}] and outlined in the book ``Matrix
Computations'' by \cite{Golub_Loan:91}. 
The method uses perturbation theorie as well as the so called Pad\'e
functions. It is possible to derive an algorithm which calculates
\begin{eqnarray}
  {\bf F} = e^{{\bf A}+{\bf E}} 
\end{eqnarray}
where 
\begin{eqnarray}
  \|{\bf E}\|_\infty \le \delta \|{\bf A}\| .
\end{eqnarray}
The accuracy of the computation given by $\delta$ can be chosen. 
The parameter q has to be the smallest non-negative integer such that
$\epsilon(q,q)\le\delta$ where
\begin{eqnarray}
  \epsilon(p,q) = 2^{3-(p+q)}\frac{p!q!}{(p+q)!(p+q+1)!}.
\end{eqnarray}
The following table shows values of epsilon for
different values of q.
\vspace{0.5cm}\\
\begin{tabular}[h]{|r|r|}
 \hline
q & $\epsilon$(q,q) \\ \hline
1 & 0.1667\\
2 & 6.9444 $\cdot$ 10$^{-4}$ \\
3 & 1.2401 $\cdot$ 10$^{-6}$ \\
4 & 1.2302 $\cdot$ 10$^{-9}$ \\
5 & 7.7667 $\cdot$ 10$^{-13}$ \\
6 & 3.3945 $\cdot$ 10$^{-16}$ \\ 
\hline
\end{tabular}
\vspace{0.5cm}\\
The algorithm is implemented in the function \artsstyle{matrix\_exp}. Input
to this function is the matrix ${\bf A}$ and the parameter $q$. As output
it gives the matrix ${\bf F}$ which is defined above.\\
The following example shows how to use the \artsstyle{matrix\_exp} function:
\begin{itemize}
\item Initialize {\bf A} and assign values:\\
  \artsstyle{Matrix A(3,3);}\\
  \artsstyle{A(1,1) = 45;}\\
  \artsstyle{A(1,2) = 3;}\\
$\cdots$ 
\item Initialize {\bf F}:\\
  \artsstyle{Matrix F(3,3);}
\item Give a paramater for the accuracy:\\
  \artsstyle{Index q=6;}
\item Call the matrix exponential function:\\
  \artsstyle{matrix\_exp(F,A,q);}
\item Print the result: \\
  \artsstyle{cout << "exp(A) = " << F;}
\end{itemize}






%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
