MatrixFunctionAtomic.h
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_MATRIX_FUNCTION_ATOMIC
26 #define EIGEN_MATRIX_FUNCTION_ATOMIC
27 
28 namespace Eigen {
29 
38 template <typename MatrixType>
40 {
41  public:
42 
43  typedef typename MatrixType::Scalar Scalar;
44  typedef typename MatrixType::Index Index;
45  typedef typename NumTraits<Scalar>::Real RealScalar;
46  typedef typename internal::stem_function<Scalar>::type StemFunction;
47  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
48 
52  MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
53 
58  MatrixType compute(const MatrixType& A);
59 
60  private:
61 
62  // Prevent copying
64  MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&);
65 
66  void computeMu();
67  bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P);
68 
70  StemFunction* m_f;
71 
73  Index m_Arows;
74 
76  Scalar m_avgEival;
77 
79  MatrixType m_Ashifted;
80 
82  RealScalar m_mu;
83 };
84 
85 template <typename MatrixType>
86 MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
87 {
88  // TODO: Use that A is upper triangular
89  m_Arows = A.rows();
90  m_avgEival = A.trace() / Scalar(RealScalar(m_Arows));
91  m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows);
92  computeMu();
93  MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows);
94  MatrixType P = m_Ashifted;
95  MatrixType Fincr;
96  for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary
97  Fincr = m_f(m_avgEival, static_cast<int>(s)) * P;
98  F += Fincr;
99  P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted;
100  if (taylorConverged(s, F, Fincr, P)) {
101  return F;
102  }
103  }
104  eigen_assert("Taylor series does not converge" && 0);
105  return F;
106 }
107 
109 template <typename MatrixType>
111 {
112  const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted;
113  VectorType e = VectorType::Ones(m_Arows);
114  N.template triangularView<Upper>().solveInPlace(e);
115  m_mu = e.cwiseAbs().maxCoeff();
116 }
117 
119 template <typename MatrixType>
120 bool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s, const MatrixType& F,
121  const MatrixType& Fincr, const MatrixType& P)
122 {
123  const Index n = F.rows();
124  const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
125  const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
126  if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
127  RealScalar delta = 0;
128  RealScalar rfactorial = 1;
129  for (Index r = 0; r < n; r++) {
130  RealScalar mx = 0;
131  for (Index i = 0; i < n; i++)
132  mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r))));
133  if (r != 0)
134  rfactorial *= RealScalar(r);
135  delta = (std::max)(delta, mx / rfactorial);
136  }
137  const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
138  if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
139  return true;
140  }
141  return false;
142 }
143 
144 } // end namespace Eigen
145 
146 #endif // EIGEN_MATRIX_FUNCTION_ATOMIC