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@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
cholesky(const MatrixType & m)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
test_eigen2_cholesky()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