uBLAS matrix

There is not such a beast like a matrix in standard C++. When we don't fallback to the raw bidimensional array, a real pain in the neck to work with and to maintain, the common approach is emulating it with a vector of vectors. It could be enough in some cases, but if our problem is not trivial, it is a good idea to look for alternative solutions provided by non standard libraries.

A popular C++ algebraic library is uBLAS, that implements the BLAS interface and it is part of the Boost Libraries.

A standard solution

If we want to manage a very simple case of a 4x4 matrix of integers that has to be just created with all its values set to zero, and then printed, we typically end up with code like this:

const int IM_SIZE_1 = 4; // 1
const int IM_SIZE_2 = 4;
IMatrix matrix(IM_SIZE_1); // 2
std::for_each(matrix.begin(), matrix.end(), [IM_SIZE_2](IVector& v) { v = IVector(IM_SIZE_2); }); // 3
dump(matrix); // 4

1. Specify the matrix dimensions, 4x4.
2. Despite its wishful name, IMatrix is not really a matrix, but a vector, see below for details.
3. Each element of the IMatrix vector has to be initialized by hand, specifying each line dimension. The IVector is passed to the lambda function by reference, so that the assignment to a new instance reflects to the actual element in IMatrix. Otherwise the assignment would have effect only on a local variable that would be discarded at the end of the lambda scope.
4. Now we can call an ad hoc function we have written to dump the matrix to the standard console.

Let's have a look of what IMatrix actually is, and how to implement the dump function:

#include <iostream>
#include <vector>
#include <algorithm>
typedef std::vector<int> IVector;
typedef std::vector<IVector> IMatrix; // 1

void dump(const IVector& vector) // 2
{
    std::cout << "( ";
    std::for_each(vector.begin(), vector.end(), [](int i){ std::cout << i << ' '; }); // 3
    std::cout << ") ";
}

void dump(const IMatrix& matrix)
{
    std::for_each(matrix.begin(), matrix.end(), [](const IVector& v){ dump(v); }); // 4
    std::cout << std::endl;
}

1. IMatrix is nothing more that a synonym for a vector of a vectors of ints.
2. To dump an IMatrix, we should know how to print its IVector components.
3. Pretty easy: for_each element in the vector, from begin() to end(), we call a tiny lambda function that does the dirty job.
4. Now printing an IMatrix is just a matter of calling (2) for each of its elements. The IVector is passed by reference just for performance reasons, we don't want to waste time to create a meaningless local copy.

uBLAS solution

The code above is sort of intuitive, and for such a simple example we won't normally ask for anything better, but let's see in any case how we can refactor it to use the uBLAS matrix concept:

// ...
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
namespace ublas = boost::numeric::ublas; // 1
//...
ublas::matrix<int> matrix(4, 4); // 2
for(unsigned int i = 0; i < matrix.size1(); ++i)
    for(unsigned int j = 0; j < matrix.size2(); ++j)
        matrix(i, j) = 0; // 3
std::cout << matrix << std::endl; // 4

1. Define an alias to the uBLAS namespace, to avoid the nuisance of using the original longish namespace.
2. Instantiate a 4x4 uBLAS matrix object.
3. For efficiency, no implicit initialization of elements is performed by the matrix ctor. The matrix dimensions are available through the getters size1() and size2().
4. Sweetly, the operator "put to" for the standard output console has been oveloaded for the uBLAS matrix.