guwi17_at_[hidden] wrote:

>BTW. I would not expect a big difference between a banded matrix with
>compile time known bandwidth compared to the current implementation. I
>still wait for someone who proves the opposite ;-)
>
>Do you already have a benchmark to measure the speed of banded_matrix?

I am doing a very simple test with the (undocumented?) diagonal_matrix in banded.hpp

A matrix product of two diagonal_matrix should return another diagonal_matrix and should have linear complexity O(size). The same goes for an element-wise vector product. Storage should not be an issue, as the data layout of the diagonal_matrix should match that of the vector.

The test performs one of both, and time measurements for a size of 1E7 double elements shows that the vector product is 10 times faster. The test also uses a matrix_vector_slice to get a view of the diagonal_matrix as a vector and does an element-wise product. This is much faster that the matrix product, but vector product is still twice as fast.

IMO compile time knowledge of the bandwith of a banded_matrix *IS* important, specially for large sizes and small bandwidth (being a large diagonal_matrix the extreme case).

I suppose that implementing a banded_matrix_c with compile time bandwith is not trivial, at least not if optimal traversal of elements is required...



test code:


#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/banded.hpp>

#ifdef NDEBUG
    const int SIZE = 10000000;
#else
    const int SIZE = 10;
#endif

void diagonal()
{
    using namespace boost::numeric::ublas;
    diagonal_matrix<double> d1 ( SIZE );
    diagonal_matrix<double> d2 ( SIZE );
   
    // how to do this with diagonal_adaptor<scalar_matrix> ?
    std::fill(d1.data().begin(), d1.data().end(), 2.0);
    std::fill(d2.data().begin(), d2.data().end(), 10.0);
   
    diagonal_matrix<double> d3 = prod(d1,d2);
}

void diagonal_as_vector()
{
    using namespace boost::numeric::ublas;
    diagonal_matrix<double> d1 ( SIZE );
    diagonal_matrix<double> d2 ( SIZE );
   
    // how to do this with diagonal_adaptor<scalar_matrix> ?
    std::fill(d1.data().begin(), d1.data().end(), 2.0);
    std::fill(d2.data().begin(), d2.data().end(), 10.0);

    matrix_vector_slice< diagonal_matrix<double> > dx1(d1, slice(0,1,SIZE), slice(0,1,SIZE));
    matrix_vector_slice< diagonal_matrix<double> > dx2(d2, slice(0,1,SIZE), slice(0,1,SIZE));
   
    // any vector_expression to diagonal_matrix ?
    vector<double> d3 = element_prod(dx1, dx2);
}

void vectors()
{
    using namespace boost::numeric::ublas;
    vector<double> v1(scalar_vector<double>(SIZE, 2.0));
    vector<double> v2(scalar_vector<double>(SIZE, 10.0));
    vector<double> v3 = element_prod(v1, v2);
}

int main(int argc, char* argv[])
{
    if (argc != 2) return 1;
    if (*argv[1] == 'd') diagonal();
    else if (*argv[1] == 'v') vectors();
    else if (*argv[1] == 'x') diagonal_as_vector();
    return 0;
}