|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: sparse linear algebra
|
|
Packit |
67cb25 |
single: linear algebra, sparse
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*********************
|
|
Packit |
67cb25 |
Sparse Linear Algebra
|
|
Packit |
67cb25 |
*********************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This chapter describes functions for solving sparse linear systems
|
|
Packit |
67cb25 |
of equations. The library provides linear algebra routines which
|
|
Packit |
67cb25 |
operate directly on the :type:`gsl_spmatrix` and :type:`gsl_vector`
|
|
Packit |
67cb25 |
objects.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this chapter are declared in the header file
|
|
Packit |
67cb25 |
:file:`gsl_splinalg.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: sparse linear algebra, overview
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Overview
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This chapter is primarily concerned with the solution of the
|
|
Packit |
67cb25 |
linear system
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A x = b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`A` is a general square :math:`n`-by-:math:`n` non-singular
|
|
Packit |
67cb25 |
sparse matrix, :math:`x` is an unknown :math:`n`-by-:math:`1` vector, and
|
|
Packit |
67cb25 |
:math:`b` is a given :math:`n`-by-1 right hand side vector. There exist
|
|
Packit |
67cb25 |
many methods for solving such sparse linear systems, which broadly
|
|
Packit |
67cb25 |
fall into either direct or iterative categories. Direct methods include
|
|
Packit |
67cb25 |
LU and QR decompositions, while iterative methods start with an
|
|
Packit |
67cb25 |
initial guess for the vector :math:`x` and update the guess through
|
|
Packit |
67cb25 |
iteration until convergence. GSL does not currently provide any
|
|
Packit |
67cb25 |
direct sparse solvers.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: sparse matrices, iterative solvers
|
|
Packit |
67cb25 |
single: sparse linear algebra, iterative solvers
|
|
Packit |
67cb25 |
single: sparse, iterative solvers
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Sparse Iterative Solvers
|
|
Packit |
67cb25 |
========================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Overview
|
|
Packit |
67cb25 |
--------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Many practical iterative methods of solving large :math:`n`-by-:math:`n`
|
|
Packit |
67cb25 |
sparse linear systems involve projecting an approximate solution for
|
|
Packit |
67cb25 |
:data:`x` onto a subspace of :math:`{\bf R}^n`. If we define a :math:`m`-dimensional
|
|
Packit |
67cb25 |
subspace :math:`{\cal K}` as the subspace of approximations to the solution
|
|
Packit |
67cb25 |
:data:`x`, then :math:`m` constraints must be imposed to determine
|
|
Packit |
67cb25 |
the next approximation. These :math:`m` constraints define another
|
|
Packit |
67cb25 |
:math:`m`-dimensional subspace denoted by :math:`{\cal L}`. The
|
|
Packit |
67cb25 |
subspace dimension :math:`m` is typically chosen to be much smaller than
|
|
Packit |
67cb25 |
:math:`n` in order to reduce the computational
|
|
Packit |
67cb25 |
effort needed to generate the next approximate solution vector.
|
|
Packit |
67cb25 |
The many iterative algorithms which exist differ mainly
|
|
Packit |
67cb25 |
in their choice of :math:`{\cal K}` and :math:`{\cal L}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Types of Sparse Iterative Solvers
|
|
Packit |
67cb25 |
---------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The sparse linear algebra library provides the following types
|
|
Packit |
67cb25 |
of iterative solvers:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_splinalg_itersolve_type
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: gmres
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_splinalg_itersolve_gmres
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This specifies the Generalized Minimum Residual Method (GMRES).
|
|
Packit |
67cb25 |
This is a projection method using :math:`{\cal K} = {\cal K}_m`
|
|
Packit |
67cb25 |
and :math:`{\cal L} = A {\cal K}_m` where :math:`{\cal K}_m` is
|
|
Packit |
67cb25 |
the :math:`m`-th Krylov subspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: {\cal K}_m = span \left\{ r_0, A r_0, ..., A^{m-1} r_0 \right\}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
K_m = span( r_0, A r_0, ..., A^(m-1) r_0)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and :math:`r_0 = b - A x_0` is the residual vector of the initial guess
|
|
Packit |
67cb25 |
:math:`x_0`. If :math:`m` is set equal to :math:`n`, then the Krylov
|
|
Packit |
67cb25 |
subspace is :math:`{\bf R}^n` and GMRES will provide the exact solution
|
|
Packit |
67cb25 |
:data:`x`. However, the goal is for the method to arrive at a very good
|
|
Packit |
67cb25 |
approximation to :data:`x` using a much smaller subspace :math:`{\cal K}_m`. By
|
|
Packit |
67cb25 |
default, the GMRES method selects :math:`m = MIN(n,10)` but the user
|
|
Packit |
67cb25 |
may specify a different value for :math:`m`. The GMRES storage
|
|
Packit |
67cb25 |
requirements grow as :math:`O(n(m+1))` and the number of flops
|
|
Packit |
67cb25 |
grow as :math:`O(4 m^2 n - 4 m^3 / 3)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In the below function :func:`gsl_splinalg_itersolve_iterate`, one
|
|
Packit |
67cb25 |
GMRES iteration is defined as projecting the approximate solution
|
|
Packit |
67cb25 |
vector onto each Krylov subspace :math:`{\cal K}_1, ..., {\cal K}_m`,
|
|
Packit |
67cb25 |
and so :math:`m` can be kept smaller by "restarting" the method
|
|
Packit |
67cb25 |
and calling :func:`gsl_splinalg_itersolve_iterate` multiple times,
|
|
Packit |
67cb25 |
providing the updated approximation :data:`x` to each new call. If
|
|
Packit |
67cb25 |
the method is not adequately converging, the user may try increasing
|
|
Packit |
67cb25 |
the parameter :math:`m`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GMRES is considered a robust general purpose iterative solver, however
|
|
Packit |
67cb25 |
there are cases where the method stagnates if the matrix is not
|
|
Packit |
67cb25 |
positive-definite and fails to reduce the residual until the very last
|
|
Packit |
67cb25 |
projection onto the subspace :math:`{\cal K}_n = {\bf R}^n`. In these
|
|
Packit |
67cb25 |
cases, preconditioning the linear system can help, but GSL does not
|
|
Packit |
67cb25 |
currently provide any preconditioners.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Iterating the Sparse Linear System
|
|
Packit |
67cb25 |
----------------------------------
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions are provided to allocate storage for the
|
|
Packit |
67cb25 |
sparse linear solvers and iterate the system to a solution.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_splinalg_itersolve * gsl_splinalg_itersolve_alloc (const gsl_splinalg_itersolve_type * T, const size_t n, const size_t m)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for the iterative solution of
|
|
Packit |
67cb25 |
:data:`n`-by-:data:`n` sparse matrix systems. The iterative solver type
|
|
Packit |
67cb25 |
is specified by :data:`T`. The argument :data:`m` specifies the size
|
|
Packit |
67cb25 |
of the solution candidate subspace :math:`{\cal K}_m`. The dimension
|
|
Packit |
67cb25 |
:data:`m` may be set to 0 in which case a reasonable default value is used.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_splinalg_itersolve_free (gsl_splinalg_itersolve * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: const char * gsl_splinalg_itersolve_name (const gsl_splinalg_itersolve * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a string pointer to the name of the solver.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_splinalg_itersolve_iterate (const gsl_spmatrix * A, const gsl_vector * b, const double tol, gsl_vector * x, gsl_splinalg_itersolve * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function performs one iteration of the iterative method for
|
|
Packit |
67cb25 |
the sparse linear system specfied by the matrix :data:`A`, right hand
|
|
Packit |
67cb25 |
side vector :data:`b` and solution vector :data:`x`. On input, :data:`x`
|
|
Packit |
67cb25 |
must be set to an initial guess for the solution. On output,
|
|
Packit |
67cb25 |
:data:`x` is updated to give the current solution estimate. The
|
|
Packit |
67cb25 |
parameter :data:`tol` specifies the relative tolerance between the residual
|
|
Packit |
67cb25 |
norm and norm of :data:`b` in order to check for convergence.
|
|
Packit |
67cb25 |
When the following condition is satisfied:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: || A x - b || \le tol \times || b ||
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|| A x - b || <= tol * || b ||
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
the method has converged, the function returns :macro:`GSL_SUCCESS` and
|
|
Packit |
67cb25 |
the final solution is provided in :data:`x`. Otherwise, the function
|
|
Packit |
67cb25 |
returns :macro:`GSL_CONTINUE` to signal that more iterations are
|
|
Packit |
67cb25 |
required. Here, :math:`|| \cdot ||` represents the Euclidean norm.
|
|
Packit |
67cb25 |
The input matrix :data:`A` may be in triplet or compressed format.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_splinalg_itersolve_normr (const gsl_splinalg_itersolve * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns the current residual norm
|
|
Packit |
67cb25 |
:math:`||r|| = ||A x - b||`, which is updated after each call to
|
|
Packit |
67cb25 |
:func:`gsl_splinalg_itersolve_iterate`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: sparse linear algebra, examples
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Examples
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This example program demonstrates the sparse linear algebra routines on
|
|
Packit |
67cb25 |
the solution of a simple 1D Poisson equation on :math:`[0,1]`:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: {d^2 u(x) \over dx^2} = f(x) = -\pi^2 \sin{(\pi x)}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
u''(x) = f(x) = -\pi^2 \sin(\pi x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with boundary conditions :math:`u(0) = u(1) = 0`. The analytic solution of
|
|
Packit |
67cb25 |
this simple problem is :math:`u(x) = \sin{\pi x}`. We will solve this
|
|
Packit |
67cb25 |
problem by finite differencing the left hand side to give
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: {1 \over h^2} \left( u_{i+1} - 2 u_i + u_{i-1} \right) = f_i
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
1/h^2 ( u_(i+1) - 2 u_i + u_(i-1) ) = f_i
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Defining a grid of :math:`N` points, :math:`h = 1/(N-1)`. In the finite
|
|
Packit |
67cb25 |
difference equation above, :math:`u_0 = u_{N-1} = 0` are known from
|
|
Packit |
67cb25 |
the boundary conditions, so we will only put the equations for
|
|
Packit |
67cb25 |
:math:`i = 1, ..., N-2` into the matrix system. The resulting
|
|
Packit |
67cb25 |
:math:`(N-2) \times (N-2)` matrix equation is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{1 \over h^2}
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{array}{cccccc}
|
|
Packit |
67cb25 |
-2 & 1 & 0 & 0 & \ldots & 0 \\
|
|
Packit |
67cb25 |
1 & -2 & 1 & 0 & \ldots & 0 \\
|
|
Packit |
67cb25 |
0 & 1 & -2 & 1 & \ldots & 0 \\
|
|
Packit |
67cb25 |
\vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\
|
|
Packit |
67cb25 |
0 & \ldots & \ldots & 1 & -2 & 1 \\
|
|
Packit |
67cb25 |
0 & \ldots & \ldots & \ldots & 1 & -2
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{array}{c}
|
|
Packit |
67cb25 |
u_1 \\
|
|
Packit |
67cb25 |
u_2 \\
|
|
Packit |
67cb25 |
u_3 \\
|
|
Packit |
67cb25 |
\vdots \\
|
|
Packit |
67cb25 |
u_{N-3} \\
|
|
Packit |
67cb25 |
u_{N-2}
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right) =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{array}{c}
|
|
Packit |
67cb25 |
f_1 \\
|
|
Packit |
67cb25 |
f_2 \\
|
|
Packit |
67cb25 |
f_3 \\
|
|
Packit |
67cb25 |
\vdots \\
|
|
Packit |
67cb25 |
f_{N-3} \\
|
|
Packit |
67cb25 |
f_{N-2}
|
|
Packit |
67cb25 |
\end{array}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
An example program which constructs and solves this system is given below.
|
|
Packit |
67cb25 |
The system is solved using the iterative GMRES solver. Here is
|
|
Packit |
67cb25 |
the output of the program::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
iter 0 residual = 4.297275996844e-11
|
|
Packit |
67cb25 |
Converged
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
showing that the method converged in a single iteration.
|
|
Packit |
67cb25 |
The calculated solution is shown in :numref:`fig_splinalg-poisson`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _fig_splinalg-poisson:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. figure:: /images/sparse_poisson.png
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solution of PDE
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/poisson.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: sparse linear algebra, references
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The implementation of the GMRES iterative solver closely follows
|
|
Packit |
67cb25 |
the publications
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* H. F. Walker, Implementation of the GMRES method using
|
|
Packit |
67cb25 |
Householder transformations, SIAM J. Sci. Stat. Comput.
|
|
Packit |
67cb25 |
9(1), 1988.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Y. Saad, Iterative methods for sparse linear systems, 2nd edition,
|
|
Packit |
67cb25 |
SIAM, 2003.
|