1
\chapter{Linear algebra}
3
\devnote{Since this chapter was written, the \dolfin{} linear algebra
4
interface has undergone some changes that are not documented here.
5
As a result, some of the information in this chapter may be
6
outdated, incomplete or simply wrong.}
8
\dolfin{} provides a high-performance linear algebra library,
9
including matrices and vectors, a set of linear solvers,
10
preconditioners, and an eigenvalue solver. The core part of the
11
functionality is provided through wrappers that provide a common
12
interface to the functionality of the linear algebra libraries
13
uBLAS~\cite{www:ublas}, PETSc~\cite{www:petsc} and
14
Trilinos~\cite{www.trilonos}.
16
%------------------------------------------------------------------------------
17
\section{Matrices and vectors}
18
\index{\texttt{Matrix}}
19
\index{\texttt{Vector}}
21
The two basic linear algebra data structures are the classes
22
\texttt{Matrix} and \texttt{Vector}, representing a sparse $M\times
23
N$ matrix and a vector of length $N$ respectively.
25
The following code demonstrates how to create a matrix and a vector:
30
Alternatively, the matrix and the vector may be created by
39
The following code demonstrates how to access the size and the
40
elements of a matrix and a vector:
48
unsigned int M = A.size(0);
49
unsigned int N = A.size(1);
54
In addition, \dolfin{} provides optimized functions for setting the
55
values of a set of entries in a matrix or vector:
57
double block[] = {2, 4, 6};
58
int rows[] = {0, 1, 2};
59
int cols[] = {0, 1, 2};
61
A.set(block, rows, cols, 3);
63
Alternatively, the set of values given by the array \texttt{block} can
64
be added to the entries given by the arrays \texttt{rows} and
67
double block[] = {2, 4, 6};
68
int rows[] = {0, 1, 2};
69
int cols[] = {0, 1, 2};
71
A.add(block, rows, cols, 3);
73
These functions are particularly useful for efficient assembly of a (sparse)
74
matrix from a bilinear form.
76
\subsection{Sparse matrices}
79
The default \dolfin{} class \texttt{Matrix} is a sparse matrix, which
80
efficiently represents the discretization of a partial differential
81
equation where most entries in the matrix are zero. Alternatively,
82
the class \texttt{SparseMatrix} may be used which is identical with
83
the class \texttt{Matrix}\footnote{The class \texttt{Matrix} is a
84
\texttt{typedef} for the class \texttt{SparseMatrix}.}.
86
If \dolfin{} has been compiled with support for PETSc, then the sparse
87
matrix is represented as a sparse PETSc matrix\footnote{By default, the
88
sparse matrix is represented as a PETSc \texttt{MATSEQAIJ} matrix, but
89
other PETSc representations are also available.}. Alternatively, the
90
class \texttt{PETScMatrix} may be used, together with the corresponding
91
class \texttt{PETScVector}.
93
If \dolfin{} has been compiled without support for PETSc, then the
94
sparse matrix is represented as a uBLAS sparse matrix. Alternatively,
95
the class \texttt{uBLASSparseMatrix} may be used, together with the
96
corresponding class \texttt{uBLASVector}.
98
\subsection{Dense matrices}
101
\dolfin{} provides the class \texttt{DenseMatrix} for representation
102
of dense matrices. A dense matrix representation is often preferable
103
when computing with matrices of small to moderate size. In particular,
104
accessing individual elements
105
is more efficient with a dense matrix representation.
107
A \texttt{DenseMatrix} is represented as uBLAS dense matrix and
108
alternatively the class \texttt{uBLASDenseMatrix} may be used,
109
together with the corresponding class \texttt{uBLASVector}.
111
\subsection{The common interface}
112
\index{\texttt{GenericMatrix}}
113
\index{\texttt{GenericVector}}
115
Although \dolfin{} differentiates between sparse and dense data
116
structures, the two classes \texttt{GenericMatrix} and
117
\texttt{GenericVector} define a common interface to all matrices and
118
vectors. Refer to the \emph{DOLFIN programmer's reference} for the
119
exact specification of these interfaces.
121
The following points should be noted:
124
If you want a specific backend like PETSc, then use
125
\texttt{PETScVector}/\texttt{PETScMatrix}.
127
If you don't care about the backend, then use \texttt{Vector}/\texttt{Matrix}.
129
If you write a \emph{function} that should be able to accept vectors
130
and matrices from any backend as input, then use
131
\texttt{GenericVector}/\texttt{GenericMatrix}.
134
%------------------------------------------------------------------------------
135
\section{Solving linear systems}
136
\index{linear systems}
138
\dolfin{} provides a set of efficient solvers for linear systems of
143
where $A$ is an $N\times N$ matrix and where $x$ and $b$ are vectors
144
of length $N$. Both iterative (Krylov subspace) solvers and direct
145
(LU) solvers are provided.
147
\subsection{Iterative methods}
148
\index{iterative methods}
149
\index{\texttt{GMRES}}
151
\index{\texttt{KrylovSolver}}
152
\index{Krylov subspace methods}
154
A linear system may be solved by the GMRES Krylov method as follows:
159
GMRES::solve(A, x, b);
161
Alternatively, the linear system may be solved by first creating an
162
object of the class \texttt{KrylovSolver}, which is more efficient for
163
repeated solution of linear systems and also allows the specification
164
of both the Krylov method and the preconditioner:
166
KrylovSolver solver(gmres, ilu);
167
solver.solve(A, x, b);
170
For uBLAS matrices and vectors, the class \texttt{uBLASKrylovSolver}
171
may be used and for PETSc matrices and vectors, the class
172
\texttt{PETScKrylovSolver} may be used.
174
\subsubsection{Krylov methods}
175
\index{conjugate gradient method}
178
\dolfin{} provides the following set of Krylov methods:
181
\begin{tabular}{|l|l|}
183
\texttt{cg} & The conjugate gradient method \\
185
\texttt{gmres} & The GMRES method \\
187
\texttt{bicgstab} & The stabilized biconjugate gradient squared
190
\texttt{default\_method} & Default choice of Krylov method \\
195
\subsubsection{Preconditioners}
196
\index{preconditioners}
198
\index{incomplete LU factorization}
200
\index{successive over-relaxation}
203
\index{algebraic multigrid}
205
\dolfin{} provides the following set of preconditioners:
208
\begin{tabular}{|l|l|}
210
\texttt{none} & No preconditioning \\
212
\texttt{jacobi} & Simple Jacobi preconditioning \\
214
\texttt{sor} & SOR, successive over-relaxation \\
216
\texttt{ilu} & Incomplete LU factorization \\
218
\texttt{icc} & Incomplete Cholesky factorization \\
220
\texttt{amg} & Algebraic multigrid (through Hypre when available) \\
222
\texttt{default\_pc} & Default choice of preconditioner \\
227
\subsubsection{Matrix-free solvers}
228
\index{matrix-free solvers}
229
\index{virtual matrix}
231
The \dolfin{} Krylov solvers may be used without direct access to a
232
matrix representation. All that is needed is to provide the
233
size of the linear system, the right-hand side, and a method
234
implementing the multiplication of the matrix with any given vector.
236
Such a ``virtual matrix'' may be defined by implementing the
237
interface defined by either the class \texttt{uBLASKrylovMatrix} of
238
\texttt{PETScKrylovMatrix}. The matrix may then be used together with
239
either the class \texttt{uBLASKrylovSolver} or \texttt{PETScKrylovSolver}.
241
\subsection{Direct methods}
243
\index{\texttt{LUSolver}}
244
\index{direct methods}
245
\index{LU factorization}
247
A linear system may be solved by a direct LU factorization as follows:
254
Alternatively, the linear system may be solved by first creating an
255
object of the class \texttt{LUSolver}, which may be more efficient for
256
repeated solution of linear systems:
259
solver.solve(A, x, b);
262
For uBLAS matrices and vectors, the class \texttt{uBLASLUSolver} may
263
be used and for PETSc matrices and vectors, the class
264
\texttt{PETScLUSolver} may be used.
266
%------------------------------------------------------------------------------
267
\section{Eigenvalue problems}
268
\index{eigenvalue problems}
271
\dolfin{} also provides a solver for eigenvalue problems\index{eigenvalue solver}.
272
The solver is only available when \dolfin{} has been compiled with support for
273
PETSc and SLEPc~\cite{www:slepc}.
275
For the basic eigenvalue problem
281
the following code demonstrates how to compute the first eigenpair:
284
SLEPcEigenvalueSolver esolver;
289
esolver.getEigenpair(lr, lc, xr, xc, 0);
292
The double and complex components of the eigenvalue are returned in \texttt{lr}
293
and \texttt{lc}, respectively, and the double and complex parts of the eigenvector
294
are returned in \texttt{xr} and \texttt{xc}, respectively.
296
For the generalized eigenvalue problem
300
the following code demonstrates how to compute the third eigenpair:
302
PETScEigenvalueSolver esolver;
307
esolver.getEigenpair(lr, lc, xr, xc, 2);
309
%------------------------------------------------------------------------------
310
\section{Linear algebra backends}
311
\index{linear algebra backends}
316
uBLAS is a C++ template library that provides BLAS level 1, 2 and 3
317
functionality (and more) for dense, packed and sparse matrices.
318
The design and implementation unify mathematical notation via operator
319
overloading and efficient code generation via expression templates.
321
\dolfin{} wrappers for uBLAS linear algebra is provided through the
322
classes \texttt{uBLASSparseMatrix}, \texttt{uBLASDenseMatrix} and
323
\texttt{uBLASVector}. These classes are implemented by subclassing the
324
corresponding uBLAS classes, which means that all standard
325
\texttt{uBLAS} operations are supported for these classes. For
326
advanced usage not covered by the common \dolfin{} interface specified
327
by the classes \texttt{GenericMatrix} and \texttt{GenericVector},
328
refer directly to the documentation of \texttt{uBLAS}.
329
%------------------------------------------------------------------------------
333
PETSc is a suite of data structures and routines for the scalable
334
(parallel) solution of scientific applications modeled by partial
335
differential equations. It employs the MPI standard for all
336
message-passing communication.
338
\dolfin{} wrappers for PETSc linear algebra is provided through the
339
classes \texttt{PETScMatrix} and \texttt{PETScVector}. Direct access
340
to the PETSc data structures is available through the member functions
341
\texttt{mat()} and \texttt{vec()}, which return the PETSc \texttt{Mat}
342
and \texttt{Vec} pointers respectively. For advanced usage not covered
343
by the common \dolfin{} interface specified by the classes
344
\texttt{GenericMatrix} and \texttt{GenericVector}, refer directly to
345
the documentation of PETSc.
347
%------------------------------------------------------------------------------
348
\subsection{Trilinos}
352
%------------------------------------------------------------------------------
353
\subsection{Matrix Template Library (MTL4)}
355
\index{Matrix Template Library}
358
%------------------------------------------------------------------------------
359
\subsection{User provided linear algebra backends}
362
\devnote{This section will explain how a user-provided linear alegra backend
363
can be used with \dolfin{}.}
added
removed
doc/manual/chapters/linearalgebra.tex
