CODEWITHSUNDEEP
c++matrix
Rastapopoulos
Rastapopoulos

Reputation: 507

Working with both dense and sparse matrices

I am writing a C++ function that operates on a matrix, passed as an argument, and would like the code to work with different types of matrices (e.g. Boost sparse matrices, std::vectors of std::vectors). My current approach is to define overloaded methods for basic matrix manipulation with the different types, to provide a unified interface to the different types of matrices, and to define my function as a template function that uses only these overloaded methods 

#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <iostream>

typedef std::vector<double> vec;
typedef std::vector<vec> mat;
typedef boost::numeric::ublas::compressed_matrix<double, boost::numeric::ublas::row_major> spmat;

namespace matrix
{
    inline void set(spmat & input, u_int i, u_int j, double val)
    {
        input(i, j) = val;
    }

    inline void set(mat & input, u_int i, u_int j, double val)
    {
        input[i][j] = val;
    }

    inline u_int size1(const mat & input)
    {
        return input.size();
    }

    inline u_int size2(const mat & input)
    {
        return input[0].size();
    }

    inline u_int size1(const spmat & input)
    {
        return input.size1();
    }

    inline u_int size2(const spmat & input)
    {
        return input.size2();
    }

    inline double get(const spmat & input, u_int i, u_int j)
    {
        return input(i, j);
    }

    inline double get(const mat & input, u_int i, u_int j)
    {
        return input[i][j];
    }
}

For simple tasks, this approach seems to be working. At the moment, however, I am trying to write a function that needs to iterate over all the entries, in the case a dense matrix, or only the nonzero entries, in the case of a sparse matrix. I know how to do this for each case individually, but would like to have only one implementation that works in both cases. What would be the standard way of achieving this?

Upvotes: 3

Views: 261

Answers (1)

Max Langhof
Max Langhof

Reputation: 23691

I have had the same problem. As I was too lazy to write iterators (that would quickly run into limitations), I decided to use a lambda approach by defining specialized functions for each matrix that call a lambda per element:

template<class Func>
void forMatrixEntries(const VecOfVecMatrix& mat, Func&& func)
{
    for (auto& row : mat.getData())
        for (auto& elem : row)
            func(elem); // You could track and pass the indices if you want.
}

template<class Func>
void forMatrixEntries(const CompressedSparseMatrix& mat, Func&& func)
{
    for (auto& elem : mat.getElements())
        func(elem); // You could track and pass the indices if you want.
}

(These could just as well be member functions so they can access internals easier - your choice). You could then use them in a unified way:

template<class Mat>
void scale(const Mat& mat, double factor)
{
    forMatrixEntries(mat, [factor](double& elem) {
        elem *= factor;
    });
}

Only drawback is that the matrix-specialized functions of course need to be in headers (because templates). But I think this approach is not only elegant but also very expressive (you get to give a name to "iterate over matrix entries" instead of convoluted loop syntax, but the loop body stays the same).

Upvotes: 2

Related Questions