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
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
test_chol_update(const MatrixType & symm)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
cholesky(const MatrixType & m)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 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
cholesky_cplx(const MatrixType & m)190 template<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
cholesky_bug241(const MatrixType & m)250 template<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
cholesky_verify_assert()265 template<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
test_cholesky()285 void 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