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@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
14static 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
29template<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
61template<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
190template<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
250template<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
265template<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
285void 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