1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <[email protected]>
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 TEST_ENABLE_TEMPORARY_TRACKING
11
12 #include "main.h"
13 #include <Eigen/Cholesky>
14 #include <Eigen/QR>
15 #include "solverbase.h"
16
17 template<typename MatrixType, int UpLo>
matrix_l1_norm(const MatrixType & m)18 typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
19 if(m.cols()==0) return typename MatrixType::RealScalar(0);
20 MatrixType symm = m.template selfadjointView<UpLo>();
21 return symm.cwiseAbs().colwise().sum().maxCoeff();
22 }
23
test_chol_update(const MatrixType & symm)24 template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
25 {
26 typedef typename MatrixType::Scalar Scalar;
27 typedef typename MatrixType::RealScalar RealScalar;
28 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
29
30 MatrixType symmLo = symm.template triangularView<Lower>();
31 MatrixType symmUp = symm.template triangularView<Upper>();
32 MatrixType symmCpy = symm;
33
34 CholType<MatrixType,Lower> chollo(symmLo);
35 CholType<MatrixType,Upper> cholup(symmUp);
36
37 for (int k=0; k<10; ++k)
38 {
39 VectorType vec = VectorType::Random(symm.rows());
40 RealScalar sigma = internal::random<RealScalar>();
41 symmCpy += sigma * vec * vec.adjoint();
42
43 // we are doing some downdates, so it might be the case that the matrix is not SPD anymore
44 CholType<MatrixType,Lower> chol(symmCpy);
45 if(chol.info()!=Success)
46 break;
47
48 chollo.rankUpdate(vec, sigma);
49 VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix());
50
51 cholup.rankUpdate(vec, sigma);
52 VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix());
53 }
54 }
55
cholesky(const MatrixType & m)56 template<typename MatrixType> void cholesky(const MatrixType& m)
57 {
58 /* this test covers the following files:
59 LLT.h LDLT.h
60 */
61 Index rows = m.rows();
62 Index cols = m.cols();
63
64 typedef typename MatrixType::Scalar Scalar;
65 typedef typename NumTraits<Scalar>::Real RealScalar;
66 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
67 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
68
69 MatrixType a0 = MatrixType::Random(rows,cols);
70 VectorType vecB = VectorType::Random(rows), vecX(rows);
71 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
72 SquareMatrixType symm = a0 * a0.adjoint();
73 // let's make sure the matrix is not singular or near singular
74 for (int k=0; k<3; ++k)
75 {
76 MatrixType a1 = MatrixType::Random(rows,cols);
77 symm += a1 * a1.adjoint();
78 }
79
80 {
81 STATIC_CHECK(( internal::is_same<typename LLT<MatrixType,Lower>::StorageIndex,int>::value ));
82 STATIC_CHECK(( internal::is_same<typename LLT<MatrixType,Upper>::StorageIndex,int>::value ));
83
84 SquareMatrixType symmUp = symm.template triangularView<Upper>();
85 SquareMatrixType symmLo = symm.template triangularView<Lower>();
86
87 LLT<SquareMatrixType,Lower> chollo(symmLo);
88 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
89
90 check_solverbase<VectorType, VectorType>(symm, chollo, rows, rows, 1);
91 check_solverbase<MatrixType, MatrixType>(symm, chollo, rows, cols, rows);
92
93 const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols));
94 RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) /
95 matrix_l1_norm<MatrixType, Lower>(symmLo_inverse);
96 RealScalar rcond_est = chollo.rcond();
97 // Verify that the estimated condition number is within a factor of 10 of the
98 // truth.
99 VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
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 // Verify that the estimated condition number is within a factor of 10 of the
110 // truth.
111 const MatrixType symmUp_inverse = cholup.solve(MatrixType::Identity(rows,cols));
112 rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
113 matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
114 rcond_est = cholup.rcond();
115 VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
116
117
118 MatrixType neg = -symmLo;
119 chollo.compute(neg);
120 VERIFY(neg.size()==0 || chollo.info()==NumericalIssue);
121
122 VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU()));
123 VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
124 VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU()));
125 VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL()));
126
127 // test some special use cases of SelfCwiseBinaryOp:
128 MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
129 m2 = m1;
130 m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB);
131 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
132 m2 = m1;
133 m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
134 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
135 m2 = m1;
136 m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
137 VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
138 m2 = m1;
139 m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
140 VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
141 }
142
143 // LDLT
144 {
145 STATIC_CHECK(( internal::is_same<typename LDLT<MatrixType,Lower>::StorageIndex,int>::value ));
146 STATIC_CHECK(( internal::is_same<typename LDLT<MatrixType,Upper>::StorageIndex,int>::value ));
147
148 int sign = internal::random<int>()%2 ? 1 : -1;
149
150 if(sign == -1)
151 {
152 symm = -symm; // test a negative matrix
153 }
154
155 SquareMatrixType symmUp = symm.template triangularView<Upper>();
156 SquareMatrixType symmLo = symm.template triangularView<Lower>();
157
158 LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
159 VERIFY(ldltlo.info()==Success);
160 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
161
162 check_solverbase<VectorType, VectorType>(symm, ldltlo, rows, rows, 1);
163 check_solverbase<MatrixType, MatrixType>(symm, ldltlo, rows, cols, rows);
164
165 const MatrixType symmLo_inverse = ldltlo.solve(MatrixType::Identity(rows,cols));
166 RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) /
167 matrix_l1_norm<MatrixType, Lower>(symmLo_inverse);
168 RealScalar rcond_est = ldltlo.rcond();
169 // Verify that the estimated condition number is within a factor of 10 of the
170 // truth.
171 VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
172
173
174 LDLT<SquareMatrixType,Upper> ldltup(symmUp);
175 VERIFY(ldltup.info()==Success);
176 VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix());
177 vecX = ldltup.solve(vecB);
178 VERIFY_IS_APPROX(symm * vecX, vecB);
179 matX = ldltup.solve(matB);
180 VERIFY_IS_APPROX(symm * matX, matB);
181
182 // Verify that the estimated condition number is within a factor of 10 of the
183 // truth.
184 const MatrixType symmUp_inverse = ldltup.solve(MatrixType::Identity(rows,cols));
185 rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
186 matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
187 rcond_est = ldltup.rcond();
188 VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
189
190 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU()));
191 VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL()));
192 VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU()));
193 VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL()));
194
195 if(MatrixType::RowsAtCompileTime==Dynamic)
196 {
197 // note : each inplace permutation requires a small temporary vector (mask)
198
199 // check inplace solve
200 matX = matB;
201 VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
202 VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
203
204
205 matX = matB;
206 VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
207 VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
208 }
209
210 // restore
211 if(sign == -1)
212 symm = -symm;
213
214 // check matrices coming from linear constraints with Lagrange multipliers
215 if(rows>=3)
216 {
217 SquareMatrixType A = symm;
218 Index c = internal::random<Index>(0,rows-2);
219 A.bottomRightCorner(c,c).setZero();
220 // Make sure a solution exists:
221 vecX.setRandom();
222 vecB = A * vecX;
223 vecX.setZero();
224 ldltlo.compute(A);
225 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
226 vecX = ldltlo.solve(vecB);
227 VERIFY_IS_APPROX(A * vecX, vecB);
228 }
229
230 // check non-full rank matrices
231 if(rows>=3)
232 {
233 Index r = internal::random<Index>(1,rows-1);
234 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r);
235 SquareMatrixType A = a * a.adjoint();
236 // Make sure a solution exists:
237 vecX.setRandom();
238 vecB = A * vecX;
239 vecX.setZero();
240 ldltlo.compute(A);
241 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
242 vecX = ldltlo.solve(vecB);
243 VERIFY_IS_APPROX(A * vecX, vecB);
244 }
245
246 // check matrices with a wide spectrum
247 if(rows>=3)
248 {
249 using std::pow;
250 using std::sqrt;
251 RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8);
252 Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows);
253 Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows);
254 for(Index k=0; k<rows; ++k)
255 d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s));
256 SquareMatrixType A = a * d.asDiagonal() * a.adjoint();
257 // Make sure a solution exists:
258 vecX.setRandom();
259 vecB = A * vecX;
260 vecX.setZero();
261 ldltlo.compute(A);
262 VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix());
263 vecX = ldltlo.solve(vecB);
264
265 if(ldltlo.vectorD().real().cwiseAbs().minCoeff()>RealScalar(0))
266 {
267 VERIFY_IS_APPROX(A * vecX,vecB);
268 }
269 else
270 {
271 RealScalar large_tol = sqrt(test_precision<RealScalar>());
272 VERIFY((A * vecX).isApprox(vecB, large_tol));
273
274 ++g_test_level;
275 VERIFY_IS_APPROX(A * vecX,vecB);
276 --g_test_level;
277 }
278 }
279 }
280
281 // update/downdate
282 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) ));
283 CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) ));
284 }
285
cholesky_cplx(const MatrixType & m)286 template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
287 {
288 // classic test
289 cholesky(m);
290
291 // test mixing real/scalar types
292
293 Index rows = m.rows();
294 Index cols = m.cols();
295
296 typedef typename MatrixType::Scalar Scalar;
297 typedef typename NumTraits<Scalar>::Real RealScalar;
298 typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType;
299 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
300
301 RealMatrixType a0 = RealMatrixType::Random(rows,cols);
302 VectorType vecB = VectorType::Random(rows), vecX(rows);
303 MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
304 RealMatrixType symm = a0 * a0.adjoint();
305 // let's make sure the matrix is not singular or near singular
306 for (int k=0; k<3; ++k)
307 {
308 RealMatrixType a1 = RealMatrixType::Random(rows,cols);
309 symm += a1 * a1.adjoint();
310 }
311
312 {
313 RealMatrixType symmLo = symm.template triangularView<Lower>();
314
315 LLT<RealMatrixType,Lower> chollo(symmLo);
316 VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
317
318 check_solverbase<VectorType, VectorType>(symm, chollo, rows, rows, 1);
319 //check_solverbase<MatrixType, MatrixType>(symm, chollo, rows, cols, rows);
320 }
321
322 // LDLT
323 {
324 int sign = internal::random<int>()%2 ? 1 : -1;
325
326 if(sign == -1)
327 {
328 symm = -symm; // test a negative matrix
329 }
330
331 RealMatrixType symmLo = symm.template triangularView<Lower>();
332
333 LDLT<RealMatrixType,Lower> ldltlo(symmLo);
334 VERIFY(ldltlo.info()==Success);
335 VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
336
337 check_solverbase<VectorType, VectorType>(symm, ldltlo, rows, rows, 1);
338 //check_solverbase<MatrixType, MatrixType>(symm, ldltlo, rows, cols, rows);
339 }
340 }
341
342 // regression test for bug 241
cholesky_bug241(const MatrixType & m)343 template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
344 {
345 eigen_assert(m.rows() == 2 && m.cols() == 2);
346
347 typedef typename MatrixType::Scalar Scalar;
348 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
349
350 MatrixType matA;
351 matA << 1, 1, 1, 1;
352 VectorType vecB;
353 vecB << 1, 1;
354 VectorType vecX = matA.ldlt().solve(vecB);
355 VERIFY_IS_APPROX(matA * vecX, vecB);
356 }
357
358 // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
359 // This test checks that LDLT reports correctly that matrix is indefinite.
360 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
cholesky_definiteness(const MatrixType & m)361 template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
362 {
363 eigen_assert(m.rows() == 2 && m.cols() == 2);
364 MatrixType mat;
365 LDLT<MatrixType> ldlt(2);
366
367 {
368 mat << 1, 0, 0, -1;
369 ldlt.compute(mat);
370 VERIFY(ldlt.info()==Success);
371 VERIFY(!ldlt.isNegative());
372 VERIFY(!ldlt.isPositive());
373 VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
374 }
375 {
376 mat << 1, 2, 2, 1;
377 ldlt.compute(mat);
378 VERIFY(ldlt.info()==Success);
379 VERIFY(!ldlt.isNegative());
380 VERIFY(!ldlt.isPositive());
381 VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
382 }
383 {
384 mat << 0, 0, 0, 0;
385 ldlt.compute(mat);
386 VERIFY(ldlt.info()==Success);
387 VERIFY(ldlt.isNegative());
388 VERIFY(ldlt.isPositive());
389 VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
390 }
391 {
392 mat << 0, 0, 0, 1;
393 ldlt.compute(mat);
394 VERIFY(ldlt.info()==Success);
395 VERIFY(!ldlt.isNegative());
396 VERIFY(ldlt.isPositive());
397 VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
398 }
399 {
400 mat << -1, 0, 0, 0;
401 ldlt.compute(mat);
402 VERIFY(ldlt.info()==Success);
403 VERIFY(ldlt.isNegative());
404 VERIFY(!ldlt.isPositive());
405 VERIFY_IS_APPROX(mat,ldlt.reconstructedMatrix());
406 }
407 }
408
409 template<typename>
cholesky_faillure_cases()410 void cholesky_faillure_cases()
411 {
412 MatrixXd mat;
413 LDLT<MatrixXd> ldlt;
414
415 {
416 mat.resize(2,2);
417 mat << 0, 1, 1, 0;
418 ldlt.compute(mat);
419 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
420 VERIFY(ldlt.info()==NumericalIssue);
421 }
422 #if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE_SSE2)
423 {
424 mat.resize(3,3);
425 mat << -1, -3, 3,
426 -3, -8.9999999999999999999, 1,
427 3, 1, 0;
428 ldlt.compute(mat);
429 VERIFY(ldlt.info()==NumericalIssue);
430 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
431 }
432 #endif
433 {
434 mat.resize(3,3);
435 mat << 1, 2, 3,
436 2, 4, 1,
437 3, 1, 0;
438 ldlt.compute(mat);
439 VERIFY(ldlt.info()==NumericalIssue);
440 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
441 }
442
443 {
444 mat.resize(8,8);
445 mat << 0.1, 0, -0.1, 0, 0, 0, 1, 0,
446 0, 4.24667, 0, 2.00333, 0, 0, 0, 0,
447 -0.1, 0, 0.2, 0, -0.1, 0, 0, 0,
448 0, 2.00333, 0, 8.49333, 0, 2.00333, 0, 0,
449 0, 0, -0.1, 0, 0.1, 0, 0, 1,
450 0, 0, 0, 2.00333, 0, 4.24667, 0, 0,
451 1, 0, 0, 0, 0, 0, 0, 0,
452 0, 0, 0, 0, 1, 0, 0, 0;
453 ldlt.compute(mat);
454 VERIFY(ldlt.info()==NumericalIssue);
455 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
456 }
457
458 // bug 1479
459 {
460 mat.resize(4,4);
461 mat << 1, 2, 0, 1,
462 2, 4, 0, 2,
463 0, 0, 0, 1,
464 1, 2, 1, 1;
465 ldlt.compute(mat);
466 VERIFY(ldlt.info()==NumericalIssue);
467 VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix());
468 }
469 }
470
cholesky_verify_assert()471 template<typename MatrixType> void cholesky_verify_assert()
472 {
473 MatrixType tmp;
474
475 LLT<MatrixType> llt;
476 VERIFY_RAISES_ASSERT(llt.matrixL())
477 VERIFY_RAISES_ASSERT(llt.matrixU())
478 VERIFY_RAISES_ASSERT(llt.solve(tmp))
479 VERIFY_RAISES_ASSERT(llt.transpose().solve(tmp))
480 VERIFY_RAISES_ASSERT(llt.adjoint().solve(tmp))
481 VERIFY_RAISES_ASSERT(llt.solveInPlace(tmp))
482
483 LDLT<MatrixType> ldlt;
484 VERIFY_RAISES_ASSERT(ldlt.matrixL())
485 VERIFY_RAISES_ASSERT(ldlt.transpositionsP())
486 VERIFY_RAISES_ASSERT(ldlt.vectorD())
487 VERIFY_RAISES_ASSERT(ldlt.isPositive())
488 VERIFY_RAISES_ASSERT(ldlt.isNegative())
489 VERIFY_RAISES_ASSERT(ldlt.solve(tmp))
490 VERIFY_RAISES_ASSERT(ldlt.transpose().solve(tmp))
491 VERIFY_RAISES_ASSERT(ldlt.adjoint().solve(tmp))
492 VERIFY_RAISES_ASSERT(ldlt.solveInPlace(tmp))
493 }
494
EIGEN_DECLARE_TEST(cholesky)495 EIGEN_DECLARE_TEST(cholesky)
496 {
497 int s = 0;
498 for(int i = 0; i < g_repeat; i++) {
499 CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
500 CALL_SUBTEST_3( cholesky(Matrix2d()) );
501 CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
502 CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
503 CALL_SUBTEST_4( cholesky(Matrix3f()) );
504 CALL_SUBTEST_5( cholesky(Matrix4d()) );
505
506 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
507 CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) );
508 TEST_SET_BUT_UNUSED_VARIABLE(s)
509
510 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
511 CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
512 TEST_SET_BUT_UNUSED_VARIABLE(s)
513 }
514 // empty matrix, regression test for Bug 785:
515 CALL_SUBTEST_2( cholesky(MatrixXd(0,0)) );
516
517 // This does not work yet:
518 // CALL_SUBTEST_2( cholesky(Matrix<double,0,0>()) );
519
520 CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
521 CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
522 CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
523 CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
524
525 // Test problem size constructors
526 CALL_SUBTEST_9( LLT<MatrixXf>(10) );
527 CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
528
529 CALL_SUBTEST_2( cholesky_faillure_cases<void>() );
530
531 TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries)
532 }
533