How to use the qr decomposition correctly?

"[email protected]" <[email protected]>
Newsgroups gmane.comp.lib.boost.ublas
Message-ID <[email protected]>
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]
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.