// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // 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 static int nb_temporaries; #define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { if(size!=0) nb_temporaries++; } #include "main.h" #include #include #define VERIFY_EVALUATION_COUNT(XPR,N) {\ nb_temporaries = 0; \ XPR; \ if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ VERIFY( (#XPR) && nb_temporaries==N ); \ } template class CholType> void test_chol_update(const MatrixType& symm) { typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix VectorType; MatrixType symmLo = symm.template triangularView(); MatrixType symmUp = symm.template triangularView(); MatrixType symmCpy = symm; CholType chollo(symmLo); CholType cholup(symmUp); for (int k=0; k<10; ++k) { VectorType vec = VectorType::Random(symm.rows()); RealScalar sigma = internal::random(); symmCpy += sigma * vec * vec.adjoint(); // we are doing some downdates, so it might be the case that the matrix is not SPD anymore CholType chol(symmCpy); if(chol.info()!=Success) break; chollo.rankUpdate(vec, sigma); VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); cholup.rankUpdate(vec, sigma); VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); } } template void cholesky(const MatrixType& m) { typedef typename MatrixType::Index Index; /* this test covers the following files: LLT.h LDLT.h */ Index rows = m.rows(); Index cols = m.cols(); typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef Matrix SquareMatrixType; typedef Matrix VectorType; MatrixType a0 = MatrixType::Random(rows,cols); VectorType vecB = VectorType::Random(rows), vecX(rows); MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); SquareMatrixType symm = a0 * a0.adjoint(); // let's make sure the matrix is not singular or near singular for (int k=0; k<3; ++k) { MatrixType a1 = MatrixType::Random(rows,cols); symm += a1 * a1.adjoint(); } SquareMatrixType symmUp = symm.template triangularView(); SquareMatrixType symmLo = symm.template triangularView(); // to test if really Cholesky only uses the upper triangular part, uncomment the following // FIXME: currently that fails !! //symm.template part().setZero(); { LLT 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); // test the upper mode LLT cholup(symmUp); VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); vecX = cholup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); matX = cholup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); MatrixType neg = -symmLo; chollo.compute(neg); VERIFY(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())); VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); } // LDLT { int sign = internal::random()%2 ? 1 : -1; if(sign == -1) { symm = -symm; // test a negative matrix } SquareMatrixType symmUp = symm.template triangularView(); SquareMatrixType symmLo = symm.template triangularView(); LDLT ldltlo(symmLo); 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); LDLT ldltup(symmUp); VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); vecX = ldltup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); matX = ldltup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL())); if(MatrixType::RowsAtCompileTime==Dynamic) { // note : each inplace permutation requires a small temporary vector (mask) // check inplace solve matX = matB; VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0); VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval()); matX = matB; VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); } // restore if(sign == -1) symm = -symm; } // test some special use cases of SelfCwiseBinaryOp: MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); m2 = m1; m2 += symmLo.template selfadjointView().llt().solve(matB); VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); m2 = m1; m2 -= symmLo.template selfadjointView().llt().solve(matB); VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); m2 = m1; m2.noalias() += symmLo.template selfadjointView().llt().solve(matB); VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView().llt().solve(matB)); m2 = m1; m2.noalias() -= symmLo.template selfadjointView().llt().solve(matB); VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); // update/downdate CALL_SUBTEST(( test_chol_update(symm) )); CALL_SUBTEST(( test_chol_update(symm) )); } template void cholesky_cplx(const MatrixType& m) { // classic test cholesky(m); // test mixing real/scalar types typedef typename MatrixType::Index Index; Index rows = m.rows(); Index cols = m.cols(); typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef Matrix RealMatrixType; typedef Matrix VectorType; RealMatrixType a0 = RealMatrixType::Random(rows,cols); VectorType vecB = VectorType::Random(rows), vecX(rows); MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); RealMatrixType symm = a0 * a0.adjoint(); // let's make sure the matrix is not singular or near singular for (int k=0; k<3; ++k) { RealMatrixType a1 = RealMatrixType::Random(rows,cols); symm += a1 * a1.adjoint(); } { RealMatrixType symmLo = symm.template triangularView(); LLT 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); } // LDLT { int sign = internal::random()%2 ? 1 : -1; if(sign == -1) { symm = -symm; // test a negative matrix } RealMatrixType symmLo = symm.template triangularView(); LDLT ldltlo(symmLo); 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); } } // regression test for bug 241 template void cholesky_bug241(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); typedef typename MatrixType::Scalar Scalar; typedef Matrix VectorType; MatrixType matA; matA << 1, 1, 1, 1; VectorType vecB; vecB << 1, 1; VectorType vecX = matA.ldlt().solve(vecB); VERIFY_IS_APPROX(matA * vecX, vecB); } // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. // This test checks that LDLT reports correctly that matrix is indefinite. // See http://forum.kde.org/viewtopic.php?f=74&t=106942 template void cholesky_indefinite(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); MatrixType mat; mat << 1, 0, 0, -1; LDLT ldlt(mat); VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } template void cholesky_verify_assert() { MatrixType tmp; LLT llt; VERIFY_RAISES_ASSERT(llt.matrixL()) VERIFY_RAISES_ASSERT(llt.matrixU()) VERIFY_RAISES_ASSERT(llt.solve(tmp)) VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) LDLT ldlt; VERIFY_RAISES_ASSERT(ldlt.matrixL()) VERIFY_RAISES_ASSERT(ldlt.permutationP()) 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)) } void test_cholesky() { int s; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( cholesky(Matrix()) ); CALL_SUBTEST_3( cholesky(Matrix2d()) ); CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); CALL_SUBTEST_3( cholesky_indefinite(Matrix2d()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_5( cholesky(Matrix4d()) ); s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); s = internal::random(1,EIGEN_TEST_MAX_SIZE/2); CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); } CALL_SUBTEST_4( cholesky_verify_assert() ); CALL_SUBTEST_7( cholesky_verify_assert() ); CALL_SUBTEST_8( cholesky_verify_assert() ); CALL_SUBTEST_2( cholesky_verify_assert() ); // Test problem size constructors CALL_SUBTEST_9( LLT(10) ); CALL_SUBTEST_9( LDLT(10) ); EIGEN_UNUSED_VARIABLE(s) }