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 // to test if really Cholesky only uses the upper triangular part, uncomment the following 87 // FIXME: currently that fails !! 88 //symm.template part<StrictlyLower>().setZero(); 89 90 { 91 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 92 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 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 // test some special use cases of SelfCwiseBinaryOp: 119 MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); 120 m2 = m1; 121 m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); 122 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 123 m2 = m1; 124 m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 125 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 126 m2 = m1; 127 m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); 128 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 129 m2 = m1; 130 m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 131 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 132 } 133 134 // LDLT 135 { 136 int sign = internal::random<int>()%2 ? 1 : -1; 137 138 if(sign == -1) 139 { 140 symm = -symm; // test a negative matrix 141 } 142 143 SquareMatrixType symmUp = symm.template triangularView<Upper>(); 144 SquareMatrixType symmLo = symm.template triangularView<Lower>(); 145 146 LDLT<SquareMatrixType,Lower> ldltlo(symmLo); 147 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 148 vecX = ldltlo.solve(vecB); 149 VERIFY_IS_APPROX(symm * vecX, vecB); 150 matX = ldltlo.solve(matB); 151 VERIFY_IS_APPROX(symm * matX, matB); 152 153 LDLT<SquareMatrixType,Upper> ldltup(symmUp); 154 VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); 155 vecX = ldltup.solve(vecB); 156 VERIFY_IS_APPROX(symm * vecX, vecB); 157 matX = ldltup.solve(matB); 158 VERIFY_IS_APPROX(symm * matX, matB); 159 160 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); 161 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); 162 VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); 163 VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL())); 164 165 if(MatrixType::RowsAtCompileTime==Dynamic) 166 { 167 // note : each inplace permutation requires a small temporary vector (mask) 168 169 // check inplace solve 170 matX = matB; 171 VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0); 172 VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval()); 173 174 175 matX = matB; 176 VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); 177 VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); 178 } 179 180 // restore 181 if(sign == -1) 182 symm = -symm; 183 184 // check matrices coming from linear constraints with Lagrange multipliers 185 if(rows>=3) 186 { 187 SquareMatrixType A = symm; 188 int c = internal::random<int>(0,rows-2); 189 A.bottomRightCorner(c,c).setZero(); 190 // Make sure a solution exists: 191 vecX.setRandom(); 192 vecB = A * vecX; 193 vecX.setZero(); 194 ldltlo.compute(A); 195 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 196 vecX = ldltlo.solve(vecB); 197 VERIFY_IS_APPROX(A * vecX, vecB); 198 } 199 200 // check non-full rank matrices 201 if(rows>=3) 202 { 203 int r = internal::random<int>(1,rows-1); 204 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r); 205 SquareMatrixType A = a * a.adjoint(); 206 // Make sure a solution exists: 207 vecX.setRandom(); 208 vecB = A * vecX; 209 vecX.setZero(); 210 ldltlo.compute(A); 211 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 212 vecX = ldltlo.solve(vecB); 213 VERIFY_IS_APPROX(A * vecX, vecB); 214 } 215 216 // check matrices with a wide spectrum 217 if(rows>=3) 218 { 219 RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); 220 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows); 221 Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows); 222 for(int k=0; k<rows; ++k) 223 d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); 224 SquareMatrixType A = a * d.asDiagonal() * a.adjoint(); 225 // Make sure a solution exists: 226 vecX.setRandom(); 227 vecB = A * vecX; 228 vecX.setZero(); 229 ldltlo.compute(A); 230 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 231 vecX = ldltlo.solve(vecB); 232 VERIFY_IS_APPROX(A * vecX, vecB); 233 } 234 } 235 236 // update/downdate 237 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) )); 238 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) )); 239 } 240 241 template<typename MatrixType> void cholesky_cplx(const MatrixType& m) 242 { 243 // classic test 244 cholesky(m); 245 246 // test mixing real/scalar types 247 248 typedef typename MatrixType::Index Index; 249 250 Index rows = m.rows(); 251 Index cols = m.cols(); 252 253 typedef typename MatrixType::Scalar Scalar; 254 typedef typename NumTraits<Scalar>::Real RealScalar; 255 typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType; 256 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 257 258 RealMatrixType a0 = RealMatrixType::Random(rows,cols); 259 VectorType vecB = VectorType::Random(rows), vecX(rows); 260 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 261 RealMatrixType symm = a0 * a0.adjoint(); 262 // let's make sure the matrix is not singular or near singular 263 for (int k=0; k<3; ++k) 264 { 265 RealMatrixType a1 = RealMatrixType::Random(rows,cols); 266 symm += a1 * a1.adjoint(); 267 } 268 269 { 270 RealMatrixType symmLo = symm.template triangularView<Lower>(); 271 272 LLT<RealMatrixType,Lower> chollo(symmLo); 273 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 274 vecX = chollo.solve(vecB); 275 VERIFY_IS_APPROX(symm * vecX, vecB); 276 // matX = chollo.solve(matB); 277 // VERIFY_IS_APPROX(symm * matX, matB); 278 } 279 280 // LDLT 281 { 282 int sign = internal::random<int>()%2 ? 1 : -1; 283 284 if(sign == -1) 285 { 286 symm = -symm; // test a negative matrix 287 } 288 289 RealMatrixType symmLo = symm.template triangularView<Lower>(); 290 291 LDLT<RealMatrixType,Lower> ldltlo(symmLo); 292 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 293 vecX = ldltlo.solve(vecB); 294 VERIFY_IS_APPROX(symm * vecX, vecB); 295 // matX = ldltlo.solve(matB); 296 // VERIFY_IS_APPROX(symm * matX, matB); 297 } 298 } 299 300 // regression test for bug 241 301 template<typename MatrixType> void cholesky_bug241(const MatrixType& m) 302 { 303 eigen_assert(m.rows() == 2 && m.cols() == 2); 304 305 typedef typename MatrixType::Scalar Scalar; 306 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 307 308 MatrixType matA; 309 matA << 1, 1, 1, 1; 310 VectorType vecB; 311 vecB << 1, 1; 312 VectorType vecX = matA.ldlt().solve(vecB); 313 VERIFY_IS_APPROX(matA * vecX, vecB); 314 } 315 316 // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. 317 // This test checks that LDLT reports correctly that matrix is indefinite. 318 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736 319 template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) 320 { 321 eigen_assert(m.rows() == 2 && m.cols() == 2); 322 MatrixType mat; 323 { 324 mat << 1, 0, 0, -1; 325 LDLT<MatrixType> ldlt(mat); 326 VERIFY(!ldlt.isNegative()); 327 VERIFY(!ldlt.isPositive()); 328 } 329 { 330 mat << 1, 2, 2, 1; 331 LDLT<MatrixType> ldlt(mat); 332 VERIFY(!ldlt.isNegative()); 333 VERIFY(!ldlt.isPositive()); 334 } 335 { 336 mat << 0, 0, 0, 0; 337 LDLT<MatrixType> ldlt(mat); 338 VERIFY(ldlt.isNegative()); 339 VERIFY(ldlt.isPositive()); 340 } 341 { 342 mat << 0, 0, 0, 1; 343 LDLT<MatrixType> ldlt(mat); 344 VERIFY(!ldlt.isNegative()); 345 VERIFY(ldlt.isPositive()); 346 } 347 { 348 mat << -1, 0, 0, 0; 349 LDLT<MatrixType> ldlt(mat); 350 VERIFY(ldlt.isNegative()); 351 VERIFY(!ldlt.isPositive()); 352 } 353 } 354 355 template<typename MatrixType> void cholesky_verify_assert() 356 { 357 MatrixType tmp; 358 359 LLT<MatrixType> llt; 360 VERIFY_RAISES_ASSERT(llt.matrixL()) 361 VERIFY_RAISES_ASSERT(llt.matrixU()) 362 VERIFY_RAISES_ASSERT(llt.solve(tmp)) 363 VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) 364 365 LDLT<MatrixType> ldlt; 366 VERIFY_RAISES_ASSERT(ldlt.matrixL()) 367 VERIFY_RAISES_ASSERT(ldlt.permutationP()) 368 VERIFY_RAISES_ASSERT(ldlt.vectorD()) 369 VERIFY_RAISES_ASSERT(ldlt.isPositive()) 370 VERIFY_RAISES_ASSERT(ldlt.isNegative()) 371 VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) 372 VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) 373 } 374 375 void test_cholesky() 376 { 377 int s = 0; 378 for(int i = 0; i < g_repeat; i++) { 379 CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); 380 CALL_SUBTEST_3( cholesky(Matrix2d()) ); 381 CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); 382 CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); 383 CALL_SUBTEST_4( cholesky(Matrix3f()) ); 384 CALL_SUBTEST_5( cholesky(Matrix4d()) ); 385 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 386 CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); 387 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 388 CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); 389 } 390 391 CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); 392 CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); 393 CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() ); 394 CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() ); 395 396 // Test problem size constructors 397 CALL_SUBTEST_9( LLT<MatrixXf>(10) ); 398 CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); 399 400 TEST_SET_BUT_UNUSED_VARIABLE(s) 401 TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) 402 } 403
