diff options
Diffstat (limited to 'test/lu.cpp')
-rw-r--r-- | test/lu.cpp | 93 |
1 files changed, 32 insertions, 61 deletions
diff --git a/test/lu.cpp b/test/lu.cpp index 9787f4d86..1bbadcbf0 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -9,6 +9,7 @@ #include "main.h" #include <Eigen/LU> +#include "solverbase.h" using namespace std; template<typename MatrixType> @@ -18,7 +19,8 @@ typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { template<typename MatrixType> void lu_non_invertible() { - typedef typename MatrixType::Index Index; + STATIC_CHECK(( internal::is_same<typename FullPivLU<MatrixType>::StorageIndex,int>::value )); + typedef typename MatrixType::RealScalar RealScalar; /* this test covers the following files: LU.h @@ -58,6 +60,10 @@ template<typename MatrixType> void lu_non_invertible() // The image of the zero matrix should consist of a single (zero) column vector VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1)); + // The kernel of the zero matrix is the entire space, and thus is an invertible matrix of dimensions cols. + KernelMatrixType kernel = MatrixType::Zero(rows,cols).fullPivLu().kernel(); + VERIFY((kernel.fullPivLu().isInvertible())); + MatrixType m1(rows, cols), m3(rows, cols2); CMatrixType m2(cols, cols2); createRandomPIMatrixOfRank(rank, rows, cols, m1); @@ -87,42 +93,24 @@ template<typename MatrixType> void lu_non_invertible() VERIFY(!lu.isInjective()); VERIFY(!lu.isInvertible()); VERIFY(!lu.isSurjective()); - VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); + VERIFY_IS_MUCH_SMALLER_THAN((m1 * m1kernel), m1); VERIFY(m1image.fullPivLu().rank() == rank); VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image); + check_solverbase<CMatrixType, MatrixType>(m1, lu, rows, cols, cols2); + m2 = CMatrixType::Random(cols,cols2); m3 = m1*m2; m2 = CMatrixType::Random(cols,cols2); // test that the code, which does resize(), may be applied to an xpr m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); - - // test solve with transposed - m3 = MatrixType::Random(rows,cols2); - m2 = m1.transpose()*m3; - m3 = MatrixType::Random(rows,cols2); - lu.template _solve_impl_transposed<false>(m2, m3); - VERIFY_IS_APPROX(m2, m1.transpose()*m3); - m3 = MatrixType::Random(rows,cols2); - m3 = lu.transpose().solve(m2); - VERIFY_IS_APPROX(m2, m1.transpose()*m3); - - // test solve with conjugate transposed - m3 = MatrixType::Random(rows,cols2); - m2 = m1.adjoint()*m3; - m3 = MatrixType::Random(rows,cols2); - lu.template _solve_impl_transposed<true>(m2, m3); - VERIFY_IS_APPROX(m2, m1.adjoint()*m3); - m3 = MatrixType::Random(rows,cols2); - m3 = lu.adjoint().solve(m2); - VERIFY_IS_APPROX(m2, m1.adjoint()*m3); } template<typename MatrixType> void lu_invertible() { /* this test covers the following files: - LU.h + FullPivLU.h */ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; Index size = MatrixType::RowsAtCompileTime; @@ -145,10 +133,12 @@ template<typename MatrixType> void lu_invertible() VERIFY(lu.isSurjective()); VERIFY(lu.isInvertible()); VERIFY(lu.image(m1).fullPivLu().isInvertible()); + + check_solverbase<MatrixType, MatrixType>(m1, lu, size, size, size); + + MatrixType m1_inverse = lu.inverse(); m3 = MatrixType::Random(size,size); m2 = lu.solve(m3); - VERIFY_IS_APPROX(m3, m1*m2); - MatrixType m1_inverse = lu.inverse(); VERIFY_IS_APPROX(m2, m1_inverse*m3); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); @@ -157,64 +147,37 @@ template<typename MatrixType> void lu_invertible() // truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); - // test solve with transposed - lu.template _solve_impl_transposed<false>(m3, m2); - VERIFY_IS_APPROX(m3, m1.transpose()*m2); - m3 = MatrixType::Random(size,size); - m3 = lu.transpose().solve(m2); - VERIFY_IS_APPROX(m2, m1.transpose()*m3); - - // test solve with conjugate transposed - lu.template _solve_impl_transposed<true>(m3, m2); - VERIFY_IS_APPROX(m3, m1.adjoint()*m2); - m3 = MatrixType::Random(size,size); - m3 = lu.adjoint().solve(m2); - VERIFY_IS_APPROX(m2, m1.adjoint()*m3); - // Regression test for Bug 302 MatrixType m4 = MatrixType::Random(size,size); VERIFY_IS_APPROX(lu.solve(m3*m4), lu.solve(m3)*m4); } -template<typename MatrixType> void lu_partial_piv() +template<typename MatrixType> void lu_partial_piv(Index size = MatrixType::ColsAtCompileTime) { /* this test covers the following files: PartialPivLU.h */ - typedef typename MatrixType::Index Index; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - Index size = internal::random<Index>(1,4); MatrixType m1(size, size), m2(size, size), m3(size, size); m1.setRandom(); PartialPivLU<MatrixType> plu(m1); + STATIC_CHECK(( internal::is_same<typename PartialPivLU<MatrixType>::StorageIndex,int>::value )); + VERIFY_IS_APPROX(m1, plu.reconstructedMatrix()); + check_solverbase<MatrixType, MatrixType>(m1, plu, size, size, size); + + MatrixType m1_inverse = plu.inverse(); m3 = MatrixType::Random(size,size); m2 = plu.solve(m3); - VERIFY_IS_APPROX(m3, m1*m2); - MatrixType m1_inverse = plu.inverse(); VERIFY_IS_APPROX(m2, m1_inverse*m3); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); const RealScalar rcond_est = plu.rcond(); // Verify that the estimate is within a factor of 10 of the truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); - - // test solve with transposed - plu.template _solve_impl_transposed<false>(m3, m2); - VERIFY_IS_APPROX(m3, m1.transpose()*m2); - m3 = MatrixType::Random(size,size); - m3 = plu.transpose().solve(m2); - VERIFY_IS_APPROX(m2, m1.transpose()*m3); - - // test solve with conjugate transposed - plu.template _solve_impl_transposed<true>(m3, m2); - VERIFY_IS_APPROX(m3, m1.adjoint()*m2); - m3 = MatrixType::Random(size,size); - m3 = plu.adjoint().solve(m2); - VERIFY_IS_APPROX(m2, m1.adjoint()*m3); } template<typename MatrixType> void lu_verify_assert() @@ -228,6 +191,8 @@ template<typename MatrixType> void lu_verify_assert() VERIFY_RAISES_ASSERT(lu.kernel()) VERIFY_RAISES_ASSERT(lu.image(tmp)) VERIFY_RAISES_ASSERT(lu.solve(tmp)) + VERIFY_RAISES_ASSERT(lu.transpose().solve(tmp)) + VERIFY_RAISES_ASSERT(lu.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(lu.determinant()) VERIFY_RAISES_ASSERT(lu.rank()) VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) @@ -240,19 +205,25 @@ template<typename MatrixType> void lu_verify_assert() VERIFY_RAISES_ASSERT(plu.matrixLU()) VERIFY_RAISES_ASSERT(plu.permutationP()) VERIFY_RAISES_ASSERT(plu.solve(tmp)) + VERIFY_RAISES_ASSERT(plu.transpose().solve(tmp)) + VERIFY_RAISES_ASSERT(plu.adjoint().solve(tmp)) VERIFY_RAISES_ASSERT(plu.determinant()) VERIFY_RAISES_ASSERT(plu.inverse()) } -void test_lu() +EIGEN_DECLARE_TEST(lu) { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() ); CALL_SUBTEST_1( lu_invertible<Matrix3f>() ); CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() ); + CALL_SUBTEST_1( lu_partial_piv<Matrix3f>() ); CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) ); CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) ); + CALL_SUBTEST_2( lu_partial_piv<Matrix2d>() ); + CALL_SUBTEST_2( lu_partial_piv<Matrix4d>() ); + CALL_SUBTEST_2( (lu_partial_piv<Matrix<double,6,6> >()) ); CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() ); CALL_SUBTEST_3( lu_invertible<MatrixXf>() ); @@ -260,7 +231,7 @@ void test_lu() CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() ); CALL_SUBTEST_4( lu_invertible<MatrixXd>() ); - CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() ); + CALL_SUBTEST_4( lu_partial_piv<MatrixXd>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ); CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() ); CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() ); @@ -269,7 +240,7 @@ void test_lu() CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() ); CALL_SUBTEST_6( lu_invertible<MatrixXcd>() ); - CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() ); + CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ); CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() ); CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() )); |