1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 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_BASIC_PRECONDITIONERS_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_BASIC_PRECONDITIONERS_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup IterativeLinearSolvers_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A preconditioner based on the digonal entries 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for: 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \code 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A.diagonal().asDiagonal() . x = b 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \endcode 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Scalar the type of the scalar. 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This preconditioner is suitable for both selfadjoint and general problems. 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The diagonal entries are pre-inverted and stored into a dense vector. 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note A variant that has yet to be implemented would attempt to preserve the norm of each column. 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename _Scalar> 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass DiagonalPreconditioner 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> Vector; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Vector::Index Index; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // this typedef is only to export the scalar type and compile-time dimensions to solve_retval 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalPreconditioner() : m_isInitialized(false) {} 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatType> 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(mat); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_invdiag.size(); } 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_invdiag.size(); } 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatType> 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalPreconditioner& analyzePattern(const MatType& ) 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatType> 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalPreconditioner& factorize(const MatType& mat) 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_invdiag.resize(mat.cols()); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int j=0; j<mat.outerSize(); ++j) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatType::InnerIterator it(mat,j); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(it && it.index()!=j) ++it; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(it && it.index()==j) 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_invdiag(j) = Scalar(1)/it.value(); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_invdiag(j) = 0; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatType> 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalPreconditioner& compute(const MatType& mat) 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return factorize(mat); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs, typename Dest> 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solve(const Rhs& b, Dest& x) const 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_invdiag.array() * b.array() ; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs> 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized."); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_invdiag.size()==b.rows() 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived()); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Vector m_invdiag; 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs> 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs> 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef DiagonalPreconditioner<_MatrixType> Dec; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup IterativeLinearSolvers_Module 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A naive preconditioner which approximates any matrix as the identity matrix 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class DiagonalPreconditioner 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass IdentityPreconditioner 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IdentityPreconditioner() {} 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IdentityPreconditioner(const MatrixType& ) {} 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IdentityPreconditioner& factorize(const MatrixType& ) { return *this; } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IdentityPreconditioner& compute(const MatrixType& ) { return *this; } 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const Rhs& solve(const Rhs& b) const { return b; } 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_BASIC_PRECONDITIONERS_H 150