1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Jitse Niesen <jitse@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
13namespace 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  */
23template <typename MatrixType>
24class 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
70template <typename MatrixType>
71MatrixType 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. */
94template <typename MatrixType>
95void 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 */
104template <typename MatrixType>
105bool 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