#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;
}