aboutsummaryrefslogtreecommitdiff
path: root/test/cholesky.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'test/cholesky.cpp')
-rw-r--r--test/cholesky.cpp87
1 files changed, 55 insertions, 32 deletions
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 8ad5ac639..0b1a7b45b 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -7,18 +7,16 @@
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef EIGEN_NO_ASSERTION_CHECKING
-#define EIGEN_NO_ASSERTION_CHECKING
-#endif
-
#define TEST_ENABLE_TEMPORARY_TRACKING
#include "main.h"
#include <Eigen/Cholesky>
#include <Eigen/QR>
+#include "solverbase.h"
template<typename MatrixType, int UpLo>
typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
+ if(m.cols()==0) return typename MatrixType::RealScalar(0);
MatrixType symm = m.template selfadjointView<UpLo>();
return symm.cwiseAbs().colwise().sum().maxCoeff();
}
@@ -57,7 +55,6 @@ template<typename MatrixType,template <typename,int> class CholType> void test_c
template<typename MatrixType> void cholesky(const MatrixType& m)
{
- typedef typename MatrixType::Index Index;
/* this test covers the following files:
LLT.h LDLT.h
*/
@@ -81,15 +78,17 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
}
{
+ STATIC_CHECK(( internal::is_same<typename LLT<MatrixType,Lower>::StorageIndex,int>::value ));
+ STATIC_CHECK(( internal::is_same<typename LLT<MatrixType,Upper>::StorageIndex,int>::value ));
+
SquareMatrixType symmUp = symm.template triangularView<Upper>();
SquareMatrixType symmLo = symm.template triangularView<Lower>();
LLT<SquareMatrixType,Lower> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
- vecX = chollo.solve(vecB);
- VERIFY_IS_APPROX(symm * vecX, vecB);
- matX = chollo.solve(matB);
- VERIFY_IS_APPROX(symm * matX, matB);
+
+ check_solverbase<VectorType, VectorType>(symm, chollo, rows, rows, 1);
+ check_solverbase<MatrixType, MatrixType>(symm, chollo, rows, cols, rows);
const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols));
RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) /
@@ -97,7 +96,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
RealScalar rcond_est = chollo.rcond();
// Verify that the estimated condition number is within a factor of 10 of the
// truth.
- VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+ VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
// test the upper mode
LLT<SquareMatrixType,Upper> cholup(symmUp);
@@ -113,12 +112,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
rcond_est = cholup.rcond();
- VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+ VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
MatrixType neg = -symmLo;
chollo.compute(neg);
- VERIFY(chollo.info()==NumericalIssue);
+ VERIFY(neg.size()==0 || chollo.info()==NumericalIssue);
VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU()));
VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
@@ -143,6 +142,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
// LDLT
{
+ STATIC_CHECK(( internal::is_same<typename LDLT<MatrixType,Lower>::StorageIndex,int>::value ));
+ STATIC_CHECK(( internal::is_same<typename LDLT<MatrixType,Upper>::StorageIndex,int>::value ));
+
int sign = internal::random<int>()%2 ? 1 : -1;
if(sign == -1)
@@ -156,10 +158,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
VERIFY(ldltlo.info()==Success);
VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
- vecX = ldltlo.solve(vecB);
- VERIFY_IS_APPROX(symm * vecX, vecB);
- matX = ldltlo.solve(matB);
- VERIFY_IS_APPROX(symm * matX, matB);
+
+ check_solverbase<VectorType, VectorType>(symm, ldltlo, rows, rows, 1);
+ check_solverbase<MatrixType, MatrixType>(symm, ldltlo, rows, cols, rows);
const MatrixType symmLo_inverse = ldltlo.solve(MatrixType::Identity(rows,cols));
RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) /
@@ -167,7 +168,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
RealScalar rcond_est = ldltlo.rcond();
// Verify that the estimated condition number is within a factor of 10 of the
// truth.
- VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+ VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
LDLT<SquareMatrixType,Upper> ldltup(symmUp);
@@ -184,7 +185,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
rcond_est = ldltup.rcond();
- VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+ VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU()));
VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL()));
@@ -289,8 +290,6 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
// test mixing real/scalar types
- typedef typename MatrixType::Index Index;
-
Index rows = m.rows();
Index cols = m.cols();
@@ -315,10 +314,9 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
LLT<RealMatrixType,Lower> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
- vecX = chollo.solve(vecB);
- VERIFY_IS_APPROX(symm * vecX, vecB);
-// matX = chollo.solve(matB);
-// VERIFY_IS_APPROX(symm * matX, matB);
+
+ check_solverbase<VectorType, VectorType>(symm, chollo, rows, rows, 1);
+ //check_solverbase<MatrixType, MatrixType>(symm, chollo, rows, cols, rows);
}
// LDLT
@@ -335,10 +333,9 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
LDLT<RealMatrixType,Lower> ldltlo(symmLo);
VERIFY(ldltlo.info()==Success);
VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
- vecX = ldltlo.solve(vecB);
- VERIFY_IS_APPROX(symm * vecX, vecB);
-// matX = ldltlo.solve(matB);
-// VERIFY_IS_APPROX(symm * matX, matB);
+
+ check_solverbase<VectorType, VectorType>(symm, ldltlo, rows, rows, 1);
+ //check_solverbase<MatrixType, MatrixType>(symm, ldltlo, rows, cols, rows);
}
}
@@ -373,6 +370,7 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
VERIFY(ldlt.info()==Success);
VERIFY(!ldlt.isNegative());
VERIFY(!ldlt.isPositive());
+ VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
}
{
mat << 1, 2, 2, 1;
@@ -380,6 +378,7 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
VERIFY(ldlt.info()==Success);
VERIFY(!ldlt.isNegative());
VERIFY(!ldlt.isPositive());
+ VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
}
{
mat << 0, 0, 0, 0;
@@ -387,6 +386,7 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
VERIFY(ldlt.info()==Success);
VERIFY(ldlt.isNegative());
VERIFY(ldlt.isPositive());
+ VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
}
{
mat << 0, 0, 0, 1;
@@ -394,6 +394,7 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
VERIFY(ldlt.info()==Success);
VERIFY(!ldlt.isNegative());
VERIFY(ldlt.isPositive());
+ VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
}
{
mat << -1, 0, 0, 0;
@@ -401,6 +402,7 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
VERIFY(ldlt.info()==Success);
VERIFY(ldlt.isNegative());
VERIFY(!ldlt.isPositive());
+ VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
}
}
@@ -452,6 +454,18 @@ void cholesky_faillure_cases()
VERIFY(ldlt.info()==NumericalIssue);
VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
}
+
+ // bug 1479
+ {
+ mat.resize(4,4);
+ mat << 1, 2, 0, 1,
+ 2, 4, 0, 2,
+ 0, 0, 0, 1,
+ 1, 2, 1, 1;
+ ldlt.compute(mat);
+ VERIFY(ldlt.info()==NumericalIssue);
+ VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
+ }
}
template<typename MatrixType> void cholesky_verify_assert()
@@ -462,19 +476,23 @@ template<typename MatrixType> void cholesky_verify_assert()
VERIFY_RAISES_ASSERT(llt.matrixL())
VERIFY_RAISES_ASSERT(llt.matrixU())
VERIFY_RAISES_ASSERT(llt.solve(tmp))
- VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))
+ VERIFY_RAISES_ASSERT(llt.transpose().solve(tmp))
+ VERIFY_RAISES_ASSERT(llt.adjoint().solve(tmp))
+ VERIFY_RAISES_ASSERT(llt.solveInPlace(tmp))
LDLT<MatrixType> ldlt;
VERIFY_RAISES_ASSERT(ldlt.matrixL())
- VERIFY_RAISES_ASSERT(ldlt.permutationP())
+ VERIFY_RAISES_ASSERT(ldlt.transpositionsP())
VERIFY_RAISES_ASSERT(ldlt.vectorD())
VERIFY_RAISES_ASSERT(ldlt.isPositive())
VERIFY_RAISES_ASSERT(ldlt.isNegative())
VERIFY_RAISES_ASSERT(ldlt.solve(tmp))
- VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp))
+ VERIFY_RAISES_ASSERT(ldlt.transpose().solve(tmp))
+ VERIFY_RAISES_ASSERT(ldlt.adjoint().solve(tmp))
+ VERIFY_RAISES_ASSERT(ldlt.solveInPlace(tmp))
}
-void test_cholesky()
+EIGEN_DECLARE_TEST(cholesky)
{
int s = 0;
for(int i = 0; i < g_repeat; i++) {
@@ -493,6 +511,11 @@ void test_cholesky()
CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
TEST_SET_BUT_UNUSED_VARIABLE(s)
}
+ // empty matrix, regression test for Bug 785:
+ CALL_SUBTEST_2( cholesky(MatrixXd(0,0)) );
+
+ // This does not work yet:
+ // CALL_SUBTEST_2( cholesky(Matrix<double,0,0>()) );
CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-20 16:08:15 +0000