\section{Solver for linear and nonlinear systems}%
\label{S:solver}

\ALBERTA provides own solvers for general linear and nonlinear
systems. The solvers use dense \code{REAL}-vectors for storing
coefficients. They are aware of \ALBERTA's DOF-vector and -matrix data
structures and work with an application provided subroutine for the
matrix-vector multiplication, and in case a preconditioner is used, a
function for preconditioning. The nonlinear solvers need subroutines
for assemblage and solution of a linearized system.

In the subsequent sections we describe the basic data structures for
the \textsf{OEM} (\textsf{O}rthogonal \textsf{E}rror \textsf{M}ethods)
module, a built-in \ALBERTA interface for solving systems involving a
\code{DOF\_MATRIX} and \code{DOF\_REAL[\_D]\_VEC[\_D]} objects, and
the access to functions for matrix-vector multiplication and
preconditioning for a direct use of the \textsf{OEM} solvers.  Then we
describe the basic data structures for multigrid solvers and for the
available solvers of nonlinear equations.  Most of the implemented
methods (and more) are described for example in
\cite{Meister:99,Saad:96}.

\subsection{Krylov-space solvers for general linear systems}%
\label{S:oem}%
\label{S:OEM}%
\idx{linear solvers|(}

Very efficient solvers for linear systems are Krylov-space solvers
(or \textsf{O}rthogonal \textsf{E}rror \textsf{M}ethods).
The \textsf{OEM} library provides such solvers for the solution of general
linear systems
\[
A\, x = b
\]
with $A \in \R^{N\times N}$ and $x,b \in \R^N$. The library solvers
work on dense flat vectors and do not need to know the storage of the
system matrix, or the matrix used for preconditioning.  Matrix-vector
multiplication and preconditioning is done by application provided
routines.

Most of the implemented \code{OEM} solvers are a \code{C}-translation
from the solvers of the FORTRAN OFM library (Orthogonale Fehler
Methoden), by D\"orfler \cite{Doerfler:95a}. \code{SymmLQ} is the
algorithm described in \cite{PaigeSaunders:75}, and \code{TfQMR} is
described in {\color{red} TO BE DETERMINED}. All solvers allow for
\emph{left} preconditioning and some also for \emph{right}
preconditioning.

The data structure (defined in \code{alberta\_util.h}) for passing information
about matrix-vector multiplication, preconditioning and tolerances,
etc. to the solvers is
\ddx{OEM_DATA@{\code{OEM\_DATA}}}
\idx{linear solver!OEM_DATA@{\code{OEM\_DATA}}}
\bv\begin{lstlisting}[label=D:OEM_DATA]
typedef int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);

typedef struct oem_data OEM_DATA;
struct oem_data
{
  OEM_MV_FCT mat_vec;
  void       *mat_vec_data;
  OEM_MV_FCT mat_vec_T;
  void       *mat_vec_T_data;
  void       (*left_precon)(void *, int, REAL *);
  void       *left_precon_data;
  void       (*right_precon)(void *, int, REAL *);
  void       *right_precon_data;

  REAL       (*scp)(void *, int, const REAL *, const REAL *);
  void       *scp_data;

  WORKSPACE  *ws;

  REAL       tolerance;
  int        restart;
  int        max_iter;
  int        info;

  REAL       initial_residual;
  REAL       residual;
};
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{mat\_vec} pointer to a function for the matrix--vector multiplication
  with the system matrix;

  \code{mat\_vec(mat\_vec\_data, dim, u, b)} applies the system matrix 
  to the input vector \code{u} and stores the product in
  \code{b}; \code{dim} is the dimension of the linear system, 
  \code{mat\_vec\_data} a pointer to application.
  %% 
  \kitem{mat\_vec\_data} pointer to application data for the
  matrix-vector multiplication, first argument to
  \code{mat\_vec()}.
  %% 
  \kitem{mat\_vec\_T} pointer to a function for the matrix--vector
  multiplication with the transposed system matrix;

  \code{mat\_vec\_T(mat\_vec\_data, dim, u, b)} applies the transposed
  system matrix to the input vector \code{u} and stores 
  the product in \code{b}; \code{dim} is the dimension of the linear
  system, \code{mat\_vec\_T\_data} a pointer to application data.
  %% 
  \kitem{mat\_vec\_T\_data} pointer to application data for the matrix-vector 
  multiplication with the transposed system matrix, first argument to 
  \code{mat\_vec\_T()}.
  %% 
  \kitem{left\_precon} pointer to function for left preconditioning;
  it may be a \nil pointer; in this case no left preconditioning is done;
  
  \code{left\_precon(left\_precon\_data, dim, r)} is the 
  implementation of the left preconditioner; \code{r} is 
  input and output vector of length \code{dim} and 
  \code{left\_precon\_data} a pointer to application data.
  %% 
  \kitem{left\_precon\_data} pointer to application data for the left
  preconditioning, first argument to \code{left\_precon()}.
  %% 
  \kitem{right\_precon} pointer to function for right preconditioning;
  it may be a \nil pointer; in this case no right preconditioning is done;

  \code{right\_precon(right\_precon\_data, dim, r)} is the
  implementation of the right preconditioner; \code{r} is input and
  output vector of length \code{dim} and \code{right\_precon\_data} a
  pointer to application data.
  %% 
  \kitem{right\_precon\_data} pointer to application data for the
  right preconditioning, first argument to \code{right\_precon()}.
  %% 
  \kitem{scp} pointer to a function for computing a problem dependent scalar
  product; it may be a \nil pointer; in this case the Euclidian
  scalar product is used;

  \code{scp(scp\_data, dim, x, y)} computes the problem dependent
  scalar product of the two vectors \code{x} and \code{y} of length
  \code{dim}; \code{scp\_data} is a pointer to application data.
  %% 
  \kitem{scp\_data} pointer to application data for computing the
  scalar product, first argument to \code{scp()}.
  %% 
  \kitem{ws} a pointer to a \code{WORKSPACE} structure for storing additional
  vectors used by a solver; if the space is not sufficient, the
  used solver will enlarge this workspace; if \code{ws} is \nil, then
  the used solver allocates memory, which is freed before exit.
  %% 
  \kitem{tolerance} tolerance for the residual; if the norm of the residual
  is less than or equal to \code{tolerance}, the solver returns
  the actual iterate as the solution of the system.
  %% 
  \kitem{restart} restart for the linear solver; used only by
  \code{oem\_gmres()} at the moment.
  %% 
  \kitem{max\_iter} maximal number of iterations to be
  performed although the tolerance may not be reached.
  %% 
  \kitem{info} the level of information produced by the solver; \code{0}
  is the lowest level of information (no information is printed)
  and \code{10} the highest level.
  %% 
  \kitem{initial\_residual} stores the norm of the initial residual on exit.
  %% 
  \kitem{residual} stores the norm of the final residual on exit.
\end{descr}

The following linear solvers are currently implemented.
Table \ref{T:OEM_methods} gives an overview over the implemented
solvers, the matrix types they apply to, and the cost of one iteration.
\begin{table}
\begin{center}
 \begin{tabular}{|c|c|l|c|} \hline
   Method & Matrix & Operations & Storage\\ \hline\hline
   BiCGstab & symmetric & 2 MV + 12 V & $5N$ \\
   CG & symmetric positive definite& 1 MV + 5 V & $3N$ \\
   GMRES & regular & k MV + ... & $(k+2)N + k(k+4)$ \\
   ODir & symmetric positive & 1 MV + 11 V & $5N$ \\
   ORes & symmetric & 1 MV + 12 V & $7N$ \\
   SymmLQ & symmetric & & $6N$\\
   TfQMR & regular & & $11N$\\
\hline
 \end{tabular}
\end{center}
\caption[Iterative solvers, storage requirements and matrix types]{OEM methods with applicable matrix types,
  numbers of operations per iteration (MV matrix-vector products, V
  vector operations), and storage requirements ($N$ number of unknowns,
  $k$ GMRES subspace dimension)}
\label{T:OEM_methods}
\end{table}
\fdx{oem_bicgstab()@{\code{oem\_bicgstab()}}}%
\idx{linear solvers!oem_bicgstab()@{\code{oem\_bicgstab()}}}%
\fdx{oem_cg()@{\code{oem\_cg()}}}%
\idx{linear solvers!oem_cg()@{\code{oem\_cg()}}}%
\fdx{oem_tfqmr()@{\code{oem\_tfqmr()}}}%
\idx{linear solvers!oem_tfqmr()@{\code{oem\_tfqmr()}}}%
\fdx{oem_gmres()@{\code{oem\_gmres()}}}%
\idx{linear solvers!oem_gmres()@{\code{oem\_gmres()}}}%
\fdx{oem_gmres_k()@{\code{oem\_gmres\_k()}}}%
\idx{linear solvers!oem_gmres_k()@{\code{oem\_gmres\_k()}}}%
\fdx{oem_odir()@{\code{oem\_odir()}}}%
\idx{linear solvers!oem_odir()@{\code{oem\_odir()}}}%
\fdx{oem_ores()@{\code{oem\_ores()}}}%
\idx{linear solvers!oem_ores()@{\code{oem\_ores()}}}%
\fdx{oem_symmlq()@{\code{oem\_symmlq()}}}%
\idx{linear solvers!oem_symmlq()@{\code{oem\_symmlq()}}}%
\bv\begin{lstlisting}
int oem_bicgstab(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_cg(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_gmres(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_gmres_k(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_odir(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_ores(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_tfqmr(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
int oem_symmlq(OEM_DATA *oem_data, int dim, const REAL *rhs, REAL *u0);
\end{lstlisting}\ev
Description:
\begin{descr}
  \kitem{oem\_bicgstab(oem\_data, dim, rhs, u0)} solves a linear system
  by a stabilized BiCG method and can be used for symmetric
  system matrices; \code{oem\_data} stores information about
  matrix vector multiplication, preconditioning, tolerances, etc.
  \code{dim} is the dimension of the linear system, \code{rhs}
  the right hand side vector, and \code{u0} the initial guess
  on input and the solution on output; \code{oem\_bicgstab()}
  needs a workspace for storing \code{5*dim} additional \code{REAL}s;
  the return value is the number of iterations;
  \code{oem\_bicgstab()} only uses left preconditioning.
  %% 
  \kitem{oem\_cg(oem\_data, dim, rhs, u0)} solves a linear system
  by the conjugate gradient method and can be used for symmetric
  positive definite system matrices; \code{oem\_data} stores
  information about matrix vector multiplication,
  preconditioning, tolerances, etc. \code{dim} is the dimension
  of the linear system, \code{rhs} the right hand side vector, and
  \code{u0} the initial guess on input and the solution on
  output; \code{oem\_cg()} needs a workspace for storing
  \code{3*dim} additional \code{REAL}s; the return value is the
  number of iterations; \code{oem\_cg()} only uses left
  preconditioning.
  %% 
  \kitem{oem\_gmres(oem\_data, dim, rhs, u0)} solves a linear system
  by the GMRes method with restart and can be used for regular
  system matrices; \code{oem\_data} stores information about
  matrix vector multiplication, preconditioning, tolerances, etc.
  \code{dim} is the dimension of the linear system, \code{rhs}
  the right hand side vector, and \code{u0} the initial guess
  on input and the solution on output; \code{oem\_data->restart}
  is the dimension of the Krylov--space for 
  the minimizing procedure; \code{oem\_data->restart} must be bigger
  than \code{0} and less or equal \code{dim}, otherwise
  \code{restart=10} will be used; \code{oem\_gmres()} needs a
  workspace for storing  
  \code{(oem\_data->restart+2)*dim +
    oem\_data->restart*(oem\_data->restart+4)} 
  additional \code{REAL}s.
  %% 
  \kitem{oem\_gmres\_k(oem\_data, dim, rhs, u0)} performs just one restart step
  (minimization on a $k$-dimensional Krylov subspace) of the GMRES
  method. This routine can be used as subroutine in other
  solvers. For example, \code{oem\_gmres()} just iterates this
  until the tolerance is met. Other applications are nonlinear
  GMRES solvers, where a new linearization is done after each
  linear GMRES restart step.
  %% 
  \kitem{oem\_odir(oem\_data, dim, rhs, u0)} solves a linear system
  by the method of orthogonal directions and can be used for
  symmetric, positive system matrices; \code{oem\_data} stores
  information about matrix vector multiplication,
  preconditioning, tolerances, etc.  \code{dim} is the dimension
  of the linear system, \code{rhs} the right hand side vector, and
  \code{u0} the initial guess on input and the solution on
  output; \code{oem\_dir()} needs a workspace for storing
  \code{5*dim} additional \code{REAL}s; the return value is the
  number of iterations; \code{oem\_odir()} only uses left
  preconditioning.
  %% 
  \kitem{oem\_ores(oem\_data, dim, rhs, u0)} solves a linear system
  by the method of orthogonal residuals and can be used for symmetric
  system matrices; 
  \code{oem\_data} stores information about
  matrix vector multiplication, preconditioning, tolerances, etc.
  \code{dim} is the dimension of the linear system, \code{rhs}
  the right hand side vector, and \code{u0} the initial guess
  on input and the solution on output; \code{oem\_res()}
  needs a workspace for storing \code{7*dim} additional \code{REAL}s;
  the return value is the number of iterations;
  \code{oem\_ores()} only uses left preconditioning.
  %% 
  \kitem{oem\_symmlq(oem\_data, dim, rhs, u0)} solves a symmetric
  linear system. \code{oem\_data} stores information about matrix
  vector multiplication, preconditioning, tolerances, etc.
  \code{dim} is the dimension of the linear system, \code{rhs}
  the right hand side vector, and \code{u0} the initial guess on
  input and the solution on output; \code{oem\_symmlq()} needs a
  workspace for storing \code{6*dim} additional \code{REAL}s; the
  return value is the number of iterations.  \code{oem\_symmlq()}
  supports uses left preconditioning.
  %% 
  \kitem{oem\_tfqmr(oem\_data, dim, rhs, u0)} solves a linear system
  using a transpose-free QMR method and can be used for regular
  system matrices; \code{oem\_data} stores
  information about matrix vector multiplication,
  preconditioning, tolerances, etc. \code{dim} is the dimension
  of the linear system, \code{rhs} the right hand side vector, and
  \code{u0} the initial guess on input and the solution on
  output; \code{oem\_tfqmr()} needs a workspace for storing
  \code{11*dim} additional \code{REAL}s; the return value is the
  number of iterations.
\end{descr}

\subsection{Krylov-space solvers for DOF matrices and vectors}%
\label{S:ALBERTA_OEM_solvers}

\begin{compatibility}
  \label{compat:precon}
  The support for additional preconditioners, as well as the
  block-matrix structure induced by the support for
  \hyperref[S:chain_impl]{direct sums of finite element spaces} (see
  \secref{S:chain_impl}) made it necessary to provide a more flexible
  and extendible interface to the implemented preconditioners.
  Additionally, some of the preconditioners need further parameters.

  Therefore, the selection of a particular preconditioner has been
  moved to separate functions
  \hyperref[S:init_oem_precon_fct]{\code{init\_oem\_precon()}},
  \hyperref[S:vinit_oem_precon_fct]{\code{vinit\_oem\_precon()}} and
  \hyperref[S:init_precon_from_type_fct]{\code{init\_precon\_from\_type()}},
  the latter requiring a special support structure
  \hyperref[S:PRECON_TYPE_struct]{\code{PRECON\_TYPE}} to pass
  parameters on to the preconditioners.

  Solver-functions, which previously accepted a mere integer to select
  a particular preconditioner, now need a pointer to a
  \hyperref[T:PRECON]{\code{PRECON}}-structure, see below
  \secref{S:precon}.
\end{compatibility}

We describe here the interface between \ALBERTA's DOF-vectors and
-matrices and the available general OEM-solvers described in the
previous \secref{S:oem}. At the highest level, there are three function,
namely \code{oem\_solve\_s()}, \code{oem\_solve\_d()} and
\code{oem\_solve\_dow()}.  The calling conventions for the three
functions are functionally identical, except for the data-type of the
DOF-vector arguments.  The function \code{oem\_solve\_s()} is used for
scalar valued problems, i.e.
\[
A\,x = b
\]
with $A \in \R^{N\times N}$ and $x,b \in \R^N$. Vector valued problems
need a closer examination, there are two cases:
\begin{enumerate}
\item \DOW-valued finite element spaces based on scalar basis functions:

\code{oem\_solve\_d()} and \code{oem\_solve\_dow()} can both either be
used for decoupled or coupled \DOW-valued problems. Decoupled problems
are of the form
\[
\left[\begin{matrix}
A      &  0      & \ldots & 0\\
0      &  A      & \ddots & \vdots\\
\vdots & \ddots  & \ddots & 0\\
0      & \ldots  &    0   & A
\end{matrix}\right]
\left[\begin{matrix} u_1\\ u_2\\ \vdots\\ u_n\end{matrix}\right]
=
\left[\begin{matrix} f_1\\ f_2\\ \vdots\\ f_n\end{matrix}\right]
\]
with $A \in \R^{N\times N}$ and $u_i,f_i \in \R^N$, $i = 1,\dots,n$,
where $n = \code{DIM\_OF\_WORLD}$. The vectors $(u_1,\dots,u_n)$ and
$(f_1,\dots,f_n)$ are stored in \code{DOF\_REAL\_D\_VEC}s, whereas
the matrix is stored as a single scalar \code{DOF\_MATRIX}.

Coupled \DOW-valued problems lead in this context to matrices of the form
\[
\left[\begin{matrix}
A^{00} &  \ldots & A^{0n}\\
\vdots &  \ddots & \vdots\\
A^{n0} &  \ldots & A^{nn}
\end{matrix}\right]
\left[\begin{matrix} u_1\\ \vdots\\ u_n\end{matrix}\right]
=
\left[\begin{matrix} f_1\\ \vdots\\ f_n\end{matrix}\right]
\]
with $A^{\mu\nu} \in \R^{N\times N}$ and $u_\nu,f_\mu \in \R^N$,
$\mu,\nu = 1,\dots,n$, where $n = \code{DIM\_OF\_WORLD}$. The vectors
$(u_1,\dots,u_n)$ and $(f_1,\dots,f_n)$ are again stored in
\code{DOF\_REAL\_D\_VEC}s. One prominent example is the discretisation
of a Stokes-problem with prescribed stresses on the boundary: in this
case the weak formulation has to be based on the deformation tensor,
which leads to matrix of above type. The matrix is still stored as a
\code{DOF\_MATRIX} structure, but its entries are $\DOW\times\DOW$
blocks: the data is stored as an $N\times N$ matrix of small $d\times
d$ blocks in analogy to \code{DOF\_REAL\_D\_VEC}s. See also
\ref{book:S:DisCoupled}.  Compare also \compatref{compat:DOWB_MATRIX}.

\item Finite element spaces based on \DOW-valued basis functions:

  In this case the DOF-vectors are scalar-valued, and the resulting
  DOF-matrix is just a scalar matrix, compare also
  \secref{book:S:DisCoupled}.
\end{enumerate}

Note that the interface routines to the OEM-solvers are aware of
direct sums of finite element spaces, as described in
\secref{S:chain_impl}, the resulting block-matrices generated by the
assemble-framework will be handled correctly, including the cases
where a standard Lagrangian finite element space is augmented by
vector-valued basis functions like face-bubbles.

An application selects a particular solver by passing one of the
following enumeration values to \code{oem\_solve\_[s|d|dow]()}:
%%
\ddx{OEM_SOLVER@{\code{OEM\_SOLVER}}}
\idx{linear solver!OEM_SOLVER@{\code{OEM\_SOLVER}}}
%%
\bv\begin{lstlisting}[label=enum:OEM_SOLVER]
typedef enum {
  NoSolver, BiCGStab, CG, GMRes, ODir, ORes, TfQMR, GMRes_k, SymmLQ
} OEM_SOLVER;
\end{lstlisting}\ev
%%
New identifiers may be added to this enumeration when new solvers are
added to \ALBERTA.
%%
In more detail, the three high-level interface function are described
below:
\begin{function}{oem\_solve\_[s|d|dow]()}
\label{S:oem_solve_fct}

\item[Prototypes] ~\hfill
%%
\fdx{oem_solve_d()@{\code{oem\_solve\_d()}}}%
\idx{linear solvers!oem_solve_d()@{\code{oem\_solve\_d()}}}%
\fdx{oem_solve_s()@{\code{oem\_solve\_s()}}}%
\idx{linear solvers!oem_solve_s()@{\code{oem\_solve\_s()}}}%
\fdx{oem_solve_dow()@{\code{oem\_solve\_dow()}}}%
\idx{linear solvers!oem_solve_dow()@{\code{oem\_solve\_dow()}}}%
%%
\bv\begin{lstlisting}
int oem_solve_s(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
                const DOF_REAL_VEC *f, DOF_REAL_VEC *u,
                OEM_SOLVER solver,
                REAL tol, const PRECON *precon,
                int restart, int max_iter, int info);
int oem_solve_d(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
                const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u,
                OEM_SOLVER solver,
                REAL tol, const PRECON *precon,
                int restart, int max_iter, int info);
int oem_solve_dow(const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
                  const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u,
                  OEM_SOLVER solver,
                  REAL tol, const PRECON *precon,
                  int restart, int max_iter, int info);
\end{lstlisting}\ev

\fdx{oem_solve_d()@{\code{oem\_solve\_d()}}}%
\idx{linear solvers!oem_solve_d()@{\code{oem\_solve\_d()}}}%
\fdx{oem_solve_s()@{\code{oem\_solve\_s()}}}%
\idx{linear solvers!oem_solve_s()@{\code{oem\_solve\_s()}}}%
\fdx{oem_solve_dowb()@{\code{oem\_solve\_dowb()}}}%
\idx{linear solvers!oem_solve_dowb()@{\code{oem\_solve\_dowb()}}}%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
iterations = oem_solve_[s|d|dow](A, mask, f, u,
                                 solver, tol, precon,
                                 restart, max_iter, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Attempt to solve the linear system defined by the
  matrix \code{A}, an optional restriction to a sub-space by masking
  out DOFs via \code{mask}, a load-vector \code{f} and an initial
  guess and storage \code{u} for the approximative solution.
\item[Parameters]~\hfill
  \begin{descr}
  \kitem{A} Pointer to a \code{DOF\_MATRIX} storing the system
    matrix.
    %% 
  \hyperitem{oem_solve:mask}{mask} Pointer to a \code{DOF\_SCHAR\_VEC}
    masking out parts of the finite element space: if
    \code{mask->vec[d] >= DIRICHLET}, then \code{A} will act as if the
    $d$-th row would be zero. Compare also the discussion in the
    section about Dirichlet boundary condition, see
    \secref{S:dirichlet_bound}
    %% 
  \kitem{f} A pointer to a \code{DOF\_REAL[\_D]\_VEC[\_D]} storing
    the load-vector of the linear system.
    %% 
  \kitem{u} A pointer to a \code{DOF\_REAL[\_D]\_VEC[\_D]} storing the
    initial guess on input and the approximative solution on output.
    In the context of interpolated Dirichlet boundary conditions
    special provisions have to be taken for the ``Dirichlet-nodes''.
    Compare the discussion in \secref{S:dirichlet_bound}.
    %% 
  \kitem{solver} Use the respective OEM-solver; see
    \hyperref[enum:OEM_SOLVER]{above} for the available keywords.
    %% 
  \kitem{tol} Tolerance for the residual; if the norm of the
    residual is less or equal \code{tol},
    \code{oem\_solve\_[s|d|dow]()} returns the actual iterate as the
    approximative solution of the system.
    %% 
  \kitem{precon} A pointer to a structure describing the
    preconditioner to use, see further below in \secref{S:precon}.
    \begin{compatibility}
      Previous versions used a simple number here, but as the
      preconditioner frame-work has become much more complicated,
      because of the support for direct sums of finite element spaces,
      the code for the selection of the preconditioner has been
      separated from the entry-point to the solvers.
    \end{compatibility}
    %% 
  \kitem{restart} Only used by \code{gmres}: the maximum dimension
    of the Krylov-space.
    %% 
  \kitem{max\_iter} Maximal number of iterations to be performed by
    the linear solver. This can be compared with the return value --
    which gives the number of iterations actually performed -- to
    determine whether the solver has achieved its goal.
    %% 
  \kitem{info} This is the level of information of the linear solver; 
    \code{0} is the lowest level of information (no information is printed)
    and \code{10} the highest level.
  \end{descr}
\item[Return Value] ~\hfill

  The number of iterations the solver needed until the norm of the
  residual was below \code{tol}, or \code{max\_iter} if the solver was
  not able to reach its goal before the prescribed maximum iteration
  count was exhausted.
\end{function}

There is also an interface to the OEM-solvers which splits the call to
the OEM-methods into an initialization part, an execution part and a
cleanup part. This is useful when the same solver applies the same
matrix to varying load-vectors. One example is the implementation of a
CG-method for Schur's complement operator of a saddle-point problem
(see \secref{S:OEM_SPCG} below). The following functions implement
this interface:

\fdx{get_oem_solver()@{\code{get\_oem\_solver()}}}%
\idx{linear solvers!get_oem_solver()@{\code{get\_oem\_solver()}}}%
%%
\fdx{init_oem_solve()@{\code{init\_oem\_solve()}}}%
\idx{linear solvers!init_oem_solve()@{\code{init\_oem\_solve()}}}%
%%
\fdx{release_oem_solve()@{\code{release\_oem\_solve()}}}%
\idx{linear solvers!release_oem_solve()@{\code{release\_oem\_solve()}}}%
%%
\fdx{call_oem_solve_s()@{\code{call\_oem\_solve\_s()}}}%
\idx{linear solvers!call_oem_solve_s()@{\code{call\_oem\_solve\_s()}}}%
%%
\fdx{call_oem_solve_d()@{\code{call\_oem\_solve\_d()}}}%
\idx{linear solvers!call_oem_solve_d()@{\code{call\_oem\_solve\_d()}}}%
%%
\fdx{call_oem_solve_dow()@{\code{call\_oem\_solve\_dow()}}}%
\idx{linear solvers!call_oem_solve_dow()@{\code{call\_oem\_solve\_dow()}}}%
%%
\bv\begin{lstlisting}
typedef int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);

OEM_MV_FCT get_oem_solver(OEM_SOLVER);
OEM_DATA *init_oem_solve(const DOF_MATRIX *A,
                         const DOF_SCHAR_VEC *mask,
                         REAL tol, const PRECON *precon,
                         int restart, int max_iter, int info);
void release_oem_solve(const OEM_DATA *oem);
int call_oem_solve_s(const OEM_DATA *oem, OEM_SOLVER solver,
                     const DOF_REAL_VEC *f, DOF_REAL_VEC *u);
int call_oem_solve_dow(const OEM_DATA *oem, OEM_SOLVER solver,
                        const DOF_REAL_VEC_D *f, DOF_REAL_VEC_D *u);
int call_oem_solve_d(const OEM_DATA *oem, OEM_SOLVER solver,
                     const DOF_REAL_D_VEC *f, DOF_REAL_D_VEC *u);
\end{lstlisting}\ev

See Example
\ref{example:oem_solve_impl}-\ref{example:get_oem_solver_direct_sums}
for short code skeletons explaining the use of these functions. The
descriptions for the individual functions are as follows:

\begin{function}{get\_oem\_solver()}
\label{S:get_oem_solver_fct}
%%
\fdx{get_oem_solver()@{\code{get\_oem\_solver()}}}%
\idx{linear solvers!get_oem_solver()@{\code{get\_oem\_solver()}}}%
%%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
solver_fct = get_oem_solver(solver_num);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Return a function pointer for the solver indicated by
  \code{solver\_num} which shuld be one of the symbols
  \code{BiCGStab}, \code{CG} \code{GMRes}, \code{ODir}, \code{ORes},
  \code{TfQMR}, \code{GMRes\_k}, \code{SymmLQ}.

\item[Parameters]~\hfill
  \begin{descr}
    \kitem{solver\_num} As explained above.
  \end{descr}
\item[Return Value]~\hfill
  
  A function pointer conforming to the type
  
\bv\begin{lstlisting}
int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);
\end{lstlisting}\ev
\end{function}

\begin{function}{init\_oem\_solve()}
\label{S:init_oem_solve_fct}

%%
\fdx{init_oem_solve()@{\code{init\_oem\_solve()}}}%
\idx{linear solvers!init_oem_solve()@{\code{init\_oem\_solve()}}}%
%%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
oem_data_handle =
  init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \code{OEM\_DATA}-handle which can be passed to the
  function pointers returned by \code{get\_oem\_solver()} (see above).
  The specific solver to use, as well as the storage for the solution
  and the load-vector, is left unspecified here; these data is given
  as parameters to \code.{call\_oem\_solve\_[s|d|dow]()}, see below.
  The data handle returned by this functions eventually should be
  deleted by a call to \code{realeas\_oem\_solve()}, which is also
  described below.
\item[Parameters]~\hfill

  The parameters have the same meaning as the respective parameters to
  \code{oem\_solve\_[s|d|dow()}; the explanations are just repeated
  here:
  \begin{descr}
    \kitem{A} Pointer to a \code{DOF\_MATRIX} storing the system
    matrix.
    %%
    \hyperitem{init_oem_solve:mask}{mask} Pointer to a
    \code{DOF\_SCHAR\_VEC} masking out parts of the finite element
    space: if \code{mask->vec[d] >= DIRICHLET}, then \code{A} will act
    as if the $d$-th row would be zero. Compare also the discussion in
    the section about Dirichlet boundary condition, see
    \secref{S:dirichlet_bound}
    %%
    \kitem{tol} Tolerance for the residual; if the norm of the
    residual is less or equal \code{tol},
    \code{oem\_solve\_[s|d|dow]()} returns the actual iterate as the
    approximative solution of the system.
    %%
    \kitem{precon} A pointer to a structure describing the
    preconditioner to use, see further below in \secref{S:precon}.
    \begin{compatibility}
      Previous versions used a simple number here, but as the
      preconditioner frame-work has become much more complicated,
      because of the support for direct sums of finite element spaces,
      the code for the selection of the preconditioner has been
      separated from the entry-point to the solvers.
    \end{compatibility}
    %%
    \kitem{restart} Only used by \code{gmres}: the maximum dimension
    of the Krylov-space.
    %%
    \kitem{max\_iter} Maximal number of iterations to be performed by
    the linear solver. This can be compared with the return value --
    which gives the number of iterations actually performed -- to
    determine whether the solver has achieved its goal.
    %%
    \kitem{info} This is the level of information of the linear solver; 
    \code{0} is the lowest level of information (no information is printed)
    and \code{10} the highest level.
  \end{descr}
\item[Return Value]~\hfill

  A pointer to an initialized \code{OEM\_DATA}-structure, see the
  source-code listing on page \pageref{D:OEM_DATA}.
\end{function}

\begin{function}{release\_oem\_solve()}
\label{S:release_oem_solve_fct}

%%
\fdx{release_oem_solve()@{\code{release\_oem\_solve()}}}%
\idx{linear solvers!release_oem_solve()@{\code{release\_oem\_solve()}}}%
%%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
release_oem_solve(oem_data_handle);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Release an \code{OEM\_DATA}-handle previously acquired by a call to
  \code{init\_oem\_solve\_[s|d|dow]()} as explained above.
\item[Parameters]~\hfill
  \begin{descr}
    \kitem{oem\_data\_handle} The \code{OEM\_DATA}-pointer to destroy.
  \end{descr}
\end{function}

\begin{function}{call\_oem\_solve\_[s|d|dow]()}
\label{S:call_oem_solve_fct}

%%
\fdx{call_oem_solve_s()@{\code{call\_oem\_solve\_s()}}}%
\idx{linear solvers!call_oem_solve_s()@{\code{call\_oem\_solve\_s()}}}%
%%
\fdx{call_oem_solve_d()@{\code{call\_oem\_solve\_d()}}}%
\idx{linear solvers!call_oem_solve_d()@{\code{call\_oem\_solve\_d()}}}%
%%
\fdx{call_oem_solve_dow()@{\code{call\_oem\_solve\_dow()}}}%
\idx{linear solvers!call_oem_solve_dow()@{\code{call\_oem\_solve\_dow()}}}%
%%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
iterations = call_oem_solve_[s|d|dow](oem_data_handle, solver, f, u);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Call an iterative solver, as indicated by \code{solver}, trying to
  solve the linear system described by \code{oem\_data\_handle} for
  the unknown \code{u}, given the load-vector \code{f}. \code{u} is at
  the same time the storage for the solution and the initial guess for
  the iterative solver.
\item[Parameters]~\hfill

  With the exception of \code{oem\_data\_handle} the parameters have
  the same meaning as the respective parameters to
  \code{oem\_solve\_[s|d|dow()}; the explanations are just repeated
  here:
  \begin{descr}
    \kitem{oem\_data\_handle} A \code{OEM\_DATA}-structure, as
    returned by a previous call to \code{init\_oem\_solve()} (or
    filled in ``by hand'').
    %%
    \kitem{f} A pointer to a \code{DOF\_REAL[\_D]\_VEC[\_D]} storing
    the load-vector of the linear system.
    %%
    \kitem{u} A pointer to a \code{DOF\_REAL[\_D]\_VEC[\_D]} storing
    the initial guess on input and the approximative solution on
    output. In the context of interpolated Dirichlet boundary
    conditions special provisions have to be taken for the
    ``Dirichlet-nodes''. Compare the discussion in
    \secref{S:dirichlet_bound}.
    %%
    \kitem{solver} Use the respective OEM-solver; see above for the
    available keywords.
  \end{descr}
\item[Return Value]~\hfill

  The number of iterations the solver needed until the norm of the
  residual was below \code{tol}, or \code{max\_iter} if the solver was
  not able to reach its goal before the prescribed maximum iteration
  count was exhausted.
\end{function}

\hrulefill

\begin{example}
  \label{example:oem_solve_impl}
The high-level function
%%
\bv\begin{lstlisting}
iterations = oem_solve_[s|d|dow](A, mask, f, u,
                                 solver, tol, precon,
                                 restart, max_iter, info);
\end{lstlisting}\ev
%%
is implemented as follows:
%%
\bv\begin{lstlisting}
int oem_solve_s(const DOF_MATRIX *A, const DOF_SCHAR_VEC *mask,
		const DOF_REAL_VEC *f,  DOF_REAL_VEC *u,
		OEM_SOLVER solver, REAL tol, const PRECON *precon,
		int restart, int max_iter, int info)
{
  const OEM_DATA *oem;
  int iter;
  
  oem = init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
  iter = call_oem_solve_s(oem, solver, f, u);
  release_oem_solve(oem);

  return iter;
}
\end{lstlisting}\ev
%%
\end{example}

\begin{example}
  \label{example:get_oem_solver}
  If it is clear which solver to use, then the call through
  \code{call\_oem\_solve\_[s|d|dow]()} in
  \exampleref{example:oem_solve_impl} can be replaced by a direct call
  to the solver-routine like follows. Note, however, that this is a
  simplified example which does not take into account that
  \code{u->fe\_space} could be a direct sum of finite element spaces,
  as explained in \secref{S:chain_impl}. Of course, it is just fine
  for application to ignore the ``direct sum'' feature if it is clear
  that it is not needed. See
  \exampleref{example:get_oem_solver_direct_sums} for an example of
  how to deal with direct sums. The reader should also remember that
  -- for simple applications -- it is sufficient to use the high-level
  routines \code{oem\_solve\_[s|d|dow]()}, see also
  \exampleref{example:oem_solve_impl} for the connection between the
  example given here and the high-level routines.
  %%
  \bv\begin{lstlisting}
const OEM_DATA *oem;
int iter;
OEM_MV_FCT solver_fct;
int dim;
  
oem        = init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
solver_fct = get_oem_solver(CG); /* e.g. */
dim        = dof_real_vec_length(u->fe_space);
FOR_ALL_FREE_DOFS(u->fe_space->admin,
  if (dof < dim) u->vec[dof] = f->vec[dof] = 0.0);
...
solver_fct(oem, dim, f->vec, u->vec); /* maybe do this multiple times ... */
...
FOR_ALL_FREE_DOFS(u->fe_space->admin,
  if (dof < dim) f_other->vec[dof] = 0.0);
solver_fct(oem, dim, f_other->vec, u->vec); /* ... with other load-vectors */
...
release_oem_solver();
  \end{lstlisting}\ev
%%
\end{example}
%%
\begin{example}
  \label{example:get_oem_solver_direct_sums}
  A similar code-skeleton, taking direct sums of finite element spaces
  into account (see \secref{S:dof_vec_skel}) would look like as quoted
  below. The interested reader maybe also wants to have a look at the
  source code \code{alberta-VERSION/alberta/src/Common/oem\_solve.c}
  in the \ALBERTA distribution.  See
  \exampleref{example:get_oem_solver} for a simpler example, ignoring
  that ``direct sum'' feature. The reader should also remember that --
  for simple applications -- it is sufficient to use the high-level
  routines \code{oem\_solve\_[s|d|dow]()}, see also
  \exampleref{example:oem_solve_impl} for the connection between the
  example given here and the high-level routines.
%%
  \bv\begin{lstlisting}
const OEM_DATA *oem;
int iter;
OEM_MV_FCT solver_fct;
int dim;
REAL *uvec, *fvec;
  
oem        = init_oem_solve(A, mask, tol, precon, restart, max_iter, info);
solver_fct = get_oem_solver(CG); /* e.g. */
dim        = dof_real_vec_length(u->fe_space);
if (!CHAIN_SINGLE(u)) {
  uvec = MEM_ALLOC(dim, REAL);
  fvec = MEM_ALLOC(dim, REAL);
  copy_from_dof_real_vec(uvec, u);
  copy_from_dof_real_vec(fvec, f);
} else {
  FOR_ALL_FREE_DOFS(u->fe_space->admin,
      if (dof < dim) u->vec[dof] = f->vec[dof] = 0.0);
  fvec = f->vec;
  uvec = u->vec;
}
...
solver_fct(oem, dim, fvec, uvec);
...
release_oem_solver();
if (!CHAIN_SINGLE(u)) {
  copy_to_dof_real_vec(u, uvec);
  MEM_FREE(uvec, dim, REAL);
  MEM_FREE(fvec, dim, REAL);
}
  \end{lstlisting}\ev
\end{example}

\subsection{SOR solvers for DOF-matrices and -vectors}
The SOR and SSOR methods are implemented directly for linear systems
defined by \code{DOF\_MATRIX} and \code{DOF\_REAL\_[D\_]VEC[\_D]}s.
\begin{remark}
  In contrast to the other solvers for linear systems, the SOR- and
  SSOR-methods described in this section do \emph{not} support direct
  sums of finite element spaces (see \secref{S:chain_impl}).
\end{remark}
\fdx{sor_d()@{\code{sor\_d()}}}%
\idx{linear solvers!sor_d()@{\code{sor\_d()}}}%
\fdx{sor_s()@{\code{sor\_s()}}}%
\idx{linear solvers!sor_s()@{\code{sor\_s()}}}%
\fdx{ssor_d()@{\code{ssor\_d()}}}%
\idx{linear solvers!ssor_d()@{\code{ssor\_d()}}}%
\fdx{ssor_s()@{\code{ssor\_s()}}}%
\idx{linear solvers!ssor_s()@{\code{ssor\_s()}}}%
\bv\begin{lstlisting}
int sor_s(DOF_MATRIX *, const DOF_REAL_VEC *, const DOF_SCHAR_VEC *, 
          DOF_REAL_VEC *, REAL, REAL, int, int);
int sor_d(DOF_MATRIX *, const DOF_REAL_D_VEC *, const DOF_SCHAR_VEC *,
          DOF_REAL_D_VEC *, REAL, REAL, int, int);
int ssor_s(DOF_MATRIX *, const DOF_REAL_VEC *, const DOF_SCHAR_VEC *, 
          DOF_REAL_VEC *, REAL, REAL, int, int);
int ssor_d(DOF_MATRIX *, const DOF_REAL_D_VEC *, const DOF_SCHAR_VEC *,
           DOF_REAL_D_VEC *, REAL, REAL, int, int);
\end{lstlisting}\ev
\begin{descr}
\kitem{[s]sor\_[s,d](matrix, f, bound, u, omega, tol, max\_iter, info)}
       solves the linear system  for a scalar or decoupled vector valued 
       problem in \ALBERTA by the [Symmetric] Successive
       Over Relaxation method; the return value is the number of used
       iterations to reach the prescribed tolerance;

       \code{matrix}: pointer to a DOF matrix storing the system
       matrix;

       \code{f}: pointer to a DOF vector storing the right hand side
       of the system;

       \code{bound}: optional pointer to a DOF vector giving Dirichlet
       boundary information;

       \code{u}: pointer to a DOF vector storing the initial 
       guess on input and the calculated solution on output;

       \code{omega}: the relaxation parameter and must be in the
       interval $(0,2]$; if it is not in this interval then
       \code{omega=1.0} is used;

       \code{tol}: tolerance for the maximum norm of the correction;
       if this norm is less than or equal to \code{tol}, then
       \code{sor\_[s,d]()} returns the actual iterate as the solution
       of the system; 

       \code{max\_iter}: maximal number of iterations to be performed
       by \code{sor\_[s,d]()} although the tolerance may not be
       reached;

       \code{info}: level of information of \code{sor\_[s,d]()};
       \code{0} is the lowest level of information (no information is printed)
       and \code{6} the highest level.
\end{descr}

\subsection{Saddle-point problems, CG solver for Schur's complement}
\label{S:OEM_SPCG}
On the linear-algebra level, a linaer saddle-point problem is of the
form
%%
\begin{equation}
  \label{eq:spproblem}
  \begin{bmatrix}
    A      & B \\
    B^\ast & 0
  \end{bmatrix}
  \,
  \begin{bmatrix}
    v \\ p
  \end{bmatrix}
  =
  \begin{bmatrix}
    f \\ g
  \end{bmatrix},
  \quad f,\,v\in X,\; g,\,p\in Y,
\end{equation}
%%
with matrices $A$ and $B$, unknown vectors $v$ and $p$ and a load
vector consisting of the vector $f$ and $g$. Usually, $A$ has its
origin in the discretization of an unconstraint minimization problem,
$B^\ast$ plays the role of a linear constraint, and $p$ is the
corresponding Lagrangian multiplier. $Y$ is the finite element space
for the Lagrangian multiplier, and $X$ a possibly different finite
element space for the principal unknown $v$:

If $A$ is invertible, then it is possible to transform
\eqref{eq:spproblem} into an equation for $p$ only:
%%
\begin{equation}
  \label{eq:spschur}
  T\, p = B^\ast\,A^{-1}\,f  - g,\quad T:= B^\ast\,A^{-1}\,B,
\end{equation}
%%
where $v$ can be reconstructed from $p$ by $v=A^{-1}(f - B\,p)$.  If
$A$ is symmetric positive definite, then so is $T$, and thus it is
possible to solve \eqref{eq:spschur} by means of a CG-method in this
case which, interestingly, even computes $v$ as a by-product.

In the same spirit as for the iterative solvers for ``ordinary''
problems, this \emph{SPCG}-method is implemented in a fairly abstract
manner, using a special data-structure to describe the saddle-point
problem. The actual CG-iteration is executed by a call to the function
\code{oem\_spcg(oem\_sp\_data,\dots)}, described below in
\secref{S:oem_spcg_fct}. It is the task of the application to fill
that \code{OEM\_SP\_DATA}-structure (see \secref{S:OEM_SP_DATA_struct}
below).  However, there are interface functions to aid the
implementation of such a saddle-point solver with \ALBERTA's
DOF-matrices and -vectors, see \secref{S:ALBERTA_SPCG} below.
%%
\fdx{oem_spcg()@{\code{oem\_spcg()}}}%
\ddx{OEM_SP_DATA()@{\code{OEM\_SP\_DATA}}}%
\bv\begin{lstlisting}
typedef struct oem_sp_data OEM_SP_DATA;
int oem_spcg(OEM_SP_DATA *data, int dimX, const REAL *f, REAL *u, int dimY,
	     const REAL *g, REAL *p);
\end{lstlisting}\ev
%%
\begin{datatype}{OEM\_SP\_DATA}

\item[Definition]~\hfill

\ddx{OEM_SP_DATA()@{\code{OEM\_SP\_DATA}}}%
\bv\begin{lstlisting}
typedef int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);
typedef void (*OEM_GEMV_FCT)(void *data,
			     REAL factor,
			     int dimX, const REAL *x, int dimY, REAL *y);

typedef struct oem_sp_data OEM_SP_DATA;
struct oem_sp_data
{
  OEM_MV_FCT   solve_Auf;
  void         *solve_Auf_data;
  OEM_GEMV_FCT B;
  void         *B_data;
  OEM_GEMV_FCT Bt;
  void         *Bt_data;
  OEM_MV_FCT   project;
  void         *project_data;
  int          (*precon)(void *ud,
			 int dimY, const REAL *g_Btu, const REAL *r, REAL *Cr);
  void         *precon_data;
  
  WORKSPACE    *ws;

  REAL         tolerance;
  int          restart;
  int          max_iter;
  int          info;

  REAL         initial_residual;
  REAL         residual;
};
\end{lstlisting}\ev

\item[Components]~\hfill

\label{S:OEM_SP_DATA_struct}
\begin{descr}
\kitem{solve\_Auf()} An application provided function for solving
  $A\,x=b$, for given initial guess and solution $x$ and load-vector
  $b$. This can, e.g. be one of the solver-functions for ordinary
  problems, see \secref{S:oem}.
  %% 
\kitem{solve\_Auf\_data} Application data passed to \code{solve\_Auf()}
  as first argument. If \code{solve\_Auf()} is one of the
  solver-functions described in \secref{S:oem} or a function pointer
  returned by \code{get\_oem\_solver()}, then this should be a
  pointer to a \code{OEM\_DATA} structure, as returned for instance by
  \code{init\_oem\_solver()}, see above in \secref{S:ALBERTA_OEM_solvers}.
  %% 
\hyperitem{OEM_SP_DATA:B}{B()} A pointer to an application provided
  function with the calling convention
  %% 
  \bv\begin{lstlisting}
B(B_data, factor, dimY, p, dimX, v);
\end{lstlisting}\ev
  %% 
  This function must implement the operation
  $\code{v} = \code{v} + \code{factor} \,\code{B},\code{p}$ In
  the abstract setting the range of the operator underlying \code{B()}
  is the same as the range of the unconstrained operator \code{A()}
  such that it makes sense to apply the inverse of \code{A()} to the
  result of \code{B()}.
  %% 
\hyperitem{OEM_SP_DATA:B_data}{B\_data} Application data passed to
  \code{B()} as first argument.
  %% 
\hyperitem{OEM_SP_DATA:Bt}{Bt()} The pendent to \code{B()}: A pointer
  to an application provided function with the calling convention
  %% 
  \bv\begin{lstlisting}
Bt(Bt_data, factor, dimX, v, dimY, p);
\end{lstlisting}\ev
  %% 
  This function must implement the operation
  $\code{p} = \code{p} + \code{factor} \,\code{B}^\ast\,\code{v}$
  %% 
  For practical reasons -- e.g. in the context of a Stokes problem --
  the range of the discrete operator \code{Bt()} need not necessarily
  be the finite element space for the Lagrangian multiplier (see
  \code{precon()} and \code{project()} below), but often is rather the
  dual of that space.
  %% 
\hyperitem{OEM_SP_DATA:Bt_data}{Bt\_data} Application data passed to
  \code{Bt()} as first argument.
  %% 
\hyperitem{OEM_SP_DATA:project}{project()} A function pointer pointing
  to an application provided function which has the task to project
  the result from \code{Bt()} to the finite element space for the
  Lagrangian multiplier. \code{project()} maybe \nil is such a
  projection is not needed.  Arguably, this could already have been
  incorporated into \code{Bt()}, however, it is sometimes more
  efficient to let the discrete operator \code{Bt()} map to the dual
  of the space for the Lagrangian multiplier. See also \code{precon()}
  below. Often \code{project()} will just be an $L^2$-projection
  involving the inversion of a mass-matrix, which can for instance be
  done by a CG-method, or maybe even more efficiently with
  mass-lumping.
  %% 
\hyperitem{OEM_SP_DATA:project_data}{project\_data} Application data
  pointer passed as first argument to \code{project()}.
  %% 
\hyperitem{OEM_SP_DATA:precon}{precon()} A function
  pointer pointing to an application provided function which should
  implement a preconditioner ``$C()$'' for the CG-method for Schur's
  complement operator. \code{precon()} may be \nil. The calling
  convention is
  %% 
  \bv\begin{lstlisting}
iterations = precon(precon_data, dim, g_Btu, r, Cr);
\end{lstlisting}\ev
  %% 
  where the non-self-explanatory arguments have the following meaning:
  %% 
  \begin{descr}
  \kitem{g\_Btu} The current value of $g - B^\ast u$, where $u$ is the
    current iterate for the principal unknown $v$ in the CG-method.
    This is the result of a call to \code{Bt()}, and most likely lives
    in the dual of the space for the Lagrangian multiplier.
    %% 
  \kitem{r} This is \code{project(g\_Btu)}, this lives in space for
    the Lagrangian multiplier.
    %% 
  \end{descr}
  %% 
  The result value of \code{precon()} must be stored in \code{Cr}.
  \code{Cr} must belong to the space for the Lagrangian multiplier.
  %% 
  As an example, it is known that for a Quasi-Stokes problem
  \[
  \mu\,u - \nu\,\Delta u +\nabla p = f,\quad \nabla\cdot u = 0,
  \]
  a good choice for a preconditioner for Schur's complement CG-method
  is
  \[
  C(\code{r}) = \nu\,\code{r} + \mu\,q,\text{ with }-\Delta q = \code{g\_Btu}.
  \]
  Note that we have omitted the boundary conditions, which, of course,
  have to be applied to close the differential equations mentioned
  above. The reader is referred to standard text-books dealing with
  the discretizations of saddle-point problems.
  %% 
\hyperitem{OEM_SP_DATA:precon_data}{precon\_data} Data pointer passed
  as first argument to \code{precon()}.
  %% 
\kitem{ws} A pointer to a work-space area. May be \nil. If supplied,
  it must point to an initialized work-space of size
  %% 
  $$2*\code{dimY}+\code{dimX}+\max(\code{dimX}, \code{dimY})$$
  %% 
  if \code{precon() == \nil} and
  %% 
  $$3*\code{dimY}+\code{dimX}+\max(\code{dimX}, \code{dimY}).$$
  %% 
  If \code{ws == \nil}, then \code{oem\_spcg()} will allocate a
  work-space area by itself.
  %% 
\kitem{tolerance} \code{oem\_spcg()} will terminate if the norm of the
  CG-residual for the Lagrangian multiplier falls below
  \code{tolerance}.
  %% 
\kitem{restart} Not used by \code{oem\_spcg()}. Could be used when
  implementing similar iterative methods for non-symmetric
  saddle-point problems, e.g. by means of applying GMRES.
  %% 
\kitem{max\_iter} \code{ome\_spcg()} will terminate after this many
  iterations of the main-loop.
  %% 
\kitem{info} An integer controlling the amount of information printed
  to the terminal the application program runs in.
  %% 
  \hyperitem{OEM_SP_DATA:initial_residual}{initial\_residual} Output.
  Upon return from \code{oem\_spcg()} this component stores the
  initial residual.
  %% 
  \hyperitem{OEM_SP_DATA:residual}{residual} Output. Upon return from
  \code{oem\_spcg()} this component stores the final residual. This
  could be used for error recovery, e.g. if the iteration terminates
  because the maximum number of iterations (as specified by
  \code{max\_iter}) was exhausted.
  %% 
\end{descr}
\end{datatype}

\begin{function}{oem\_spcg()}
\label{S:oem_spcg_fct}

\fdx{oem_spcg()@{\code{oem\_spcg()}}}%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
iterations = oem_spcg(sp_data, dimX, f, u, dimY, g, p);
\end{lstlisting}\ev
\item[Description] ~\hfill

  This function implement a CG-method for the inversion of Schur's
  complement operator for a linear symmetric saddle-point problem.
\item[Parameters] ~\hfill
\begin{descr}
\kitem{sp\_data} A pointer to a correctly filled
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}} structure, as
  explained above. Upon return from \code{oem\_spcg()}, the fields
  \hyperlink{OEM_SP_DATA:initial_residual}{\code{initial\_residual}}
  and \hyperlink{OEM_SP_DATA:residual}{\code{residual}} will contain
  the initial and the final residual of the CG-iterations.
%%
\kitem{dimX} The dimension of the space for the principal unknown \code{u}.
%%
\kitem{f} Load-vector for the principal equation.
%%
\kitem{u} Storage for the principal unknown, and start-value for the
  principal unknown for the CG-method.
%%
\kitem{dimY} Dimension of the space for the Lagrangian multiplier.
%%
\kitem{g} Load-vector for the constraint equation.
%%
\kitem{p} Storage for the Lagrangian multiplier and start-value for
  the CG-method.
\end{descr}
\item[Return Value] ~\hfill

  The number of times the main-loop of the CG-iteration was executed.
  If this is equal to \code{sp\_data->max\_iter}, then the application
  should also inspect \code{sp\_data->residual} to determine whether
  the approximative solution is still acceptable.
\end{function}

\subsection{Saddle-pointer solvers for DOF-matrices and -vectors}
\label{S:ALBERTA_SPCG}

Similar to the functions explained in \secref{S:ALBERTA_OEM_solvers}
there are also interface functions to mediate between the more
low-level \hyperref[S:oem_spcg_fct]{oem\_spcg()} function described in
the previous \secref{S:OEM_SPCG} and the DOF-vectors and -matrices
generated by \ALBERTA's assemble frame-work, as described in
\secref{S:ass_tools}. The functions below have the slight disadvantage
that they take too many arguments. The interface functions support
direct sums of finite element spaces (see \secref{S:chain_impl}) which
is of some importance in the context of mixed discretizations for the
Stokes-problem.

There are two interfaces available: one for a saddle-point problem
with a single linear constraint, and one for a saddle-point problem
with multiple linear constraints, with the restriction that the
constraints are decoupled. We start with the single-constraint version
\hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}} in
\secref{S:oem_sp_solve_fct} and continue with the multiple-constraint
functions
\hyperref[S:init_sp_constraint_fct]{\code{init\_sp\_constraint()}},
\hyperref[S:release_sp_constraint_fct]{\code{release\_sp\_constraint()}}
and
\hyperref[S:oem_sp_schur_solve_fct]{\code{oem\_sp\_schur\_solve()}} in
the Sections
\ref{S:init_sp_constraint_fct}-\ref{S:oem_sp_schur_solve_fct}. There
is one additional support function
\hyperref[S:sp_dirichlet_bound_fct]{\code{sp\_dirichlet\_bound()}}
which deals with compatibility conditions in the context of a
divergence constraint and Dirichlet boundary conditions, see
\secref{S:sp_dirichlet_bound_fct}.

The suite of demo-programs contains example programs for the
discretization of Stokes and Quasi-Stokes problems, the interested
reader is referred to the programs
\bv\begin{verbatim}
alberta-VERSION-demo/src/Common/stokes.c
\end{verbatim}\ev
and
\bv\begin{verbatim}
alberta-VERSION-demo/src/Common/quasi-stokes.c.
\end{verbatim}\ev

The prototypes for the available functions read as follows:
%%
\fdx{oem_sp_solve_dow_scl()@{\code{oem\_sp\_solve\_dow\_scl()}}}%
\idx{linear solvers!oem_sp_solve_dow_scl()@{\code{oem\_sp\_solve\_dow\_scl()}}}%
\fdx{oem_sp_solve_ds()@{\code{oem\_sp\_solve\_ds()}}}%
\idx{linear solvers!oem_sp_solve_ds()@{\code{oem\_sp\_solve\_ds()}}}%
%%
\fdx{sp_dirichlet_bound_dow_scl()@{\code{sp\_dirichlet\_bound\_dow\_scl()}}}%
\idx{linear solvers!sp_dirichlet_bound_dow_scl()@{\code{sp\_dirichlet\_bound\_dow\_scl()}}}%
%%
\fdx{sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}}%
\idx{linear solvers!sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}}%
%%
\ddx{SP_CONSTRAINT@{\code{SP\_CONSTRAINT}}}%
\idx{linear solvers!SP_CONSTRAINT@{\code{SP\_CONSTRAINT}}}%
%%
\fdx{oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}}%
\idx{linear solvers!oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}}%
%%
\fdx{init_sp_constraint()@{\code{init\_sp\_constraint()}}}%
\idx{linear solvers!init_sp_constraint()@{\code{init\_sp\_constraint()}}}%
%%
\fdx{release_sp_constraint()@{\code{release\_sp\_constraint()}}}%
\idx{linear solvers!release_sp_constraint()@{\code{release\_sp\_constraint()}}}%
%%
\bv\begin{lstlisting}
int oem_sp_solve_dow_scl(OEM_SOLVER sp_solver,
                         REAL sp_tol, REAL tol_incr,
                         int sp_max_iter, int sp_info,
                         const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
                         OEM_SOLVER A_solver,
                         int A_max_iter, const PRECON *A_precon,
                         DOF_MATRIX *B,
                         DOF_MATRIX *Bt,
                         DOF_MATRIX *Yproj,
                         OEM_SOLVER Yproj_solver,
                         int Yproj_max_iter, const PRECON *Yproj_precon,
                         DOF_MATRIX *Yprec,
                         OEM_SOLVER Yprec_solver,
                         int Yprec_max_iter, const PRECON *Yprec_precon,
                         REAL Yproj_frac, REAL Ypre_frac,
                         const DOF_REAL_VEC_D *f,
                         const DOF_REAL_VEC *g,
                         DOF_REAL_VEC_D *x,
                         DOF_REAL_VEC *y);
int oem_sp_solve_ds(OEM_SOLVER sp_solver,
                    REAL sp_tol, REAL tol_incr,
                    int sp_max_iter, int sp_info,
                    const DOF_MATRIX *A, const DOF_SCHAR_VEC *bound,
                    OEM_SOLVER A_solver,
                    int A_max_iter, const PRECON *A_precon,
                    DOF_MATRIX *B,
                    DOF_MATRIX *Bt,
                    DOF_MATRIX *Yproj,
                    OEM_SOLVER Yproj_solver,
                    int Yproj_max_iter, const PRECON *Yproj_precon,
                    DOF_MATRIX *Yprec,
                    OEM_SOLVER Yprec_solver,
                    int Yprec_max_iter, const PRECON *Yprec_precon,
                    REAL Yproj_frac, REAL Ypre_frac,
                    const DOF_REAL_D_VEC *f,
                    const DOF_REAL_VEC *g,
                    DOF_REAL_D_VEC *x,
                    DOF_REAL_VEC *y);

REAL sp_dirichlet_bound_dow_scl(MatrixTranspose transpose,
                                const DOF_MATRIX *Bt,
                                const DOF_SCHAR_VEC *bound,
                                const DOF_REAL_VEC_D *u_h,
                                DOF_REAL_VEC *g_h);
REAL sp_dirichlet_bound_ds(MatrixTranspose transpose,
                           const DOF_MATRIX *Bt,
                           const DOF_SCHAR_VEC *bound,
                           const DOF_REAL_D_VEC *u_h,
                           DOF_REAL_VEC *g_h);

typedef struct sp_constraint
{
  const DOF_MATRIX    *B, *Bt;
  const DOF_SCHAR_VEC *bound;
  OEM_MV_FCT          project;
  OEM_DATA            *project_data;
  OEM_MV_FCT          precon;
  OEM_DATA            *precon_data;
  REAL                proj_factor, prec_factor;
} SP_CONSTRAINT;

SP_CONSTRAINT *init_sp_constraint(const DOF_MATRIX *B,
                                  const DOF_MATRIX *Bt,
                                  const DOF_SCHAR_VEC *bound,
                                  REAL tol, int info,
                                  const DOF_MATRIX *Yproj,
                                  OEM_SOLVER Yproj_solver,
                                  int Yproj_max_iter,
                                  const PRECON *Yproj_prec,
                                  const DOF_MATRIX *Yprec,
                                  OEM_SOLVER Yprec_solver,
                                  int Yprec_max_iter,
                                  const PRECON *Yprec_prec,
				  void (*Yprec_bndry)(void *data,
						      const DOF_REAL_VEC *r,
						      DOF_REAL_VEC *mod_r,
						      DOF_REAL_VEC *Cr),
				  void *Yprec_bndry_data,
                                  REAL Yproj_frac, REAL Yprec_frac);
void release_sp_constraint(SP_CONSTRAINT *constraint_data);
int oem_sp_schur_solve(OEM_SOLVER sp_solver,
                       REAL sp_tol, int sp_max_iter, int sp_info,
                       OEM_MV_FCT principal_inverse,
                       OEM_DATA *principal_data,
                       const DOF_REAL_VEC_D *f,
                       DOF_REAL_VEC_D *u,
                       SP_CONSTRAINT *constraint,
                       const DOF_REAL_VEC *g,
                       DOF_REAL_VEC *p,
                       ...);
\end{lstlisting}\ev

\begin{function}{oem\_sp\_solve\_[dow\_scl|ds]()}
\label{S:oem_sp_solve_fct}

%%
\fdx{oem_sp_solve_dow_scl()@{\code{oem\_sp\_solve\_dow\_scl()}}}%
\idx{linear solvers!oem_sp_solve_dow_scl()@{\code{oem\_sp\_solve\_dow\_scl()}}}%
\fdx{oem_sp_solve_ds()@{\code{oem\_sp\_solve\_ds()}}}%
\idx{linear solvers!oem_sp_solve_ds()@{\code{oem\_sp\_solve\_ds()}}}%
%%
\item[Synopsis]~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
iterations = oem_sp_solve_[dow_scl|ds](
  sp_solver,
  sp_tol, tol_incr, sp_max_iter, sp_info,
  A, mask, A_solver, A_max_iter, A_precon,
  B, Bt,
  Yproj, Yproj_solver, Yproj_max_iter, Yproj_precon,
  Yprec, Yprec_solver, Yprec_max_iter, Yprec_precon,
  Yproj_frac, Yprec_frac,
  f, g, x, y);
\end{lstlisting}\ev

\item[Description]~\hfill

  This function implements an interface between the DOF-vector and
  -matrix level and the low-level functions described in
  \secref{S:OEM_SPCG} above. Internally, \code{oem\_sp\_solve()} emits
  calls to \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solve()}}
  and initializes the support data-structure
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}}. Then finally
  the function \hyperref[S:oem_spcg_fct]{\code{oem\_spcg()}} is
  called, see also \secref{S:OEM_SPCG}.
  
  \code{oem\_sp\_solve()} implements a preconditioner $C$ of the form
  \begin{equation}
    \label{eq:sppreconimpl}
    C(r) = \code{Yproj\_frac}\ast\code{Yproj}(r) +
    \code{Yprec\_frac}\ast\code{Yprec}^{-1}(r),
  \end{equation}
  which has the form of the usual preconditioner for a Quasi-Stokes
  problem, which was already mentioned in the explanation for the
  parameter \hyperlink{oem_spcg_fct:precon}{\code{precon()}} for the
  function \hyperref[S:oem_spcg_fct]{\code{oem\_spcg()}}, see
  \secref{S:oem_spcg_fct}.

\item[Parameters]~\hfill

\hyperitem{oem_sp_solve:sp_solver}{sp\_solver} The solver used for the
  \emph{outer} iteration.  Currently, only a CG-method for a symmetric
  and positive (semi-) definite Schur's complement operator is
  implemented, so \code{sp\_solver} must equal the symbol \code{CG}.
  %%
\hyperitem{oem_sp_solve:sp_tol}{sp\_tol} The tolerance for the
  \emph{outer} CG-loop.
  %%
\hyperitem{oem_sp_solve:tol_incr}{tol\_incr} A decrease in tolerance
  for the iterative solvers for the sub-problems, like inverting the
  principal part \code{A} of the operator. The tolerances for the
  solvers for the sub-problems will be \code{sp\_tol / tol\_incr}.
  %%
\kitem{sp\_max\_iter} The maximum number of iterations for the outer
  CG-loop.
  %%
\hyperitem{oem_sp_solve:sp_info}{sp\_info} The verbosity level. The
  solvers for the sub-problems will inherit a decreased verbosity
  level of $\max(0,\,\code{sp\-info}-3)$.
  %%
\kitem{A} The matrix for the principal part of the saddle-point problem.
  %%
\hyperitem{oem_sp_solve:bound}{bound} A \code{DOF\_SCHAR\_VEC} used to
  exclude DOFs from the operation of the matrix-vector routines. See
  semantics are as explained in the explanations for the argument
  \hyperlink{init_oem_solve:mask}{\code{mask}} to the function
  \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solve()}}, see
  \secref{S:init_oem_solve_fct}.
  %%
\kitem{A\_solver} The solver to use to invert \code{A}, compare with
  the explanations for
  \hyperref[S:get_oem_solver_fct]{\code{get\_oem\_solver()}} in
  \secref{S:get_oem_solver_fct} and the parameter
  \hyperlink{oem_solve_fct:solver}{\code{solver}} to
  \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solver()}}.
  %%
\kitem{A\_max\_iter} The maximum number of iterations for the linear
  solver used to invert \code{A}.
  %%
\kitem{A\_precon} A pointer to the descriptor for the preconditioner
  to use for the inversion of \code{A}, see \secref{S:precon} below.
  %%
\hyperitem{oem_sp_solve:B}{B} A pointer to the matrix implementing $B$, see
  \eqref{eq:spproblem}.
  %%
\hyperitem{oem_sp_solve:Bt}{Bt} A pointer to the matrix implementing
  $B^\ast$, see \eqref{eq:spproblem}. \code{Bt} may be \nil, in which
  case the matrix \code{B} is used, passing the
  \hyperref[enum:MatrixTranspose]{\code{Transpose}} flag to the
  matrix-vector routines, see \secref{S:DOF_BLAS}. An application
  calling \code{oem\_sp\_solve()} with \code{Bt == \nil} most likely
  will want to make use of the optional parameter \code{mask} above in
  order to implement Dirichlet boundary conditions.
  %%
\kitem{Yproj} The matrix for the back-projection of the result from
  applying \code{Bt} to the finite element space for the constraint.
  Compare the remarks in the explanation of the component
  \hyperlink{OEM_SP_DATA:project}{\code{project()}} of the
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}} structure.
  %%
\hyperitem{oem_sp_solve:Yproj_solver}{Yproj\_solver} The solver to use
  for inverting \code{Yproj}.
  %%
\hyperitem{oem_sp_solve:Yproj_max_iter}{Yproj\_max\_iter} The maximum
  number of iterations for inverting \code{Yproj}.
  %%
\hyperitem{oem_sp_solve:Yproj_precon}{Yproj\_precon} The
  preconditioner for the iterative solver for the inversion of
  \code{Yproj}. See \secref{S:precon} below.
  %%
\hyperitem{oem_sp_solve:Yprec}{Yprec} A part defining one part of the
  preconditioner as explained in equation \eqref{eq:sppreconimpl}.
  Maybe \nil, in which case no preconditioner will be applied in the
  outer CG-loop for inverting Schur's complement.
  %%
\hyperitem{oem_sp_solve:Yprec_solver}{Yprec\_solver} The solver to use
  for inverting \code{Yprec}.
  %%
\hyperitem{oem_sp_solve:Yprec_max_iter}{Yprec\_max\_iter} The maximum
  number of iterations for inverting \code{Yprec}.
  %%
\hyperitem{oem_sp_solve:Yprec_precon}{Yprec\_precon} The
  preconditioner for the iterative solver for the inversion of
  \code{Yprec}. See \secref{S:precon} below.
  %%
\hyperitem{oem_sp_solve:Yproj_frac}{Yproj\_frac} See equation
  \eqref{eq:sppreconimpl} above.
  %%
\hyperitem{oem_sp_solve:Yprec_frac}{Yprec\_frac} See equation
  \eqref{eq:sppreconimpl} above.
  %%
\hyperitem{oem_sp_solve:f}{f} The load vector for the principal
  unknown.
  %% 
\hyperitem{oem_sp_solve:g}{g} The load vector for the linear
  constraint. Even in the case when the non-discrete problem is
  subject to a homogeneous constraint, it can be necessary to impose a
  slightly inhomogeneous constraint in the discrete setting. One
  notable example is the implementation of Dirichlet boundary
  conditions in the context of a divergence constraint. In this case
  interpolated Dirichlet boundary values will in general fail to
  fulfill the compatibility condition the discrete divergence
  constraint imposes on the discrete boundary values. Compare with the
  explanations for
  \hyperref[S:sp_dirichlet_bound_fct]{\code{sp\_dirichlet\_bound()}}
  below.
  %% 
\hyperitem{oem_sp_solve:x}{x} Storage for the principal of the
  solution, and initial guess for the CG-method.
  %%
\hyperitem{oem_sp_solve:y}{y} Storage for the Lagrangian multiplier,
  and initial guess for the CG-method for Schur's complement.
\end{function}

\begin{datatype}{SP\_CONSTRAINT}
\label{S:SP_CONSTRAINT_struct}

\item[Description]~\hfill

In the multi-constraint case, each single constraint is described by a
\code{SP\_CONSTRAINT} structure, in order to reduce the number of
parameters which have to be passed to the saddle-point solver. Such a
structure can be obtained by a call to
\hyperref[S:init_sp_constraint_fct]{init\_sp\_constraint()}, see below
\secref{S:init_sp_constraint_fct}.

The meaning of the individual structure components is identical to the
meaning of the respective component of the
\hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}} or parameter of
the \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}} function,
the reader is therefore referred to \secref{S:OEM_SP_DATA_struct} and
\secref{S:oem_sp_solve_fct} for a detailed discussion.

\item[Definition]~\hfill

\ddx{SP_CONSTRAINT@{\code{SP\_CONSTRAINT}}}%
\idx{linear solvers!SP_CONSTRAINT@{\code{SP\_CONSTRAINT}}}%
%%
\bv\begin{lstlisting}
typedef struct sp_constraint
{
  const DOF_MATRIX    *B, *Bt;
  const DOF_SCHAR_VEC *bound;
  OEM_MV_FCT          project;
  void                *project_data;
  OEM_MV_FCT          precon;
  void                *precon_data;
  REAL                proj_factor, prec_factor;
} SP_CONSTRAINT;
\end{lstlisting}\ev

\item[Components]~\hfill
%%
\begin{descr}
\hyperitem{SP_CONSTRAINT:B}{B} See parameter
  \hyperlink{oem_sp_solve:B}{\code{B}} of
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
  %% 
\hyperitem{SP_CONSTRAINT:Bt}{Bt} See parameter
  \hyperlink{oem_sp_solve:Bt}{\code{Bt}} of
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
  %%
\hyperitem{SP_CONSTRAINT:bound}{bound} See parameter
  \hyperlink{oem_sp_solve:bound}{\code{bound}} of
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
  %%
\hyperitem{SP_CONSTRAINT:project}{project()} See component
  \hyperlink{OEM_SP_DATA:project}{\code{project()}} of
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}}.
  %%
\hyperitem{SP_CONSTRAINT:project_data}{project\_data} See component
  \hyperlink{OEM_SP_DATA:project_data}{\code{project\_data}} of
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}}.
  %%
\hyperitem{SP_CONSTRAINT:precon}{precon()} See component
  \hyperlink{OEM_SP_DATA:precon}{\code{precon()}} of
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}}.
  %%
\hyperitem{SP_CONSTRAINT:precon_data}{precon\_data} See component
  \hyperlink{OEM_SP_DATA:precon_data}{\code{precon\_data}} of
  \hyperref[S:OEM_SP_DATA_struct]{\code{OEM\_SP\_DATA}}.
  %%
\hyperitem{SP_CONSTRAINT:proj_factor}{proj\_factor} See parameter
  \hyperlink{oem_sp_solve:Yproj_frac}{\code{Yproj\_frac}} of
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
  %%
\hyperitem{SP_CONSTRAINT:prec_factor}{prec\_factor} See parameter
  \hyperlink{oem_sp_solve:Yprec_frac}{\code{Yprec\_frac}} of
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
  %%
\end{descr}
\end{datatype}

\begin{function}{init\_sp\_constraint()}
\label{S:init_sp_constraint_fct}

%%
\fdx{init_sp_constraint()@{\code{init\_sp\_constraint()}}}%
\idx{linear solvers!init_sp_constraint()@{\code{init\_sp\_constraint()}}}%
%%
\item[Synopsis]~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
constraint_data =
  init_sp_constraint(B, Bt, bound, tol, info,
                     Yproj, Yproj_solver, Yproj_max_iter,
                     Yproj_prec,
                     Yprec, Yprec_solver, Yprec_max_iter,
                     Yprec_prec,
                     Yprec_bndry, Yprec_bndry_data,
                     Yproj_frac, Yprec_frac);
\end{lstlisting}\ev

\item[Description]~\hfill

  Allocate and initialize a
  \hyperref[S:SP_CONSTRAINT_struct]{\code{SP\_CONSTRAINT}} structure,
  for later use with
  \hyperref[S:oem_sp_schur_solve_fct]{\code{oem\_sp\_schur\_solve()}},
  see \secref{S:oem_sp_schur_solve_fct} below. The meaning of the
  parameters is almost identical to the corresponding parameters to
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}, see
  \secref{S:oem_sp_solve_fct} above.

\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{init_sp_constraint:B}{B} See parameter
    \hyperlink{oem_sp_solve:B}{\code{B}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Bt}{Bt} See parameter
    \hyperlink{oem_sp_solve:Bt}{\code{Bt}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:bound}{bound} See parameter
    \hyperlink{oem_sp_solve:bound}{\code{bound}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:tol}{tol} The tolerance for the
    sub-solvers used to invert \code{Yproj} and \code{Yprec} (if
    present). Compare parameter
    \hyperlink{oem_sp_solve:tol_incr}{\code{tol\_incr}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:info}{info} Control the amount of
    messages printed to the terminal the application has been started
    from. Compare parameter
    \hyperlink{oem_sp_solve:sp_info}{\code{sp\_info}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yproj}{Yproj} See parameter
    \hyperlink{oem_sp_solve:Yproj}{\code{Yproj}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yproj_solver}{Yproj\_solver} See
    parameter
    \hyperlink{oem_sp_solve:Yproj_solver}{\code{Yproj\_solver}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yproj_max_iter}{Yproj\_max\_iter} See
    parameter
    \hyperlink{oem_sp_solve:Yproj_max_iter}{\code{Yproj\_max\_iter}}
    of \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yproj_prec}{Yproj\_prec} See parameter
    \hyperlink{oem_sp_solve:Yproj_prec}{\code{Yproj\_prec}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec}{Yprec} See parameter
    \hyperlink{oem_sp_solve:Yprec}{\code{Yprec}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_solver}{Yprec\_solver} See
    parameter
    \hyperlink{oem_sp_solve:Yprec_solver}{\code{Yprec\_solver}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_max_iter}{Yprec\_max\_iter} See
    parameter
    \hyperlink{oem_sp_solve:Yprec_max_iter}{\code{Yprec\_max\_iter}}
    of \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_prec}{Yprec\_prec} See parameter
    \hyperlink{oem_sp_solve:Yprec_prec}{\code{Yprec\_prec}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_frac}{Yprec\_frac} See parameter
    \hyperlink{oem_sp_solve:Yprec_frac}{\code{Yprec\_frac}}
    of \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_bndry}{Yprec\_bndry(data, r,
      mod\_r, Cr)}
    A callback for cases where the constraint has to fulfil special
    boundary conditions. \code{Yprec\_bndry} may be \nil. The first
    argument to the call-back is the application provided
    \code{Yprec\_bndry\_data}-pointer specified by the following
    argument. \code{r} is the current residual which normally serves
    as load-vector for the preconditioner (see equation
    \eqref{eq:sppreconimpl}), \code{mod\_r} is a modifiable copy of
    \code{r}, and \code{Cr} is the preconditioned residual which is
    solved for when inverting
    \hyperlink{oem_sp_constraint:Yprec}{\code{Yprec}}.
    %% 
  \hyperitem{init_sp_constraint:Yprec_bndry_data}{Yprec\_bndry\_data}
    See the description for \code{Yprec\_bndry()} above;
    \code{Yprec\_bndry\_data} is the application-data pointer for that
    callback.
    %% 
  \hyperitem{init_sp_constraint:Yprec_frac}{Yprec\_frac} See parameter
    \hyperlink{oem_sp_solve:Yprec_frac}{\code{Yprec\_frac}} of
    \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}}.
    %% 
  \end{descr}
\item[Return Value]~\hfill

  A pointer to an initialized
  \hyperref[S:SP_CONSTRAINT_struct]{\code{SP\_CONSTRAINT}} structure,
  which can be passed as argument to
  \hyperref[S:oem_sp_schur_solve_fct]{\code{oem\_sp\_schur\_solve()}}
  described in \secref{S:oem_sp_schur_solve_fct} below. The return
  structure should be deleted by a call to
  \hyperref[S:release_sp_constraint_fct]{\code{release\_sp\_constraint()}},
  see below.
\end{function}

\begin{function}{release\_sp\_constraint()}
\label{S:release_sp_constraint_fct}

%%
\fdx{release_sp_constraint()@{\code{release\_sp\_constraint()}}|(}
\idx{linear solvers!release_sp_constraint()@{\code{release\_sp\_constraint()}}|(}
%%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
release_sp_constraint(constraint_data);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Release the resources associated with a
  \hyperlink{S:SP_CONSTRAINT_stuct}{\code{SP\_CONSTRAINT}} structure
  as returned by
  \hyperlink{S:init_sp_constraint_fct}{\code{init\_sp\_constraint()}}.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{release_sp_constraint:constraint_data}{constraint\_data}
    A pointer to a
    \hyperlink{S:SP_CONSTRAINT_struct}{\code{SP\_CONSTRAINT}} structure
    previously acquired by a call to
    \hyperlink{S:init_sp_constraint_fct}{\code{init\_sp\_constraint()}},
    see \secref{S:init_sp_constraint_fct}.
    %% 
  \end{descr}
\end{function}
%%
\fdx{release_sp_constraint()@{\code{release\_sp\_constraint()}}|)}
\idx{linear solvers!release_sp_constraint()@{\code{release\_sp\_constraint()}}|)}
%%

\begin{function}{oem\_sp\_schur\_solve()}
\label{S:oem_sp_schur_solve_fct}

%%
\fdx{oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}|(}
\idx{linear solvers!oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}|(}
%%
\item[Prototype] ~\hfill
\bv\begin{lstlisting}
int oem_sp_schur_solve(OEM_SOLVER sp_solver,
                       REAL sp_tol, int sp_max_iter, int sp_info,
                       OEM_MV_FCT principal_inverse,
                       OEM_DATA *principal_data,
                       const DOF_REAL_VEC_D *f,
                       DOF_REAL_VEC_D *u,
                       SP_CONSTRAINT *constraint,
                       const DOF_REAL_VEC *g,
                       DOF_REAL_VEC *p,
                       ...);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
iterations =
  oem_sp_schur_solve(sp_solver,
                     sp_tol, sp_max_iter, sp_info,
                     A_inverse, A_data, f, u,
                     constraint, g, p,
                     ...);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Solve a saddle-point problem with possibly multiple, decoupled
  linear constraints by inverting the associated Schur's complement
  operator by means of an iterative method. Currently, only a
  CG-method is implemented, so the principal operator \code{A} has to
  be symmetric and positive (semi-) definite.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{oem_sp_solve:sp_solver}{sp\_solver} The solver used for
    the \emph{outer} iteration.  Currently, only a CG-method for a
    symmetric and positive (semi-) definite Schur's complement
    operator is implemented, so \code{sp\_solver} must equal the
    symbol \code{CG}.
    %% 
  \hyperitem{oem_sp_solve:sp_tol}{sp\_tol} The tolerance for the
    \emph{outer} CG-loop.
    %% 
  \hyperitem{oem_sp_solve:sp_max_iter}{sp\_max\_iter} The maximum
    number of iterations for the outer CG-loop.
    %% 
  \hyperitem{oem_sp_solve:sp_info}{sp\_info} A ``verbosity-level''
    controlling the amount of information printed to the terminal the
    application is running from.
    %% 
  \hyperitem{oem_sp_solve:A_inverse}{A\_inverse()} Pointer to a
    solver-function, for instance as returned by
    \hyperref[S:get_oem_solver_fct]{\code{get\_oem\_solver()}}.
    %% 
  \hyperitem{oem_sp_solve:A_data}{A\_data} A pointer to a data
    structure needed by \code{A\_inverse()}, the pointer is passed as
    first argument to \code{A\_inverse()}. See also
    \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solver()}} in
    \secref{S:init_oem_solve_fct}.
    %% 
  \hyperitem{oem_sp_solve:f}{f} The load-vector for the principal equation.
    %% 
  \hyperitem{oem_sp_solve:u}{u} Storage for the principal unknown
    (solution), and initial guess for the CG-method.
    %% 
  \hyperitem{oem_sp_solve:constraint}{constraint} A
    \hyperref[S:SP_CONSTRAINT_struct]{\code{SP\_CONSTRAINT}}
    structure, for instance as generated by a call to
    \hyperref[S:init_sp_constraint_fct]{\code{init\_sp\_constraint()}},
    see \secref{S:init_sp_constraint_fct}, see also
    \hyperref[S:release_sp_constraint_fct]{\code{release\_sp\_constraint()}},
    \secref{S:release_sp_constraint_fct}.
    %% 
  \hyperitem{oem_sp_solve:g}{g} The load vector for the possibly
    inhomogeneous linear constraint described by the parameter
    \code{constraint}. Note that only \emph{scalar} constraints are
    supported by this function, consequently \code{g} is a scalar
    \code{DOF\_REAL\_VEC}.
    %% 
  \hyperitem{oem_sp_solve:p}{p} Storage for the Lagrangian multiplier
    associated with \code{constraint} and initial guess for the
    Lagrangian multiplier in the outer CG-loop.
    %% 
  \hyperitem{oem_sp_solve:va_args}{\dots} More constraints may be
    added after the parameter \code{p}, each as a triple 
\bv\begin{lstlisting}
..., constraint_data, load_vector, lagrangian_multiplier, ...
\end{lstlisting}\ev
    All constraints must be decoupled from each other. After the final
    constraint a \nil-pointer must be passed to
    \code{oem\_sp\_schur\_solve()}, if only a single constraint is
    needed, then the first argument after the parameter \code{p} must
    already be a \nil-pointer.
    %% 
  \end{descr}
\item[Return Value] ~\hfill

  The number of iterations of the outer CG-loop for the inversion of
  Schur's complement.
\item[Examples] ~\hfill

  The single-constraint
  \hyperref[S:oem_sp_solve_fct]{\code{oem\_sp\_solve()}} functions are
  implemented on top of \code{oem\_sp\_schur\_solve()}. The interested
  reader may want to have a look at
  \code{alberta-VERSION/alberta/src/Common/oem\_sp\_solve.c}. See also
  \exampleref{example:oem_sp_schur_solve} below.

\end{function}
%%
\fdx{oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}|)}
\idx{linear solvers!oem_sp_schur_solve()@{\code{oem\_sp\_schur\_solve()}}|)}
%%

\begin{example}
  \label{example:oem_sp_schur_solve}
  A brief demonstration of how
  \hyperref[S:oem_sp_schur_solve_fct]{\code{oem\_sp\_schur\_solve()}}
  could be used in the single constraint case is given below. The
  reader is referred to \secref{S:precon} below for the documentation
  of the functions related to preconditioning.
  %%
\bv\begin{lstlisting}
... /* other stuff */

A_prec = init_precon_from_type(A, NULL /* bound */, sub_info, &A_prec_type);
A_oem  = init_oem_solve(A, NULL, tol, A_prec, -1, A_miter, sub_info);

Yproj_prec = init_precon_from_type(Yproj, NULL /* bound */, sub_info,
				   Yproj_prec_type);
Yprec_prec = init_precon_from_type(Yprec, NULL /* bound */, sub_info,
				   Yprec_prec_type);
SP_CONSTRAINT *div_constraint =
  init_sp_constraint(B, Bt, NULL, tol / 100.0, MAX(0, info - 3),
                     Yproj, Yproj_solver, Yproj_miter, Yproj_prec,
                     Yprec, Yprec_solver, Yprec_miter, Yprec_prec,
                     nu, mu);

oem_sp_schur_solve(solver, tol, miter, info,
                   get_oem_solver(A_solver), A_oem,
                   f_h, u_h,
                   div_constraint,
                   g_h, p_h,
                   NULL);

release_sp_constraint(div_constraint);
release_oem_solve(A_oem);

... /* other stuff */
\end{lstlisting}\ev
\end{example}

\begin{function}{sp\_dirichlet\_bound\_[dow\_scl|ds]()}
\label{S:sp_dirichlet_bound_fct}

\fdx{sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}|(}
\idx{linear solvers!sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}|(}
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
REAL sp_dirichlet_bound_dow_scl(MatrixTranspose transpose,
                                const DOF_MATRIX *Bt,
                                const DOF_SCHAR_VEC *bound,
                                const DOF_REAL_VEC_D *u,
                                DOF_REAL_VEC *g);
REAL sp_dirichlet_bound_ds(MatrixTranspose transpose,
                           const DOF_MATRIX *Bt,
                           const DOF_SCHAR_VEC *bound,
                           const DOF_REAL_D_VEC *u,
                           DOF_REAL_VEC *g);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
flux_excess = sp_dirichlet_bound_[dow_scl|ds](
  transpose, Bt, bound, u, g);
\end{lstlisting}\ev
\item[Description] ~\hfill

  If a flow field $u$ is subject to a divergence constraint and has to
  satisfy Dirichlet boundary values $h$ on the entire boundary of a
  domain $\Omega$, and if the test-space for the Lagrangian multiplier
  contains the function which is constant and equal to $1$ on the
  entire domain, then the Dirichlet boundary values have to satisfy
  the compatibility condition
  %%
  \begin{equation}
    \label{eq:divdirichletcompat}
    0 = \int_\Omega 1\,\div u = -\int_{\partial\Omega} u\cdot\nu =
    -\int_{\partial\Omega} h\cdot\nu.
  \end{equation}
  %%
  This compatibility conditions has also to be satisfied in the
  discrete setting, however, if one simply uses Lagrange-interpolation
  to implement Dirichlet boundary values, then the discrete Dirichlet
  boundary values in general violate this condition, and consequently
  the discrete saddle point problem does not have a solution in this
  case. One way to cope with this problem is to solve a slightly
  inhomogeneous saddle-point problem, where a load-vector for the
  Lagrangian multiplier compensates for the ``flux-excess'' of the
  interpolated Dirichlet boundary conditions (another way would be to
  modify the boundary values, of course).
  
  \code{sp\_dirichlet\_bound()} computes a load-vector for the
  Lagrangian multiplier by applying the $B^\ast$ operator to the
  boundary values. Of course, this makes only sense if the discrete
  boundary values asymptotically approximate the compatibility
  condition in the limit $h\rightarrow\infty$.

\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{sp_dirichlet_bound:transpose}{transpose} If equal to
    \hyperref[enum:MatrixTranspose]{\code{Tranpose}}, then the
    following parameter \code{Bt} actually is not $B^\ast$, but $B$.
    \code{sp\_dirichlet\_bound()} internally uses the transposed
    matrix for computing the load-vector \code{g}. If the parameter
    \code{Bt} is actually $B^\ast$, the \code{transpose} should be set to 
    \hyperref[enum:MatrixTranspose]{\code{NoTranpose}}.
    %% 
  \hyperitem{sp_dirichlet_bound:Bt}{Bt} A pointer to the DOF-matrix
    implementing the $B^\ast$ matrix from equation \eqref{eq:spschur},
    or the $B$-matrix if \code{transpose == Transpose}.
    %% 
  \hyperitem{sp_dirichlet_bound:bound}{bound} A
    \code{DOF\_SCHAR\_VEC}, if \code{bound->vec[dof] >= DIRICHLET},
    then the corresponding DOF belongs to a Dirichlet boundary.
    \code{bound} \emph{must not} be \nil,
    \code{sp\_dirichlet\_bound()} just works on the linear algebra
    level and does not loop over the mesh-elements. A suitable
    boundary-flag vector can be obtained by a call to the function
    \hyperref[S:dirichlet_bound]{\code{dirichlet\_bound()}}, see also
    \secref{S:dirichlet_bound}.

    If \code{sp\_dirichlet\_bound()} encounters DOFs with
    \code{bound->vec[dof] <= NEUMANN}, then it returns immediately to
    the caller and does not modify the load-vector \code{g}. See also
    \secref{S:boundary}.
    %% 
  \hyperitem{sp_dirichlet_bound:u}{u} The initial value for the
    principle unknown, \code{sp\_dirichlet\_bound()} expects that
    \code{u} already carries the Dirichlet boundary values.
    %% 
  \hyperitem{sp_dirichlet_bound:g}{g} Storage for the load-vector to
    compensate for the flux-excess. Note that the application has to
    initialize \code{g} prior to calling
    \code{sp\_dirichlet\_bound()}, which works also in the case of an
    inhomogeneous divergence constraint. In that case the
    compatibility condition has to be modified in the obvious manner.
    Anyhow, \code{sp\_dirichlet\_bound()} works additive.
    %% 
  \end{descr}
\item[Return Value] ~\hfill
  
  The total flux excess over the boundary segments of the domain, or
  \code{0.0} if for any DOF with \code{bound->vec[DOF] <= NEUMANN} was
  encountered.
\item[Examples] ~\hfill

  The interested read is referred to the program
\bv\begin{verbatim}
alberta-VERSION-demo/src/Common/stokes.c
\end{verbatim}\ev
\end{function}

\fdx{sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}|)}
\idx{linear solvers!sp_dirichlet_bound_ds()@{\code{sp\_dirichlet\_bound\_ds()}}|)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{OEM matrix-vector functions for DOF-matrices and -vectors}%
\label{S:matvec}

The general \code{oem\_...()} solvers all need pointers to
matrix-vector multiplication routines which do not accept arguments of
type \code{DOF\_REAL\_[D\_]VEC[\_D]} and a \code{DOF\_MATRIX} but work
directly on flat \code{REAL}-arrays. For the application to a scalar
or vector-valued linear system described by a \code{DOF\_MATRIX} (and
an optional \code{DOF\_SCHAR\_VEC} which can be used to honour
Dirichlet boundary conditions, see \secref{S:dirichlet_bound}), the
following routines are provided:
%%
\fdx{init_oem_mat_vec()@{\code{init\_oem\_mat\_vec()}}}%
\fdx{exit_oem_mat_vec()@{\code{exit\_oem\_mat\_vec()}}}%
%%
\bv\begin{lstlisting}
typedef int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);

OEM_MV_FCT oem_init_mat_vec(void **dataptrptr,
			    MatrixTranspose transpose, const DOF_MATRIX *A, 
			    const DOF_SCHAR_VEC *mask);
void exit_oem_mat_vec(void *dataptr)
\end{lstlisting}\ev

\begin{example}
  \label{E:init_oem_solve_impl}
  A short example demonstrating the function listed above. These are
  stripped-down versions of \code{init/release\_oem\_solve()}
  explained in \secref{S:ALBERTA_OEM_solvers}. The interested reader is
  referred to \code{alberta-VERSION/alberta/src/Common/oem\_solve.c}
  for the full source code.
  %% 
  \bv\begin{lstlisting}
OEM_DATA *simple_init_oem_solve(const DOF_MATRIX *A,
                                const DOF_SCHAR_VEC *mask,
                                REAL tol, int max_iter, int info)
{
  OEM_DATA     *oem;  
  const MatrixTranspose transpose = NoTranspose;

  oem = MEM_CALLOC(1, OEM_DATA); 
  oem->mat_vec   = init_oem_mat_vec(&oem->mat_vec_data, transpose, A, mask);
  oem->ws        = NULL; /* work-space,
                          * let the solvers handle this point for themselves.
                          */
  oem->tolerance = tol;
  oem->max_iter  = max_iter;
  oem->info      = MAX(0, info);

  return oem;
}

void simple_release_oem_solve(const OEM_DATA *_oem)
{
  OEM_DATA *oem = (OEM_DATA *)_oem;
  
  exit_oem_mat_vec(oem->mat_vec_data);
  MEM_FREE(oem, 1, OEM_DATA);
}
\end{lstlisting}\ev  
\end{example}

\begin{function}{init\_oem\_mat\_vec()}
\label{S:init_oem_mat_vec_fct}
%%
\fdx{init_oem_mat_vec()@{\code{init\_oem\_mat\_vec()}}}%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
mat_vec_fct =
  oem_init_mat_vec(&mv_data_ptr, transpose, A, mask);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Return a pointer to a function implementing the matrix-vector
  operation of the matrix \code{A} with a
  \code{DOF\_REAL[\_D]\_VEC[\_D]}. Of course, a matrix-vector product
  between a $\DOW\times\DOW$ block-matrix and a scalar
  \code{DOF\_REAL\_VEC} does not make sense. This function is fully
  aware of \ALBERTA's implementation of direct sums of finite element
  spaces, as described in \secref{S:chain_impl}.

\item[Parameters]~\hfill
  \begin{descr}
  \kitem{mv\_data\_ptr} After calling this function,
    \code{mv\_data\_ptr} will point to a control structure which must
    be passed as first argument to the function returned by
    \code{init\_oem\_mat\_vec()}. The application can call
    \code{exit\_oem\_mat\_vec()} to release the memory resources
    allocated by \code{init\_oem\_mat\_vec()}.
    %%
  \kitem{transpose} One of \code{Transpose} or \code{NoTranspose},
    indicating the matrix-vector operation should be performed with
    either the transposed or non-transposed matrix.
    %%
  \kitem{A} A pointer to a \code{DOF\_MATRIX}.
    %%
  \kitem{mask} A pointer to a \code{DOF\_SCHAR\_VEC} which can be used
    to exclude DOFs from the matrix-vector product. \code{mask} can be
    \nil. See \secref{S:dirichlet_bound} for further explanations.
  \end{descr}
\item[Return Value] ~\hfill

  A function pointer, pointing to the function actually implementing
  the matrix-vector operation. This function obeys the calling
  convention for the matrix-vector routines in the \code{OEM\_DATA}
  structure, see \secref{S:OEM} above.
\item[Examples]
  See \exampleref{E:init_oem_solve_impl}.
\end{function}

\begin{function}{exit\_oem\_mat\_vec()}
\label{S:exit_oem_mat_vec_fct}

\fdx{exit_oem_mat_vec()@{\code{exit\_oem\_mat\_vec()}}}%
\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
exit_oem_mat_vec(mv_data_ptr);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Release the resources previously allocated by a call to
  \code{init\_oem\_mat\_vec()}.

\item[Parameters]~\hfill
  \begin{descr}
  \kitem{mv\_data\_ptr} The data-pointer allocated by
    \code{init\_oem\_mat\_vec()}.
  \end{descr}

\item[Examples]
  See \exampleref{E:init_oem_solve_impl}.
\end{function}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Preconditioners}
\label{S:precon}

\begin{compatibility}
  The \code{get\_XXX\_precon()} functions no longer carry a
  \code{\dots\_[s|d|dow|}-suffix. This has been dropped, because the
  \hyperref[T:DOF_MATRIX]{\code{DOF\_MATRIX}} structure now carries
  its own block-type, and the finite element spaces described by the
  \hyperref[T:FE_SPACE]{\code{FE\_SPACE}} structure now know about the
  dimension of the range their elements are mapping to.

  See also \compatref{compat:precon} above for further remarks.
\end{compatibility}

The interface functions described in \secref{S:ALBERTA_OEM_solvers}
and \secref{S:ALBERTA_SPCG} which call the iterative solvers described
in \secref{S:OEM} and \secref{S:OEM_SPCG} all need a pointer to a
\code{PRECON} structure. Such a structure can either be initialized by
calls to one of the \code{get\_XXX\_precon()} functions described in
the Sections \ref{S:get_diag_precon_fct}-\ref{S:get_ILUk_precon_fct}:
%%
\fdx{get_diag_precon()@{\code{get\_diag\_precon()}}}%
\fdx{get_HB_precon()@{\code{get\_HB\_precon()}}}%
\fdx{get_BPX_precon()@{\code{get\_BPX\_precon()}}}%
\fdx{get_SSOR_precon()@{\code{get\_SSOR\_precon()}}}%
\fdx{get_ILUk_precon()@{\code{get\_ILUk\_precon()}}}%
%%
\bv\begin{lstlisting}
const PRECON *get_diag_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *bound);
const PRECON *get_HB_precon(const DOF_MATRIX *matrix,
                            const DOF_SCHAR_VEC *bound,
                            int info);
const PRECON *get_BPX_precon(const DOF_MATRIX *matrix,
                             const DOF_SCHAR_VEC *bound,
                             int info);
const PRECON *get_SSOR_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *bound,
                              REAL omega,
                              int n_iter);
const PRECON *get_ILUk_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *mask,
                              int ilu_level, int info);
\end{lstlisting}\ev
%%
These functions implement a diagonal and an SSOR preconditioner and
two hierarchical basis preconditioners (classical Yserentant
\cite{Yserentant:86} and Bramble-Pasciak-Xu \cite{BPX:90} types). The
\hyperref[S:get_ILUk_precon_fct]{$ILU(k)$ preconditioner} is the one
described in \cite{templates:94}.

Another possibility to get access to preconditioners are calls to the
following functions (see Sections
\ref{S:init_oem_precon_fct}-\ref{S:init_precon_from_type_fct}), which
also implement preconditioners for the block-matrices which arise in
the context of \hyperref[S:chain_impl]{direct sums of finite element
  spaces} (see \secref{S:chain_impl}):
%%
\bv
\begin{lstlisting}
const PRECON *init_oem_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *mask,
                              int info, OEM_PRECON precon,
                              ... /* ssor_omega, ssor_n_iter etc. */);
const PRECON *vinit_oem_precon(const DOF_MATRIX *A,
                               const DOF_SCHAR_VEC *mask,
                               int info, OEM_PRECON,
                               va_list ap);
const PRECON *init_precon_from_type(const DOF_MATRIX *A,
                                    const DOF_SCHAR_VEC *mask,
                                    int info,
                                    const PRECON_TYPE *prec_type);
\end{lstlisting}\ev

\begin{datatype}{PRECON}
\label{S:PRECON_struct}

\item[Description]~\hfill

%%
A preconditioner may need some initialization phase, which depends on
the matrix of the linear system, but is independent of the actual
application of the preconditioner to a vector. Thus, a preconditioner
is described by three functions for initialization, application, and a
final exit routine which may free memory which was allocated
during initialization, e.g. All three functions are collected in the
structure

\item[Definition]~\hfill

\ddx{PRECON@{\code{PRECON}}}%
\bv\begin{lstlisting}[label=T:PRECON]
typedef struct precon PRECON;
struct precon
{
  void    *precon_data;

  bool    (*init_precon)(void *precon_data);
  void    (*precon)(void *precon_data, int n, REAL *vec);
  void    (*exit_precon)(void *precon_data);
};
\end{lstlisting}\ev

\item[Components]~\hfill

\begin{descr}
\kitem{precon\_data} data for the preconditioner; always the first
argument to the functions \code{init\_precon()}, \code{precon()}, 
and \code{exit\_precon()}.
\kitem{init\_precon(precon\_data)} pointer to a function for initializing 
the preconditioning method; the return value is \code{false}
if initialization fails, otherwise \code{true}.
\kitem{precon(precon\_data)} pointer to a function for executing 
the preconditioning method; 

\code{precon} can be used as the entry \code{left\_precon} or
\code{right\_precon} in an \code{OEM\_DATA} structure together with
\code{precon\_data} as the corresponding pointer
\code{left\_precon\_data} respectively \code{right\_precon\_data}.

\kitem{exit\_precon(precon\_data)} frees all data used by the
preconditioning method.
\end{descr}
\end{datatype}

\begin{function}{get\_diag\_precon()}
\label{S:get_diag_precon_fct}
%%
\fdx{get_diag_precon()@{\code{get\_diag\_precon()}}|(}
\idx{preconditioner!get_diag_precon()@{\code{get\_diag\_precon()}}|(}

\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *get_diag_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *bound);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon_ptr = get_diag_precon(A, bound);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure
  describing a diagonal preconditioner for \code{A}. The application
  should call \code{precon\_ptr->exit\_precon(precon\_ptr)} to release
  the resources associated with \code{precon\_ptr} ones the
  preconditioner is no longer needed. But note that the solver
  interface-functions \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}
  and \hyperref[S:release_oem_solve_fct]{\code{release\_oem\_solve()}}
  call \code{exit\_precon()} on their own.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_diag_precon:A}{A} The matrix to compute the diagonal
    preconditioner for.
    %% 
  \hyperitem{get_diag_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %% 
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
\fdx{get_diag_precon()@{\code{get\_diag\_precon()}}|)}
\idx{preconditioner!get_diag_precon()@{\code{get\_diag\_precon()}}|)}
%%

\begin{function}{get\_HB\_precon()}
\label{S:get_HB_precon_fct}
%%
\fdx{get_HB_precon()@{\code{get\_HB\_precon()}}|(}
\idx{preconditioner!get_HB_precon()@{\code{get\_HB\_precon()}}|(}
\idx{preconditioner!hierarchical basis}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *get_HB_precon(const DOF_MATRIX *matrix,
                            const DOF_SCHAR_VEC *bound,
                            int info);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon_ptr = get_HB_precon(A, bound, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure
  describing a hierarchical preconditioner, as described in
  \cite{Yserentant:86}. The application should call
  \code{precon\_ptr->exit\_precon(precon\_ptr)} to release the
  resources associated with \code{precon\_ptr} once the preconditioner
  is no longer needed. But note that the solver interface-functions
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} and
  \hyperref[S:release_oem_solve_fct]{\code{release\_oem\_solve()}}
  call \code{exit\_precon()} on their own.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_HB_precon:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{get_HB_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %% 
  \hyperitem{get_HB_precon:info}{info} An integer controlling the
    amount of information printed to the terminal the application is
    running in (larger values mean more ``noise'').
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
\fdx{get_HB_precon()@{\code{get\_HB\_precon()}}|)}
\idx{preconditioner!get_HB_precon()@{\code{get\_HB\_precon()}}|)}
%%

\begin{function}{get\_BPX\_precon()}
\label{S:get_BPX_precon_fct}
%%
\fdx{get_BPX_precon()@{\code{get\_BPX\_precon()}}|(}
\idx{preconditioner!get_BPX_precon()@{\code{get\_BPX\_precon()}}|(}
\idx{preconditioner!BPX}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *get_BPX_precon(const DOF_MATRIX *matrix,
                             const DOF_SCHAR_VEC *bound,
                             int info);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon_ptr = get_BPX_precon(A, bound, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure
  describing the BPX-preconditioner, as described in \cite{BPX:90}.
  The application should call
  \code{precon\_ptr->exit\_precon(precon\_ptr)} to release the
  resources associated with \code{precon\_ptr} once the preconditioner
  is no longer needed. But note that the solver interface-functions
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} and
  \hyperref[S:release_oem_solve_fct]{\code{release\_oem\_solve()}}
  call \code{exit\_precon()} on their own.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_BPX_precon:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{get_BPX_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %% 
  \hyperitem{get_BPX_precon:info}{info} An integer controlling the
    amount of information printed to the terminal the application is
    running in (larger values mean more ``noise'').
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
\fdx{get_BPX_precon()@{\code{get\_BPX\_precon()}}|)}
\idx{preconditioner!get_BPX_precon()@{\code{get\_BPX\_precon()}}|)}
%%

\begin{function}{get\_SSOR\_precon()}
\label{S:get_SSOR_precon_fct}
%%
\fdx{get_SSOR_precon()@{\code{get\_SSOR\_precon()}}|(}
\idx{preconditioner!get_SSOR_precon()@{\code{get\_SSOR\_precon()}}|(}
\idx{preconditioner!SSOR}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *get_SSOR_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *bound,
                              REAL omega,
                              int n_iter);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon_ptr = get_SSOR_precon(A, bound, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure
  describing an SSOR-preconditioner. The application should call
  \code{precon\_ptr->exit\_precon(precon\_ptr)} to release the
  resources associated with \code{precon\_ptr} once the preconditioner
  is no longer needed. But note that the solver interface-functions
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} and
  \hyperref[S:release_oem_solve_fct]{\code{release\_oem\_solve()}}
  call \code{exit\_precon()} on their own.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_SSOR_precon:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{get_SSOR_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %% 
  \hyperitem{get_SSOR_precon:omega}{omega} The relaxation parameter.
    %% 
  \hyperitem{get_SSOR_precon:n_iter}{n\_iter} The number of
    SSOR-iterations to perform.
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
\fdx{get_SSOR_precon()@{\code{get\_SSOR\_precon()}}|)}
\idx{preconditioner!get_SSOR_precon()@{\code{get\_SSOR\_precon()}}|)}
%%

\begin{function}{get\_ILUk\_precon()}
\label{S:get_ILUk_precon_fct}
%%
\fdx{get_ILUk_precon()@{\code{get\_ILUk\_precon()}}|(}
\idx{preconditioner!get_ILUk_precon()@{\code{get\_ILUk\_precon()}}|(}
\idx{preconditioner!ILUk}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *get_ILUk_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *mask,
                              int ilu_level, int info);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon_ptr = get_ILUk_precon(A, bound, info);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure
  describing an $ILU(k)$-preconditioner as described
  \cite{templates:94}. This preconditioner uses a combinatorical,
  ``level''-based strategy to control the amount of fill-in generated
  by the incomplete $LU$-factorization. The preconditioner can benefit
  from re-ordering the DOFs in a way that the amount of fill-in
  generated by a complete $LU$-factorization would be minimized.
  Currently, \ALBERTA searches for a library \code{libgpskca} and
  expects that this library contains the functions of the
  \code{GPSKCA} package from
  \href{http://www.netlib.org}{www.netlib.org}, \cite{GPSKCA:82}.

  Note the level-based fill-in control has the disadvantage that the
  generated preconditioner may not even be positive definite, even if
  \code{A} is spd. On the other hand, $ILU(k)$ may still be spd even
  if \code{A} is not.

  The application should call
  \code{precon\_ptr->exit\_precon(precon\_ptr)} to release the
  resources associated with \code{precon\_ptr} once the preconditioner
  is no longer needed. But note that the solver interface-functions
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} and
  \hyperref[S:release_oem_solve_fct]{\code{release\_oem\_solve()}}
  call \code{exit\_precon()} on their own.
\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{get_ILUk_precon:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{get_ILUk_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %% 
  \hyperitem{get_ILUk_precon:level}{level} The control parameter for
    the amount of fill-in, see \cite{templates:94}.
    %% 
  \hyperitem{get_ILUk_precon:info}{info} An integer controlling the
    amount of information printed to the terminal the application is
    running in (larger values mean more ``noise'').
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
\fdx{get_ILUk_precon()@{\code{get\_ILUk\_precon()}}|)}
\idx{preconditioner!get_ILUk_precon()@{\code{get\_ILUk\_precon()}}|)}
%%

\begin{datatype}{OEM\_PRECON}
\label{S:OEM_PRECON_enum}

\item[Definition]~\hfill

\ddx{OEM_PRECON@{\code{OEM\_PRECON}}}
\idx{preconditioner!OEM\_PRECON@{\code{OEM\_PRECON}}}
\bv\begin{lstlisting}
typedef enum {
  PreconEnd  = -1,
  PreconRepeat = PreconEnd,
  NoPrecon   = 0,
  DiagPrecon = 1,
  HBPrecon   = 2,
  BPXPrecon  = 3,
  SSORPrecon = 4,
  __SSORPrecon = 5,
  ILUkPrecon = 6,
  BlkDiagPrecon = 512,
  BlkSSORPrecon = 513,
} OEM_PRECON;
\end{lstlisting}\ev

\item[Symbols]~\hfill
\begin{descr}
\hyperitem{OEM_PRECON:PreconEnd}{PreconEnd}
\hyperitem{OEM_PRECON:PreconRepeat}{PreconRepeat} Terminate the
  variable argument list of
  \hyperref[S:init_oem_precon_fct]{\code{init\_oem\_precon()}}, see
  \secref{S:init_oem_precon_fct} in the context of block-matrix
  preconditioners for block-matrices having their origin in direct-sum
  structure of the underlying finite element spaces (see
  \secref{S:chain_impl}).
  %%
\hyperitem{OEM_PRECON:NoPrecon}{NoPrecon}
\hyperitem{OEM_PRECON:DiagPrecon}{DiagPrecon}
\hyperitem{OEM_PRECON:HBPrecon}{HBPrecon}
\hyperitem{OEM_PRECON:BPXPrecon}{BPXPrecon} Self-explanatory, select
  the respective preconditioner.
  %%
\hyperitem{OEM_PRECON:SSORPrecon}{SSORPrecon} Select an
  SSOR-preconditioner with
  \hyperlink{get_SSOR_precon:omega}{\code{omega == 1.0}} and
  \hyperlink{get_SSOR_precon:n_iter}{\code{n\_iter == 2}}.
  %%
\hyperitem{OEM_PRECON:__SSORPrecon}{\_\_SSORPrecon}
  Select an  SSOR-preconditioner with control over
  \hyperlink{get_SSOR_precon:omega}{\code{omega}} and
  \hyperlink{get_SSOR_precon:n_iter}{\code{n\_iter}}.
  %%
\hyperitem{OEM_PRECON:ILUkPrecon}{ILUkPrecon}
  Self explanatory.
  %%
\hyperitem{OEM_PRECON:BlkDiagPrecon}{BlkDiagPrecon} Select a
  preconditioner which acts on a block-matrix structure induced by a
  finite element space with is composed of several components as a
  direct sum (see \secref{S:chain_impl}).
  %%
\hyperitem{OEM_PRECON:BlkSSORPrecon}{BlkSSORPrecon}
  Currently not supported.
\end{descr}
\end{datatype}

\begin{function}{init\_oem\_precon()}
\label{S:init_oem_precon_fct}
\label{S:vinit_oem_precon_fct}
%%
\fdx{init_oem_precon()@{\code{init\_oem\_precon()}}|(}
\idx{preconditioner!init_oem_precon()@{\code{init\_oem\_precon()}}|(}
\fdx{vinit_oem_precon()@{\code{vinit\_oem\_precon()}}|(}
\idx{preconditioner!vinit_oem_precon()@{\code{vinit\_oem\_precon()}}|(}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *init_oem_precon(const DOF_MATRIX *A,
                              const DOF_SCHAR_VEC *bound,
                              int info, OEM_PRECON precon_enum,
                              ... /* ssor_omega, ssor_n_iter etc. */);
const PRECON *vinit_oem_precon(const DOF_MATRIX *A,
                               const DOF_SCHAR_VEC *bound,
                               int info, OEM_PRECON precon_enum,
                               va_list ap);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon = init_oem_precon(A, bound, info, precon_enum, ...);
precon = vinit_oem_precon(A, bound, info, precon_enum, ap);
\end{lstlisting}\ev
\item[Description] ~\hfill

  These two function initialize a
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure, based on the
  value of a descriptive enumeration symbol. The returned structure
  can then be passed to
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} or
  \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solve()}}, as
  described in \secref{S:ALBERTA_OEM_solvers}.  In contrast to the
  \code{get\_XXX\_precon()} functions described above these two
  functions support matrices with the block-matrix structure implied
  by using
  \hyperref[S:chain_impl]{direct sums of finite element spaces}, see
  \secref{S:chain_impl} for further explanations.

  For the difference between the \code{\dots} ``argument'' and the
  \code{ap} argument the reader is referred to any text-book dealing
  with the \code{C}-programming language.

\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{init_oem_precon:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{init_oem_precon:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
  \hyperitem{init_oem_precon:info}{info} An integer controlling the
    amount of information printed to the terminal the application is
    running in (larger values mean more ``noise'').
    %%
  \hyperitem{init_oem_precon:precon_enum}{precon\_enum} An enumeration
    value as defined by
    \hyperlink{S:OEM_PRECON_enum}{\code{OEM\_PRECON}}, see
    \secref{S:OEM_PRECON_enum}, selecting the respective
    preconditioner to use.
   %% 
  \hyperitem{init_oem_precon:va_list}{\dots, ap} A variable-length
    argument list, providing additional parameters needed by some of
    the preconditioners, as explained below:
    \begin{descr}
    \kitem{\_\_SSORPrecon} The two arguments following
      \code{precon\_enum} must specify the relaxation parameter
      \hyperlink{get_SSOR_precon:omega}{\code{omega}} and the number
      of iterations
      \hyperlink{get_SSOR_precon:n_iter}{\code{n\_iter}}
      to perform.
      %% 
    \kitem{ILUkPrecon} The argument following \code{precon\_enum} must
      specify the control-parameter
      \hyperlink{get_ILUk_precon:level}{$k$} controlling the amount of
      fill-in.
      %%
    \kitem{BlkDiagPrecon} The parameters following \code{precon\_enum}
      must specify the type and parameters for the preconditioners for
      the diagonal blocks. It is the responsibility of the calling
      application to ensure that enough preconditioners are defined.
      An example to generate a block-diagonal preconditioner for a
      $3\times 3$ block-matrix (e.g. in the context of a
      ``Crouzeix-Raviart'' discretization for the Stokes-problem in
      3d) would be
      %%
      \bv\begin{lstlisting}
precon = init_oem_precon(A, NULL, 3 /* info */, BlkDiagPrecon,
                         __SSORPrecon, 1.5, 2, DiagPrecon, DiagPrecon);
\end{lstlisting}\ev
      %%
      The symbol
      \hyperlink{OEM_PRECON:PreconRepeat}{\code{PreconRepeat}} has a
      special meaning: it indicates that the last specified
      preconditioner should also be used for all other blocks. In the
      $3\times 3$ example given above, the following code-fragment
      would select diagonal preconditioning for all blocks>
      %%
      \bv\begin{lstlisting}
precon = init_oem_precon(A, NULL, 3 /* info */, BlkDiagPrecon,
                         DiagPrecon, PreconRepeat);
\end{lstlisting}\ev
      %%
    \end{descr}
  \end{descr}
\item[Return Value] ~\hfill

  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\end{function}
%%
%%
\fdx{init_oem_precon()@{\code{init\_oem\_precon()}}|)}
\idx{preconditioner!init_oem_precon()@{\code{init\_oem\_precon()}}|)}
\fdx{vinit_oem_precon()@{\code{vinit\_oem\_precon()}}|)}
\idx{preconditioner!vinit_oem_precon()@{\code{vinit\_oem\_precon()}}|)}
%%

\begin{datatype}{PRECON\_TYPE}
\label{S:PRECON_TYPE_struct}

\item[Description]~\hfill

A data structure which can be use to define more complex
preconditioners. The purpose of this structure is to avoid defining
functions with an endless number of arguments. This
``parameter-transport-structure'' can be passed to
\hyperref[S:init_precon_from_type_fct]{\code{init\_precon\_from\_type()}},
instead of calling
\hyperref[S:init_oem_precon_fct]{\code{init\_oem\_precon()}}. The
actual definition looks somewhat complicated and maybe ugly, but using
this structure is more or less straight-forward, have a look at
\exampleref{example:PRECON_TYPE} below.

\item[Definition]~\hfill

\bv\begin{lstlisting}
#define N_BLOCK_PRECON_MAX 10

struct __precon_type {
  OEM_PRECON type;
  union {
    struct {
      REAL omega;
      int n_iter;  
    } __SSOR;
    struct {
      int level;
    } ILUk;
  } param;
};    

typedef struct precon_type
{
  OEM_PRECON type;
  union {
    struct {
      REAL omega;
      int n_iter;  
    } __SSOR;
    struct {
      int level;
    } ILUk;
    struct {
      struct __precon_type precon[N_BLOCK_PRECON_MAX];
    } BlkDiag;
    struct {
      struct __precon_type precon[N_BLOCK_PRECON_MAX];
      REAL omega;
      int  n_iter;
    } BlkSSOR;
  } param;
} PRECON_TYPE;
\end{lstlisting}\ev

\item[Components]~\hfill
\begin{descr}
\hyperitem{PRECON_TYPE:type}{type} One of the symbolic constants
  defined by the \hyperref[S:OEM_PRECON_enum]{\code{OEM\_PRECON}}
  enumeration type. See \secref{S:OEM_PRECON_enum}.
  %% 
\hyperitem{PRECON_TYPE:param}{param} If the preconditioner defined by
  \code{type} needs additional parameters, then the corresponding
  section in the \code{param} component has to be filled. The names of
  the structure components correspond to the parameters for the
  \code{get\_XXX\_precon()} functions described above, currently, only
  \hyperref[S:get_SSOR_precon_fct]{\code{\_\_SSORPrecon}},
  \hyperref[S:get_ILUk_precon_fct]{\code{ILUkPrecon}} and, of course,
  \code{BlkDiagPrecon} need additional parameters. For the latter, the
  \code{param} component contains an array of
  \code{N\_BLOCK\_PRECON\_MAX} many \code{struct \_\_precon\_type}
  sub-structures for storing additional parameters possibly needed by
  the sub-preconditioners.
\end{descr}
\end{datatype}

\begin{example}
  \label{example:PRECON_TYPE}

  Two short examples demonstrating the use of the
  \hyperref[S:PRECON_TYPE_struct]{\code{PRECON\_TYPE}} structure
  defined above.
  \begin{itemize}
  \item Defining an SSOR preconditioner with control over the
    relaxation parameter and the number of iterations:
    %% 
    \bv\begin{lstlisting}
PRECON_TYPE prec;
prec.type = __SSORPrecon;
prec.param.__SSOR.omega  = 1.5;
prec.param.__SSOR.n_iter = 2;
\end{lstlisting}\ev
    
  \item Defining a preconditioner for a block-matrix resulting from
    using a
    \hyperref[S:chain_impl]{direct sum} of finite element spaces
    %%
    \bv\begin{lstlisting}
PRECON_TYPE prec;
prec.type = BlkDiagPrecon;
prec.param.BlkDiag.precon[0].type = __SSOR;
prec.param.BlkDiag.param.precon[0].__SSOR.omega  = 1.0;
prec.param.BlkDiag.param.precon[0].__SSOR.n_iter = 1;
for (i = 1; i < 3; i++) {
  prec.param.BlkDiag.precon[i].type = DiagPrecon;
}
\end{lstlisting}\ev
  \end{itemize}
\end{example}

\begin{function}{init\_precon\_from\_type()}
\label{S:init_precon_from_type_fct}
%%
\fdx{init_precon_from_type()@{\code{init\_precon\_from\_type()}}|(}
\idx{preconditioner!init_precon_from_type()@{\code{init\_\_precon\_from\_type()}}|(}
%%
\item[Prototype] ~\hfill
%%
\bv\begin{lstlisting}
const PRECON *init_precon_from_type(const DOF_MATRIX *A,
                                    const DOF_SCHAR_VEC *bound,
                                    int info,
                                    const PRECON_TYPE *prec_type);
\end{lstlisting}\ev

\item[Synopsis] ~\hfill

\bv\begin{lstlisting}[basicstyle=\normalsize]
precon = init_precon_from_type(A, bound, info, prec_type);
\end{lstlisting}\ev
\item[Description] ~\hfill

  Initialize a \hyperref[S:PRECON_struct]{\code{PRECON}} structure,
  based on contents of the
  \hyperref[S:PRECON_TYPE_struct]{\code{prec\_type}} parameter.  The
  returned structure can then be passed to
  \hyperref[S:oem_solve_fct]{\code{oem\_solve()}} or
  \hyperref[S:init_oem_solve_fct]{\code{init\_oem\_solve()}}, as
  described in \secref{S:ALBERTA_OEM_solvers}.  In contrast to the
  \code{get\_XXX\_precon()} functions described above these two
  functions support matrices with the block-matrix structure implied
  by using
  \hyperref[S:chain_impl]{direct sums of finite element spaces}, see
  \secref{S:chain_impl} for further explanations.

\item[Parameters]~\hfill
  \begin{descr}
  \hyperitem{init_precon_from_type:A}{A} The matrix to compute the
    preconditioner for.
    %% 
  \hyperitem{init_precon_from_type:bound}{bound} A flag-vector, masking out
    specific DOFs, compare the explanations for the
    \hyperlink{oem_solve:mask}{\code{mask}} parameter to
    \hyperref[S:oem_solve_fct]{\code{oem\_solve()}}, see
    \secref{S:oem_solve_fct}. \code{bound} may be \nil.
    %%
  \hyperitem{init_precon_from_type:info}{info} An integer controlling the
    amount of information printed to the terminal the application is
    running in (larger values mean more ``noise'').
    %%
  \hyperitem{init_precon_from_type:prec_type}{prec\_type} A pointer to
    a structure of type
    \hyperref[S:PRECON_TYPE_struct]{\code{PRECON\_TYPE}}, as described
    in \secref{S:PRECON_TYPE_struct} above, describing the
    preconditioner to generate.
  \end{descr}
\item[Return Value] ~\hfill
  
  A pointer to an initialized
  \hyperref[S:PRECON_struct]{\code{PRECON}} structure implementing the
  preconditioner, see \secref{S:PRECON_struct}.
\item[Examples] ~\hfill

  The function
  \hyperref[S:init_oem_precon_fct]{\code{init\_oem\_precon()}} (see
  \secref{S:init_oem_precon_fct}) is implemented on top of
  \code{init\_precon\_from\_type()}. The interested reader is referred
  to the source code in
  \code{alberta-VERSION/alberta/src/Common/oem\_solver.c}
\end{function}
%%
\fdx{init_oem_precon()@{\code{init\_oem\_precon()}}|)}
\idx{preconditioner!init_oem_precon()@{\code{init\_oem\_precon()}}|)}
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Multigrid solvers}\label{S:multigrid_solver}

A abstract framework for multigrid solvers is available. The main 
data structure for the multigrid solver \code{MG()} is
\ddx{MULTI_GRID_INFO@\code{MULTI\_GRID\_INFO}}%
\bv\begin{lstlisting}
typedef struct multi_grid_info MULTI_GRID_INFO;
struct multi_grid_info
{
  REAL             tolerance;                 /* tol. for resid         */
  REAL             exact_tolerance;           /* tol. for exact_solver  */

  int              cycle;                     /* 1=V-cycle, 2=W-cycle   */
  int              n_pre_smooth, n_in_smooth; /* no of smoothing loops  */
  int              n_post_smooth;             /* no of smoothing loops  */
  int              mg_levels;                 /* current no. of levels  */
  int              exact_level;               /* level for exact_solver */
  int              max_iter;                  /* max. no of MG iter's   */
  int              info;

  int              (*init_multi_grid)(MULTI_GRID_INFO *mg_info);
  void             (*pre_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*in_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*post_smooth)(MULTI_GRID_INFO *mg_info, int level, int n);
  void             (*mg_restrict)(MULTI_GRID_INFO *mg_info, int level);
  void             (*mg_prolongate)(MULTI_GRID_INFO *mg_info, int level);
  void             (*exact_solver)(MULTI_GRID_INFO *mg_info, int level);
  REAL             (*mg_resid)(MULTI_GRID_INFO *mg_info, int level);
  void             (*exit_multi_grid)(MULTI_GRID_INFO *mg_info);

  void             *data;                         /* application dep. data */
};
\end{lstlisting}\ev
The entries yield following information:
\begin{descr}
\kitem{tolerance} tolerance for norm of residual.
\kitem{exact\_tolerance} tolerance for ``exact solver'' on coarsest level.
\kitem{cycle} selection of multigrid cycle type: $1=$V-cycle, $2=$W-cycle, ....
\kitem{n\_pre\_smooth} number of smoothing steps on each level before (first)
      coarse level correction.
\kitem{n\_in\_smooth} number of smoothing steps on each level between
      coarse level corrections (for \code{cycle} $\geq 2$).
\kitem{n\_post\_smooth} number of smoothing steps on each level after (last)
      coarse level correction.
\kitem{mg\_levels} number of levels.
\kitem{exact\_level} selection of grid level where the ``exact'' solver is used
      (and no further coarse grid correction), usually \code{exact\_level=0}.
\kitem{max\_iter} maximal number of multigrid iterations.
\kitem{info} level of information produced by the multigrid method.
\kitem{init\_multi\_grid} pointer to a function for initializing the
  multigrid method; may be \nil;

  if not \nil, \code{init\_multi\_grid(mg\_info)} initializes data
  needed by the multigrid method, returns \code{true} if an error occurs.

\kitem{pre\_smooth} pointer to a function for performing the smoothing
  step before coarse grid corrections;

 \code{pre\_smooth(mg\_info, level, n)} performs \code{n} smoothing iterations
      on grid \code{level}.

\kitem{in\_smooth} pointer to a function for performing the smoothing
  step between coarse grid corrections;

 \code{in\_smooth(mg\_info, level, n)} performs \code{n} smoothing iterations
      on grid \code{level}.

\kitem{post\_smooth} pointer to a function for performing the smoothing
  step after coarse grid corrections;

\code{{post\_smooth(mg\_info, level, n)}} performs \code{n} smoothing 
  iterations on grid \code{level}.

\kitem{mg\_restrict} pointer to a function for computing and restricting the
  residual to a coarser level;

 \code{mg\_restrict(mg\_info, level)} computes and restricts the
   residual from grid \code{level} to next coarser grid (\code{level-1}).

\kitem{mg\_prolongate} pointer to a function for prolongating and adding
  coarse grid corrections to the fine grid solution;

  \code{mg\_prolongate(mg\_info, level)} prolongates and adds the coarse grid 
      (\code{level-1}) correction to the fine grid solution on grid 
      \code{level}.

\kitem{exact\_solver} pointer to a function for the ``exact'' solver;

\code{exact\_solver(mg\_info, level)} computes the ``exact'' solution of 
  the problem on grid \code{level} with tolerance 
  \code{mg\_info->exact\_tolerance}.

\kitem{mg\_resid} pointer to a function for computing the norm of
 the actual residual;

 \code{mg\_resid(mg\_info, level)} returns the norm of residual on grid
      \code{level}.

\kitem{exit\_multi\_grid} a pointer to a cleanup routine, may be \nil;

  if not \nil \code{exit\_multi\_grid(mg\_info)} is called after
  termination of the multigrid method for freeing used data.

\kitem{data} pointer to application dependent data, holding information on
      or about different grid levels, e.g.
\end{descr}
%
The abstract multigrid solver is implemented in the routine
\fdx{MG()@\code{MG()}}
\bv\begin{lstlisting}
int MG(MULTI_GRID_INFO *)
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{MG(mg\_info)} 
based upon information given in the data structure \code{mg\_info}, the
subroutine \verb|MG()| iterates until the prescribed tolerance is met or
the prescribed number of multigrid cycles is performed.
\end{descr}
%
Main parts of the \code{MG()} routine are:
\bv\begin{lstlisting}
{
  int iter;
  REAL resid;

  if (mg_info->init_multi_grid) 
    if (mg_info->init_multi_grid(mg_info))
      return(-1);

  resid = mg_info->resid(mg_info, mg_info->mg_levels-1);
  if (resid <= mg_info->tolerance)
    return(0);

  for (iter = 0; iter < mg_info->max_iter; iter++)
  {
    recursive_MG_iteration(mg_info, mg_info->mg_levels-1);
    resid = mg_info->resid(mg_info, mg_info->mg_levels-1);
    if (resid <= mg_info->tolerance)
      break;
  }
  if (mg_info->exit_multi_grid)
    mg_info->exit_multi_grid(mg_info);

  return(iter+1);
}
\end{lstlisting}\ev
The subroutine \code{recursive\_MG\_iteration()} performs smoothing,
restriction of the residual and prolongation of the coarse grid correction:
\bv\begin{lstlisting}
static void recursive_MG_iteration(MULTI_GRID_INFO *mg_info, int level)
{
  int  cycle;

  if (level <= mg_info->exact_level) {
    mg_info->exact_solver(mg_info, level);
  }
  else {
    if (mg_info->pre_smooth)
      mg_info->pre_smooth(mg_info, level, mg_info->n_pre_smooth);

    for (cycle = 0; cycle < mg_info->cycle; cycle++) {
      if ((cycle > 0) && mg_info->in_smooth)
        mg_info->in_smooth(mg_info, level, mg_info->n_in_smooth);

      mg_info->mg_restrict(mg_info, level);
      recursive_MG_iteration(mg_info, level-1);
      mg_info->prolongate(mg_info, level);
    }

    if (mg_info->post_smooth)
      mg_info->post_smooth(mg_info, level, mg_info->n_post_smooth);
  }
}
\end{lstlisting}\ev

For multigrid solution of a scalar linear system
\[
  A u = f
\]
given by a \code{DOF\_MATRIX A} and a \code{DOF\_REAL\_VEC f},
the following subroutine is available:
\fdx{mg_s()@\code{mg\_s()}}
\bv\begin{lstlisting}
int mg_s(DOF_MATRIX *, DOF_REAL_VEC *, const DOF_REAL_VEC *, 
         const DOF_SCHAR_VEC *, REAL, int, int, char *);
\end{lstlisting}\ev
%
Description:
\begin{descr}
\kitem{mg\_s(matrix, u, f, bound, tol, max\_iter, info, prefix)} 
solves the linear system  for a scalar valued problem by a
multigrid method; the return value is the number of performed
iterations; 

\code{matrix} is a pointer to a DOF matrix storing the system matrix,
\code{u} is a pointer to a DOF vector for the solution,
holding an initial guess on input;
\code{f} is a pointer to a DOF vector storing the right hand side and
\code{bound} a pointer to a DOF vector with information about boundary 
DOFs; \code{bound} must not be \nil if Dirichlet DOFs are used;

\code{tol} is the tolerance for multigrid solver, \code{max\_iter}
the maximal number of multigrid iterations and
\code{info} gives the level of information for the solver;

\code{prefix} is a parameter key prefix for the initialization of additional
      data via \code{GET\_PARAMETER}, see Table \ref{T:mg_parms}, may be \nil;
\begin{table}
\begin{center}
 \begin{tabular}{|l|c|l|} \hline
   member & default & key\\ \hline\hline
   \code{mg\_info->cycle} & \code{1} &\code{prefix->cycle}\\
   \code{mg\_info->n\_pre\_smooth} & \code{1} &\code{prefix->n\_pre\_smooth}\\
   \code{mg\_info->n\_in\_smooth} & \code{1} & \code{prefix->n\_in\_smooth}\\
   \code{mg\_info->n\_post\_smooth} & \code{1} & \code{prefix->n\_post\_smooth}\\
   \code{mg\_info->exact\_level} & \code{0} & \code{prefix->exact\_level}\\
   \code{mg\_info->info} & \code{info} & \code{prefix->info}\\
  \hline
   \code{mg\_s\_info->smoother} & \code{1} & \code{prefix->smoother}\\
   \code{mg\_s\_info->smooth\_omega} & \code{1.0} 
   &\code{prefix->smooth\_omega} \\
   \code{mg\_s\_info->exact\_solver} & \code{1} &\code{prefix->exact\_solver}\\
   \code{mg\_s\_info->exact\_omega} & \code{1.0} &\code{prefix->exact\_omega}\\
\hline
 \end{tabular}
\end{center}
\caption{Parameters read by \code{mg\_s()} and \code{mg\_s\_init()}}
\label{T:mg_parms}
\end{table}
%
an SOR smoother (\code{mg\_s\_info->smoother=1}) and an SSOR smoother
(\code{smoother=2}) are available; under-- or over relaxation
parameter can be specified by \code{mg\_s\_info->smooth\_omega}. These
SOR/SSOR smoothers are used for \code{exact\_solver}, too.
\end{descr}

For applications, where several systems with the same matrix have to
be solved, computing time can be saved by doing all initializations
like setup of grid levels and restriction of matrices only once.
For such cases, three subroutines are available:
\fdx{mg_s_init()@\code{mg\_s\_init()}}%
\fdx{mg_s_solve()@\code{mg\_s\_solve()}}%
\fdx{mg_s_exit()@\code{mg\_s\_exit()}}%
\bv\begin{lstlisting}
MG_S_INFO *mg_s_init(DOF_MATRIX *, const DOF_SCHAR_VEC *, int, char *);
int mg_s_solve(MG_S_INFO *, DOF_REAL_VEC *, const DOF_REAL_VEC *, REAL, int);
void mg_s_exit(MG_S_INFO *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{mg\_s\_init(matrix, bound, info, prefix)}
function for initializing a multigrid method for solving
a scalar valued problem by \code{mg\_s\_solve()}; the return
value is a pointer to data used by \code{mg\_s\_solve()} and 
is the first argument to this function; the structure \code{MG\_S\_INFO} 
contains matrices and vectors for linear problems on all used grid levels.

\code{matrix} is a pointer to a DOF matrix storing the system matrix,
\code{bound} a pointer to a DOF vector with information about boundary 
DOFs; \code{bound} must not be \nil if Dirichlet DOFs are used;

\code{info} gives the level of information for \code{mg\_s\_solve()};
\code{prefix} is a parameter key prefix for the initialization of additional
data via \code{GET\_PARAMETER}, see Table \ref{T:mg_parms}, may be \nil.

\kitem{mg\_s\_solve(mg\_s\_info, u, f, tol, max\_iter)} 
solves the linear system  for a scalar valued problem by a
multigrid method; the routine has to be initialize by
\code{mg\_s\_init()} and the return value \code{mg\_s\_info} of
\code{mg\_s\_init()} is the first argument; the return value 
of \code{mg\_s\_solve()} is the number of performed iterations; 
  
\code{u} is a pointer to a DOF vector for the solution, holding an
initial guess on input; \code{f} is a pointer to a DOF vector storing
the right hand side; \code{tol} is the tolerance for multigrid solver,
\code{max\_iter} the maximal number of multigrid iterations;

the function may be called several times with different right
hand sides \code{f}.

\kitem{mg\_s\_exit(mg\_s\_info)} frees data needed for the
multigrid method and which is allocated by \code{mg\_s\_init()}.
\end{descr}

\begin{remark}
The multigrid solver is currently available only for Lagrange finite
elements of first order (\code{lagrange1}). An implementation for
higher order elements is future work.
\end{remark}%
\idx{linear solvers|)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Nonlinear solvers}%
\label{S:nls}%
\idx{nonlinear solvers|(}
For the solution of a nonlinear equation
\begin{equation}\label{E:nonlin}
u \in \R^N: \qquad F(u) = 0 \qquad\mbox{in } \R^N
\end{equation}
several Newton methods are provided. For testing the convergence
a (problem dependent) norm of either the correction $d_k$ in
the $k$th step, i.e.
\[
\|d_k\| = \|u_{k+1} - u_k\|,
\]
or the residual, i.e.
\[
\|F(u_{k+1})\|,
\] 
is used. 

The data structure (defined in \code{alberta\_util.h}) for passing information
about assembling and solving a linearized equation, tolerances,
etc. to the solvers is
\ddx{NLS_DATA@{\code{NLS\_DATA}}}
\idx{nonlinear solvers!NLS_DATA@{\code{NLS\_DATA}}}
\bv\begin{lstlisting}
typedef struct nls_data NLS_DATA;
struct nls_data
{
  void       (*update)(void *, int, const REAL *, int, REAL *);
  void       *update_data;
  int        (*solve)(void *, int, const REAL *, REAL *);
  void       *solve_data;
  REAL       (*norm)(void *, int, const REAL *);
  void       *norm_data;

  WORKSPACE  *ws;

  REAL       tolerance;
  int        restart;
  int        max_iter;
  int        info;

  REAL       initial_residual;
  REAL       residual;
};
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{update} subroutine for computing a linearized system;

      \code{update(update\_data, dim, uk, update\_matrix, F)} computes
      a linearization of the system matrix, if \code{update\_matrix}
      is not zero, and the right hand side \code{F}, if \code{F} is
      not \nil, around the actual iterate \code{uk}; \code{dim} is the
      dimension of the nonlinear system, and \code{update\_data} a
      pointer to user data.
\kitem{update\_data} pointer to user data for the update of a linearized
      equation, first argument to \code{update()}.
\kitem{solve} function for solving a linearized system for the new 
      correction; the return value is the number of iterations used by
      an iterative solver or zero; this number is printed, if
      information about the solution process should be produced;

      \code{solve(solve\_data, dim, F, d)} solves the linearized equation
      of dimension \code{dim} with right hand side \code{F} for a
      correction \code{d} of the actual iterate; \code{d} is initialized
      with zeros and \code{update\_data} is a pointer to user
      data.
\kitem{solve\_data} pointer to user data for solution of the linearized
      equation, first argument to \code{solve()};

      the nonlinear solver does not know how the system matrix is stored;
      such information can be passed from \code{update()} to \code{solve()}
      by using pointers to the same \code{DOF} matrix in both 
      \code{update\_data} and \code{solve\_data}, e.g.
\kitem{norm} function for computing a problem dependent norm $\|.\|$;
      if \code{norm} is \nil, the Euclidian norm is used;

      \code{norm(norm\_data, dim, x)} returns the norm
      of the vector \code{x}; \code{dim} is the dimension of the nonlinear
      system, and \code{norm\_data} pointer to user data.
\kitem{norm\_data} pointer to user data for the calculation of the 
      problem dependent norm, first argument to \code{norm()}.
\kitem{ws} a pointer to a \code{WORKSPACE} structure for storing additional
       vectors used by a solver; if the space is not sufficient, the
       used solver will enlarge this workspace; if \code{ws} is \nil, then
       the used solver allocates memory, which is freed before exit.
\kitem{tolerance} tolerance for the nonlinear solver; if the norm of the 
       correction/residual is less or equal \code{tolerance}, the
       solver returns the actual iterate as the solution of the nonlinear
       system.
\kitem{restart} restart for the nonlinear solver.
\kitem{max\_iter} is a maximal number of iterations to be
       performed, even if the tolerance may not be reached.
\kitem{info} the level of information produced by the solver; \code{0}
       is the lowest level of information (no information is printed)
       and \code{4} the highest level.
\kitem{initial\_residual} stores the norm of the initial correction/residual
       on exit.
\kitem{residual} stores the norm of the last correction/residual on exit.
\end{descr}
The following Newton methods for solving \mathref{E:nonlin}
are currently implemented:
\fdx{nls_newton()@{\code{nls\_newton()}}}%
\idx{nonlinear solvers!nls_newton()@{\code{nls\_newton()}}}%
\fdx{nls_newton_ds()@{\code{nls\_newton\_ds()}}}%
\idx{nonlinear solvers!nls_newton_ds()@{\code{nls\_newton\_ds()}}}%
\fdx{nls_newton_fs()@{\code{nls\_newton\_fs()}}}%
\idx{nonlinear solvers!nls_newton_fs()@{\code{nls\_newton\_fs()}}}%
\fdx{nls_newton_br()@{\code{nls\_newton\_br()}}}%
\idx{nonlinear solvers!nls_newton_br()@{\code{nls\_newton\_br()}}}%
\bv\begin{lstlisting}
int nls_newton(NLS_DATA *, int, REAL *);
int nls_newton_ds(NLS_DATA *, int, REAL *);
int nls_newton_fs(NLS_DATA *, int, REAL *);
int nls_newton_br(NLS_DATA *, REAL, int, REAL *);
\end{lstlisting}\ev
Description:
\begin{descr}
\kitem{nls\_newton(nls\_data, dim, u0)} solves a nonlinear system by the
       classical Newton method; the return value is the number of iterations;

       \code{nls\_data} stores information about functions for the
       assemblage and solution of $DF(u_k)$, $F(u_k)$, calculation of
       a norm, tolerances, etc.  \code{dim} is the dimension of the
       nonlinear system, and \code{u0} the initial guess on input and
       the solution on output; \code{nls\_newton()} stops if the norm
       of the \textbf{correction} is less or equal
       \code{nls\_data->tolerance}; it needs a workspace for storing
       \code{2*dim} additional \code{REAL}s.
\kitem{nls\_newton\_ds(nls\_data, dim, u0)} solves a nonlinear system by a
       Newton method with step size control; the return value is the number 
       of iterations;

       \code{nls\_data} stores information about functions for the
       assembling and solving of $DF(u_k)$, $F(u_k)$, calculation of a
       norm, tolerances, etc.  \code{dim} is the dimension of the
       nonlinear system, and \code{u0} the initial guess on input and
       the solution on output; \code{nls\_newton\_ds()} stops if the norm
       of the \textbf{correction} is less or equal \code{nls\_data->tolerance};
       in each iteration at most \code{nls\_data->restart} steps for
       controlling the step size $\tau$ are performed; the aim is to
       choose $\tau$ such that
\[
\|DF(u_k)^{-1} F(u_k+\tau d_k)\| \le (1-{\textstyle\frac12}\tau) \|d_k\|
\]
       holds, where $\|.\|$ is the problem dependent norm, if 
       \code{nls\_data->norm} is not \nil, otherwise the Euclidian norm;
       each step needs the update of $F$, the solution of one
       linearized problem (the system matrix for the linearized
       system does not change during step size control) and the
       calculation of a norm;

       \code{nls\_newton\_ds()}
       needs a workspace for storing \code{4*dim} additional \code{REAL}s.
\kitem{nls\_newton\_fs(nls\_data, dim, u0)} solves a nonlinear system by a
       Newton method with step size control; the return value is the number 
       of iterations;

       \code{nls\_data} stores information about functions for the
       assembling and solving of $DF(u_k)$, $F(u_k)$, calculation of a
       norm, tolerances, etc.  \code{dim} is the dimension of the
       nonlinear system, and \code{u0} the initial guess on input and
       the solution on output; \code{nls\_newton\_fs()} stops if the norm
       of the \textbf{residual} is less or equal \code{nls\_data->tolerance};
       in each iteration at most \code{nls\_data->restart} steps for
       controlling the step size $\tau$ are performed; the aim is to
       choose $\tau$ such that
\[
\|F(u_k+\tau d_k)\| \le (1-{\textstyle\frac12}\tau) \|F(u_k)\|
\]
       holds, where $\|.\|$ is the problem dependent norm, if 
       \code{nls\_data->norm} is not \nil, otherwise the Euclidian norm;
       the step size control is not expensive, since in each step only
       an update of $F$ and the calculation of $\|F\|$ are 
       involved;

       \code{nls\_newton\_fs()}
       needs a workspace for storing \code{3*dim} additional \code{REAL}s.

\kitem{nls\_newton\_br(nls\_data, delta, dim, u0)} solves a nonlinear 
       system by a global Newton method by Bank and Rose \cite{BankRose:81};
       the return value is the number of iterations;
       
       \code{nls\_data} stores information about functions for the
       assembling and solving of $DF(u_k)$, $F(u_k)$, calculation of a
       norm, tolerances, etc.  \code{delta} is a parameter with
       $\delta\in(0,1-\alpha_0)$, where $\alpha_0 = \|DF(u_0)\, u_0 +
       F(u_0)\|/\|F(u_0)\|$; \code{dim} is the dimension of the
       nonlinear system, and \code{u0} the initial guess on input and
       the solution on output; \code{nls\_newton\_br()} stops if the
       norm of the \textbf{residual} is less or equal
       \code{nls\_data->tolerance}; in each iteration at most
       \code{nls\_data->restart} steps for controlling the step size
       by the method of Bank and Rose are performed; the step size
       control is not expensive, since in each step only an update of
       $F$ and the calculation of $\|F\|$ are involved;

       \code{nls\_newton\_br()}
       needs a workspace for storing \code{3*dim} additional \code{REAL}s.
\end{descr}
\idx{nonlinear solvers|)}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "alberta-man"
%%% End: 

% LocalWords:  Krylov DOF