[ublas] lapack bindings

Subject: [ublas] lapack bindings
From: David Mebane (mebane_at_[hidden])
Date: 2011-12-13 15:09:28

Hi everyone,

Trying to get the bindings up and working; having trouble.  Wondering if
there's anything wrong with the following code, which was adapted from the
examples:

// solving A * X = B
// A symmetric positive definite
// factor (potrf()) and solve (potrs())

#include <cstddef>
#include <iostream>
#include <boost/numeric/bindings/lapack/computational/potrf.hpp>
#include <boost/numeric/bindings/lapack/computational/potrs.hpp>
#include <boost/numeric/bindings/ublas/matrix.hpp>
#include <boost/numeric/bindings/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>

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

using std::size_t;
using std::cout;
using std::endl;

typedef float real_t;

typedef ublas::matrix<real_t, ublas::column_major> m_t;

typedef ublas::symmetric_adaptor<m_t, ublas::lower> symml_t;

typedef ublas::symmetric_adaptor<m_t, ublas::upper> symmu_t;

int main() {

  // for more descriptive comments see ublas_posv.cc
  cout << endl;

  // symmetric
  cout << "real symmetric\n" << endl;

  size_t n = 5;
  size_t nrhs = 2;
  m_t al (n, n), au (n, n);  // matrices (storage)
  symml_t sal (al);   // symmetric adaptor
  symmu_t sau (au);   // symmetric adaptor
  m_t x (n, nrhs);
  m_t bl (n, nrhs), bu (n, nrhs);  // RHS matrices

  for (unsigned i = 0; i < sal.size1 (); ++ i) {
      for (unsigned j = 0; j <= i; ++ j)
          sal (i, j) = 3 * i + j + 1;
  }

  for (unsigned i = 0; i < sau.size1 (); ++ i) {
      for (unsigned j = 0; j <= i; ++ j)
          sau (i, j) = 3 * i + j + 1;
  }

  cout << "sal = " << sal << "\n";
  cout << endl;
  cout << "al = " << al << "\n";
  cout << endl;

  cout << "sau = " << sau << "\n";
  cout << endl;
  cout << "au = " << au << "\n";
  cout << endl;

  for (int i = 0; i < x.size1(); ++i) {
    x (i, 0) = 1.;
    x (i, 1) = 2.;
  }
  bl = prod (sal, x);
  bu = prod (sau, x);

  cout << "bl = " << bl << "\n";
  cout << endl;
  cout << "bu = " << bu << "\n";
  cout << endl;

  int ierr = lapack::potrf (sal);
  if (!ierr) {
    lapack::potrs (sal, bl);
    cout << "xl = " << bl << "\n";
  }
  cout << endl;

  ierr = lapack::potrf (sau);
  if (!ierr) {
    lapack::potrs (sau, bu);
    cout << "xu = " << bu << "\n";
  }
  cout << endl;

  return 0;

}

It compiles, but the lapack routines fail.  Could it have something to do
with my version of lapack?

Thanks,
-David