mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-11 19:29:02 +08:00
Clean up MatrixFunction and MatrixLogarithm.
This commit is contained in:
Jitse Niesen
parent
084dc63b4c
commit
a3a55357db
@ -1,7 +1,7 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2009-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
||||
// Copyright (C) 2009-2011, 2013 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
|
||||
@ -11,394 +11,244 @@
|
||||
#define EIGEN_MATRIX_FUNCTION
|
||||
|
||||
#include "StemFunction.h"
|
||||
#include "MatrixFunctionAtomic.h"
|
||||
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
namespace internal {
|
||||
|
||||
/** \brief Maximum distance allowed between eigenvalues to be considered "close". */
|
||||
static const float matrix_function_separation = 0.1;
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
* \brief Class for computing matrix functions.
|
||||
* \tparam MatrixType type of the argument of the matrix function,
|
||||
* expected to be an instantiation of the Matrix class template.
|
||||
* \tparam AtomicType type for computing matrix function of atomic blocks.
|
||||
* \tparam IsComplex used internally to select correct specialization.
|
||||
* \class MatrixFunctionAtomic
|
||||
* \brief Helper class for computing matrix functions of atomic matrices.
|
||||
*
|
||||
* This class implements the Schur-Parlett algorithm for computing matrix functions. The spectrum of the
|
||||
* matrix is divided in clustered of eigenvalues that lies close together. This class delegates the
|
||||
* computation of the matrix function on every block corresponding to these clusters to an object of type
|
||||
* \p AtomicType and uses these results to compute the matrix function of the whole matrix. The class
|
||||
* \p AtomicType should have a \p compute() member function for computing the matrix function of a block.
|
||||
*
|
||||
* \sa class MatrixFunctionAtomic, class MatrixLogarithmAtomic
|
||||
* Here, an atomic matrix is a triangular matrix whose diagonal entries are close to each other.
|
||||
*/
|
||||
template <typename MatrixType,
|
||||
typename AtomicType,
|
||||
int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
|
||||
class MatrixFunction : internal::noncopyable
|
||||
{
|
||||
public:
|
||||
|
||||
/** \brief Constructor.
|
||||
*
|
||||
* \param[in] A argument of matrix function, should be a square matrix.
|
||||
* \param[in] atomic class for computing matrix function of atomic blocks.
|
||||
*
|
||||
* The class stores references to \p A and \p atomic, so they should not be
|
||||
* changed (or destroyed) before compute() is called.
|
||||
*/
|
||||
MatrixFunction(const MatrixType& A, AtomicType& atomic);
|
||||
|
||||
/** \brief Compute the matrix function.
|
||||
*
|
||||
* \param[out] result the function \p f applied to \p A, as
|
||||
* specified in the constructor.
|
||||
*
|
||||
* See MatrixBase::matrixFunction() for details on how this computation
|
||||
* is implemented.
|
||||
*/
|
||||
template <typename ResultType>
|
||||
void compute(ResultType &result);
|
||||
};
|
||||
|
||||
|
||||
/** \internal \ingroup MatrixFunctions_Module
|
||||
* \brief Partial specialization of MatrixFunction for real matrices
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
class MatrixFunction<MatrixType, AtomicType, 0> : internal::noncopyable
|
||||
{
|
||||
private:
|
||||
|
||||
typedef internal::traits<MatrixType> Traits;
|
||||
typedef typename Traits::Scalar Scalar;
|
||||
static const int Rows = Traits::RowsAtCompileTime;
|
||||
static const int Cols = Traits::ColsAtCompileTime;
|
||||
static const int Options = MatrixType::Options;
|
||||
static const int MaxRows = Traits::MaxRowsAtCompileTime;
|
||||
static const int MaxCols = Traits::MaxColsAtCompileTime;
|
||||
|
||||
typedef std::complex<Scalar> ComplexScalar;
|
||||
typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
|
||||
|
||||
public:
|
||||
|
||||
/** \brief Constructor.
|
||||
*
|
||||
* \param[in] A argument of matrix function, should be a square matrix.
|
||||
* \param[in] atomic class for computing matrix function of atomic blocks.
|
||||
*/
|
||||
MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { }
|
||||
|
||||
/** \brief Compute the matrix function.
|
||||
*
|
||||
* \param[out] result the function \p f applied to \p A, as
|
||||
* specified in the constructor.
|
||||
*
|
||||
* This function converts the real matrix \c A to a complex matrix,
|
||||
* uses MatrixFunction<MatrixType,1> and then converts the result back to
|
||||
* a real matrix.
|
||||
*/
|
||||
template <typename ResultType>
|
||||
void compute(ResultType& result)
|
||||
{
|
||||
ComplexMatrix CA = m_A.template cast<ComplexScalar>();
|
||||
ComplexMatrix Cresult;
|
||||
MatrixFunction<ComplexMatrix, AtomicType> mf(CA, m_atomic);
|
||||
mf.compute(Cresult);
|
||||
result = Cresult.real();
|
||||
}
|
||||
|
||||
private:
|
||||
typename internal::nested<MatrixType>::type m_A; /**< \brief Reference to argument of matrix function. */
|
||||
AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */
|
||||
};
|
||||
|
||||
|
||||
/** \internal \ingroup MatrixFunctions_Module
|
||||
* \brief Partial specialization of MatrixFunction for complex matrices
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
class MatrixFunction<MatrixType, AtomicType, 1> : internal::noncopyable
|
||||
template <typename MatrixType>
|
||||
class MatrixFunctionAtomic
|
||||
{
|
||||
private:
|
||||
public:
|
||||
|
||||
typedef internal::traits<MatrixType> Traits;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::Index Index;
|
||||
static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
|
||||
static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
|
||||
static const int Options = MatrixType::Options;
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> VectorType;
|
||||
typedef Matrix<Index, Traits::RowsAtCompileTime, 1> IntVectorType;
|
||||
typedef Matrix<Index, Dynamic, 1> DynamicIntVectorType;
|
||||
typedef std::list<Scalar> Cluster;
|
||||
typedef std::list<Cluster> ListOfClusters;
|
||||
typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
|
||||
typedef typename stem_function<Scalar>::type StemFunction;
|
||||
|
||||
public:
|
||||
/** \brief Constructor
|
||||
* \param[in] f matrix function to compute.
|
||||
*/
|
||||
MatrixFunctionAtomic(StemFunction f) : m_f(f) { }
|
||||
|
||||
MatrixFunction(const MatrixType& A, AtomicType& atomic);
|
||||
template <typename ResultType> void compute(ResultType& result);
|
||||
/** \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:
|
||||
|
||||
void computeSchurDecomposition();
|
||||
void partitionEigenvalues();
|
||||
typename ListOfClusters::iterator findCluster(Scalar key);
|
||||
void computeClusterSize();
|
||||
void computeBlockStart();
|
||||
void constructPermutation();
|
||||
void permuteSchur();
|
||||
void swapEntriesInSchur(Index index);
|
||||
void computeBlockAtomic();
|
||||
Block<MatrixType> block(MatrixType& A, Index i, Index j);
|
||||
void computeOffDiagonal();
|
||||
DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C);
|
||||
|
||||
typename internal::nested<MatrixType>::type m_A; /**< \brief Reference to argument of matrix function. */
|
||||
AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */
|
||||
MatrixType m_T; /**< \brief Triangular part of Schur decomposition */
|
||||
MatrixType m_U; /**< \brief Unitary part of Schur decomposition */
|
||||
MatrixType m_fT; /**< \brief %Matrix function applied to #m_T */
|
||||
ListOfClusters m_clusters; /**< \brief Partition of eigenvalues into clusters of ei'vals "close" to each other */
|
||||
DynamicIntVectorType m_eivalToCluster; /**< \brief m_eivalToCluster[i] = j means i-th ei'val is in j-th cluster */
|
||||
DynamicIntVectorType m_clusterSize; /**< \brief Number of eigenvalues in each clusters */
|
||||
DynamicIntVectorType m_blockStart; /**< \brief Row index at which block corresponding to i-th cluster starts */
|
||||
IntVectorType m_permutation; /**< \brief Permutation which groups ei'vals in the same cluster together */
|
||||
|
||||
/** \brief Maximum distance allowed between eigenvalues to be considered "close".
|
||||
*
|
||||
* This is morally a \c static \c const \c Scalar, but only
|
||||
* integers can be static constant class members in C++. The
|
||||
* separation constant is set to 0.1, a value taken from the
|
||||
* paper by Davies and Higham. */
|
||||
static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
|
||||
StemFunction* m_f;
|
||||
};
|
||||
|
||||
/** \brief Constructor.
|
||||
*
|
||||
* \param[in] A argument of matrix function, should be a square matrix.
|
||||
* \param[in] atomic class for computing matrix function of atomic blocks.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
MatrixFunction<MatrixType,AtomicType,1>::MatrixFunction(const MatrixType& A, AtomicType& atomic)
|
||||
: m_A(A), m_atomic(atomic)
|
||||
template <typename MatrixType>
|
||||
typename NumTraits<typename MatrixType::Scalar>::Real matrix_function_compute_mu(const MatrixType& A)
|
||||
{
|
||||
/* empty body */
|
||||
typedef typename plain_col_type<MatrixType>::type VectorType;
|
||||
typename MatrixType::Index rows = A.rows();
|
||||
const MatrixType N = MatrixType::Identity(rows, rows) - A;
|
||||
VectorType e = VectorType::Ones(rows);
|
||||
N.template triangularView<Upper>().solveInPlace(e);
|
||||
return e.cwiseAbs().maxCoeff();
|
||||
}
|
||||
|
||||
/** \brief Compute the matrix function.
|
||||
*
|
||||
* \param[out] result the function \p f applied to \p A, as
|
||||
* specified in the constructor.
|
||||
template <typename MatrixType>
|
||||
MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
|
||||
{
|
||||
// TODO: Use that A is upper triangular
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
typedef typename MatrixType::Index Index;
|
||||
Index rows = A.rows();
|
||||
Scalar avgEival = A.trace() / Scalar(RealScalar(rows));
|
||||
MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows);
|
||||
RealScalar mu = matrix_function_compute_mu(Ashifted);
|
||||
MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows);
|
||||
MatrixType P = Ashifted;
|
||||
MatrixType Fincr;
|
||||
for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary
|
||||
Fincr = m_f(avgEival, static_cast<int>(s)) * P;
|
||||
F += Fincr;
|
||||
P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted;
|
||||
|
||||
// test whether Taylor series converged
|
||||
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 < rows; r++) {
|
||||
RealScalar mx = 0;
|
||||
for (Index i = 0; i < rows; i++)
|
||||
mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + 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 (mu * delta * P_norm < NumTraits<Scalar>::epsilon() * F_norm) // series converged
|
||||
break;
|
||||
}
|
||||
}
|
||||
return F;
|
||||
}
|
||||
|
||||
/** \brief Find cluster in \p clusters containing some value
|
||||
* \param[in] key Value to find
|
||||
* \returns Iterator to cluster containing \p key, or \c clusters.end() if no cluster in \p m_clusters
|
||||
* contains \p key.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
template <typename ResultType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::compute(ResultType& result)
|
||||
template <typename Scalar, typename ListOfClusters>
|
||||
typename ListOfClusters::iterator matrix_function_find_cluster(Scalar key, ListOfClusters& clusters)
|
||||
{
|
||||
computeSchurDecomposition();
|
||||
partitionEigenvalues();
|
||||
computeClusterSize();
|
||||
computeBlockStart();
|
||||
constructPermutation();
|
||||
permuteSchur();
|
||||
computeBlockAtomic();
|
||||
computeOffDiagonal();
|
||||
result = m_U * (m_fT.template triangularView<Upper>() * m_U.adjoint());
|
||||
}
|
||||
|
||||
/** \brief Store the Schur decomposition of #m_A in #m_T and #m_U */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::computeSchurDecomposition()
|
||||
{
|
||||
const ComplexSchur<MatrixType> schurOfA(m_A);
|
||||
m_T = schurOfA.matrixT();
|
||||
m_U = schurOfA.matrixU();
|
||||
typename std::list<Scalar>::iterator j;
|
||||
for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) {
|
||||
j = std::find(i->begin(), i->end(), key);
|
||||
if (j != i->end())
|
||||
return i;
|
||||
}
|
||||
return clusters.end();
|
||||
}
|
||||
|
||||
/** \brief Partition eigenvalues in clusters of ei'vals close to each other
|
||||
*
|
||||
* This function computes #m_clusters. This is a partition of the
|
||||
* eigenvalues of #m_T in clusters, such that
|
||||
* # Any eigenvalue in a certain cluster is at most separation() away
|
||||
* from another eigenvalue in the same cluster.
|
||||
* # The distance between two eigenvalues in different clusters is
|
||||
* more than separation().
|
||||
* The implementation follows Algorithm 4.1 in the paper of Davies
|
||||
* and Higham.
|
||||
* \param[in] eivals Eigenvalues
|
||||
* \param[out] clusters Resulting partition of eigenvalues
|
||||
*
|
||||
* The partition satisfies the following two properties:
|
||||
* # Any eigenvalue in a certain cluster is at most matrix_function_separation() away from another eigenvalue
|
||||
* in the same cluster.
|
||||
* # The distance between two eigenvalues in different clusters is more than matrix_function_separation().
|
||||
* The implementation follows Algorithm 4.1 in the paper of Davies and Higham.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::partitionEigenvalues()
|
||||
template <typename EivalsType, typename Cluster>
|
||||
void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters)
|
||||
{
|
||||
using std::abs;
|
||||
const Index rows = m_T.rows();
|
||||
VectorType diag = m_T.diagonal(); // contains eigenvalues of A
|
||||
|
||||
for (Index i=0; i<rows; ++i) {
|
||||
// Find set containing diag(i), adding a new set if necessary
|
||||
typename ListOfClusters::iterator qi = findCluster(diag(i));
|
||||
if (qi == m_clusters.end()) {
|
||||
typedef typename EivalsType::Index Index;
|
||||
for (Index i=0; i<eivals.rows(); ++i) {
|
||||
// Find set containing eivals(i), adding a new set if necessary
|
||||
typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(eivals(i), clusters);
|
||||
if (qi == clusters.end()) {
|
||||
Cluster l;
|
||||
l.push_back(diag(i));
|
||||
m_clusters.push_back(l);
|
||||
qi = m_clusters.end();
|
||||
l.push_back(eivals(i));
|
||||
clusters.push_back(l);
|
||||
qi = clusters.end();
|
||||
--qi;
|
||||
}
|
||||
|
||||
// Look for other element to add to the set
|
||||
for (Index j=i+1; j<rows; ++j) {
|
||||
if (abs(diag(j) - diag(i)) <= separation() && std::find(qi->begin(), qi->end(), diag(j)) == qi->end()) {
|
||||
typename ListOfClusters::iterator qj = findCluster(diag(j));
|
||||
if (qj == m_clusters.end()) {
|
||||
qi->push_back(diag(j));
|
||||
for (Index j=i+1; j<eivals.rows(); ++j) {
|
||||
if (abs(eivals(j) - eivals(i)) <= matrix_function_separation
|
||||
&& std::find(qi->begin(), qi->end(), eivals(j)) == qi->end()) {
|
||||
typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(eivals(j), clusters);
|
||||
if (qj == clusters.end()) {
|
||||
qi->push_back(eivals(j));
|
||||
} else {
|
||||
qi->insert(qi->end(), qj->begin(), qj->end());
|
||||
m_clusters.erase(qj);
|
||||
clusters.erase(qj);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Find cluster in #m_clusters containing some value
|
||||
* \param[in] key Value to find
|
||||
* \returns Iterator to cluster containing \c key, or
|
||||
* \c m_clusters.end() if no cluster in m_clusters contains \c key.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
typename MatrixFunction<MatrixType,AtomicType,1>::ListOfClusters::iterator MatrixFunction<MatrixType,AtomicType,1>::findCluster(Scalar key)
|
||||
/** \brief Compute size of each cluster given a partitioning */
|
||||
template <typename ListOfClusters, typename Index>
|
||||
void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix<Index, Dynamic, 1>& clusterSize)
|
||||
{
|
||||
typename Cluster::iterator j;
|
||||
for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) {
|
||||
j = std::find(i->begin(), i->end(), key);
|
||||
if (j != i->end())
|
||||
return i;
|
||||
const Index numClusters = static_cast<Index>(clusters.size());
|
||||
clusterSize.setZero(numClusters);
|
||||
Index clusterIndex = 0;
|
||||
for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
|
||||
clusterSize[clusterIndex] = cluster->size();
|
||||
++clusterIndex;
|
||||
}
|
||||
return m_clusters.end();
|
||||
}
|
||||
|
||||
/** \brief Compute #m_clusterSize and #m_eivalToCluster using #m_clusters */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::computeClusterSize()
|
||||
/** \brief Compute start of each block using clusterSize */
|
||||
template <typename VectorType>
|
||||
void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart)
|
||||
{
|
||||
const Index rows = m_T.rows();
|
||||
VectorType diag = m_T.diagonal();
|
||||
const Index numClusters = static_cast<Index>(m_clusters.size());
|
||||
blockStart.resize(clusterSize.rows());
|
||||
blockStart(0) = 0;
|
||||
for (typename VectorType::Index i = 1; i < clusterSize.rows(); i++) {
|
||||
blockStart(i) = blockStart(i-1) + clusterSize(i-1);
|
||||
}
|
||||
}
|
||||
|
||||
m_clusterSize.setZero(numClusters);
|
||||
m_eivalToCluster.resize(rows);
|
||||
/** \brief Compute mapping of eigenvalue indices to cluster indices */
|
||||
template <typename EivalsType, typename ListOfClusters, typename VectorType>
|
||||
void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster)
|
||||
{
|
||||
typedef typename EivalsType::Index Index;
|
||||
eivalToCluster.resize(eivals.rows());
|
||||
Index clusterIndex = 0;
|
||||
for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) {
|
||||
for (Index i = 0; i < diag.rows(); ++i) {
|
||||
if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) {
|
||||
++m_clusterSize[clusterIndex];
|
||||
m_eivalToCluster[i] = clusterIndex;
|
||||
for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
|
||||
for (Index i = 0; i < eivals.rows(); ++i) {
|
||||
if (std::find(cluster->begin(), cluster->end(), eivals(i)) != cluster->end()) {
|
||||
eivalToCluster[i] = clusterIndex;
|
||||
}
|
||||
}
|
||||
++clusterIndex;
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Compute #m_blockStart using #m_clusterSize */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::computeBlockStart()
|
||||
/** \brief Compute permutation which groups ei'vals in same cluster together */
|
||||
template <typename DynVectorType, typename VectorType>
|
||||
void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation)
|
||||
{
|
||||
m_blockStart.resize(m_clusterSize.rows());
|
||||
m_blockStart(0) = 0;
|
||||
for (Index i = 1; i < m_clusterSize.rows(); i++) {
|
||||
m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1);
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Compute #m_permutation using #m_eivalToCluster and #m_blockStart */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::constructPermutation()
|
||||
{
|
||||
DynamicIntVectorType indexNextEntry = m_blockStart;
|
||||
m_permutation.resize(m_T.rows());
|
||||
for (Index i = 0; i < m_T.rows(); i++) {
|
||||
Index cluster = m_eivalToCluster[i];
|
||||
m_permutation[i] = indexNextEntry[cluster];
|
||||
typedef typename VectorType::Index Index;
|
||||
DynVectorType indexNextEntry = blockStart;
|
||||
permutation.resize(eivalToCluster.rows());
|
||||
for (Index i = 0; i < eivalToCluster.rows(); i++) {
|
||||
Index cluster = eivalToCluster[i];
|
||||
permutation[i] = indexNextEntry[cluster];
|
||||
++indexNextEntry[cluster];
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Permute Schur decomposition in #m_U and #m_T according to #m_permutation */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::permuteSchur()
|
||||
/** \brief Permute Schur decomposition in U and T according to permutation */
|
||||
template <typename VectorType, typename MatrixType>
|
||||
void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T)
|
||||
{
|
||||
IntVectorType p = m_permutation;
|
||||
for (Index i = 0; i < p.rows() - 1; i++) {
|
||||
typedef typename VectorType::Index Index;
|
||||
for (Index i = 0; i < permutation.rows() - 1; i++) {
|
||||
Index j;
|
||||
for (j = i; j < p.rows(); j++) {
|
||||
if (p(j) == i) break;
|
||||
for (j = i; j < permutation.rows(); j++) {
|
||||
if (permutation(j) == i) break;
|
||||
}
|
||||
eigen_assert(p(j) == i);
|
||||
eigen_assert(permutation(j) == i);
|
||||
for (Index k = j-1; k >= i; k--) {
|
||||
swapEntriesInSchur(k);
|
||||
std::swap(p.coeffRef(k), p.coeffRef(k+1));
|
||||
JacobiRotation<typename MatrixType::Scalar> rotation;
|
||||
rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k));
|
||||
T.applyOnTheLeft(k, k+1, rotation.adjoint());
|
||||
T.applyOnTheRight(k, k+1, rotation);
|
||||
U.applyOnTheRight(k, k+1, rotation);
|
||||
std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Swap rows \a index and \a index+1 in Schur decomposition in #m_U and #m_T */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::swapEntriesInSchur(Index index)
|
||||
{
|
||||
JacobiRotation<Scalar> rotation;
|
||||
rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index));
|
||||
m_T.applyOnTheLeft(index, index+1, rotation.adjoint());
|
||||
m_T.applyOnTheRight(index, index+1, rotation);
|
||||
m_U.applyOnTheRight(index, index+1, rotation);
|
||||
}
|
||||
|
||||
/** \brief Compute block diagonal part of #m_fT.
|
||||
/** \brief Compute block diagonal part of matrix function.
|
||||
*
|
||||
* This routine computes the matrix function applied to the block diagonal part of #m_T, with the blocking
|
||||
* given by #m_blockStart. The matrix function of each diagonal block is computed by #m_atomic. The
|
||||
* off-diagonal parts of #m_fT are set to zero.
|
||||
* This routine computes the matrix function applied to the block diagonal part of \p T (which should be
|
||||
* upper triangular), with the blocking given by \p blockStart and \p clusterSize. The matrix function of
|
||||
* each diagonal block is computed by \p atomic. The off-diagonal parts of \p fT are set to zero.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::computeBlockAtomic()
|
||||
template <typename MatrixType, typename AtomicType, typename VectorType>
|
||||
void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
|
||||
{
|
||||
m_fT.resize(m_T.rows(), m_T.cols());
|
||||
m_fT.setZero();
|
||||
for (Index i = 0; i < m_clusterSize.rows(); ++i) {
|
||||
block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i));
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Return block of matrix according to blocking given by #m_blockStart */
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
Block<MatrixType> MatrixFunction<MatrixType,AtomicType,1>::block(MatrixType& A, Index i, Index j)
|
||||
{
|
||||
return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j));
|
||||
}
|
||||
|
||||
/** \brief Compute part of #m_fT above block diagonal.
|
||||
*
|
||||
* This routine assumes that the block diagonal part of #m_fT (which
|
||||
* equals the matrix function applied to #m_T) has already been computed and computes
|
||||
* the part above the block diagonal. The part below the diagonal is
|
||||
* zero, because #m_T is upper triangular.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
|
||||
{
|
||||
for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) {
|
||||
for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) {
|
||||
// compute (blockIndex, blockIndex+diagIndex) block
|
||||
DynMatrixType A = block(m_T, blockIndex, blockIndex);
|
||||
DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex);
|
||||
DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex);
|
||||
C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex);
|
||||
for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) {
|
||||
C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex);
|
||||
C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex);
|
||||
}
|
||||
block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C);
|
||||
}
|
||||
fT.setZero(T.rows(), T.cols());
|
||||
for (typename VectorType::Index i = 0; i < clusterSize.rows(); ++i) {
|
||||
fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
|
||||
= atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)));
|
||||
}
|
||||
}
|
||||
|
||||
@ -410,8 +260,8 @@ void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
|
||||
*
|
||||
* \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
|
||||
* 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
|
||||
* \f[
|
||||
* \sum_{k=i}^m A_{ik} X_{kj} + \sum_{k=1}^j X_{ik} B_{kj} = C_{ij}.
|
||||
* \f]
|
||||
@ -420,16 +270,12 @@ void MatrixFunction<MatrixType,AtomicType,1>::computeOffDiagonal()
|
||||
* X_{ij} = \frac{1}{A_{ii} + B_{jj}} \Bigl( C_{ij}
|
||||
* - \sum_{k=i+1}^m A_{ik} X_{kj} - \sum_{k=1}^{j-1} X_{ik} B_{kj} \Bigr).
|
||||
* \f]
|
||||
* It is assumed that A and B are such that the numerator is never
|
||||
* zero (otherwise the Sylvester equation does not have a unique
|
||||
* solution). In that case, these equations can be evaluated in the
|
||||
* order \f$ i=m,\ldots,1 \f$ and \f$ j=1,\ldots,n \f$.
|
||||
* It is assumed that A and B are such that the numerator is never zero (otherwise the Sylvester equation
|
||||
* does not have a unique solution). In that case, these equations can be evaluated in the order
|
||||
* \f$ i=m,\ldots,1 \f$ and \f$ j=1,\ldots,n \f$.
|
||||
*/
|
||||
template <typename MatrixType, typename AtomicType>
|
||||
typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<MatrixType,AtomicType,1>::solveTriangularSylvester(
|
||||
const DynMatrixType& A,
|
||||
const DynMatrixType& B,
|
||||
const DynMatrixType& C)
|
||||
template <typename MatrixType>
|
||||
MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C)
|
||||
{
|
||||
eigen_assert(A.rows() == A.cols());
|
||||
eigen_assert(A.isUpperTriangular());
|
||||
@ -438,9 +284,12 @@ typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<M
|
||||
eigen_assert(C.rows() == A.rows());
|
||||
eigen_assert(C.cols() == B.rows());
|
||||
|
||||
typedef typename MatrixType::Index Index;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
|
||||
Index m = A.rows();
|
||||
Index n = B.rows();
|
||||
DynMatrixType X(m, n);
|
||||
MatrixType X(m, n);
|
||||
|
||||
for (Index i = m - 1; i >= 0; --i) {
|
||||
for (Index j = 0; j < n; ++j) {
|
||||
@ -469,17 +318,164 @@ typename MatrixFunction<MatrixType,AtomicType,1>::DynMatrixType MatrixFunction<M
|
||||
return X;
|
||||
}
|
||||
|
||||
/** \brief Compute part of matrix function above block diagonal.
|
||||
*
|
||||
* This routine completes the computation of \p fT, denoting a matrix function applied to the triangular
|
||||
* matrix \p T. It assumes that the block diagonal part of \p fT has already been computed. The part below
|
||||
* the diagonal is zero, because \p T is upper triangular.
|
||||
*/
|
||||
template <typename MatrixType, typename VectorType>
|
||||
void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT)
|
||||
{
|
||||
typedef internal::traits<MatrixType> Traits;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::Index Index;
|
||||
static const int RowsAtCompileTime = Traits::RowsAtCompileTime;
|
||||
static const int ColsAtCompileTime = Traits::ColsAtCompileTime;
|
||||
static const int Options = MatrixType::Options;
|
||||
typedef Matrix<Scalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
|
||||
|
||||
for (Index k = 1; k < clusterSize.rows(); k++) {
|
||||
for (Index i = 0; i < clusterSize.rows() - k; i++) {
|
||||
// compute (i, i+k) block
|
||||
DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i));
|
||||
DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
|
||||
DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))
|
||||
* T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k));
|
||||
C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
|
||||
* fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k));
|
||||
for (Index m = i + 1; m < i + k; m++) {
|
||||
C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
|
||||
* T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
|
||||
C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m))
|
||||
* fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k));
|
||||
}
|
||||
fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k))
|
||||
= matrix_function_solve_triangular_sylvester(A, B, C);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
* \brief Class for computing matrix functions.
|
||||
* \tparam MatrixType type of the argument of the matrix function,
|
||||
* expected to be an instantiation of the Matrix class template.
|
||||
* \tparam AtomicType type for computing matrix function of atomic blocks.
|
||||
* \tparam IsComplex used internally to select correct specialization.
|
||||
*
|
||||
* This class implements the Schur-Parlett algorithm for computing matrix functions. The spectrum of the
|
||||
* matrix is divided in clustered of eigenvalues that lies close together. This class delegates the
|
||||
* computation of the matrix function on every block corresponding to these clusters to an object of type
|
||||
* \p AtomicType and uses these results to compute the matrix function of the whole matrix. The class
|
||||
* \p AtomicType should have a \p compute() member function for computing the matrix function of a block.
|
||||
*
|
||||
* \sa class MatrixFunctionAtomic, class MatrixLogarithmAtomic
|
||||
*/
|
||||
template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
|
||||
struct matrix_function_compute
|
||||
{
|
||||
/** \brief Compute the matrix function.
|
||||
*
|
||||
* \param[in] A argument of matrix function, should be a square matrix.
|
||||
* \param[in] atomic class for computing matrix function of atomic blocks.
|
||||
* \param[out] result the function \p f applied to \p A, as
|
||||
* specified in the constructor.
|
||||
*
|
||||
* See MatrixBase::matrixFunction() for details on how this computation
|
||||
* is implemented.
|
||||
*/
|
||||
template <typename AtomicType, typename ResultType>
|
||||
static void run(const MatrixType& A, AtomicType& atomic, ResultType &result);
|
||||
};
|
||||
|
||||
/** \internal \ingroup MatrixFunctions_Module
|
||||
* \brief Partial specialization of MatrixFunction for real matrices
|
||||
*
|
||||
* This converts the real matrix to a complex matrix, compute the matrix function of that matrix, and then
|
||||
* converts the result back to a real matrix.
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
struct matrix_function_compute<MatrixType, 0>
|
||||
{
|
||||
template <typename AtomicType, typename ResultType>
|
||||
static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
|
||||
{
|
||||
typedef internal::traits<MatrixType> Traits;
|
||||
typedef typename Traits::Scalar Scalar;
|
||||
static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime;
|
||||
static const int Options = MatrixType::Options;
|
||||
static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime;
|
||||
|
||||
typedef std::complex<Scalar> ComplexScalar;
|
||||
typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix;
|
||||
|
||||
ComplexMatrix CA = A.template cast<ComplexScalar>();
|
||||
ComplexMatrix Cresult;
|
||||
matrix_function_compute<ComplexMatrix>::run(CA, atomic, Cresult);
|
||||
result = Cresult.real();
|
||||
}
|
||||
};
|
||||
|
||||
/** \internal \ingroup MatrixFunctions_Module
|
||||
* \brief Partial specialization of MatrixFunction for complex matrices
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
struct matrix_function_compute<MatrixType, 1>
|
||||
{
|
||||
template <typename AtomicType, typename ResultType>
|
||||
static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
|
||||
{
|
||||
typedef internal::traits<MatrixType> Traits;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::Index Index;
|
||||
|
||||
// compute Schur decomposition of A
|
||||
const ComplexSchur<MatrixType> schurOfA(A);
|
||||
MatrixType T = schurOfA.matrixT();
|
||||
MatrixType U = schurOfA.matrixU();
|
||||
|
||||
// partition eigenvalues into clusters of ei'vals "close" to each other
|
||||
std::list<std::list<Scalar> > clusters;
|
||||
matrix_function_partition_eigenvalues(T.diagonal(), clusters);
|
||||
|
||||
// compute size of each cluster
|
||||
Matrix<Index, Dynamic, 1> clusterSize;
|
||||
matrix_function_compute_cluster_size(clusters, clusterSize);
|
||||
|
||||
// blockStart[i] is row index at which block corresponding to i-th cluster starts
|
||||
Matrix<Index, Dynamic, 1> blockStart;
|
||||
matrix_function_compute_block_start(clusterSize, blockStart);
|
||||
|
||||
// compute map so that eivalToCluster[i] = j means that ei'val T(i,i) is in j-th cluster
|
||||
Matrix<Index, Dynamic, 1> eivalToCluster;
|
||||
matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster);
|
||||
|
||||
// compute permutation which groups ei'vals in same cluster together
|
||||
Matrix<Index, Traits::RowsAtCompileTime, 1> permutation;
|
||||
matrix_function_compute_permutation(blockStart, eivalToCluster, permutation);
|
||||
|
||||
// permute Schur decomposition
|
||||
matrix_function_permute_schur(permutation, U, T);
|
||||
|
||||
// compute result
|
||||
MatrixType fT; // matrix function applied to T
|
||||
matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT);
|
||||
matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT);
|
||||
result = U * (fT.template triangularView<Upper>() * U.adjoint());
|
||||
}
|
||||
};
|
||||
|
||||
} // end of namespace internal
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
*
|
||||
* \brief Proxy for the matrix function of some matrix (expression).
|
||||
*
|
||||
* \tparam Derived Type of the argument to the matrix function.
|
||||
*
|
||||
* This class holds the argument to the matrix function until it is
|
||||
* assigned or evaluated for some other reason (so the argument
|
||||
* should not be changed in the meantime). It is the return type of
|
||||
* matrixBase::matrixFunction() and related functions and most of the
|
||||
* time this is the only way it is used.
|
||||
* This class holds the argument to the matrix function until it is assigned or evaluated for some other
|
||||
* reason (so the argument should not be changed in the meantime). It is the return type of
|
||||
* matrixBase::matrixFunction() and related functions and most of the time this is the only way it is used.
|
||||
*/
|
||||
template<typename Derived> class MatrixFunctionReturnValue
|
||||
: public ReturnByValue<MatrixFunctionReturnValue<Derived> >
|
||||
@ -490,18 +486,16 @@ template<typename Derived> class MatrixFunctionReturnValue
|
||||
typedef typename Derived::Index Index;
|
||||
typedef typename internal::stem_function<Scalar>::type StemFunction;
|
||||
|
||||
/** \brief Constructor.
|
||||
/** \brief Constructor.
|
||||
*
|
||||
* \param[in] A %Matrix (expression) forming the argument of the
|
||||
* matrix function.
|
||||
* \param[in] A %Matrix (expression) forming the argument of the matrix function.
|
||||
* \param[in] f Stem function for matrix function under consideration.
|
||||
*/
|
||||
MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { }
|
||||
|
||||
/** \brief Compute the matrix function.
|
||||
*
|
||||
* \param[out] result \p f applied to \p A, where \p f and \p A
|
||||
* are as in the constructor.
|
||||
* \param[out] result \p f applied to \p A, where \p f and \p A are as in the constructor.
|
||||
*/
|
||||
template <typename ResultType>
|
||||
inline void evalTo(ResultType& result) const
|
||||
@ -513,12 +507,12 @@ template<typename Derived> class MatrixFunctionReturnValue
|
||||
static const int Options = PlainObject::Options;
|
||||
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
|
||||
typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
|
||||
typedef MatrixFunctionAtomic<DynMatrixType> AtomicType;
|
||||
|
||||
typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
|
||||
AtomicType atomic(m_f);
|
||||
|
||||
const PlainObject Aevaluated = m_A.eval();
|
||||
MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
|
||||
mf.compute(result);
|
||||
internal::matrix_function_compute<PlainObject>::run(Aevaluated, atomic, result);
|
||||
}
|
||||
|
||||
Index rows() const { return m_A.rows(); }
|
||||
|
||||
@ -1,127 +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 : internal::noncopyable
|
||||
{
|
||||
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:
|
||||
|
||||
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
|
||||
@ -17,81 +17,30 @@
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
* \class MatrixLogarithmAtomic
|
||||
* \brief Helper class for computing matrix logarithm of atomic matrices.
|
||||
*
|
||||
* \internal
|
||||
* Here, an atomic matrix is a triangular matrix whose diagonal
|
||||
* entries are close to each other.
|
||||
*
|
||||
* \sa class MatrixFunctionAtomic, MatrixBase::log()
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
class MatrixLogarithmAtomic : internal::noncopyable
|
||||
namespace internal {
|
||||
|
||||
template <typename Scalar>
|
||||
struct matrix_log_min_pade_degree
|
||||
{
|
||||
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. */
|
||||
MatrixLogarithmAtomic() { }
|
||||
|
||||
/** \brief Compute matrix logarithm of atomic matrix
|
||||
* \param[in] A argument of matrix logarithm, should be upper triangular and atomic
|
||||
* \returns The logarithm of \p A.
|
||||
*/
|
||||
MatrixType compute(const MatrixType& A);
|
||||
|
||||
private:
|
||||
|
||||
void compute2x2(const MatrixType& A, MatrixType& result);
|
||||
void computeBig(const MatrixType& A, MatrixType& result);
|
||||
int getPadeDegree(float normTminusI);
|
||||
int getPadeDegree(double normTminusI);
|
||||
int getPadeDegree(long double normTminusI);
|
||||
void computePade(MatrixType& result, const MatrixType& T, int degree);
|
||||
void computePade3(MatrixType& result, const MatrixType& T);
|
||||
void computePade4(MatrixType& result, const MatrixType& T);
|
||||
void computePade5(MatrixType& result, const MatrixType& T);
|
||||
void computePade6(MatrixType& result, const MatrixType& T);
|
||||
void computePade7(MatrixType& result, const MatrixType& T);
|
||||
void computePade8(MatrixType& result, const MatrixType& T);
|
||||
void computePade9(MatrixType& result, const MatrixType& T);
|
||||
void computePade10(MatrixType& result, const MatrixType& T);
|
||||
void computePade11(MatrixType& result, const MatrixType& T);
|
||||
|
||||
static const int minPadeDegree = 3;
|
||||
static const int maxPadeDegree = std::numeric_limits<RealScalar>::digits<= 24? 5: // single precision
|
||||
std::numeric_limits<RealScalar>::digits<= 53? 7: // double precision
|
||||
std::numeric_limits<RealScalar>::digits<= 64? 8: // extended precision
|
||||
std::numeric_limits<RealScalar>::digits<=106? 10: // double-double
|
||||
11; // quadruple precision
|
||||
static const int value = 3;
|
||||
};
|
||||
|
||||
/** \brief Compute logarithm of triangular matrix with clustered eigenvalues. */
|
||||
template <typename MatrixType>
|
||||
MatrixType MatrixLogarithmAtomic<MatrixType>::compute(const MatrixType& A)
|
||||
template <typename Scalar>
|
||||
struct matrix_log_max_pade_degree
|
||||
{
|
||||
using std::log;
|
||||
MatrixType result(A.rows(), A.rows());
|
||||
if (A.rows() == 1)
|
||||
result(0,0) = log(A(0,0));
|
||||
else if (A.rows() == 2)
|
||||
compute2x2(A, result);
|
||||
else
|
||||
computeBig(A, result);
|
||||
return result;
|
||||
}
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
static const int value = std::numeric_limits<RealScalar>::digits<= 24? 5: // single precision
|
||||
std::numeric_limits<RealScalar>::digits<= 53? 7: // double precision
|
||||
std::numeric_limits<RealScalar>::digits<= 64? 8: // extended precision
|
||||
std::numeric_limits<RealScalar>::digits<=106? 10: // double-double
|
||||
11; // quadruple precision
|
||||
};
|
||||
|
||||
/** \brief Compute logarithm of 2x2 triangular matrix. */
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixType& result)
|
||||
void matrix_log_compute_2x2(const MatrixType& A, MatrixType& result)
|
||||
{
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
using std::abs;
|
||||
using std::ceil;
|
||||
using std::imag;
|
||||
@ -116,47 +65,14 @@ void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixTy
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Compute logarithm of triangular matrices with size > 2.
|
||||
* \details This uses a inverse scale-and-square algorithm. */
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computeBig(const MatrixType& A, MatrixType& result)
|
||||
{
|
||||
using std::pow;
|
||||
int numberOfSquareRoots = 0;
|
||||
int numberOfExtraSquareRoots = 0;
|
||||
int degree;
|
||||
MatrixType T = A, sqrtT;
|
||||
const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1: // single precision
|
||||
maxPadeDegree<= 7? 2.6429608311114350e-1: // double precision
|
||||
maxPadeDegree<= 8? 2.32777776523703892094e-1L: // extended precision
|
||||
maxPadeDegree<=10? 1.05026503471351080481093652651105e-1L: // double-double
|
||||
1.1880960220216759245467951592883642e-1L; // quadruple precision
|
||||
|
||||
while (true) {
|
||||
RealScalar normTminusI = (T - MatrixType::Identity(T.rows(), T.rows())).cwiseAbs().colwise().sum().maxCoeff();
|
||||
if (normTminusI < maxNormForPade) {
|
||||
degree = getPadeDegree(normTminusI);
|
||||
int degree2 = getPadeDegree(normTminusI / RealScalar(2));
|
||||
if ((degree - degree2 <= 1) || (numberOfExtraSquareRoots == 1))
|
||||
break;
|
||||
++numberOfExtraSquareRoots;
|
||||
}
|
||||
matrix_sqrt_triangular(T, sqrtT);
|
||||
T = sqrtT.template triangularView<Upper>();
|
||||
++numberOfSquareRoots;
|
||||
}
|
||||
|
||||
computePade(result, T, degree);
|
||||
result *= pow(RealScalar(2), numberOfSquareRoots);
|
||||
}
|
||||
|
||||
/* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = float) */
|
||||
template <typename MatrixType>
|
||||
int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(float normTminusI)
|
||||
int matrix_log_get_pade_degree(float normTminusI)
|
||||
{
|
||||
const float maxNormForPade[] = { 2.5111573934555054e-1 /* degree = 3 */ , 4.0535837411880493e-1,
|
||||
5.3149729967117310e-1 };
|
||||
int degree = 3;
|
||||
const int minPadeDegree = matrix_log_min_pade_degree<float>::value;
|
||||
const int maxPadeDegree = matrix_log_max_pade_degree<float>::value;
|
||||
int degree = minPadeDegree;
|
||||
for (; degree <= maxPadeDegree; ++degree)
|
||||
if (normTminusI <= maxNormForPade[degree - minPadeDegree])
|
||||
break;
|
||||
@ -164,12 +80,13 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(float normTminusI)
|
||||
}
|
||||
|
||||
/* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */
|
||||
template <typename MatrixType>
|
||||
int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(double normTminusI)
|
||||
int matrix_log_get_pade_degree(double normTminusI)
|
||||
{
|
||||
const double maxNormForPade[] = { 1.6206284795015624e-2 /* degree = 3 */ , 5.3873532631381171e-2,
|
||||
1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 };
|
||||
int degree = 3;
|
||||
const int minPadeDegree = matrix_log_min_pade_degree<double>::value;
|
||||
const int maxPadeDegree = matrix_log_max_pade_degree<double>::value;
|
||||
int degree = minPadeDegree;
|
||||
for (; degree <= maxPadeDegree; ++degree)
|
||||
if (normTminusI <= maxNormForPade[degree - minPadeDegree])
|
||||
break;
|
||||
@ -177,8 +94,7 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(double normTminusI)
|
||||
}
|
||||
|
||||
/* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */
|
||||
template <typename MatrixType>
|
||||
int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(long double normTminusI)
|
||||
int matrix_log_get_pade_degree(long double normTminusI)
|
||||
{
|
||||
#if LDBL_MANT_DIG == 53 // double precision
|
||||
const long double maxNormForPade[] = { 1.6206284795015624e-2L /* degree = 3 */ , 5.3873532631381171e-2L,
|
||||
@ -200,7 +116,9 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(long double normTminusI)
|
||||
3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L,
|
||||
8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L };
|
||||
#endif
|
||||
int degree = 3;
|
||||
const int minPadeDegree = matrix_log_min_pade_degree<long double>::value;
|
||||
const int maxPadeDegree = matrix_log_max_pade_degree<long double>::value;
|
||||
int degree = minPadeDegree;
|
||||
for (; degree <= maxPadeDegree; ++degree)
|
||||
if (normTminusI <= maxNormForPade[degree - minPadeDegree])
|
||||
break;
|
||||
@ -209,197 +127,168 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(long double normTminusI)
|
||||
|
||||
/* \brief Compute Pade approximation to matrix logarithm */
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade(MatrixType& result, const MatrixType& T, int degree)
|
||||
void matrix_log_compute_pade(MatrixType& result, const MatrixType& T, int degree)
|
||||
{
|
||||
switch (degree) {
|
||||
case 3: computePade3(result, T); break;
|
||||
case 4: computePade4(result, T); break;
|
||||
case 5: computePade5(result, T); break;
|
||||
case 6: computePade6(result, T); break;
|
||||
case 7: computePade7(result, T); break;
|
||||
case 8: computePade8(result, T); break;
|
||||
case 9: computePade9(result, T); break;
|
||||
case 10: computePade10(result, T); break;
|
||||
case 11: computePade11(result, T); break;
|
||||
default: assert(false); // should never happen
|
||||
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
|
||||
const int minPadeDegree = 3;
|
||||
const int maxPadeDegree = 11;
|
||||
assert(degree >= minPadeDegree && degree <= maxPadeDegree);
|
||||
|
||||
const RealScalar nodes[][maxPadeDegree] = {
|
||||
{ 0.1127016653792583114820734600217600L, 0.5000000000000000000000000000000000L, // degree 3
|
||||
0.8872983346207416885179265399782400L },
|
||||
{ 0.0694318442029737123880267555535953L, 0.3300094782075718675986671204483777L, // degree 4
|
||||
0.6699905217924281324013328795516223L, 0.9305681557970262876119732444464048L },
|
||||
{ 0.0469100770306680036011865608503035L, 0.2307653449471584544818427896498956L, // degree 5
|
||||
0.5000000000000000000000000000000000L, 0.7692346550528415455181572103501044L,
|
||||
0.9530899229693319963988134391496965L },
|
||||
{ 0.0337652428984239860938492227530027L, 0.1693953067668677431693002024900473L, // degree 6
|
||||
0.3806904069584015456847491391596440L, 0.6193095930415984543152508608403560L,
|
||||
0.8306046932331322568306997975099527L, 0.9662347571015760139061507772469973L },
|
||||
{ 0.0254460438286207377369051579760744L, 0.1292344072003027800680676133596058L, // degree 7
|
||||
0.2970774243113014165466967939615193L, 0.5000000000000000000000000000000000L,
|
||||
0.7029225756886985834533032060384807L, 0.8707655927996972199319323866403942L,
|
||||
0.9745539561713792622630948420239256L },
|
||||
{ 0.0198550717512318841582195657152635L, 0.1016667612931866302042230317620848L, // degree 8
|
||||
0.2372337950418355070911304754053768L, 0.4082826787521750975302619288199080L,
|
||||
0.5917173212478249024697380711800920L, 0.7627662049581644929088695245946232L,
|
||||
0.8983332387068133697957769682379152L, 0.9801449282487681158417804342847365L },
|
||||
{ 0.0159198802461869550822118985481636L, 0.0819844463366821028502851059651326L, // degree 9
|
||||
0.1933142836497048013456489803292629L, 0.3378732882980955354807309926783317L,
|
||||
0.5000000000000000000000000000000000L, 0.6621267117019044645192690073216683L,
|
||||
0.8066857163502951986543510196707371L, 0.9180155536633178971497148940348674L,
|
||||
0.9840801197538130449177881014518364L },
|
||||
{ 0.0130467357414141399610179939577740L, 0.0674683166555077446339516557882535L, // degree 10
|
||||
0.1602952158504877968828363174425632L, 0.2833023029353764046003670284171079L,
|
||||
0.4255628305091843945575869994351400L, 0.5744371694908156054424130005648600L,
|
||||
0.7166976970646235953996329715828921L, 0.8397047841495122031171636825574368L,
|
||||
0.9325316833444922553660483442117465L, 0.9869532642585858600389820060422260L },
|
||||
{ 0.0108856709269715035980309994385713L, 0.0564687001159523504624211153480364L, // degree 11
|
||||
0.1349239972129753379532918739844233L, 0.2404519353965940920371371652706952L,
|
||||
0.3652284220238275138342340072995692L, 0.5000000000000000000000000000000000L,
|
||||
0.6347715779761724861657659927004308L, 0.7595480646034059079628628347293048L,
|
||||
0.8650760027870246620467081260155767L, 0.9435312998840476495375788846519636L,
|
||||
0.9891143290730284964019690005614287L } };
|
||||
|
||||
const RealScalar weights[][maxPadeDegree] = {
|
||||
{ 0.2777777777777777777777777777777778L, 0.4444444444444444444444444444444444L, // degree 3
|
||||
0.2777777777777777777777777777777778L },
|
||||
{ 0.1739274225687269286865319746109997L, 0.3260725774312730713134680253890003L, // degree 4
|
||||
0.3260725774312730713134680253890003L, 0.1739274225687269286865319746109997L },
|
||||
{ 0.1184634425280945437571320203599587L, 0.2393143352496832340206457574178191L, // degree 5
|
||||
0.2844444444444444444444444444444444L, 0.2393143352496832340206457574178191L,
|
||||
0.1184634425280945437571320203599587L },
|
||||
{ 0.0856622461895851725201480710863665L, 0.1803807865240693037849167569188581L, // degree 6
|
||||
0.2339569672863455236949351719947755L, 0.2339569672863455236949351719947755L,
|
||||
0.1803807865240693037849167569188581L, 0.0856622461895851725201480710863665L },
|
||||
{ 0.0647424830844348466353057163395410L, 0.1398526957446383339507338857118898L, // degree 7
|
||||
0.1909150252525594724751848877444876L, 0.2089795918367346938775510204081633L,
|
||||
0.1909150252525594724751848877444876L, 0.1398526957446383339507338857118898L,
|
||||
0.0647424830844348466353057163395410L },
|
||||
{ 0.0506142681451881295762656771549811L, 0.1111905172266872352721779972131204L, // degree 8
|
||||
0.1568533229389436436689811009933007L, 0.1813418916891809914825752246385978L,
|
||||
0.1813418916891809914825752246385978L, 0.1568533229389436436689811009933007L,
|
||||
0.1111905172266872352721779972131204L, 0.0506142681451881295762656771549811L },
|
||||
{ 0.0406371941807872059859460790552618L, 0.0903240803474287020292360156214564L, // degree 9
|
||||
0.1303053482014677311593714347093164L, 0.1561735385200014200343152032922218L,
|
||||
0.1651196775006298815822625346434870L, 0.1561735385200014200343152032922218L,
|
||||
0.1303053482014677311593714347093164L, 0.0903240803474287020292360156214564L,
|
||||
0.0406371941807872059859460790552618L },
|
||||
{ 0.0333356721543440687967844049466659L, 0.0747256745752902965728881698288487L, // degree 10
|
||||
0.1095431812579910219977674671140816L, 0.1346333596549981775456134607847347L,
|
||||
0.1477621123573764350869464973256692L, 0.1477621123573764350869464973256692L,
|
||||
0.1346333596549981775456134607847347L, 0.1095431812579910219977674671140816L,
|
||||
0.0747256745752902965728881698288487L, 0.0333356721543440687967844049466659L },
|
||||
{ 0.0278342835580868332413768602212743L, 0.0627901847324523123173471496119701L, // degree 11
|
||||
0.0931451054638671257130488207158280L, 0.1165968822959952399592618524215876L,
|
||||
0.1314022722551233310903444349452546L, 0.1364625433889503153572417641681711L,
|
||||
0.1314022722551233310903444349452546L, 0.1165968822959952399592618524215876L,
|
||||
0.0931451054638671257130488207158280L, 0.0627901847324523123173471496119701L,
|
||||
0.0278342835580868332413768602212743L } };
|
||||
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k) {
|
||||
RealScalar weight = weights[degree-minPadeDegree][k];
|
||||
RealScalar node = nodes[degree-minPadeDegree][k];
|
||||
result += weight * (MatrixType::Identity(T.rows(), T.rows()) + node * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
}
|
||||
|
||||
/** \brief Compute logarithm of triangular matrices with size > 2.
|
||||
* \details This uses a inverse scale-and-square algorithm. */
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade3(MatrixType& result, const MatrixType& T)
|
||||
void matrix_log_compute_big(const MatrixType& A, MatrixType& result)
|
||||
{
|
||||
const int degree = 3;
|
||||
const RealScalar nodes[] = { 0.1127016653792583114820734600217600L, 0.5000000000000000000000000000000000L,
|
||||
0.8872983346207416885179265399782400L };
|
||||
const RealScalar weights[] = { 0.2777777777777777777777777777777778L, 0.4444444444444444444444444444444444L,
|
||||
0.2777777777777777777777777777777778L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
using std::pow;
|
||||
|
||||
int numberOfSquareRoots = 0;
|
||||
int numberOfExtraSquareRoots = 0;
|
||||
int degree;
|
||||
MatrixType T = A, sqrtT;
|
||||
|
||||
int maxPadeDegree = matrix_log_max_pade_degree<Scalar>::value;
|
||||
const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1: // single precision
|
||||
maxPadeDegree<= 7? 2.6429608311114350e-1: // double precision
|
||||
maxPadeDegree<= 8? 2.32777776523703892094e-1L: // extended precision
|
||||
maxPadeDegree<=10? 1.05026503471351080481093652651105e-1L: // double-double
|
||||
1.1880960220216759245467951592883642e-1L; // quadruple precision
|
||||
|
||||
while (true) {
|
||||
RealScalar normTminusI = (T - MatrixType::Identity(T.rows(), T.rows())).cwiseAbs().colwise().sum().maxCoeff();
|
||||
if (normTminusI < maxNormForPade) {
|
||||
degree = matrix_log_get_pade_degree(normTminusI);
|
||||
int degree2 = matrix_log_get_pade_degree(normTminusI / RealScalar(2));
|
||||
if ((degree - degree2 <= 1) || (numberOfExtraSquareRoots == 1))
|
||||
break;
|
||||
++numberOfExtraSquareRoots;
|
||||
}
|
||||
matrix_sqrt_triangular(T, sqrtT);
|
||||
T = sqrtT.template triangularView<Upper>();
|
||||
++numberOfSquareRoots;
|
||||
}
|
||||
|
||||
matrix_log_compute_pade(result, T, degree);
|
||||
result *= pow(RealScalar(2), numberOfSquareRoots);
|
||||
}
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
* \class MatrixLogarithmAtomic
|
||||
* \brief Helper class for computing matrix logarithm of atomic matrices.
|
||||
*
|
||||
* Here, an atomic matrix is a triangular matrix whose diagonal entries are close to each other.
|
||||
*
|
||||
* \sa class MatrixFunctionAtomic, MatrixBase::log()
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade4(MatrixType& result, const MatrixType& T)
|
||||
class MatrixLogarithmAtomic
|
||||
{
|
||||
const int degree = 4;
|
||||
const RealScalar nodes[] = { 0.0694318442029737123880267555535953L, 0.3300094782075718675986671204483777L,
|
||||
0.6699905217924281324013328795516223L, 0.9305681557970262876119732444464048L };
|
||||
const RealScalar weights[] = { 0.1739274225687269286865319746109997L, 0.3260725774312730713134680253890003L,
|
||||
0.3260725774312730713134680253890003L, 0.1739274225687269286865319746109997L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
public:
|
||||
/** \brief Compute matrix logarithm of atomic matrix
|
||||
* \param[in] A argument of matrix logarithm, should be upper triangular and atomic
|
||||
* \returns The logarithm of \p A.
|
||||
*/
|
||||
MatrixType compute(const MatrixType& A);
|
||||
};
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade5(MatrixType& result, const MatrixType& T)
|
||||
MatrixType MatrixLogarithmAtomic<MatrixType>::compute(const MatrixType& A)
|
||||
{
|
||||
const int degree = 5;
|
||||
const RealScalar nodes[] = { 0.0469100770306680036011865608503035L, 0.2307653449471584544818427896498956L,
|
||||
0.5000000000000000000000000000000000L, 0.7692346550528415455181572103501044L,
|
||||
0.9530899229693319963988134391496965L };
|
||||
const RealScalar weights[] = { 0.1184634425280945437571320203599587L, 0.2393143352496832340206457574178191L,
|
||||
0.2844444444444444444444444444444444L, 0.2393143352496832340206457574178191L,
|
||||
0.1184634425280945437571320203599587L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
using std::log;
|
||||
MatrixType result(A.rows(), A.rows());
|
||||
if (A.rows() == 1)
|
||||
result(0,0) = log(A(0,0));
|
||||
else if (A.rows() == 2)
|
||||
matrix_log_compute_2x2(A, result);
|
||||
else
|
||||
matrix_log_compute_big(A, result);
|
||||
return result;
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade6(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 6;
|
||||
const RealScalar nodes[] = { 0.0337652428984239860938492227530027L, 0.1693953067668677431693002024900473L,
|
||||
0.3806904069584015456847491391596440L, 0.6193095930415984543152508608403560L,
|
||||
0.8306046932331322568306997975099527L, 0.9662347571015760139061507772469973L };
|
||||
const RealScalar weights[] = { 0.0856622461895851725201480710863665L, 0.1803807865240693037849167569188581L,
|
||||
0.2339569672863455236949351719947755L, 0.2339569672863455236949351719947755L,
|
||||
0.1803807865240693037849167569188581L, 0.0856622461895851725201480710863665L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade7(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 7;
|
||||
const RealScalar nodes[] = { 0.0254460438286207377369051579760744L, 0.1292344072003027800680676133596058L,
|
||||
0.2970774243113014165466967939615193L, 0.5000000000000000000000000000000000L,
|
||||
0.7029225756886985834533032060384807L, 0.8707655927996972199319323866403942L,
|
||||
0.9745539561713792622630948420239256L };
|
||||
const RealScalar weights[] = { 0.0647424830844348466353057163395410L, 0.1398526957446383339507338857118898L,
|
||||
0.1909150252525594724751848877444876L, 0.2089795918367346938775510204081633L,
|
||||
0.1909150252525594724751848877444876L, 0.1398526957446383339507338857118898L,
|
||||
0.0647424830844348466353057163395410L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade8(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 8;
|
||||
const RealScalar nodes[] = { 0.0198550717512318841582195657152635L, 0.1016667612931866302042230317620848L,
|
||||
0.2372337950418355070911304754053768L, 0.4082826787521750975302619288199080L,
|
||||
0.5917173212478249024697380711800920L, 0.7627662049581644929088695245946232L,
|
||||
0.8983332387068133697957769682379152L, 0.9801449282487681158417804342847365L };
|
||||
const RealScalar weights[] = { 0.0506142681451881295762656771549811L, 0.1111905172266872352721779972131204L,
|
||||
0.1568533229389436436689811009933007L, 0.1813418916891809914825752246385978L,
|
||||
0.1813418916891809914825752246385978L, 0.1568533229389436436689811009933007L,
|
||||
0.1111905172266872352721779972131204L, 0.0506142681451881295762656771549811L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade9(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 9;
|
||||
const RealScalar nodes[] = { 0.0159198802461869550822118985481636L, 0.0819844463366821028502851059651326L,
|
||||
0.1933142836497048013456489803292629L, 0.3378732882980955354807309926783317L,
|
||||
0.5000000000000000000000000000000000L, 0.6621267117019044645192690073216683L,
|
||||
0.8066857163502951986543510196707371L, 0.9180155536633178971497148940348674L,
|
||||
0.9840801197538130449177881014518364L };
|
||||
const RealScalar weights[] = { 0.0406371941807872059859460790552618L, 0.0903240803474287020292360156214564L,
|
||||
0.1303053482014677311593714347093164L, 0.1561735385200014200343152032922218L,
|
||||
0.1651196775006298815822625346434870L, 0.1561735385200014200343152032922218L,
|
||||
0.1303053482014677311593714347093164L, 0.0903240803474287020292360156214564L,
|
||||
0.0406371941807872059859460790552618L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade10(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 10;
|
||||
const RealScalar nodes[] = { 0.0130467357414141399610179939577740L, 0.0674683166555077446339516557882535L,
|
||||
0.1602952158504877968828363174425632L, 0.2833023029353764046003670284171079L,
|
||||
0.4255628305091843945575869994351400L, 0.5744371694908156054424130005648600L,
|
||||
0.7166976970646235953996329715828921L, 0.8397047841495122031171636825574368L,
|
||||
0.9325316833444922553660483442117465L, 0.9869532642585858600389820060422260L };
|
||||
const RealScalar weights[] = { 0.0333356721543440687967844049466659L, 0.0747256745752902965728881698288487L,
|
||||
0.1095431812579910219977674671140816L, 0.1346333596549981775456134607847347L,
|
||||
0.1477621123573764350869464973256692L, 0.1477621123573764350869464973256692L,
|
||||
0.1346333596549981775456134607847347L, 0.1095431812579910219977674671140816L,
|
||||
0.0747256745752902965728881698288487L, 0.0333356721543440687967844049466659L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
|
||||
template <typename MatrixType>
|
||||
void MatrixLogarithmAtomic<MatrixType>::computePade11(MatrixType& result, const MatrixType& T)
|
||||
{
|
||||
const int degree = 11;
|
||||
const RealScalar nodes[] = { 0.0108856709269715035980309994385713L, 0.0564687001159523504624211153480364L,
|
||||
0.1349239972129753379532918739844233L, 0.2404519353965940920371371652706952L,
|
||||
0.3652284220238275138342340072995692L, 0.5000000000000000000000000000000000L,
|
||||
0.6347715779761724861657659927004308L, 0.7595480646034059079628628347293048L,
|
||||
0.8650760027870246620467081260155767L, 0.9435312998840476495375788846519636L,
|
||||
0.9891143290730284964019690005614287L };
|
||||
const RealScalar weights[] = { 0.0278342835580868332413768602212743L, 0.0627901847324523123173471496119701L,
|
||||
0.0931451054638671257130488207158280L, 0.1165968822959952399592618524215876L,
|
||||
0.1314022722551233310903444349452546L, 0.1364625433889503153572417641681711L,
|
||||
0.1314022722551233310903444349452546L, 0.1165968822959952399592618524215876L,
|
||||
0.0931451054638671257130488207158280L, 0.0627901847324523123173471496119701L,
|
||||
0.0278342835580868332413768602212743L };
|
||||
eigen_assert(degree <= maxPadeDegree);
|
||||
MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows());
|
||||
result.setZero(T.rows(), T.rows());
|
||||
for (int k = 0; k < degree; ++k)
|
||||
result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI)
|
||||
.template triangularView<Upper>().solve(TminusI);
|
||||
}
|
||||
} // end of namespace internal
|
||||
|
||||
/** \ingroup MatrixFunctions_Module
|
||||
*
|
||||
@ -441,12 +330,11 @@ public:
|
||||
static const int Options = PlainObject::Options;
|
||||
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
|
||||
typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType;
|
||||
typedef MatrixLogarithmAtomic<DynMatrixType> AtomicType;
|
||||
typedef internal::MatrixLogarithmAtomic<DynMatrixType> AtomicType;
|
||||
AtomicType atomic;
|
||||
|
||||
const PlainObject Aevaluated = m_A.eval();
|
||||
MatrixFunction<PlainObject, AtomicType> mf(Aevaluated, atomic);
|
||||
mf.compute(result);
|
||||
internal::matrix_function_compute<PlainObject>::run(Aevaluated, atomic, result);
|
||||
}
|
||||
|
||||
Index rows() const { return m_A.rows(); }
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
||||
// Copyright (C) 2011, 2013 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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user