Using uBLAS

From NA-Wiki

Jump to: navigation, search

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

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.
Personal tools