[ublas] How to use the qr decomposition correctly?

Subject: [ublas] How to use the qr decomposition correctly?
From: raymond03152410_at_[hidden]
Date: 2015-03-27 02:05:06

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

raymond03152410_at_[hidden]