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

Subject: Re: [ublas] How to use the qr decomposition correctly?
From: Thomas Klimpel (thomas.klimpel_at_[hidden])
Date: 2015-03-27 05:33:10

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: "[e-mail]" <[e-mail]>
An: ublas <[e-mail]>
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
>
>  
>
> [e-mail]
_______________________________________________ ublas mailing list [e-mail] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [e-mail]