Re: [ublas] tridiagonal_matrix

From: dariomt_at_[hidden]
Date: 2008-01-18 06:00:07

*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;
}