Eigen::MatrixPower< MatrixType > Class Template Reference

Class for computing matrix powers. More...

#include <MatrixPower.h>

Inheritance diagram for Eigen::MatrixPower< MatrixType >:

Public Member Functions
	MatrixPower (const MatrixType &A)
	Constructor. More...
const MatrixPowerParenthesesReturnValue< MatrixType >	operator() (RealScalar p)
	Returns the matrix power. More...
template<typename ResultType >
void	compute (ResultType &res, RealScalar p)
	Compute the matrix power. More...
Index	rows () const
Index	cols () const

Private Types
typedef MatrixType::Scalar	Scalar
typedef MatrixType::RealScalar	RealScalar
typedef std::complex< RealScalar >	ComplexScalar
typedef Matrix< ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime >	ComplexMatrix

Private Member Functions
void	split (RealScalar &p, RealScalar &intpart)
	Split p into integral part and fractional part. More...
void	initialize ()
	Perform Schur decomposition for fractional power. More...
template<typename ResultType >
void	computeIntPower (ResultType &res, RealScalar p)
template<typename ResultType >
void	computeFracPower (ResultType &res, RealScalar p)
Private Member Functions inherited from Eigen::internal::noncopyable
EIGEN_DEVICE_FUNC	noncopyable ()
EIGEN_DEVICE_FUNC	~noncopyable ()

Static Private Member Functions
template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
static void	revertSchur (Matrix< ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols > &res, const ComplexMatrix &T, const ComplexMatrix &U)
template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
static void	revertSchur (Matrix< RealScalar, Rows, Cols, Options, MaxRows, MaxCols > &res, const ComplexMatrix &T, const ComplexMatrix &U)

Private Attributes
MatrixType::Nested	m_A
	Reference to the base of matrix power. More...
MatrixType	m_tmp
	Temporary storage. More...
ComplexMatrix	m_T
	Store the result of Schur decomposition. More...
ComplexMatrix	m_U
ComplexMatrix	m_fT
	Store fractional power of m_T. More...
RealScalar	m_conditionNumber
	Condition number of m_A. More...
Index	m_rank
	Rank of m_A. More...
Index	m_nulls
	Rank deficiency of m_A. More...

Detailed Description

template<typename MatrixType>
class Eigen::MatrixPower< MatrixType >

Class for computing matrix powers.

Template Parameters

MatrixType

type of the base, expected to be an instantiation of the Matrix class template.

This class is capable of computing real/complex matrices raised to an arbitrary real power. Meanwhile, it saves the result of Schur decomposition if an non-integral power has even been calculated. Therefore, if you want to compute multiple (>= 2) matrix powers for the same matrix, using the class directly is more efficient than calling MatrixBase::pow().

Example:

#include <unsupported/Eigen/MatrixFunctions>
#include <iostream>
 
using namespace Eigen;
 
int main() {
  Matrix4cd A = Matrix4cd::Random();
  MatrixPower<Matrix4cd> Apow(A);
 
  std::cout << "The matrix A is:\n"
            << A
            << "\n\n"
               "A^3.1 is:\n"
            << Apow(3.1)
            << "\n\n"
               "A^3.3 is:\n"
            << Apow(3.3)
            << "\n\n"
               "A^3.7 is:\n"
            << Apow(3.7)
            << "\n\n"
               "A^3.9 is:\n"
            << Apow(3.9) << std::endl;
  return 0;
}

Output:

Member Typedef Documentation

ComplexMatrix

template<typename MatrixType >

typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> Eigen::MatrixPower< MatrixType >::ComplexMatrix

private

ComplexScalar

template<typename MatrixType >

typedef std::complex<RealScalar> Eigen::MatrixPower< MatrixType >::ComplexScalar

private

RealScalar

template<typename MatrixType >

typedef MatrixType::RealScalar Eigen::MatrixPower< MatrixType >::RealScalar

private

Scalar

template<typename MatrixType >

typedef MatrixType::Scalar Eigen::MatrixPower< MatrixType >::Scalar

private

Constructor & Destructor Documentation

MatrixPower()

template<typename MatrixType >

Eigen::MatrixPower< MatrixType >::MatrixPower ( const MatrixType &  A )

inlineexplicit

Constructor.

Parameters

[in]

A

the base of the matrix power.

The class stores a reference to A, so it should not be changed (or destroyed) before evaluation.

                                            : m_A(A), m_conditionNumber(0), m_rank(A.cols()), m_nulls(0) {
    eigen_assert(A.rows() == A.cols());
  }

References Eigen::PlainObjectBase< Derived >::cols(), eigen_assert, and Eigen::PlainObjectBase< Derived >::rows().

Member Function Documentation

cols()

template<typename MatrixType >

Index Eigen::MatrixPower< MatrixType >::cols ( ) const

inline

{ return m_A.cols(); }

References Eigen::MatrixPower< MatrixType >::m_A.

Referenced by gdb.printers._MatrixEntryIterator::__next__(), gdb.printers.EigenMatrixPrinter::children(), gdb.printers.EigenSparseMatrixPrinter::children(), gdb.printers.EigenMatrixPrinter::to_string(), and gdb.printers.EigenSparseMatrixPrinter::to_string().

compute()

template<typename MatrixType >

template<typename ResultType >

void Eigen::MatrixPower< MatrixType >::compute	(	ResultType &	res,
		RealScalar	p
	)

Compute the matrix power.

Parameters

[in]	p	exponent, a real scalar.
[out]	res	\( A^p \) where A is specified in the constructor.

                                                                   {
  using std::pow;
  switch (cols()) {
    case 0:
      break;
    case 1:
      res(0, 0) = pow(m_A.coeff(0, 0), p);
      break;
    default:
      RealScalar intpart;
      split(p, intpart);
 
      res = MatrixType::Identity(rows(), cols());
      computeIntPower(res, intpart);
      if (p) computeFracPower(res, p);
  }
}

References cols, p, Eigen::bfloat16_impl::pow(), res, rows, and oomph::DoubleVectorHelpers::split().

Referenced by Eigen::MatrixPowerReturnValue< Derived >::evalTo().

computeFracPower()

template<typename MatrixType >

template<typename ResultType >

void Eigen::MatrixPower< MatrixType >::computeFracPower ( ResultType &  res, RealScalar  p  )

private

                                                                            {
  Block<ComplexMatrix, Dynamic, Dynamic> blockTp(m_fT, 0, 0, m_rank, m_rank);
  eigen_assert(m_conditionNumber);
  eigen_assert(m_rank + m_nulls == rows());
 
  MatrixPowerAtomic<ComplexMatrix>(m_T.topLeftCorner(m_rank, m_rank), p).compute(blockTp);
  if (m_nulls) {
    m_fT.topRightCorner(m_rank, m_nulls) = m_T.topLeftCorner(m_rank, m_rank)
                                               .template triangularView<Upper>()
                                               .solve(blockTp * m_T.topRightCorner(m_rank, m_nulls));
  }
  revertSchur(m_tmp, m_fT, m_U);
  res = m_tmp * res;
}

References Eigen::MatrixPowerAtomic< MatrixType >::compute(), eigen_assert, p, res, and rows.

computeIntPower()

template<typename MatrixType >

template<typename ResultType >

void Eigen::MatrixPower< MatrixType >::computeIntPower ( ResultType &  res, RealScalar  p  )

private

                                                                           {
  using std::abs;
  using std::fmod;
  RealScalar pp = abs(p);
 
  if (p < 0)
    m_tmp = m_A.inverse();
  else
    m_tmp = m_A;
 
  while (true) {
    if (fmod(pp, 2) >= 1) res = m_tmp * res;
    pp /= 2;
    if (pp < 1) break;
    m_tmp *= m_tmp;
  }
}

References abs(), Eigen::bfloat16_impl::fmod(), p, and res.

initialize()

template<typename MatrixType >

void Eigen::MatrixPower< MatrixType >::initialize

private

Perform Schur decomposition for fractional power.

                                         {
  const ComplexSchur<MatrixType> schurOfA(m_A);
  JacobiRotation<ComplexScalar> rot;
  ComplexScalar eigenvalue;
 
  m_fT.resizeLike(m_A);
  m_T = schurOfA.matrixT();
  m_U = schurOfA.matrixU();
  m_conditionNumber = m_T.diagonal().array().abs().maxCoeff() / m_T.diagonal().array().abs().minCoeff();
 
  // Move zero eigenvalues to the bottom right corner.
  for (Index i = cols() - 1; i >= 0; --i) {
    if (m_rank <= 2) return;
    if (m_T.coeff(i, i) == RealScalar(0)) {
      for (Index j = i + 1; j < m_rank; ++j) {
        eigenvalue = m_T.coeff(j, j);
        rot.makeGivens(m_T.coeff(j - 1, j), eigenvalue);
        m_T.applyOnTheRight(j - 1, j, rot);
        m_T.applyOnTheLeft(j - 1, j, rot.adjoint());
        m_T.coeffRef(j - 1, j - 1) = eigenvalue;
        m_T.coeffRef(j, j) = RealScalar(0);
        m_U.applyOnTheRight(j - 1, j, rot);
      }
      --m_rank;
    }
  }
 
  m_nulls = rows() - m_rank;
  if (m_nulls) {
    eigen_assert(m_T.bottomRightCorner(m_nulls, m_nulls).isZero() &&
                 "Base of matrix power should be invertible or with a semisimple zero eigenvalue.");
    m_fT.bottomRows(m_nulls).fill(RealScalar(0));
  }
}

References cols, eigen_assert, i, j, rot(), rows, and schurOfA().

operator()()

template<typename MatrixType >

const MatrixPowerParenthesesReturnValue<MatrixType> Eigen::MatrixPower< MatrixType >::operator() ( RealScalar  p )

inline

Returns the matrix power.

Parameters

[in]

p

exponent, a real scalar.

Returns

The expression \( A^p \), where A is specified in the constructor.

                                                                               {
    return MatrixPowerParenthesesReturnValue<MatrixType>(*this, p);
  }

References p.

revertSchur() [1/2]

template<typename MatrixType >

template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>

void Eigen::MatrixPower< MatrixType >::revertSchur ( Matrix< ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols > &  res, const ComplexMatrix &  T, const ComplexMatrix &  U  )

inlinestaticprivate

                                                                                                 {
  res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint());
}

References res, and RachelsAdvectionDiffusion::U.

revertSchur() [2/2]

template<typename MatrixType >

template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>

void Eigen::MatrixPower< MatrixType >::revertSchur ( Matrix< RealScalar, Rows, Cols, Options, MaxRows, MaxCols > &  res, const ComplexMatrix &  T, const ComplexMatrix &  U  )

inlinestaticprivate

                                                                                                 {
  res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real();
}

References res, and RachelsAdvectionDiffusion::U.

rows()

template<typename MatrixType >

Index Eigen::MatrixPower< MatrixType >::rows ( ) const

inline

{ return m_A.rows(); }

References Eigen::MatrixPower< MatrixType >::m_A.

Referenced by gdb.printers._MatrixEntryIterator::__next__(), gdb.printers.EigenMatrixPrinter::children(), gdb.printers.EigenSparseMatrixPrinter::children(), gdb.printers.EigenMatrixPrinter::to_string(), and gdb.printers.EigenSparseMatrixPrinter::to_string().

split()

template<typename MatrixType >

void Eigen::MatrixPower< MatrixType >::split ( RealScalar &  p, RealScalar &  intpart  )

private

Split p into integral part and fractional part.

Parameters

[in]	p	The exponent.
[out]	p	The fractional part ranging in \( (-1, 1) \).
[out]	intpart	The integral part.

Only if the fractional part is nonzero, it calls initialize().

                                                                      {
  using std::floor;
  using std::pow;
 
  intpart = floor(p);
  p -= intpart;
 
  // Perform Schur decomposition if it is not yet performed and the power is
  // not an integer.
  if (!m_conditionNumber && p) initialize();
 
  // Choose the more stable of intpart = floor(p) and intpart = ceil(p).
  if (p > RealScalar(0.5) && p > (1 - p) * pow(m_conditionNumber, p)) {
    --p;
    ++intpart;
  }
}

References Eigen::bfloat16_impl::floor(), p, and Eigen::bfloat16_impl::pow().

Member Data Documentation

m_A

template<typename MatrixType >

MatrixType::Nested Eigen::MatrixPower< MatrixType >::m_A

private

Reference to the base of matrix power.

Referenced by Eigen::MatrixPower< MatrixType >::cols(), and Eigen::MatrixPower< MatrixType >::rows().

m_conditionNumber

template<typename MatrixType >

RealScalar Eigen::MatrixPower< MatrixType >::m_conditionNumber

private

Condition number of m_A.

It is initialized as 0 to avoid performing unnecessary Schur decomposition, which is the bottleneck.

m_fT

template<typename MatrixType >

ComplexMatrix Eigen::MatrixPower< MatrixType >::m_fT

private

Store fractional power of m_T.

m_nulls

template<typename MatrixType >

Index Eigen::MatrixPower< MatrixType >::m_nulls

private

Rank deficiency of m_A.

m_rank

template<typename MatrixType >

Index Eigen::MatrixPower< MatrixType >::m_rank

private

Rank of m_A.

m_T

template<typename MatrixType >

ComplexMatrix Eigen::MatrixPower< MatrixType >::m_T

private

Store the result of Schur decomposition.

m_tmp

template<typename MatrixType >

MatrixType Eigen::MatrixPower< MatrixType >::m_tmp

private

Temporary storage.

m_U

template<typename MatrixType >

ComplexMatrix Eigen::MatrixPower< MatrixType >::m_U

private

---

The documentation for this class was generated from the following file:

• MatrixPower.h