aboutsummaryrefslogtreecommitdiff
path: root/test/lu.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'test/lu.cpp')
-rw-r--r--test/lu.cpp93
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> >() ));
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-24 20:46:56 +0000