1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <g.gael (at) free.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 #define EIGEN_NO_ASSERTION_CHECKING 11 #include "main.h" 12 #include <Eigen/Cholesky> 13 #include <Eigen/LU> 14 15 #ifdef HAS_GSL 16 #include "gsl_helper.h" 17 #endif 18 19 template<typename MatrixType> void cholesky(const MatrixType& m) 20 { 21 /* this test covers the following files: 22 LLT.h LDLT.h 23 */ 24 int rows = m.rows(); 25 int cols = m.cols(); 26 27 typedef typename MatrixType::Scalar Scalar; 28 typedef typename NumTraits<Scalar>::Real RealScalar; 29 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; 30 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 31 32 MatrixType a0 = MatrixType::Random(rows,cols); 33 VectorType vecB = VectorType::Random(rows), vecX(rows); 34 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 35 SquareMatrixType symm = a0 * a0.adjoint(); 36 // let's make sure the matrix is not singular or near singular 37 MatrixType a1 = MatrixType::Random(rows,cols); 38 symm += a1 * a1.adjoint(); 39 40 #ifdef HAS_GSL 41 if (ei_is_same_type<RealScalar,double>::ret) 42 { 43 typedef GslTraits<Scalar> Gsl; 44 typename Gsl::Matrix gMatA=0, gSymm=0; 45 typename Gsl::Vector gVecB=0, gVecX=0; 46 convert<MatrixType>(symm, gSymm); 47 convert<MatrixType>(symm, gMatA); 48 convert<VectorType>(vecB, gVecB); 49 convert<VectorType>(vecB, gVecX); 50 Gsl::cholesky(gMatA); 51 Gsl::cholesky_solve(gMatA, gVecB, gVecX); 52 VectorType vecX(rows), _vecX, _vecB; 53 convert(gVecX, _vecX); 54 symm.llt().solve(vecB, &vecX); 55 Gsl::prod(gSymm, gVecX, gVecB); 56 convert(gVecB, _vecB); 57 // test gsl itself ! 58 VERIFY_IS_APPROX(vecB, _vecB); 59 VERIFY_IS_APPROX(vecX, _vecX); 60 61 Gsl::free(gMatA); 62 Gsl::free(gSymm); 63 Gsl::free(gVecB); 64 Gsl::free(gVecX); 65 } 66 #endif 67 68 { 69 LDLT<SquareMatrixType> ldlt(symm); 70 VERIFY(ldlt.isPositiveDefinite()); 71 // in eigen3, LDLT is pivoting 72 //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); 73 ldlt.solve(vecB, &vecX); 74 VERIFY_IS_APPROX(symm * vecX, vecB); 75 ldlt.solve(matB, &matX); 76 VERIFY_IS_APPROX(symm * matX, matB); 77 } 78 79 { 80 LLT<SquareMatrixType> chol(symm); 81 VERIFY(chol.isPositiveDefinite()); 82 VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint()); 83 chol.solve(vecB, &vecX); 84 VERIFY_IS_APPROX(symm * vecX, vecB); 85 chol.solve(matB, &matX); 86 VERIFY_IS_APPROX(symm * matX, matB); 87 } 88 89 #if 0 // cholesky is not rank-revealing anyway 90 // test isPositiveDefinite on non definite matrix 91 if (rows>4) 92 { 93 SquareMatrixType symm = a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint(); 94 LLT<SquareMatrixType> chol(symm); 95 VERIFY(!chol.isPositiveDefinite()); 96 LDLT<SquareMatrixType> cholnosqrt(symm); 97 VERIFY(!cholnosqrt.isPositiveDefinite()); 98 } 99 #endif 100 } 101 102 void test_eigen2_cholesky() 103 { 104 for(int i = 0; i < g_repeat; i++) { 105 CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); 106 CALL_SUBTEST_2( cholesky(Matrix2d()) ); 107 CALL_SUBTEST_3( cholesky(Matrix3f()) ); 108 CALL_SUBTEST_4( cholesky(Matrix4d()) ); 109 CALL_SUBTEST_5( cholesky(MatrixXcd(7,7)) ); 110 CALL_SUBTEST_6( cholesky(MatrixXf(17,17)) ); 111 CALL_SUBTEST_7( cholesky(MatrixXd(33,33)) ); 112 } 113 } 114
