# Using uBLAS

### From NA-Wiki

## Contents |

## Competing with Fortran - finding a good matrix/vector library for C++

Fortran does vectors and matrices well. It also has a very well established set of subroutines for matrix/vector operations (BLAS) and for dense linear systems (LAPACK). So in order to get serious about C++ as an alternative for numerics, we really have to choose carefully the libraries we use for matrices and linear systems. What are we looking for?

- An alternative to the classical C pointers-to-pointers representation
- A library that is actually used in serious projects (it must be robust and dependable in the long run).
- A library that provides a standardized syntax (e.g. STL-compliant)
- High performance matrix/vector operations
- Classes with a reasonable level of abstraction (e.g. things like vector.length())

Alas, to get all of this is unrealistic. For instance, there is only one alternative to Fortran-BLAS in terms of performance: CBLAS in ATLAS or Intel MKL. But this means abandoning all abstraction and going back to the eighties. For projects where extreme performance is valuable - but not critical - some abstraction will probably make the coding more enjoyable.

One very attractive candidate is uBLAS. It is a part of BOOST, a set of high quality libraries for C++.

Another interesting candidate is MTL4

## Why uBLAS?

- uBLAS is a part of BOOST, which is referred to as "...one of the most highly regarded and expertly designed C++ library projects in the world" by Sutter and Alexandrescu in "C++ Coding Standards".
- uBLAS is used by serious and complicated software projects, such as DOLFIN.
- uBLAS extends the concepts of the STL vector class in a suitable way for numerics.
- uBLAS has clean and good interfaces.
- Good computational performance for operations on matrices and vectors (but generally not competitive with highly tuned BLAS such as Goto, ATLAS or MKL).
- Bindings to ATLAS, when performance really counts.
- Sparse storage layouts: banded, compressed, coordinate and mapped.
- Expressions on matrices like in Fortran or MATLAB.
- Efficient fixed size arrays, as well as dynamic arrays.

## Links

- uBLAS Documentation
- Overview
- Overview of operations (don't miss reading this)

## Examples

There are some examples online, at the Effective UBLAS wiki. Also, there are some examples at Gunter Winkler's homepage.

### Expressions at work

Download code

uBLAS provides interfaces for manipulating matrices and vectors directly, much as in MATLAB. However, beware of nested multiplications, since the complexity will grow unless a temporary is used. The following code snippet illustrates a small subset of the possible expressions and member functions:

#include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/io.hpp> int main () { using namespace boost::numeric::ublas; unsigned i; matrix<double> mat(6,6); vector<double> v (6); for (i = 0; i < v.size(); ++ i) v (i) = i; vector<double> w = 3.0*v; // multiplication by scalar vector<double> x = w + 2.543*w; // addition double c = norm_2(x); // euclidean norm c = norm_inf(x); // max-norm c = inner_prod(x,w); // inner product for(i = 0; i<mat.size1(); i++) mat(i,i) = 2.0; vector<double> ax = prod(mat,x); // matrix/vector product matrix<double> mat2 = prod(mat,mat) // matrix/matrix product x += v*2.0; // x <- x + 2*v vector<double> xx = x; // copy in memory x -= v*2.0; // x <- x - 2*v v = w + 2.0*v; }

### CG with sparse matrices

Download code

This little example implements the conjugate gradient method for solving linear systems. There are a few notable things here:

- The actual matrix class is hidden in the matrix
`typedef`. This means that the CG-routine works without modification for all storage classes in uBLAS. For instance, get compressed row storage by including`matrix_sparse.hpp`and defining`Matrix`as`compressed_matrix<double, ublas::row_major>`or get banded storage by including`banded.hpp`and defining`Matrix`as`banded_matrix<double>`. - All operations in the CG loop are written as operations on vectors directly.

#include <boost/numeric/ublas/vector.hpp> //#include <boost/numeric/ublas/matrix.hpp> //#include <boost/numeric/ublas/matrix_sparse.hpp> #include <boost/numeric/ublas/banded.hpp> #include <boost/numeric/ublas/io.hpp> using namespace std; namespace ublas = boost::numeric::ublas; // dense matrix storage (don't want this): //typedef ublas::matrix<double> Matrix; // sparse matrix storage: typedef ublas::banded_matrix<double> Matrix; //typedef ublas::compressed_matrix<double, ublas::row_major> Matrix; //typedef ublas::compressed_matrix<double, ublas::column_major> Matrix; // dense vector storage typedef ublas::vector<double> Vector; bool cg_solve(const Matrix &A, const Vector &b, Vector &x, double tol) { int niter; bool conv = false; double alpha, beta, residn; Vector resid = b - prod(A,x); // residual Vector d = resid; // search direction Vector resid_old; Vector temp; // CG loop for(niter = 1; niter <= N; niter++) { temp = prod(A,d); alpha = inner_prod(resid,resid)/inner_prod(d,temp); x += (d*alpha); resid_old = resid; resid -= (temp*alpha); residn = norm_2(resid); if(residn < tol){ conv = true; break; } beta = inner_prod(resid,resid)/inner_prod(resid_old,resid_old); d = resid + d*beta; } return conv; }

### Fixed size vectors and matrices

(soon)

Bounded

### Subrange indexing and assignment - FD with no loops

Using `project`, `row`, `column`, `range`, `slice`, `subrange` and `subslice` it is possible to both access and assign to index ranges of a matrix directly (as in MATLAB)! This is documented here. As an example, consider this function that numerically approximates the derivative of a function (defined on an equidistant grid):

#include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/vector_proxy.hpp> using namespace std; namespace ublas = boost::numeric::ublas; typedef ublas::vector<double> Vector; Vector dfdx(const Vector& f, double h) { const unsigned N = f.size(); Vector d(N); // central diff in interior subrange(d,1,N-1) = (subrange(f,2,N)-subrange(f,0,N-2))/(2*h); // unsymetric at the ends d(0) = (-f(2)+4*f(1)-3*f(0))/(2*h); d(N-1) = (3*f(N-1)-4*f(N-2)+f(N-3))/(2*h); return d; }

Note that the subrange convention is a somewhat counterintuitive: The range `range r(0,N)` gives the indices `0,1,2,...,N-1`.

## Performance tips

(soon)

- Use compiler optimization and define
`NDEBUG`!

- There are performance-related tips scattered around the official documentation site, notably on the Effective UBLAS wiki.

- There are performance tips on Gunter Winkler's homepage.

- Use ATLAS bindings when performance is critical.