1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_NO_ASSERTION_CHECKING 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_NO_ASSERTION_CHECKING 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstatic int nb_temporaries; 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { if(size!=0) nb_temporaries++; } 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Cholesky> 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/QR> 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define VERIFY_EVALUATION_COUNT(XPR,N) {\ 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nb_temporaries = 0; \ 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath XPR; \ 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY( (#XPR) && nb_temporaries==N ); \ 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType symmLo = symm.template triangularView<Lower>(); 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType symmUp = symm.template triangularView<Upper>(); 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType symmCpy = symm; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CholType<MatrixType,Lower> chollo(symmLo); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CholType<MatrixType,Upper> cholup(symmUp); 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k<10; ++k) 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType vec = VectorType::Random(symm.rows()); 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar sigma = internal::random<RealScalar>(); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symmCpy += sigma * vec * vec.adjoint(); 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we are doing some downdates, so it might be the case that the matrix is not SPD anymore 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CholType<MatrixType,Lower> chol(symmCpy); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(chol.info()!=Success) 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chollo.rankUpdate(vec, sigma); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholup.rankUpdate(vec, sigma); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void cholesky(const MatrixType& m) 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* this test covers the following files: 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT.h LDLT.h 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a0 = MatrixType::Random(rows,cols); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType vecB = VectorType::Random(rows), vecX(rows); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SquareMatrixType symm = a0 * a0.adjoint(); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // let's make sure the matrix is not singular or near singular 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k<3; ++k) 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a1 = MatrixType::Random(rows,cols); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symm += a1 * a1.adjoint(); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // to test if really Cholesky only uses the upper triangular part, uncomment the following 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // FIXME: currently that fails !! 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //symm.template part<StrictlyLower>().setZero(); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SquareMatrixType symmUp = symm.template triangularView<Upper>(); 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SquareMatrixType symmLo = symm.template triangularView<Lower>(); 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<SquareMatrixType,Lower> chollo(symmLo); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = chollo.solve(vecB); 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = chollo.solve(matB); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * matX, matB); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test the upper mode 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<SquareMatrixType,Upper> cholup(symmUp); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = cholup.solve(vecB); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = cholup.solve(matB); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * matX, matB); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType neg = -symmLo; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chollo.compute(neg); 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(chollo.info()==NumericalIssue); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU())); 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // test some special use cases of SelfCwiseBinaryOp: 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 = m1; 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB); 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 = m1; 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 = m1; 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB); 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB)); 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2 = m1; 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // LDLT 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int sign = internal::random<int>()%2 ? 1 : -1; 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(sign == -1) 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symm = -symm; // test a negative matrix 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SquareMatrixType symmUp = symm.template triangularView<Upper>(); 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SquareMatrixType symmLo = symm.template triangularView<Lower>(); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT<SquareMatrixType,Lower> ldltlo(symmLo); 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = ldltlo.solve(vecB); 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = ldltlo.solve(matB); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * matX, matB); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT<SquareMatrixType,Upper> ldltup(symmUp); 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = ldltup.solve(vecB); 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = ldltup.solve(matB); 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * matX, matB); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(MatrixType(ldltup.matrixU().transpose().conjugate()), MatrixType(ldltup.matrixL())); 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(MatrixType::RowsAtCompileTime==Dynamic) 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // note : each inplace permutation requires a small temporary vector (mask) 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check inplace solve 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = matB; 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0); 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval()); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matX = matB; 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // restore 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(sign == -1) 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symm = -symm; 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check matrices coming from linear constraints with Lagrange multipliers 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(rows>=3) 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SquareMatrixType A = symm; 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int c = internal::random<int>(0,rows-2); 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.bottomRightCorner(c,c).setZero(); 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Make sure a solution exists: 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setRandom(); 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecB = A * vecX; 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setZero(); 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ldltlo.compute(A); 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX = ldltlo.solve(vecB); 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A * vecX, vecB); 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check non-full rank matrices 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(rows>=3) 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int r = internal::random<int>(1,rows-1); 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r); 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SquareMatrixType A = a * a.adjoint(); 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Make sure a solution exists: 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setRandom(); 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecB = A * vecX; 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setZero(); 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ldltlo.compute(A); 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX = ldltlo.solve(vecB); 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A * vecX, vecB); 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check matrices with a wide spectrum 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(rows>=3) 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows); 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows); 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int k=0; k<rows; ++k) 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SquareMatrixType A = a * d.asDiagonal() * a.adjoint(); 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Make sure a solution exists: 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setRandom(); 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecB = A * vecX; 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX.setZero(); 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ldltlo.compute(A); 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vecX = ldltlo.solve(vecB); 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A * vecX, vecB); 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // update/downdate 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) )); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) )); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void cholesky_cplx(const MatrixType& m) 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // classic test 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholesky(m); 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test mixing real/scalar types 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealMatrixType a0 = RealMatrixType::Random(rows,cols); 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType vecB = VectorType::Random(rows), vecX(rows); 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealMatrixType symm = a0 * a0.adjoint(); 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // let's make sure the matrix is not singular or near singular 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k<3; ++k) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealMatrixType a1 = RealMatrixType::Random(rows,cols); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symm += a1 * a1.adjoint(); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealMatrixType symmLo = symm.template triangularView<Lower>(); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<RealMatrixType,Lower> chollo(symmLo); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = chollo.solve(vecB); 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// matX = chollo.solve(matB); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// VERIFY_IS_APPROX(symm * matX, matB); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // LDLT 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int sign = internal::random<int>()%2 ? 1 : -1; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(sign == -1) 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath symm = -symm; // test a negative matrix 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealMatrixType symmLo = symm.template triangularView<Lower>(); 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT<RealMatrixType,Lower> ldltlo(symmLo); 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecX = ldltlo.solve(vecB); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(symm * vecX, vecB); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// matX = ldltlo.solve(matB); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// VERIFY_IS_APPROX(symm * matX, matB); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// regression test for bug 241 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void cholesky_bug241(const MatrixType& m) 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m.rows() == 2 && m.cols() == 2); 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matA; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA << 1, 1, 1, 1; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType vecB; 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vecB << 1, 1; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType vecX = matA.ldlt().solve(vecB); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(matA * vecX, vecB); 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This test checks that LDLT reports correctly that matrix is indefinite. 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> void cholesky_definiteness(const MatrixType& m) 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m.rows() == 2 && m.cols() == 2); 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType mat; 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat << 1, 0, 0, -1; 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT<MatrixType> ldlt(mat); 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isNegative()); 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isPositive()); 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat << 1, 2, 2, 1; 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT<MatrixType> ldlt(mat); 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isNegative()); 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isPositive()); 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat << 0, 0, 0, 0; 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT<MatrixType> ldlt(mat); 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(ldlt.isNegative()); 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(ldlt.isPositive()); 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat << 0, 0, 0, 1; 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT<MatrixType> ldlt(mat); 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isNegative()); 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(ldlt.isPositive()); 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat << -1, 0, 0, 0; 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT<MatrixType> ldlt(mat); 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(ldlt.isNegative()); 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(!ldlt.isPositive()); 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void cholesky_verify_assert() 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType tmp; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<MatrixType> llt; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(llt.matrixL()) 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(llt.matrixU()) 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(llt.solve(tmp)) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT<MatrixType> ldlt; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.matrixL()) 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.permutationP()) 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.vectorD()) 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.isPositive()) 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.isNegative()) 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.solve(tmp)) 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp)) 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_cholesky() 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int s = 0; 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( cholesky(Matrix2d()) ); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( cholesky(Matrix3f()) ); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_5( cholesky(Matrix4d()) ); 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() ); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() ); 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test problem size constructors 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_9( LLT<MatrixXf>(10) ); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(s) 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 403