|
|
Boost.Numeric.Bindings.UMFPACK |
Ax=b and related systems, where A is
sparse and unsymmetric coefficient matrix, b is right-hand side vector and
x is vector of unknowns.
Matrix A can be square or rectangular, regular or
singular.
Rectangular matrices can only be factorised; to solve the system
Ax=b, matrix A must be square.
(If solving system with singular A, divide by zero will
occur and solutions will contain Infs and/or
NaNs, but parts of the solution may be valid.)
Matrix elements can be real (C/C++ type double) or
complex (real and imaginary part are represented by doubles).
UMFPACK is based on the
Unsymmetric-pattern MultiFrontal method.
It factorises PAQ, PRAQ or
PR-1AQ into LU, where:
P and Q are permutation matrices --
column ordering Q is selected to give a good a
priori upper bound on fill-in and is refined during numerical
factorisation, while the row ordering P has the dual
role of preserving sparsity and maintaining numerical accuracy via
relaxed partial pivoting and row interchanges;
R is a diagonal matrix of row scaling factors;
L and U are lower and upper triangular matrices.
UMFPACK (version 4.1) is written in ANSI C. Interface of the so-called primary routines, with a simple example, is described in: Timothy A. Davis, UMFPACK Version 4.1 Quick Start Guide, TR-03-009 (PDF); see also: Timothy A. Davis, Algorithm 8xx: UMFPACK V4.1, an unsymmetric-pattern multifrontal method, TR-03-007 (PDF or PS). Detailed description of the C interface is given in: Timothy A. Davis, UMFPACK Version 4.1 User Guide, TR-03-008 (PDF).
UMFPACK Bindings Library provides generic, higher level interface to UMFPACK functions. By `generic' we mean that bindings are based on traits and therefore can be used with various matrix and vector classes, and by `higher level' that some complexities of the C interface are hidden from the user.
boost/numeric/bindings/umfpack.
C++ program should include umfpack.hpp.
All `public'
functions are in the namespace boost::numeric::bindings::umfpack.
The location of UMFPACK header files may depend on your installation, so
probably you should edit umfpack_incl.hpp providing your path
to UMFPACK's header umfpack.h.
#include <iostream>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/bindings/traits/ublas_sparse.hpp>
#include <boost/numeric/bindings/umfpack/umfpack.hpp>
#include <boost/numeric/ublas/io.hpp>
namespace ublas = boost::numeric::ublas;
namespace umf = boost::numeric::bindings::umfpack;
int main() {
ublas::compressed_matrix<double, ublas::column_major, 0,
ublas::unbounded_array<int>, ublas::unbounded_array<double> > A (5,5,12);
ublas::vector<double> B (5), X (5);
A(0,0) = 2.; A(0,1) = 3;
A(1,0) = 3.; A(1,2) = 4.; A(1,4) = 6;
A(2,1) = -1.; A(2,2) = -3.; A(2,3) = 2.;
A(3,2) = 1.;
A(4,1) = 4.; A(4,2) = 2.; A(4,4) = 1.;
B(0) = 8.; B(1) = 45.; B(2) = -3.; B(3) = 3.; B(4) = 19.;
umf::symbolic_type<double> Symbolic;
umf::numeric_type<double> Numeric;
umf::symbolic (A, Symbolic);
umf::numeric (A, Symbolic, Numeric);
umf::solve (A, X, B, Numeric);
std::cout << X << std::endl; // output: [5](1,2,3,4,5)
}
Solution of the linear system AX = B consists of three
steps: symbolic analysis of the matrix A, numerical
factorisation and backward substition. In more details:
symbolic()
pre-orders the columns of A to reduce fill-in (using COLAMD, AMD
or user-provided ordering), finds the supernodal column elimination
tree and post-orders the tree.
(Matrix A can be rectangular or square.)
Symbolic analysis details are returned in an object of type symbolic_type<> (see sec. 4) which is passed by reference.
numeric()
performs the numerical factorisation PAQ = LU,
PRAQ = LU or PR-1AQ = LU,
using the symbolic analysis previously computed by symbolic().
Numerical factorisation is returned in an object of type numeric_type<> (see sec. 4)
which is passed by reference.
solve()
solves a linear system for the solution X,
using numerical factorisation previously computed by numeric().
Only square systems are handled. Optionally, function performs
iterative refinement.
Functions symbolic(), numeric() and
solve() correspond to UMFPACK's functions
umfpack_*_symbolic()/umfpack_*_qsymbolic(),
umfpack_*_numeric() and umfpack_*_solve(),
where * can be di or zi,
that is, matrix and vector elements can be doubles or
complex numbers represented by two doubles, while
indices, sizes etc. have type int.
Bindings provide two additional functions:
factor() combines symbolic analysis and
numerical factorisation (symbolic_type<> object is
used only internally):
// ...
umf::numeric_type Numeric;
umf::factor (A, Numeric);
umf::solve (A, X, B, Numeric);
umf_solve() combines all three steps
(both symbolic_type<> and numeric_type<> objects are used internally):
// ...
umf::umf_solve (A, X, B);
| Symbol | Type |
|---|---|
status, sys | int
|
A | matrix (see sec. 3) |
B, X | vector |
Qinit | vector of ints
|
Symbolic | symbolic_type<>
(see sec. 4)
|
Numeric | numeric_type<>
(see sec. 4)
|
Control | control_type<>
(see sec. 4)
|
Info | info_type<>
(see sec. 4).
|
For possible values assigned to status
see UMFPACK User Guide.
symbolic()umfpack_*_symbolic() or
umfpack_*_qsymbolic().)
Syntax:
| 1 | status = symbolic (A, Symbolic);
|
| 2 | status = symbolic (A, Symbolic, Control);
|
| 3 | status = symbolic (A, Symbolic, Control, Info);
|
| 4 | status = symbolic (A, Qinit, Symbolic);
|
| 5 | status = symbolic (A, Qinit, Symbolic, Control);
|
| 6 | status = symbolic (A, Qinit, Symbolic, Control, Info);
|
Symbolic analysis details are returned in Symbolic.
1, 2, 3 use COLAMD or AMD for column pre-ordering, while 4, 5, 6
use pre-ordering given in Qinit (for more details
see UMFPACK documentation).
1, 4 use default control settings; 2, 3, 5, 6 use control settings
given in Control.
3, 6 return statistics about ordering in Info.
numeric()Symbolic
computed previously by symbolic(). (Calls umfpack_*_numeric()).
Syntax:
| 1 | status = numeric (A, Symbolic, Numeric);
|
| 2 | status = numeric (A, Symbolic, Numeric, Control);
|
| 3 | status = numeric (A, Symbolic, Numeric, Control, Info);
|
Numerical factorisation is returned in Numeric.
1 uses default control settings; 2, 3 use control settings
given in Control.
3 returns statistics about factorisation in Info.
factor()umfpack_*_symbolic() or
umfpack_*_qsymbolic() and then umfpack_*_numeric().)
Syntax:
| 1 | status = factor (A, Numeric);
|
| 2 | status = factor (A, Numeric, Control);
|
| 3 | status = factor (A, Numeric, Control, Info);
|
Numerical factorisation is returned in Numeric.
1 uses default control settings; 2, 3 use control settings
given in Control.
3 returns statistics about factorisation in Info.
solve()AX = B or related systems (see
below) using Numeric computed previously by
numeric() or factor().
(Calls umfpack_*_solve()).
Syntax:
| 1 | status = solve (A, X, B, Numeric);
|
| 2 | status = solve (A, X, B, Numeric, Control);
|
| 3 | status = solve (A, X, B, Numeric, Control, Info);
|
| 4 | status = solve (sys, A, X, B, Numeric);
|
| 5 | status = solve (sys, A, X, B, Numeric, Control);
|
| 6 | status = solve (sys, A, X, B, Numeric, Control, Info);
|
1, 2, 3 solve system AX = B; 4, 5, 6 solve system denoted
by sys:
sys | System |
|---|---|
UMFPACK_A | AX = B
|
UMFPACK_At | AHX = B
|
UMFPACK_Aat | ATX = B
|
UMFPACK_Pt_L | PTLX = B
|
UMFPACK_L | LX = B
|
UMFPACK_Lt_P | LHPX = B
|
UMFPACK_Lat_P | LTPX = B
|
UMFPACK_Lt | LHX = B
|
UMFPACK_Lat | LTX = B
|
UMFPACK_U_Qt | UQTX = B
|
UMFPACK_U | UX = B
|
UMFPACK_Q_Ut | QUHX = B
|
UMFPACK_Q_Uat | QUTX = B
|
UMFPACK_Ut | UHX = B
|
UMFPACK_Uat | UTX = B
|
1, 4 use default control settings; 2, 3, 5, 6 use control settings
given in Control.
3, 6 return statistics about ordering in Info.
umf_solve()AX = B.
(Calls
umfpack_*_symbolic()/umfpack_*_qsymbolic(),
umfpack_*_numeric() and umfpack_*_solve().)
Syntax:
| 1 | status = umf_solve (A, X, B);
|
| 2 | status = umf_solve (A, X, B, Control);
|
| 3 | status = umf_solve (A, X, B, Control, Info);
|
1 uses default control settings; 2, 3 use control settings
given in Control.
3 returns statistics about factorisation in Info.
free()symbolic_type<> or
numeric_type<> objects.
(Calls umfpack_*_free_symbolic() or
umfpack_*_free_numeric().)
Syntax:
free (Symbolic);
| |
Symbolic.free();
| |
free (Numeric);
| |
Numeric.free();
|
defaults()control_type<> object's array.
(Calls umfpack_*_defaults().
Syntax:
defaults (Control);
| |
Control.defaults();
|
scale()X = B, X = RB
or X = R-1B, as appropriate.
(Calls umfpack_*_scale().)
Syntax:
status = scale (X, B, Numeric);
|
report_status()umfpack_*_report_status().)
Syntax:
report_status (Control, status);
|
For control parameters used in various report_* functions
see UMFPACK User Guide.
report_control()umfpack_*_report_control().)
Syntax:
report_control (Control);
|
report_info()Info by
symbolic(), numeric(),
factor(), solve() and umf_solve().
(Calls umfpack_*_report_info().)
Syntax:
report_info (Control, Info);
|
report_matrix()umfpack_*_report_matrix() or umfpack_*_report_triplet().)
Syntax:
status = report_matrix (A, Control);
|
report_vector()umfpack_*_report_vector().)
Syntax:
status = report_vector (X, Control);
|
report_symbolic()symbolic_type<> object.
(Calls umfpack_*_report_symbolic().)
Syntax:
status = report_symbolic (Symbolic, Control);
|
report_numeric()numeric_type<> object.
(Calls umfpack_*_report_numeric().)
Syntax:
status = report_numeric (Numeric, Control);
|
report_permutation()umfpack_*_report_perm().)
Syntax:
status = report_permutation (Qinit, Control);
|
UMFPACK functions manipulate matrices stored in compressed column format. Let
A be m-by-n with nnz nonzero double entries.
Its compressed column storage consists of three arrays:
int Ap[n+1];
int Ai[nnz];
double Ax[nnz];
Ax contains nonzero entries: elements of the matrix column
k, sorted in ascending order of corresponding row indices
and without duplicates, are in Ax[Ap[k]] to (and including)
Ax[Ap[k+1]-1] and their row indices are in
Ai[Ap[k]] to Ai[Ap[k+1]-1].
Thus, Ap is the array of column start locations in
Ai and Ax;
Ap[0] == 0 and Ap[n] == nnz.
This is the storage scheme of
ublas::compressed_matrix<double, ublas::column_major, 0>.
For example, the matrix
[2 3 0 0 0]
[3 0 4 0 6]
[0 -1 -3 2 0]
[0 0 1 0 0]
[0 4 2 0 1]
is represented by
int Ap[] = { 0, 2, 5, 9, 10, 12 };
int Ai[] = { 0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4 };
double Ax[] = { 2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1. };
For complex matrices, array Ax contains the real parts of
complex numbers and there is an additional array for imaginary parts:
double Az[nnz];
Note that this is not the storage scheme of
ublas::compressed_matrix<std::complex<double>, ublas::column_major, 0>
which contains one array of complex values.
But, UMFPACK bindings internally use function traits::detail::disentangle(),
which converts array of complex numbers into two real arrays holding
real and imaginary parts, and
function traits::detail::interlace(), which performs
opposite conversion, so you can still write:
ublas::compressed_matrix<std::complex<double>, ublas::column_major> C (5,5,12);
std::vector<std::complex<double> > B (5), X (5);
// ...
umf::umf_solve (C, X, B);
UMFPACK provide functions umfpack_*_triplet_to_col()
which convert coordinate (AKA triplet) storage format of a sparse matrix to
compressed column format. The coordinate format consists of three arrays:
int Ti[nnz];
int Tj[nnz];
double Tx[nnz];
Matrix entry aij is stored in kth triplet such
that Ti[k] == i, Tj[k] == j and
Tx[k] == aij. In general, triplets can be in any order
(moreover, duplicate entries can exist). If column indices are stored
in ascending order and if, within each column, row indices are sorted
in ascending order (and if there is no duplicates), that is, for
previous example
int Ti[] = { 0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4 };
int Tj[] = { 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4 };
double Tx[] = { 2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1. };
UMFPACK bindings can be used with matrices in coordinate format
(internally, bindings functions use umfpack_*_triplet_to_col()
to create just the pattern of the compressed column matrix
(arrays Ap and Ai), and array
Tx is used instead of Ax).
Instead of ublas::compressed_matrix<> one can use
ublas::coordinate_matrix<double, ublas::column_major, 0> A (5,5,12);
in the example in section 2.
(Conversion will be performed three times -- in symbolic(),
numeric() and solve(). But if you use
umf_solve(), conversion will be performed only once.)
symbolic_type<> and numeric_type<>
wrap void pointers which UMFPACK uses to pass around
symbolic analysis and numerical factorisation information.
Their destructors call umfpack_*_free_symbolic() and
umfpack_*_free_numeric(), respectively;
see also free().
Classes control_type<> and info_type<>
wrap UMFPACK's control and info arrays.
For arrays' entries (e.g. control parameters and error codes)
see UMFPACK User Guide.
Both classes provide operator[] for accessing individual entries; e.g.:
control_type<double> Control;
Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
control_type<>'s constructor sets default
parameters in the control array (i.e. calls umfpack_*_defaults());
see also defaults().
Template argument when constructing objects of these classes must be the type of matrix and vector elements:
compressed_matrix<std::complex<double>, column_major> C (5,5,12);
control_type<std::complex<double> > Control;
info_type<std::complex<double> > Info;
symbolic_type<std::complex<double> > Symbolic;
numeric_type<std::complex<double> > Numeric;
UMFPACK License:
Your use or distribution of UMFPACK or any modified version of UMFPACK implies that you agree to this License.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use or copy this program, provided that the Copyright, this License, and the Availability of the original version is retained on all copies. User documentation of any code that uses UMFPACK or any modified version of UMFPACK code must cite the Copyright, this License, the Availability note, and "Used by permission." Permission to modify the code and to distribute modified code is granted, provided the Copyright, this License, and the Availability note are retained, and a notice that the code was modified is included. This software was developed with support from the National Science Foundation, and is provided to you free of charge.
Availability:
http://www.cise.ufl.edu/research/sparse/umfpack
Used by permission.