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
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 */
MatrixFunctionAtomic(StemFunction f)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>
compute(const MatrixType & A)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>
computeMu()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>
taylorConverged(Index s,const MatrixType & F,const MatrixType & Fincr,const MatrixType & P)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