How to use the qr decomposition correctly?

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

How to use the qr decomposition correctly?

raymond03152410@163.com
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


_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: How to use the qr decomposition correctly?

PB-19

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) );


My guess is you're hitting the problem where MSVC's Windows include files have a #define for max, so the compiler is treating the token 'max' as a macro. One solution is to put parethenses around std::max:

(std::max)(1,n_w*32)

Pete

_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: How to use the qr decomposition correctly?

Thomas Klimpel-2
In reply to this post by raymond03152410@163.com
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: "[hidden email]" <[hidden email]>
An: ublas <[hidden email]>
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
 
_______________________________________________ ublas mailing list [hidden email] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [hidden email]

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