1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_NO_ASSERTION_CHECKING 11 #define EIGEN_NO_ASSERTION_CHECKING 12 #endif 13 14 static int nb_temporaries; 15 16 #define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { if(size!=0) nb_temporaries++; } 17 18 #include "main.h" 19 #include <Eigen/Cholesky> 20 #include <Eigen/QR> 21 22 #define VERIFY_EVALUATION_COUNT(XPR,N) {\ 23 nb_temporaries = 0; \ 24 XPR; \ 25 if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ 26 VERIFY( (#XPR) && nb_temporaries==N ); \ 27 } 28 29 template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm) 30 { 31 typedef typename MatrixType::Scalar Scalar; 32 typedef typename MatrixType::RealScalar RealScalar; 33 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 34 35 MatrixType symmLo = symm.template triangularView<Lower>(); 36 MatrixType symmUp = symm.template triangularView<Upper>(); 37 MatrixType symmCpy = symm; 38 39 CholType<MatrixType,Lower> chollo(symmLo); 40 CholType<MatrixType,Upper> cholup(symmUp); 41 42 for (int k=0; k<10; ++k) 43 { 44 VectorType vec = VectorType::Random(symm.rows()); 45 RealScalar sigma = internal::random<RealScalar>(); 46 symmCpy += sigma * vec * vec.adjoint(); 47 48 // we are doing some downdates, so it might be the case that the matrix is not SPD anymore 49 CholType<MatrixType,Lower> chol(symmCpy); 50 if(chol.info()!=Success) 51 break; 52 53 chollo.rankUpdate(vec, sigma); 54 VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); 55 56 cholup.rankUpdate(vec, sigma); 57 VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); 58 } 59 } 60 61 template<typename MatrixType> void cholesky(const MatrixType& m) 62 { 63 typedef typename MatrixType::Index Index; 64 /* this test covers the following files: 65 LLT.h LDLT.h 66 */ 67 Index rows = m.rows(); 68 Index cols = m.cols(); 69 70 typedef typename MatrixType::Scalar Scalar; 71 typedef typename NumTraits<Scalar>::Real RealScalar; 72 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; 73 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 74 75 MatrixType a0 = MatrixType::Random(rows,cols); 76 VectorType vecB = VectorType::Random(rows), vecX(rows); 77 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 78 SquareMatrixType symm = a0 * a0.adjoint(); 79 // let's make sure the matrix is not singular or near singular 80 for (int k=0; k<3; ++k) 81 { 82 MatrixType a1 = MatrixType::Random(rows,cols); 83 symm += a1 * a1.adjoint(); 84 } 85 86 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 87 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 88 89 // to test if really Cholesky only uses the upper triangular part, uncomment the following 90 // FIXME: currently that fails !! 91 //symm.template part<StrictlyLower>().setZero(); 92 93 { 94 LLT<SquareMatrixType,Lower> chollo(symmLo); 95 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 96 vecX = chollo.solve(vecB); 97 VERIFY_IS_APPROX(symm * vecX, vecB); 98 matX = chollo.solve(matB); 99 VERIFY_IS_APPROX(symm * matX, matB); 100 101 // test the upper mode 102 LLT<SquareMatrixType,Upper> cholup(symmUp); 103 VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); 104 vecX = cholup.solve(vecB); 105 VERIFY_IS_APPROX(symm * vecX, vecB); 106 matX = cholup.solve(matB); 107 VERIFY_IS_APPROX(symm * matX, matB); 108 109 MatrixType neg = -symmLo; 110 chollo.compute(neg); 111 VERIFY(chollo.info()==NumericalIssue); 112 113 VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU())); 114 VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); 115 VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); 116 VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); 117 } 118 119 // LDLT 120 { 121 int sign = internal::random<int>()%2 ? 1 : -1; 122 123 if(sign == -1) 124 { 125 symm = -symm; // test a negative matrix 126 } 127 128 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 129 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 130 131 LDLT<SquareMatrixType,Lower> ldltlo(symmLo); 132 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 133 vecX = ldltlo.solve(vecB); 134 VERIFY_IS_APPROX(symm * vecX, vecB); 135 matX = ldltlo.solve(matB); 136 VERIFY_IS_APPROX(symm * matX, matB); 137 138 LDLT<SquareMatrixType,Upper> ldltup(symmUp); 139 VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); 140 vecX = ldltup.solve(vecB); 141 VERIFY_IS_APPROX(symm * vecX, vecB); 142 matX = ldltup.solve(matB); 143 VERIFY_IS_APPROX(symm * matX, matB); 144 145 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); 146 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); 147 VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); 148 VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL())); 149 150 if(MatrixType::RowsAtCompileTime==Dynamic) 151 { 152 // note : each inplace permutation requires a small temporary vector (mask) 153 154 // check inplace solve 155 matX = matB; 156 VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0); 157 VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval()); 158 159 160 matX = matB; 161 VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); 162 VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); 163 } 164 165 // restore 166 if(sign == -1) 167 symm = -symm; 168 } 169 170 // test some special use cases of SelfCwiseBinaryOp: 171 MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); 172 m2 = m1; 173 m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); 174 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 175 m2 = m1; 176 m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 177 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 178 m2 = m1; 179 m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); 180 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 181 m2 = m1; 182 m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 183 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 184 185 // update/downdate 186 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) )); 187 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) )); 188 } 189 190 template<typename MatrixType> void cholesky_cplx(const MatrixType& m) 191 { 192 // classic test 193 cholesky(m); 194 195 // test mixing real/scalar types 196 197 typedef typename MatrixType::Index Index; 198 199 Index rows = m.rows(); 200 Index cols = m.cols(); 201 202 typedef typename MatrixType::Scalar Scalar; 203 typedef typename NumTraits<Scalar>::Real RealScalar; 204 typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType; 205 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 206 207 RealMatrixType a0 = RealMatrixType::Random(rows,cols); 208 VectorType vecB = VectorType::Random(rows), vecX(rows); 209 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 210 RealMatrixType symm = a0 * a0.adjoint(); 211 // let's make sure the matrix is not singular or near singular 212 for (int k=0; k<3; ++k) 213 { 214 RealMatrixType a1 = RealMatrixType::Random(rows,cols); 215 symm += a1 * a1.adjoint(); 216 } 217 218 { 219 RealMatrixType symmLo = symm.template triangularView<Lower>(); 220 221 LLT<RealMatrixType,Lower> chollo(symmLo); 222 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 223 vecX = chollo.solve(vecB); 224 VERIFY_IS_APPROX(symm * vecX, vecB); 225 // matX = chollo.solve(matB); 226 // VERIFY_IS_APPROX(symm * matX, matB); 227 } 228 229 // LDLT 230 { 231 int sign = internal::random<int>()%2 ? 1 : -1; 232 233 if(sign == -1) 234 { 235 symm = -symm; // test a negative matrix 236 } 237 238 RealMatrixType symmLo = symm.template triangularView<Lower>(); 239 240 LDLT<RealMatrixType,Lower> ldltlo(symmLo); 241 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 242 vecX = ldltlo.solve(vecB); 243 VERIFY_IS_APPROX(symm * vecX, vecB); 244 // matX = ldltlo.solve(matB); 245 // VERIFY_IS_APPROX(symm * matX, matB); 246 } 247 } 248 249 // regression test for bug 241 250 template<typename MatrixType> void cholesky_bug241(const MatrixType& m) 251 { 252 eigen_assert(m.rows() == 2 && m.cols() == 2); 253 254 typedef typename MatrixType::Scalar Scalar; 255 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 256 257 MatrixType matA; 258 matA << 1, 1, 1, 1; 259 VectorType vecB; 260 vecB << 1, 1; 261 VectorType vecX = matA.ldlt().solve(vecB); 262 VERIFY_IS_APPROX(matA * vecX, vecB); 263 } 264 265 template<typename MatrixType> void cholesky_verify_assert() 266 { 267 MatrixType tmp; 268 269 LLT<MatrixType> llt; 270 VERIFY_RAISES_ASSERT(llt.matrixL()) 271 VERIFY_RAISES_ASSERT(llt.matrixU()) 272 VERIFY_RAISES_ASSERT(llt.solve(tmp)) 273 VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) 274 275 LDLT<MatrixType> ldlt; 276 VERIFY_RAISES_ASSERT(ldlt.matrixL()) 277 VERIFY_RAISES_ASSERT(ldlt.permutationP()) 278 VERIFY_RAISES_ASSERT(ldlt.vectorD()) 279 VERIFY_RAISES_ASSERT(ldlt.isPositive()) 280 VERIFY_RAISES_ASSERT(ldlt.isNegative()) 281 VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) 282 VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) 283 } 284 285 void test_cholesky() 286 { 287 int s; 288 for(int i = 0; i < g_repeat; i++) { 289 CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); 290 CALL_SUBTEST_3( cholesky(Matrix2d()) ); 291 CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); 292 CALL_SUBTEST_4( cholesky(Matrix3f()) ); 293 CALL_SUBTEST_5( cholesky(Matrix4d()) ); 294 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 295 CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); 296 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 297 CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); 298 } 299 300 CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); 301 CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); 302 CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() ); 303 CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() ); 304 305 // Test problem size constructors 306 CALL_SUBTEST_9( LLT<MatrixXf>(10) ); 307 CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); 308 309 EIGEN_UNUSED_VARIABLE(s) 310 } 311
