eigen2_cholesky.cpp revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
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
19template<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
102void 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