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 | */
|
---|
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
|
---|