Search code examples
c++optimizationboostmatrixsparse-matrix

boost compressed matrix basics


I am confused on how the boost::compressed_matrix works. Suppose I declare the compressed_matrix like this:

boost::numeric::ublas::compressed_matrix<double> T(1000, 1000, 3*1000);

This allocates space for 3*1000 elements in a 1000x1000 matrix. Now how do I give it the locations which are the non-zero elements? When and how are the non-zero elements set? Is it each time I assign an element in the matrix, e.g. B(4,4)=4, it would mark that element as non-zero?

I would really appreciate if you could help me learn this using an example if possible. Some insight into the internal implementation would be great. I want to make sure I don't write programs that are sub-optimal by guess work.

thank you!


Solution

  • compressed matrix has an underlying linear container (unbounded_array by default, but you can make it bounded_array or std::vector if you want), which contains all non-zero elements of the matrix, in row-major (by default) order. That means that whenever you write a new non-zero element to compressed matrix, it is inserted into that underlying array. If you're not filling the matrix in (row-major) order, every insert will be O(n). When you're changing an existing non-zero element, it is simply changed in the underlying array.

    Here's a simple test to see what the underlying structure looks like:

    #include <boost/numeric/ublas/matrix_sparse.hpp>
    #include <boost/numeric/ublas/storage.hpp>
    namespace ublas = boost::numeric::ublas;
    void show_array(const ublas::unbounded_array<double>& a)
    {
        for(size_t i=0; i<a.size(); ++i)
                std::cout << a[i] << ' ';
        std::cout << '\n';
    }
    int main()
    {
        ublas::compressed_matrix<double> m (10, 10, 3 * 10);
        m(0, 5) = 1; // underlying array is {1, 0, 0, 0, ...}
        show_array(m.value_data());
        m(0, 6) = 2; // underlying array is {1, 2, 0, 0, ...}
        show_array(m.value_data());
        m(0, 4) = 3;  // underlying array is {3, 1, 2, 0, ...}
        show_array(m.value_data());
        m(0, 4) = 7;  // underlying array is {7, 1, 2, 0, ...}
        show_array(m.value_data());
    }