17#include <boost/numeric/ublas/matrix.hpp>
18#include <boost/numeric/ublas/vector.hpp>
19#include <boost/numeric/ublas/io.hpp>
28using namespace boost::numeric::ublas;
38template<
class matrix_type>
41 typedef typename matrix_type::size_type size_type;
44 for (size_type i=1; i<
A.size2() && upper; i++)
46 for (size_type j=0; j<i && upper; j++)
48 upper = fabs(
A(i,j)) < eps;
58 if (
os_) *
os_ <<
"testReflector() begin" << endl;
79 if (
os_) *
os_ <<
"testReflector() end" << endl;
84 if (
os_) *
os_ <<
"testQR() begin" << endl;
90 for (dmatrix::size_type i=0; i<
A.size1(); i++)
92 for (dmatrix::size_type j=0; j<
A.size2(); j++)
94 A(i,j) = ((i==j) || (j == 0));
107 identity_matrix<dmatrix::value_type> eye(Q.size1());
108 diff = prod(Q, herm(Q)) - eye;
112 catch (boost::numeric::ublas::bad_argument ba)
114 if (
os_) *
os_ <<
"exception: " << ba.what() << endl;
117 if (
os_) *
os_ <<
"testQR() end" << endl;
122 if (
os_) *
os_ <<
"testRectangularQR() begin" << endl;
128 for (dmatrix::size_type i=0; i<
A.size1(); i++)
130 for (dmatrix::size_type j=0; j<
A.size2(); j++)
132 A(i,j) = ((i==j) || (j == 0));
137 qr<dmatrix>(
A, Q, R);
145 identity_matrix<dmatrix::value_type> eye(Q.size1());
146 diff = prod(Q, herm(Q)) - eye;
150 catch (boost::numeric::ublas::bad_argument ba)
152 if (
os_) *
os_ <<
"exception: " << ba.what() << endl;
155 if (
os_) *
os_ <<
"testRectangularQR() end" << endl;
160 if (argc>1 && !strcmp(argv[1],
"-v"))
os_ = &cout;
161 if (
os_) *
os_ <<
"qrTest\n";
void diff(const string &filename1, const string &filename2)
KernelTraitsBase< Kernel >::space_type::abscissa_type x
boost::numeric::ublas::vector< double > dvector
boost::numeric::ublas::matrix< double > dmatrix
void Reflector(const vector_type &x, matrix_type &F)
int main(int argc, char **argv)
bool isUpperTriangular(const matrix_type &A, double eps)
#define unit_assert_equal(x, y, epsilon)