Re: [ublas] lapack bindings

Subject: Re: [ublas] lapack bindings
From: David Mebane (mebane_at_[hidden])
Date: 2011-12-14 17:50:44

(Sorry, forgot to turn off digest before I posted.)

Rutger:

What you're getting there is exactly what I'm coming up with. If the
routines were succeeding, the code would produce output for xl and xu.
This would seem to mean (assuming lapack::potrf works for you normally)
that there's something amiss with the test code itself.

-David

(For reference, the original posts:)

Message: 1
Date: Tue, 13 Dec 2011 15:09:28 -0500
From: David Mebane <mebane_at_[hidden]>
To: ublas_at_[hidden]
Subject: [ublas] lapack bindings
Message-ID:
       <CAA7ZD+9C-
0AMpn9YPgCMKX380S3C=0EKJc6Gm0Sktfbh9AH=Mg_at_[hidden]>
Content-Type: text/plain; charset="iso-8859-1"

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
-------------- next part --------------
HTML attachment scrubbed and removed

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

Message: 2
Date: Wed, 14 Dec 2011 09:57:49 +0100
From: Rutger ter Borg <rutger_at_[hidden]>
To: ublas_at_[hidden]
Subject: Re: [ublas] lapack bindings
Message-ID: <jc9oed$8vj$1_at_[hidden]>
Content-Type: text/plain; charset=UTF-8; format=flowed

On 2011-12-13 21:09, David Mebane wrote:
>
>
> It compiles, but the lapack routines fail.  Could it have something to
> do with my version of lapack?
>
>
> Thanks,
> -David
>

Hello David,

on my machine this code produces the output shown below. I guess it has
something to do with your toolchain/library environment. Could you
provide more information on that?

Cheers,

Rutger

$ g++ test.cpp -llapack
$ ./a.out

real symmetric

sal =
[5,5]((1,4,7,10,13),(4,5,8,11,14),(7,8,9,12,15),(10,11,12,13,16),(13,14,15,16,17))

al =
[5,5]((1,0,0,0,0),(4,5,0,0,0),(7,8,9,0,0),(10,11,12,13,0),(13,14,15,16,17))

sau =
[5,5]((1,4,7,10,13),(4,5,8,11,14),(7,8,9,12,15),(10,11,12,13,16),(13,14,15,16,17))

au =
[5,5]((1,4,7,10,13),(0,5,8,11,14),(0,0,9,12,15),(0,0,0,13,16),(0,0,0,0,17))

bl = [5,2]((35,70),(42,84),(51,102),(62,124),(75,150))

bu = [5,2]((35,70),(42,84),(51,102),(62,124),(75,150))