aboutsummaryrefslogtreecommitdiff
path: root/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h
diff options
context:
space:
mode:
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h')
-rw-r--r--unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h131
1 files changed, 0 insertions, 131 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h
deleted file mode 100644
index efe332c48..000000000
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h
+++ /dev/null
@@ -1,131 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_MATRIX_FUNCTION_ATOMIC
-#define EIGEN_MATRIX_FUNCTION_ATOMIC
-
-namespace Eigen {
-
-/** \ingroup MatrixFunctions_Module
- * \class MatrixFunctionAtomic
- * \brief Helper class for computing matrix functions of atomic matrices.
- *
- * \internal
- * Here, an atomic matrix is a triangular matrix whose diagonal
- * entries are close to each other.
- */
-template <typename MatrixType>
-class MatrixFunctionAtomic
-{
- public:
-
- typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::Index Index;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef typename internal::stem_function<Scalar>::type StemFunction;
- typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
-
- /** \brief Constructor
- * \param[in] f matrix function to compute.
- */
- MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
-
- /** \brief Compute matrix function of atomic matrix
- * \param[in] A argument of matrix function, should be upper triangular and atomic
- * \returns f(A), the matrix function evaluated at the given matrix
- */
- MatrixType compute(const MatrixType& A);
-
- private:
-
- // Prevent copying
- MatrixFunctionAtomic(const MatrixFunctionAtomic&);
- MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&);
-
- void computeMu();
- bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P);
-
- /** \brief Pointer to scalar function */
- StemFunction* m_f;
-
- /** \brief Size of matrix function */
- Index m_Arows;
-
- /** \brief Mean of eigenvalues */
- Scalar m_avgEival;
-
- /** \brief Argument shifted by mean of eigenvalues */
- MatrixType m_Ashifted;
-
- /** \brief Constant used to determine whether Taylor series has converged */
- RealScalar m_mu;
-};
-
-template <typename MatrixType>
-MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
-{
- // TODO: Use that A is upper triangular
- m_Arows = A.rows();
- m_avgEival = A.trace() / Scalar(RealScalar(m_Arows));
- m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows);
- computeMu();
- MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows);
- MatrixType P = m_Ashifted;
- MatrixType Fincr;
- for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary
- Fincr = m_f(m_avgEival, static_cast<int>(s)) * P;
- F += Fincr;
- P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted;
- if (taylorConverged(s, F, Fincr, P)) {
- return F;
- }
- }
- eigen_assert("Taylor series does not converge" && 0);
- return F;
-}
-
-/** \brief Compute \c m_mu. */
-template <typename MatrixType>
-void MatrixFunctionAtomic<MatrixType>::computeMu()
-{
- const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted;
- VectorType e = VectorType::Ones(m_Arows);
- N.template triangularView<Upper>().solveInPlace(e);
- m_mu = e.cwiseAbs().maxCoeff();
-}
-
-/** \brief Determine whether Taylor series has converged */
-template <typename MatrixType>
-bool MatrixFunctionAtomic<MatrixType>::taylorConverged(Index s, const MatrixType& F,
- const MatrixType& Fincr, const MatrixType& P)
-{
- const Index n = F.rows();
- const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff();
- const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff();
- if (Fincr_norm < NumTraits<Scalar>::epsilon() * F_norm) {
- RealScalar delta = 0;
- RealScalar rfactorial = 1;
- for (Index r = 0; r < n; r++) {
- RealScalar mx = 0;
- for (Index i = 0; i < n; i++)
- mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast<int>(s+r))));
- if (r != 0)
- rfactorial *= RealScalar(r);
- delta = (std::max)(delta, mx / rfactorial);
- }
- const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff();
- if (m_mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm)
- return true;
- }
- return false;
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_MATRIX_FUNCTION_ATOMIC