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