mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-01-18 14:34:17 +08:00
Refactor MatrixFunction class: Split new class MatrixFunctionAtomic off.
This commit is contained in:
Jitse Niesen
parent
a25c9b1e46
commit
d35cc381fe
@ -78,6 +78,7 @@ EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase<Derived>& M,
|
||||
typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f,
|
||||
typename MatrixBase<Derived>::PlainMatrixType* result);
|
||||
|
||||
#include "MatrixFunctionAtomic.h"
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
* \class MatrixFunction
|
||||
@ -169,15 +170,10 @@ class MatrixFunction<MatrixType, 1, 1>
|
||||
void computeTriangular(const MatrixType& T, MatrixType& result, const IntVectorType& blockSize);
|
||||
void computeBlockAtomic(const MatrixType& T, MatrixType& result, const IntVectorType& blockSize);
|
||||
MatrixType solveTriangularSylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C);
|
||||
MatrixType computeAtomic(const MatrixType& T);
|
||||
void divideInBlocks(const VectorType& v, listOfLists* result);
|
||||
void constructPermutation(const VectorType& diag, const listOfLists& blocks,
|
||||
IntVectorType& blockSize, IntVectorType& permutation);
|
||||
|
||||
RealScalar computeMu(const MatrixType& M);
|
||||
bool taylorConverged(const MatrixType& T, int s, const MatrixType& F,
|
||||
const MatrixType& Fincr, const MatrixType& P, RealScalar mu);
|
||||
|
||||
static const RealScalar separation() { return static_cast<RealScalar>(0.01); }
|
||||
StemFunction *m_f;
|
||||
};
|
||||
@ -271,11 +267,11 @@ void MatrixFunction<MatrixType,1,1>::computeTriangular(const MatrixType& T, Matr
|
||||
|
||||
/** \brief Solve a triangular Sylvester equation AX + XB = C
|
||||
*
|
||||
* \param[in] A The matrix A; should be square and upper triangular
|
||||
* \param[in] B The matrix B; should be square and upper triangular
|
||||
* \param[in] C The matrix C; should have correct size.
|
||||
* \param[in] A the matrix A; should be square and upper triangular
|
||||
* \param[in] B the matrix B; should be square and upper triangular
|
||||
* \param[in] C the matrix C; should have correct size.
|
||||
*
|
||||
* \returns The solution X.
|
||||
* \returns the solution X.
|
||||
*
|
||||
* If A is m-by-m and B is n-by-n, then both C and X are m-by-n.
|
||||
* The (i,j)-th component of the Sylvester equation is
|
||||
@ -346,8 +342,9 @@ void MatrixFunction<MatrixType,1,1>::computeBlockAtomic(const MatrixType& T, Mat
|
||||
result.resize(T.rows(), T.cols());
|
||||
result.setZero();
|
||||
for (int i = 0; i < blockSize.rows(); i++) {
|
||||
MatrixFunctionAtomic<MatrixType> mfa(m_f);
|
||||
result.block(blockStart, blockStart, blockSize(i), blockSize(i))
|
||||
= computeAtomic(T.block(blockStart, blockStart, blockSize(i), blockSize(i)));
|
||||
= mfa.compute(T.block(blockStart, blockStart, blockSize(i), blockSize(i)));
|
||||
blockStart += blockSize(i);
|
||||
}
|
||||
}
|
||||
@ -434,64 +431,6 @@ void MatrixFunction<MatrixType,1,1>::constructPermutation(const VectorType& diag
|
||||
}
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
MatrixType MatrixFunction<MatrixType,1,1>::computeAtomic(const MatrixType& T)
|
||||
{
|
||||
// TODO: Use that T is upper triangular
|
||||
const int n = T.rows();
|
||||
const Scalar sigma = T.trace() / Scalar(n);
|
||||
const MatrixType M = T - sigma * MatrixType::Identity(n, n);
|
||||
const RealScalar mu = computeMu(M);
|
||||
MatrixType F = m_f(sigma, 0) * MatrixType::Identity(n, n);
|
||||
MatrixType P = M;
|
||||
MatrixType Fincr;
|
||||
for (int s = 1; s < 1.1*n + 10; s++) { // upper limit is fairly arbitrary
|
||||
Fincr = m_f(sigma, s) * P;
|
||||
F += Fincr;
|
||||
P = (1/(s + 1.0)) * P * M;
|
||||
if (taylorConverged(T, s, F, Fincr, P, mu)) {
|
||||
return F;
|
||||
}
|
||||
}
|
||||
ei_assert("Taylor series does not converge" && 0);
|
||||
return F;
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
typename MatrixFunction<MatrixType,1,1>::RealScalar MatrixFunction<MatrixType,1,1>::computeMu(const MatrixType& M)
|
||||
{
|
||||
const int n = M.rows();
|
||||
const MatrixType N = MatrixType::Identity(n, n) - M;
|
||||
VectorType e = VectorType::Ones(n);
|
||||
N.template triangularView<UpperTriangular>().solveInPlace(e);
|
||||
return e.cwise().abs().maxCoeff();
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
bool MatrixFunction<MatrixType,1,1>::taylorConverged(const MatrixType& T, int s, const MatrixType& F,
|
||||
const MatrixType& Fincr, const MatrixType& P, RealScalar mu)
|
||||
{
|
||||
const int n = F.rows();
|
||||
const RealScalar F_norm = F.cwise().abs().rowwise().sum().maxCoeff();
|
||||
const RealScalar Fincr_norm = Fincr.cwise().abs().rowwise().sum().maxCoeff();
|
||||
if (Fincr_norm < epsilon<Scalar>() * F_norm) {
|
||||
RealScalar delta = 0;
|
||||
RealScalar rfactorial = 1;
|
||||
for (int r = 0; r < n; r++) {
|
||||
RealScalar mx = 0;
|
||||
for (int i = 0; i < n; i++)
|
||||
mx = std::max(mx, std::abs(m_f(T(i, i), s+r)));
|
||||
if (r != 0)
|
||||
rfactorial *= r;
|
||||
delta = std::max(delta, mx / rfactorial);
|
||||
}
|
||||
const RealScalar P_norm = P.cwise().abs().rowwise().sum().maxCoeff();
|
||||
if (mu * delta * P_norm < epsilon<Scalar>() * F_norm)
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
template <typename Derived>
|
||||
EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase<Derived>& M,
|
||||
typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f,
|
||||
|
||||
142
unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h
Normal file
142
unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h
Normal file
@ -0,0 +1,142 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
||||
//
|
||||
// Eigen is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU Lesser General Public
|
||||
// License as published by the Free Software Foundation; either
|
||||
// version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Alternatively, you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of
|
||||
// the License, or (at your option) any later version.
|
||||
//
|
||||
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU Lesser General Public
|
||||
// License and a copy of the GNU General Public License along with
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#ifndef EIGEN_MATRIX_FUNCTION_ATOMIC
|
||||
#define EIGEN_MATRIX_FUNCTION_ATOMIC
|
||||
|
||||
/** \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 ei_traits<MatrixType> Traits;
|
||||
typedef typename Traits::Scalar Scalar;
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
typedef typename ei_stem_function<Scalar>::type StemFunction;
|
||||
typedef Matrix<Scalar, Traits::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(int s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P);
|
||||
|
||||
/** \brief Pointer to scalar function */
|
||||
StemFunction* m_f;
|
||||
|
||||
/** \brief Size of matrix function */
|
||||
int 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(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 (int s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary
|
||||
Fincr = m_f(m_avgEival, s) * P;
|
||||
F += Fincr;
|
||||
P = (1/(s + 1.0)) * P * m_Ashifted;
|
||||
if (taylorConverged(s, F, Fincr, P)) {
|
||||
return F;
|
||||
}
|
||||
}
|
||||
ei_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<UpperTriangular>().solveInPlace(e);
|
||||
m_mu = e.cwise().abs().maxCoeff();
|
||||
}
|
||||
|
||||
/** \brief Determine whether Taylor series has converged */
|
||||
template <typename MatrixType>
|
||||
bool MatrixFunctionAtomic<MatrixType>::taylorConverged(int s, const MatrixType& F,
|
||||
const MatrixType& Fincr, const MatrixType& P)
|
||||
{
|
||||
const int n = F.rows();
|
||||
const RealScalar F_norm = F.cwise().abs().rowwise().sum().maxCoeff();
|
||||
const RealScalar Fincr_norm = Fincr.cwise().abs().rowwise().sum().maxCoeff();
|
||||
if (Fincr_norm < epsilon<Scalar>() * F_norm) {
|
||||
RealScalar delta = 0;
|
||||
RealScalar rfactorial = 1;
|
||||
for (int r = 0; r < n; r++) {
|
||||
RealScalar mx = 0;
|
||||
for (int i = 0; i < n; i++)
|
||||
mx = std::max(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, s+r)));
|
||||
if (r != 0)
|
||||
rfactorial *= r;
|
||||
delta = std::max(delta, mx / rfactorial);
|
||||
}
|
||||
const RealScalar P_norm = P.cwise().abs().rowwise().sum().maxCoeff();
|
||||
if (m_mu * delta * P_norm < epsilon<Scalar>() * F_norm)
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
#endif // EIGEN_MATRIX_FUNCTION_ATOMIC
|
||||
Loading…
Reference in New Issue
Block a user