Re: [ublas] lapack bindings for sparse matrices?

From: Vardan Akopian (vakopian_at_[hidden])
Date: 2006-10-09 00:22:08

On 10/6/06, Vardan Akopian <vakopian_at_[hidden]> wrote:
> On 10/6/06, Matthias Troyer <troyer_at_[hidden]> wrote:
> >
> > On Oct 5, 2006, at 11:51 AM, Markus Werle wrote:
> >
> > > Hi!
> > >
> > > It seems that all bindings only work for dense matrices.
> > > But maybe I obverlooked something, so please tell me if it is easy
> > > to get
> > > code similar to the one below compiled with a few changes ...
> > >
> > > regards, Markus
> >
> > Lapack is a library for dense matrices and hence only dense matrix
> > bindings can exist.
> >
>
> Well, banded_matrix is not really dense, but it can be used with
> lapack bindings (I've been using it for long time). I can post some
> code if anyone is interested.
>

This is for Jit (and others who might be interested).
Here is the code I use for solving equations for banded matrices.
Attached is a file called "gbsv.hpp", which is very similar to the
other files in bindings/lapack. I only did the double precision, but
it would be very easy to add single precision (and then maybe it could
be added to the lapack bindings package). Here is how I use it:

#include <boost/numeric/bindings/traits/ublas_banded.hpp>
#include "gbsv.hpp"
#include <vector>

static const char NORMAL = 'N';
static const char TRANSPOSE = 'T';

// solves the equation Ax = B, and puts the solution in B
// A is mutated by this routine
template <typename MatrA, typename MatrB> inline
void InPlaceSolve(MatrA& a, MatrB& b)
{
        using namespace boost::numeric::bindings::lapack;
        std::vector<int> piv(a.size1());
        int ret = gbtrf(a, piv);
        if (ret < 0) {
                CStdString err;
                err.Format("banded::Solve: argument %d in DGBTRF had an illigal value", -ret);
                throw RuntimeError(err);
        }
        if (ret > 0) {
                CStdString err;
                err.Format("banded::Solve: the (%d,%d) diagonal element is 0 after
DGBTRF", ret, ret);
                throw RuntimeError(err);
        }
        
        ret = gbtrs(NORMAL, a, piv, b);
        if (ret < 0) {
                CStdString err;
                err.Format("banded::Solve: argument %d in DGBTRS had an illigal value", -ret);
                throw RuntimeError(err);
        }
}

void test()
{
        // if the matrix has kl lower and ku upper diagonals, then we should
        // allocate kl lower and kl+ku upper diagonals
        size_t sz = 1000, kl = 1, ku = 1;
        ublas::banded_matrix<double> a(sz, sz, kl, kl+ku);
        ublas::vector<double> b(sz);
        // fill values in a and b
        // ...
        InPlaceSolve(a, b);
}

Hope it's useful.
-Vardan