1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2008-2012 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_SIMPLICIAL_CHOLESKY_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SIMPLICIAL_CHOLESKY_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathenum SimplicialCholeskyMode { 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholeskyLLT, 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholeskyLDLT 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename CholMatrixType, typename InputMatrixType> 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang struct simplicial_cholesky_grab_input { 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef CholMatrixType const * ConstCholMatrixPtr; 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp) 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang tmp = input; 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang pmat = &tmp; 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixType> 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang struct simplicial_cholesky_grab_input<MatrixType,MatrixType> { 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef MatrixType const * ConstMatrixPtr; 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/) 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang pmat = &input; 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SparseCholesky_Module 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief A base class for direct sparse Cholesky factorizations 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * selfadjoint and positive definite. These factorizations allow for solving A.X = B where 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * X and B can be either dense or sparse. 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * such that the factorized matrix is P A P^-1. 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \tparam Derived the type of the derived class, that is the actual factorization type. 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass SimplicialCholeskyBase : public SparseSolverBase<Derived> 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseSolverBase<Derived> Base; 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::m_isInitialized; 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::traits<Derived>::MatrixType MatrixType; 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::traits<Derived>::OrderingType OrderingType; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = internal::traits<Derived>::UpLo }; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef CholMatrixType const * ConstCholMatrixPtr; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<StorageIndex,Dynamic,1> VectorI; 712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = MatrixType::ColsAtCompileTime, 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::derived; 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor */ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholeskyBase() 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_info(Success), m_shiftOffset(0), m_shiftScale(1) 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit SimplicialCholeskyBase(const MatrixType& matrix) 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_info(Success), m_shiftOffset(0), m_shiftScale(1) 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().compute(matrix); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~SimplicialCholeskyBase() 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& derived() { return *static_cast<Derived*>(this); } 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Derived& derived() const { return *static_cast<const Derived*>(this); } 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the matrix.appears to be negative. 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Decomposition is not initialized."); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the permutation P 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa permutationPinv() */ 1152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return m_P; } 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the inverse P^-1 of the permutation P 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa permutationP() */ 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return m_Pinv; } 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization. 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c d_ii = \a offset + \a scale * \c d_ii 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default is the identity transformation with \a offset=0, and \a scale=1. 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to \c *this. 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_shiftOffset = offset; 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_shiftScale = scale; 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_PARSED_BY_DOXYGEN 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Stream> 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void dumpMemory(Stream& s) 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int total = 0; 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==b.rows()); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_info!=Success) 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_P.size()>0) 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = m_P * b; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = b; 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_matrix.nonZeros()>0) // otherwise L==I 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().matrixL().solveInPlace(dest); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_diag.size()>0) 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = m_diag.asDiagonal().inverse() * dest; 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_matrix.nonZeros()>0) // otherwise U==I 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().matrixU().solveInPlace(dest); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_P.size()>0) 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = m_Pinv * dest; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Rhs,typename Dest> 1832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 1842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::solve_sparse_through_dense_panels(derived(), b, dest); 1862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_PARSED_BY_DOXYGEN 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the sparse Cholesky decomposition of \a matrix */ 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<bool DoLDLT> 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute(const MatrixType& matrix) 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(matrix.rows()==matrix.cols()); 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = matrix.cols(); 1982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CholMatrixType tmp(size,size); 1992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ConstCholMatrixPtr pmat; 2002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ordering(matrix, pmat, tmp); 2012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang analyzePattern_preordered(*pmat, DoLDLT); 2022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang factorize_preordered<DoLDLT>(*pmat); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<bool DoLDLT> 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& a) 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 2092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index size = a.cols(); 2102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CholMatrixType tmp(size,size); 2112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ConstCholMatrixPtr pmat; 2122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_P.size()==0 && (UpLo&Upper)==Upper) 2142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // If there is no ordering, try to directly use the input matrix without any copy 2162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp); 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 2192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 2212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang pmat = &tmp; 2222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang factorize_preordered<DoLDLT>(*pmat); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<bool DoLDLT> 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize_preordered(const CholMatrixType& a); 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& a, bool doLDLT) 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 2332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index size = a.cols(); 2342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CholMatrixType tmp(size,size); 2352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ConstCholMatrixPtr pmat; 2362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ordering(a, pmat, tmp); 2372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang analyzePattern_preordered(*pmat,doLDLT); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** keeps off-diagonal entries; drops diagonal entries */ 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct keep_diag { 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool operator() (const Index& row, const Index& col, const Scalar&) const 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return row!=col; 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable ComputationInfo m_info; 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_factorizationIsOk; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_analysisIsOk; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CholMatrixType m_matrix; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType m_diag; // the diagonal coefficients (LDLT mode) 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorI m_parent; // elimination tree 2582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorI m_nonZerosPerCol; 2592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation 2602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_shiftOffset; 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_shiftScale; 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT; 2672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT; 2682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky; 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _Ordering OrderingType; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 2782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 2792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL; 2812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; 2822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 2832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _Ordering OrderingType; 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 2932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 2942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL; 2952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 2972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _Ordering OrderingType; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SparseCholesky_Module 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SimplicialLLT 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A direct sparse LLT Cholesky factorizations 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class provides a LL^T Cholesky factorizations of sparse matrices that are 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * selfadjoint and positive definite. The factorization allows for solving A.X = B where 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * X and B can be either dense or sparse. 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * such that the factorized matrix is P A P^-1. 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or Upper. Default is Lower. 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int _UpLo, typename _Ordering> 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SimplicialCholeskyBase<SimplicialLLT> Base; 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 3382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::traits<SimplicialLLT> Traits; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::MatrixL MatrixL; 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::MatrixU MatrixU; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor */ 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialLLT() : Base() {} 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Constructs and performs the LLT factorization of \a matrix */ 3482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit SimplicialLLT(const MatrixType& matrix) 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(matrix) {} 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of the factor L */ 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixL matrixL() const { 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getL(Base::m_matrix); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of the factor U (= L^*) */ 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixU matrixU() const { 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getU(Base::m_matrix); 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the sparse Cholesky decomposition of \a matrix */ 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialLLT& compute(const MatrixType& matrix) 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template compute<false>(matrix); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a symbolic decomposition on the sparcity of \a matrix. 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is particularly useful when solving for several problems having the same structure. 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& a) 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(a, false); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a numeric decomposition of \a matrix 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& a) 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template factorize<false>(a); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the determinant of the underlying matrix from the current factorization */ 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar determinant() const 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar detL = Base::m_matrix.diagonal().prod(); 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return numext::abs2(detL); 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SparseCholesky_Module 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SimplicialLDLT 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A direct sparse LDLT Cholesky factorizations without square root. 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * selfadjoint and positive definite. The factorization allows for solving A.X = B where 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * X and B can be either dense or sparse. 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * such that the factorized matrix is P A P^-1. 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or Upper. Default is Lower. 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 4162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 4172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int _UpLo, typename _Ordering> 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SimplicialCholeskyBase<SimplicialLDLT> Base; 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 4302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::traits<SimplicialLDLT> Traits; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::MatrixL MatrixL; 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Traits::MatrixU MatrixU; 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor */ 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialLDLT() : Base() {} 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Constructs and performs the LLT factorization of \a matrix */ 4402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit SimplicialLDLT(const MatrixType& matrix) 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(matrix) {} 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a vector expression of the diagonal D */ 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const VectorType vectorD() const { 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Base::m_diag; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of the factor L */ 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixL matrixL() const { 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getL(Base::m_matrix); 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns an expression of the factor U (= L^*) */ 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixU matrixU() const { 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getU(Base::m_matrix); 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the sparse Cholesky decomposition of \a matrix */ 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialLDLT& compute(const MatrixType& matrix) 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template compute<true>(matrix); 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a symbolic decomposition on the sparcity of \a matrix. 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is particularly useful when solving for several problems having the same structure. 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& a) 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(a, true); 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a numeric decomposition of \a matrix 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& a) 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template factorize<true>(a); 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the determinant of the underlying matrix from the current factorization */ 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar determinant() const 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Base::m_diag.prod(); 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \deprecated use SimplicialLDLT or class SimplicialLLT 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \ingroup SparseCholesky_Module 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SimplicialCholesky 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class SimplicialLDLT, class SimplicialLLT 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int _UpLo, typename _Ordering> 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SimplicialCholeskyBase<SimplicialCholesky> Base; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 5112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 5122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::traits<SimplicialCholesky> Traits; 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholesky() : Base(), m_LDLT(true) {} 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit SimplicialCholesky(const MatrixType& matrix) 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(), m_LDLT(true) 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholesky& setMode(SimplicialCholeskyMode mode) 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath switch(mode) 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case SimplicialCholeskyLLT: 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_LDLT = false; 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case SimplicialCholeskyLDLT: 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_LDLT = true; 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath default: 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const VectorType vectorD() const { 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Base::m_diag; 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const CholMatrixType rawMatrix() const { 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Base::m_matrix; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the sparse Cholesky decomposition of \a matrix */ 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SimplicialCholesky& compute(const MatrixType& matrix) 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_LDLT) 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template compute<true>(matrix); 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template compute<false>(matrix); 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a symbolic decomposition on the sparcity of \a matrix. 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is particularly useful when solving for several problems having the same structure. 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& a) 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(a, m_LDLT); 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Performs a numeric decomposition of \a matrix 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& a) 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_LDLT) 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template factorize<true>(a); 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::template factorize<false>(a); 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 5892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::m_matrix.rows()==b.rows()); 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(Base::m_info!=Success) 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(Base::m_P.size()>0) 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = Base::m_P * b; 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = b; 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(Base::m_matrix.nonZeros()>0) // otherwise L==I 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_LDLT) 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLTTraits::getL(Base::m_matrix).solveInPlace(dest); 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(Base::m_diag.size()>0) 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = Base::m_diag.asDiagonal().inverse() * dest; 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Base::m_matrix.nonZeros()>0) // otherwise I==I 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_LDLT) 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLTTraits::getU(Base::m_matrix).solveInPlace(dest); 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(Base::m_P.size()>0) 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest = Base::m_Pinv * dest; 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 6252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \internal */ 6262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Rhs,typename Dest> 6272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 6282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::solve_sparse_through_dense_panels(*this, b, dest); 6302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar determinant() const 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_LDLT) 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Base::m_diag.prod(); 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); 6417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return numext::abs2(detL); 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_LDLT; 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 6502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap) 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = a.rows(); 6542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang pmat = ≈ 6552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Note that ordering methods compute the inverse permutation 6562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 6582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CholMatrixType C; 6602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang C = a.template selfadjointView<UpLo>(); 6612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang OrderingType ordering; 6632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ordering(C,m_Pinv); 6642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_Pinv.size()>0) m_P = m_Pinv.inverse(); 6672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else m_P.resize(0); 6687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ap.resize(size,size); 6702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 6732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_Pinv.resize(0); 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_P.resize(0); 6762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor) 6772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // we have to transpose the lower part to to the upper one 6792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ap.resize(size,size); 6802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>(); 6812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 6832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap); 6842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SIMPLICIAL_CHOLESKY_H 690