1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk> 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_MATRIX_FUNCTION_ATOMIC 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_MATRIX_FUNCTION_ATOMIC 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class MatrixFunctionAtomic 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Helper class for computing matrix functions of atomic matrices. 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \internal 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Here, an atomic matrix is a triangular matrix whose diagonal 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * entries are close to each other. 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixFunctionAtomic 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::stem_function<Scalar>::type StemFunction; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] f matrix function to compute. 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixFunctionAtomic(StemFunction f) : m_f(f) { } 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute matrix function of atomic matrix 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] A argument of matrix function, should be upper triangular and atomic 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns f(A), the matrix function evaluated at the given matrix 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType compute(const MatrixType& A); 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Prevent copying 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixFunctionAtomic(const MatrixFunctionAtomic&); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeMu(); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P); 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Pointer to scalar function */ 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StemFunction* m_f; 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Size of matrix function */ 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index m_Arows; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Mean of eigenvalues */ 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar m_avgEival; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Argument shifted by mean of eigenvalues */ 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_Ashifted; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constant used to determine whether Taylor series has converged */ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_mu; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A) 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: Use that A is upper triangular 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_Arows = A.rows(); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_avgEival = A.trace() / Scalar(RealScalar(m_Arows)); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath computeMu(); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType P = m_Ashifted; 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType Fincr; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Fincr = m_f(m_avgEival, static_cast<int>(s)) * P; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath F += Fincr; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (taylorConverged(s, F, Fincr, P)) { 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return F; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert("Taylor series does not converge" && 0); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return F; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \brief Compute \c m_mu. */ 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixFunctionAtomic<MatrixType>::computeMu() 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType e = VectorType::Ones(m_Arows); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N.template triangularView<Upper>().solveInPlace(e); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_mu = e.cwiseAbs().maxCoeff(); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \brief Determine whether Taylor series has converged */ 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s, const MatrixType& F, 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& Fincr, const MatrixType& P) 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index n = F.rows(); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) { 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar delta = 0; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar rfactorial = 1; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index r = 0; r < n; r++) { 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar mx = 0; 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < n; i++) 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r)))); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (r != 0) 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rfactorial *= RealScalar(r); 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath delta = (std::max)(delta, mx / rfactorial); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_MATRIX_FUNCTION_ATOMIC 132