I just downloaded the latest boost numeric bindings from the boost sandbox, and tried to compile and run the sample program ublas_gesv.cpp, adding some lines to print out the solutions obtained to the original equation A x = b.

If I make the matrix small, e.g., size below 10, things seem to work fine – when I take the original A matrix and multiply by the proposed solution, I do indeed get the original b vector (note – I make copies of both A and b before testing this). However, at size 20, for example, here’s what I get for the double version of the code:

x =   [20](-1.36885e+14,-0.015625,1.36885e+14,-2.73771e+14,4.10656e+14,-5.47541e+14,6.84426e+14,-8.21312e+14,9.58197e+14,-1.09508e+15,1.23197e+15,-1.36885e+15,1.50574e+15,-1.64262e+15,1.77951e+15,-1.91639e+15,2.05324e+15,-2.18962e+15,2.31842e+15,-2.31842e+15)

B =   [20](0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)

Prod = [20](-0.015625,1,2,3.0625,4,5,6.125,7,7.625,9.5,10,10.25,10.25,13,13.25,15.75,15,21.5,11.5,19.5)

The last line should just repeat the line above, if everything went according to plan.

I tried this on two different machines (Windows/cygwin and OS X Snow Leopard), both running gcc, with the same results. In case this matters, I linked the program (on my Windows machine) with the command

g++ -o ublas_gbsv ublas_gbsv.o -L/RHS/local/lib -L../Clib/lib -llapack -lf77blas -lcblas -latlas -lgfortran  

In case I did something stupid in printing out the results, below is my (slightly) edited version of ublas_gesv.cpp, including extra lines to print the results. (Note that I also added one extra line, line 71, to set the uninitialized superdiagonal to zero).

Thanks.

Richard Stanton

--------------------

#include <iostream>

#include <boost/numeric/bindings/std/vector.hpp>

#include <boost/numeric/bindings/ublas/banded.hpp>

#include <boost/numeric/bindings/ublas/vector.hpp>

#include <boost/numeric/bindings/bandwidth.hpp>

#include <boost/numeric/bindings/lapack/computational/gbtrf.hpp>

#include <boost/numeric/bindings/lapack/computational/gbtrs.hpp>

#include <boost/numeric/ublas/io.hpp>

#include <vector>

#include <stdexcept>

namespace ublas = boost::numeric::ublas;

namespace lapack = boost::numeric::bindings::lapack;

namespace bindings = boost::numeric::bindings;

using namespace std;

// solves the equation Ax = B, and puts the solution in B

// A is mutated by this routine

template <typename MatrA, typename MatrB>

void InPlaceSolve(MatrA& a, MatrB& b)

{

  // if the matrix has kl lower and ku upper diagonals, then we should have

  // allocated kl lower and kl+ku upper diagonals

  fortran_int_t const kl = bindings::bandwidth_lower (a);

  fortran_int_t const ku = bindings::bandwidth_upper (a) - kl;

  std::vector<fortran_int_t> piv(a.size1());

  int ret = lapack::gbtrf(/*kl, ku, */a, piv);

 if (ret < 0) {

    //CStdString err;

    //err.Format("banded::Solve: argument %d in DGBTRF had an illegal value", -ret);

    //throw RuntimeError(err);

    throw std::runtime_error("banded::Solve: argument %d in DGBTRF had an illegal value");

  }

  if (ret > 0) {

    //CStdString err;

    //err.Format("banded::Solve: the (%d,%d) diagonal element is 0 after DGBTRF", ret, ret);

    //throw RuntimeError(err);

    throw std::runtime_error("banded::Solve: the (%d,%d) diagonal element is 0 after DGBTRF");

  }

  ret = lapack::gbtrs(/*kl, ku, */a, piv, b);

  if (ret < 0) {

    //CStdString err;

    //err.Format("banded::Solve: argument %d in DGBTRS had an illegal value", -ret);

    //throw RuntimeError(err);

    throw std::runtime_error("banded::Solve: argument %d in DGBTRS had an illegal value");

  }

}

template<typename T>

void do_typename()

{

  using namespace boost::numeric::ublas;

  // if the matrix has kl lower and ku upper diagonals, then we should

  // allocate kl lower and kl+ku upper diagonals

  size_t sz = 20, kl = 1, ku = 1;

  ublas::banded_matrix<T> a(sz, sz, kl, kl+ku);

  ublas::vector<T> b(sz);

  // fill values in a and b

  for (size_t i = 0; i < sz; ++i) {

    a(i,i) = i;

    b(i) = i;

  }

  for (size_t i = 1; i < sz; ++i) {

    a(i,i-1) = i;

    a(i-1,i) = 1;

    // RHS addition, in case the uninitialized elements of a cause a problem

    if (i<sz-1) a(i-1,i+1)=0;

  }

  ublas::banded_matrix<T> a1(a);

  ublas::vector<T> b1(b);

  InPlaceSolve(a, b);

  ublas::vector<T> prodVec(prod(a1, b));

  cout << " A =   " << a1 << endl;

  cout << " x =   " << b << endl;

  cout << " B =   " << b1 << endl;

  cout << "Prod = " << prodVec << endl;

}

int main()

{

  do_typename<float>();

  do_typename<double>();

  do_typename<std::complex<float> >();

  do_typename<std::complex<double> >();

  return 0;

}