
Ublas : 
From: Vardan Akopian (vakopian_at_[hidden])
Date: 20061009 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