Re: How to use the qr decomposition correctly?

Hi Raymond,

This looks like an error in boost\numeric\bindings\lapack\ormqr.hpp to me, but I can see no error like this in the version of that file that I have. I have no idea from where you got your version of numeric-bindings. I see that you are using numeric_bindings-v1. A maintained source for that (old) variant is

https://svn.boost.org/svn/boost/sandbox/numeric_bindings-v1

You can grab the sources via

svn co https://svn.boost.org/svn/boost/sandbox/numeric_bindings-v1

The test of ormqr.hpp is at

numeric_bindings-v1/libs/numeric/bindings/lapack/test/ublas_geqrf.cpp

This test is part of the test suite, maybe I should run it with your architecture/compiler to be sure that it should work. If all else fails, you can tell me that information and I will try.

(The newer variant is https://svn.boost.org/svn/boost/sandbox/numeric_bindings, but I guess right now you just want to get your existing code working. Later you may consider upgrading...)

Regards,

Thomas Klimpel

Arbeit: +49 - (0)89 - 330 91 97 78

Mobil : 01523 - 17 157 13

Privat: +49 - (0)89 - 66 00 89 44

Gesendet: Freitag, 27. März 2015 um 07:05 Uhr

Von: "[email protected]" <[email protected]>

An: ublas <[email protected]>

Betreff: [ublas] How to use the qr decomposition correctly?

Hi,everybody.

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.

-- Begin code sample --

//qrf.h

#ifndef QRF_H

#define QRF_H

#if defined(_MSC_VER) && (_MSC_VER >= 1020)

# pragma once

#endif

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

void qr_factorize(

const boost::numeric::ublas::matrix<double>& A,

boost::numeric::ublas::matrix<double>& Q,

boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& R);

#endif

//qrf.cpp

#include "stdafx.h"

#include "qrf.h"

#include <vector>

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

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

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

#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK

#include <boost/numeric/bindings/lapack/geqrf.hpp>

#include <boost/numeric/bindings/lapack/ormqr.hpp>

#include <boost/numeric/bindings/traits/std_vector.hpp>

#include <boost/numeric/bindings/traits/ublas_matrix.hpp>

#undef BOOST_NUMERIC_BINDINGS_USE_CLAPACK

void qr_factorize(

const boost::numeric::ublas::matrix<double>& A,

boost::numeric::ublas::matrix<double>& Q,

boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& R)

{

namespace ublas = boost::numeric::ublas;

BOOST_UBLAS_CHECK(A.size1() >= A.size2(), ublas::external_logic());

std::vector<double> tau(A.size2());

ublas::matrix<double, ublas::column_major> CQ(A);

int info;

info = boost::numeric::bindings::lapack::geqrf(CQ, tau);

BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

ublas::triangular_matrix<double, ublas::upper> CR =

ublas::triangular_adaptor<const ublas::matrix_range<ublas::matrix<double, ublas::column_major> >, ublas::upper>(

ublas::project(CQ, ublas::range(0, A.size2()), ublas::range(0, A.size2()))

);

ublas::identity_matrix<double,ublas::column_major> I(A.size1());

ublas::matrix<double, ublas::column_major> tempQ(I);

info = boost::numeric::bindings::lapack::ormqr('L', 'N', CQ, tau, tempQ, boost::numeric::bindings::lapack::optimal_workspace());

BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

#if BOOST_UBLAS_TYPE_CHECK

BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(CQ, CR), A), ublas::internal_logic());

#endif

ublas::matrix<double> CCQ(CQ);

Q.assign_temporary(CCQ);

R.assign_temporary(CR);

}

// test_LAPACK_QRF.cpp : Defines the entry point for the console application.

//

#include "stdafx.h"

#include "test_LAPACK_QRF.h"

#include "qrf.h"

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

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

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

#ifdef _DEBUG

#ifdef _MSC_VER

# pragma comment(lib, "libf2cd.lib")

# pragma comment(lib, "BLASd.lib")

# pragma comment(lib, "clapackd.lib")

#endif

#else

#ifdef _MSC_VER

# pragma comment(lib, "libf2c.lib")

# pragma comment(lib, "BLAS.lib")

# pragma comment(lib, "clapack.lib")

#endif

#endif

#ifdef _DEBUG

#define new DEBUG_NEW

#endif

// The one and only application object

CWinApp theApp;

using namespace std;

namespace ublas = boost::numeric::ublas;

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])

{

int nRetCode = 0;

HMODULE hModule = ::GetModuleHandle(NULL);

if (hModule != NULL)

{

// initialize MFC and print and error on failure

if (!AfxWinInit(hModule, NULL, ::GetCommandLine(), 0))

{

// TODO: change error code to suit your needs

_tprintf(_T("Fatal Error: MFC initialization failed\n"));

nRetCode = 1;

}

else

{

// TODO: code your application's behavior here.

//ublas::matrix<double,ublas::column_major> m(3, 3);

ublas::matrix<double> m(3, 3);

ublas::matrix<double> Q;

ublas::triangular_matrix<double, ublas::upper> R;

m(0, 0) = 12; m(0, 1) = -51; m(0, 2) = 4;

m(1, 0) = 6; m(1, 1) = 167; m(1, 2) = -68;

m(2, 0) = -4; m(2, 1) = 24; m(2, 2) = -41;

qr_factorize(m, Q, R);

std::cout << Q << std::endl;

std::cout << R << std::endl;

std::cout << ublas::prod(Q, R) << std::endl;

return 0;

}

}

else

{

// TODO: change error code to suit your needs

_tprintf(_T("Fatal Error: GetModuleHandle failed\n"));

nRetCode = 1;

}

return nRetCode;

}

-- End code sample --

Compiling this code gives an error from lapack:

Error 1 error C2589: '(' : illegal token on right side of '::' d:\3rd\include\boost\numeric\bindings\lapack\ormqr.hpp 189

traits::detail::array<value_type> work( std::max(1,n_w*32) );

----------

raymond Chen Shanghai

[email protected]

_______________________________________________ ublas mailing list [email protected] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [email protected]