Re: [ublas] Error in ublas_gbsv.cpp?

Subject: Re: [ublas] Error in ublas_gbsv.cpp?
From: Richard Stanton (stanton_at_[hidden])
Date: 2011-03-25 11:51:06

I apologize for writing the wrong file name earlier . I'm talking about ublas_gbsv.cpp, which uses banded matrices.

From: Richard Stanton
Sent: Friday, March 25, 2011 8:44 AM
To: 'ublas_at_[hidden]'
Subject: Error in ublas_gesv.cpp?

I just downloaded the latest boost numeric bindings from the boost sandbox, and tried to compile and run the sample program ublas_gbsv.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_gbsv.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;
}