#include <algorithm>
#include<boost/numeric/ublas/io.hpp>
#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/vector.hpp>
#include<boost/numeric/ublas/traits.hpp>
#include<boost/numeric/bindings/ublas.hpp>
#include<boost/numeric/bindings/lapack.hpp>
#include <complex>
#include <cstddef>
#include <iostream>


namespace ublas  = boost::numeric::ublas;
namespace lapack = boost::numeric::bindings::lapack;


namespace detail { namespace /*<unnamed>*/ {

template <typename T>
void test()
{
	typedef T value_type;
	typedef typename ublas::type_traits<value_type>::real_type real_type;
	typedef ublas::matrix<value_type,ublas::column_major> matrix_type;
	typedef ublas::vector<real_type> vector_type;

	const ::std::size_t m(1);
	const ::std::size_t n(3);
	const ::std::size_t k(::std::min(m,n));
	const char jobu('N');
	const char jobvt('N');

	matrix_type A(m,n);
	A(0,0) = value_type(-1); A(0,1) = value_type(0); A(0,2) = value_type(1);

	vector_type s(k);
	matrix_type U(1,1); // JOBU=='N' --> not referenced
	matrix_type VT(1,1); // JOBVT=='N' --> not referenced
	lapack::gesvd(jobu, jobvt, A, s, U, VT);

	::std::cout << "A = " << A << ::std::endl;
	::std::cout << "s = " << s << ::std::endl;
}

}} // Namespace detail::<unnamed>


int main()
{
	std::cout << "Testing sgesvd..." << std::endl;
	detail::test<float>();

	std::cout << "Testing dgesvd..." << std::endl;
	detail::test<double>();

	std::cout << "Testing cgesvd..." << std::endl;
	detail::test< ::std::complex<float> >();

	std::cout << "Testing zgesvd..." << std::endl;
	detail::test< ::std::complex<double> >();

    return 0;
}