1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2011-2014 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_ITERATIVE_SOLVER_BASE_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_ITERATIVE_SOLVER_BASE_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType> 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct is_ref_compatible_impl 192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprivate: 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template <typename T0> 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang struct any_conversion 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template <typename T> any_conversion(const volatile T&); 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template <typename T> any_conversion(T&); 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang struct yes {int a[1];}; 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang struct no {int a[2];}; 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename T> 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static yes test(const Ref<const T>&, int); 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename T> 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static no test(any_conversion<T>, ...); 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpublic: 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static MatrixType ms_from; 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) }; 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType> 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct is_ref_compatible 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value }; 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value> 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass generic_matrix_wrapper; 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// We have an explicit matrix at hand, compatible with Ref<> 502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType> 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass generic_matrix_wrapper<MatrixType,false> 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpublic: 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<const MatrixType> ActualMatrixType; 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<int UpLo> struct ConstSelfAdjointViewReturnType { 562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType<UpLo>::Type Type; 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixFree = false 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_matrix_wrapper() 642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_dummy(0,0), m_matrix(m_dummy) 652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_matrix_wrapper(const InputType &mat) 692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_matrix(mat) 702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const ActualMatrixType& matrix() const 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_matrix; 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void grab(const EigenBase<MatrixDerived> &mat) 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.~Ref<const MatrixType>(); 812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ::new (&m_matrix) Ref<const MatrixType>(mat.derived()); 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void grab(const Ref<const MatrixType> &mat) 852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(&(mat.derived()) != &m_matrix) 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix.~Ref<const MatrixType>(); 892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ::new (&m_matrix) Ref<const MatrixType>(mat); 902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprotected: 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixType m_dummy; // used to default initialize the Ref<> object 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ActualMatrixType m_matrix; 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// MatrixType is not compatible with Ref<> -> matrix-free wrapper 992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType> 1002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass generic_matrix_wrapper<MatrixType,true> 1012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 1022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpublic: 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef MatrixType ActualMatrixType; 1042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<int UpLo> struct ConstSelfAdjointViewReturnType 1052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef ActualMatrixType Type; 1072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 1102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixFree = true 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_matrix_wrapper() 1142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : mp_matrix(0) 1152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 1162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang generic_matrix_wrapper(const MatrixType &mat) 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : mp_matrix(&mat) 1192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const ActualMatrixType& matrix() const 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *mp_matrix; 1242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void grab(const MatrixType &mat) 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang mp_matrix = &mat; 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprotected: 1322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const ActualMatrixType *mp_matrix; 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 1342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 1362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup IterativeLinearSolvers_Module 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Base class for linear iterative solvers 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename Derived> 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass IterativeSolverBase : public SparseSolverBase<Derived> 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangprotected: 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseSolverBase<Derived> Base; 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::m_isInitialized; 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::traits<Derived>::MatrixType MatrixType; 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::traits<Derived>::Preconditioner Preconditioner; 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::StorageIndex StorageIndex; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = MatrixType::ColsAtCompileTime, 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::derived; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor. */ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IterativeSolverBase() 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initialize the solver with matrix \a A for further \c Ax=b solving. 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor is a shortcut for the default constructor followed 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by a call to compute(). 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \warning this class stores a reference to the matrix A as well as some 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * precomputed values that depend on it. Therefore, if \a A is changed 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class becomes invalid. Call compute() to update it with the new 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A, or modify a copy of A. 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit IterativeSolverBase(const EigenBase<MatrixDerived>& A) 1832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_matrixWrapper(A.derived()) 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 1862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang compute(matrix()); 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~IterativeSolverBase() {} 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** Initializes the iterative solver for the sparsity pattern of the matrix \a A for further solving \c Ax=b problems. 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 1932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Currently, this function mostly calls analyzePattern on the preconditioner. In the future 1942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * we might, for instance, implement column reordering for faster matrix vector products. 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 1972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Derived& analyzePattern(const EigenBase<MatrixDerived>& A) 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang grab(A.derived()); 2002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_preconditioner.analyzePattern(matrix()); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 2032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = m_preconditioner.info(); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems. 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 2092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Currently, this function mostly calls factorize on the preconditioner. 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \warning this class stores a reference to the matrix A as well as some 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * precomputed values that depend on it. Therefore, if \a A is changed 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class becomes invalid. Call compute() to update it with the new 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A, or modify a copy of A. 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Derived& factorize(const EigenBase<MatrixDerived>& A) 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 2202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang grab(A.derived()); 2212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_preconditioner.factorize(matrix()); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = m_preconditioner.info(); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems. 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 2292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Currently, this function mostly initializes/computes the preconditioner. In the future 2302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * we might, for instance, implement column reordering for faster matrix vector products. 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \warning this class stores a reference to the matrix A as well as some 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * precomputed values that depend on it. Therefore, if \a A is changed 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class becomes invalid. Call compute() to update it with the new 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A, or modify a copy of A. 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 2382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Derived& compute(const EigenBase<MatrixDerived>& A) 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang grab(A.derived()); 2412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_preconditioner.compute(matrix()); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 2452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = m_preconditioner.info(); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 2502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rows() const { return matrix().rows(); } 2512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 2532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index cols() const { return matrix().cols(); } 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the tolerance threshold used by the stopping criteria. 2562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa setTolerance() 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar tolerance() const { return m_tolerance; } 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** Sets the tolerance threshold used by the stopping criteria. 2612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 2622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. 2632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * The default value is the machine precision given by NumTraits<Scalar>::epsilon() 2642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Derived& setTolerance(const RealScalar& tolerance) 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tolerance = tolerance; 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a read-write reference to the preconditioner for custom configuration. */ 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Preconditioner& preconditioner() { return m_preconditioner; } 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a read-only reference to the preconditioner. */ 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Preconditioner& preconditioner() const { return m_preconditioner; } 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the max number of iterations. 2782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * It is either the value setted by setMaxIterations or, by default, 2792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * twice the number of columns of the matrix. 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 2812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index maxIterations() const 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** Sets the max number of iterations. 2872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Default is twice the number of columns of the matrix. 2882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 2892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Derived& setMaxIterations(Index maxIters) 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_maxIterations = maxIters; 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return derived(); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the number of iterations performed during the last solve */ 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index iterations() const 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iterations; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the tolerance error reached during the last solve. 3032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * It is a close approximation of the true relative residual error |Ax-b|/|b|. 3042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar error() const 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_error; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * and \a x0 as an initial solution. 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 3142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa solve(), compute() 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 3162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Rhs,typename Guess> 3172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline const SolveWithGuess<Derived, Rhs, Guess> 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_isInitialized && "Solver is not initialized."); 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); 3222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns Success if the iterations converged, and NoConvergence otherwise. */ 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 3332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Rhs, typename DestDerived> 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==b.rows()); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rhsCols = b.cols(); 3392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index size = b.rows(); 3402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang DestDerived& dest(aDest.derived()); 3412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename DestDerived::Scalar DestScalar; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Eigen::Matrix<DestScalar,Dynamic,1> tb(size); 3432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Eigen::Matrix<DestScalar,Dynamic,1> tx(cols()); 3442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // We do not directly fill dest because sparse expressions have to be free of aliasing issue. 3452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other. 3462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename DestDerived::PlainObject tmp(cols(),rhsCols); 3472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(Index k=0; k<rhsCols; ++k) 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tb = b.col(k); 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tx = derived().solve(tb); 3512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang tmp.col(k) = tx.sparseView(0); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 3532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dest.swap(tmp); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = false; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = false; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_maxIterations = -1; 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tolerance = NumTraits<Scalar>::epsilon(); 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 3652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper; 3672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; 3682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const ActualMatrixType& matrix() const 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_matrixWrapper.matrix(); 3722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 3752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void grab(const InputType &A) 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrixWrapper.grab(A); 3782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixWrapper m_matrixWrapper; 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Preconditioner m_preconditioner; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_maxIterations; 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_tolerance; 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable RealScalar m_error; 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang mutable Index m_iterations; 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable ComputationInfo m_info; 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang mutable bool m_analysisIsOk, m_factorizationIsOk; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_ITERATIVE_SOLVER_BASE_H 395