Boost logo

Ublas :

From: wedekind (wedekind_at_[hidden])
Date: 2007-07-30 11:32:23


Hello all,

I'm quite new to numerical computing, especially using boost's
lapack-bindings. So I hope to find some generous numerical computing guru
among you to set me on the right path :)

The following example does not work for me. Seemingly because of wrong
dimensions of one of the matrices involved... Could you please check this
example and tell me what's wrong with it? I've played around with it for a
bit and can't find any examples on the web:

-- Begin code sample --

#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/bindings/lapack/geqrf.hpp>
#include <boost/numeric/bindings/lapack/ormqr.hpp>

int main(int argc, char * argv[])
{
   unsigned int m = 4;
   unsigned int n = 3;

   matrix<double> A(m, n);
   vector<double> tau(n);

   A(0,0) = 1; A(0,1) = 2; A(0,2) = 3;
   A(1,0) = 4; A(1,1) = 5; A(1,2) = 6;
   A(2,0) = 7; A(2,1) = 8; A(2,2) = 9;
   A(3,0) = 10; A(3,1) = 11; A(3,2) = 12;

   identity_matrix<double> I(m);
   matrix<double> Q(I);

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

   std::cout << "matrix A is:\n" << A << std::endl << std::endl;
   std::cout << "Q is:\n" << Q << std::endl << std::endl;

   lapack::geqrf(A, tau);

   std::cout << "matrix A is now:\n" << A << std::endl << std::endl;

   lapack::ormqr( 'L', 'T', A, tau, Q, lapack::optimal_workspace());

   std::cout << "matrix A is:\n" << A << std::endl << std::endl;
   std::cout << "Q is:\n" << Q << std::endl << std::endl;

   return 0;
}

-- End of code sample --

Running this code gives an error from lapack:

** On entry to DORMQR parameter number 7 had an illegal value

The parameter 7 to DORMQR is the "leading dimension" of matrix A: in
ormqr.hpp line 168: "traits::leading_dimension (a)".

The matrix A given to ormqr() has a leading dimension of 3, Q has a leading
dimension of 4. Maybe this difference leads to the error message?

So in general, I want to know, how to solve Ax = b using QR decomposition
with ublast and lapack.

Cheers!

Marco