Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
Eigen::MatrixPowerAtomic< MatrixType > Class Template Reference
Matrix functions module

Class for computing matrix powers. More...

Inherits internal::noncopyable.

Public Member Functions

void compute (ResultType &res) const
 Compute the matrix power. More...
 
 MatrixPowerAtomic (const MatrixType &T, RealScalar p)
 Constructor. More...
 

Private Types

enum  {
  RowsAtCompileTime ,
  MaxRowsAtCompileTime
}
 
typedef std::complex< RealScalar > ComplexScalar
 
typedef MatrixType::RealScalar RealScalar
 
typedef Block< MatrixType, Dynamic, Dynamic > ResultType
 
typedef MatrixType::Scalar Scalar
 

Private Member Functions

void compute2x2 (ResultType &res, RealScalar p) const
 
void computeBig (ResultType &res) const
 
void computePade (int degree, const MatrixType &IminusT, ResultType &res) const
 

Static Private Member Functions

static ComplexScalar computeSuperDiag (const ComplexScalar &, const ComplexScalar &, RealScalar p)
 
static RealScalar computeSuperDiag (RealScalar, RealScalar, RealScalar p)
 
static int getPadeDegree (double normIminusT)
 
static int getPadeDegree (float normIminusT)
 
static int getPadeDegree (long double normIminusT)
 

Private Attributes

const MatrixType & m_A
 
RealScalar m_p
 

Detailed Description

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

Class for computing matrix powers.

Template Parameters
MatrixTypetype of the base, expected to be an instantiation of the Matrix class template.

This class is capable of computing triangular real/complex matrices raised to a power in the interval \( (-1, 1) \).

Note
Currently this class is only used by MatrixPower. One may insist that this be nested into MatrixPower. This class is here to facilitate future development of triangular matrix functions.

Definition at line 88 of file MatrixPower.h.

Member Typedef Documentation

◆ ComplexScalar

template<typename MatrixType >
typedef std::complex<RealScalar> Eigen::MatrixPowerAtomic< MatrixType >::ComplexScalar
private

Definition at line 97 of file MatrixPower.h.

◆ RealScalar

template<typename MatrixType >
typedef MatrixType::RealScalar Eigen::MatrixPowerAtomic< MatrixType >::RealScalar
private

Definition at line 96 of file MatrixPower.h.

◆ ResultType

template<typename MatrixType >
typedef Block<MatrixType,Dynamic,Dynamic> Eigen::MatrixPowerAtomic< MatrixType >::ResultType
private

Definition at line 98 of file MatrixPower.h.

◆ Scalar

template<typename MatrixType >
typedef MatrixType::Scalar Eigen::MatrixPowerAtomic< MatrixType >::Scalar
private

Definition at line 95 of file MatrixPower.h.

Member Enumeration Documentation

◆ anonymous enum

template<typename MatrixType >
anonymous enum
private
Enumerator
RowsAtCompileTime 
MaxRowsAtCompileTime 

Definition at line 91 of file MatrixPower.h.

91  {
92  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
93  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
94  };

Constructor & Destructor Documentation

◆ MatrixPowerAtomic()

template<typename MatrixType >
Eigen::MatrixPowerAtomic< MatrixType >::MatrixPowerAtomic ( const MatrixType &  T,
RealScalar  p 
)

Constructor.

Parameters
[in]Tthe base of the matrix power.
[in]pthe exponent of the matrix power, should be in \( (-1, 1) \).

The class stores a reference to T, so it should not be changed (or destroyed) before evaluation. Only the upper triangular part of T is read.

Definition at line 136 of file MatrixPower.h.

136  :
137  m_A(T), m_p(p)
138 {
139  eigen_assert(T.rows() == T.cols());
140  eigen_assert(p > -1 && p < 1);
141 }
eigen_assert
#define eigen_assert(x)
Eigen::MatrixPowerAtomic::m_A
const MatrixType & m_A
Definition: MatrixPower.h:100

Member Function Documentation

◆ compute()

template<typename MatrixType >
void Eigen::MatrixPowerAtomic< MatrixType >::compute ( ResultType &  res) const

Compute the matrix power.

Parameters
[out]res\( A^p \) where A and p are specified in the constructor.

Definition at line 144 of file MatrixPower.h.

145 {
146  using std::pow;
147  switch (m_A.rows()) {
148  case 0:
149  break;
150  case 1:
151  res(0,0) = pow(m_A(0,0), m_p);
152  break;
153  case 2:
154  compute2x2(res, m_p);
155  break;
156  default:
157  computeBig(res);
158  }
159 }
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Eigen::MatrixPowerAtomic::compute2x2
void compute2x2(ResultType &res, RealScalar p) const
Definition: MatrixPower.h:176
Eigen::MatrixPowerAtomic::computeBig
void computeBig(ResultType &res) const
Definition: MatrixPower.h:195
Eigen::pow
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(internal::remove_all_t< DerType >, typename internal::traits< internal::remove_all_t< DerType >>::Scalar, product) > pow(const Eigen::AutoDiffScalar< DerType > &x, const typename internal::traits< internal::remove_all_t< DerType >>::Scalar &y)
Definition: AutoDiffScalar.h:646

◆ compute2x2()

template<typename MatrixType >
void Eigen::MatrixPowerAtomic< MatrixType >::compute2x2 ( ResultType &  res,
RealScalar  p 
) const
private

Definition at line 176 of file MatrixPower.h.

177 {
178  using std::abs;
179  using std::pow;
180  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
181 
182  for (Index i=1; i < m_A.cols(); ++i) {
183  res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
184  if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
185  res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
186  else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1)))
187  res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
188  else
189  res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
190  res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
191  }
192 }
i
int i
p
float * p
Eigen::MatrixPowerAtomic::computeSuperDiag
static ComplexScalar computeSuperDiag(const ComplexScalar &, const ComplexScalar &, RealScalar p)
Definition: MatrixPower.h:293
Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
adtl::abs
adouble abs(const adouble &x)
Definition: AdolcForward:74

◆ computeBig()

template<typename MatrixType >
void Eigen::MatrixPowerAtomic< MatrixType >::computeBig ( ResultType &  res) const
private

Definition at line 195 of file MatrixPower.h.

196 {
197  using std::ldexp;
198  const int digits = std::numeric_limits<RealScalar>::digits;
199  const RealScalar maxNormForPade = RealScalar(
200  digits <= 24? 4.3386528e-1L // single precision
201  : digits <= 53? 2.789358995219730e-1L // double precision
202  : digits <= 64? 2.4471944416607995472e-1L // extended precision
203  : digits <= 106? 1.1016843812851143391275867258512e-1L // double-double
204  : 9.134603732914548552537150753385375e-2L); // quadruple precision
205  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
206  RealScalar normIminusT;
207  int degree, degree2, numberOfSquareRoots = 0;
208  bool hasExtraSquareRoot = false;
209 
210  for (Index i=0; i < m_A.cols(); ++i)
211  eigen_assert(m_A(i,i) != RealScalar(0));
212 
213  while (true) {
214  IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
215  normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
216  if (normIminusT < maxNormForPade) {
217  degree = getPadeDegree(normIminusT);
218  degree2 = getPadeDegree(normIminusT/2);
219  if (degree - degree2 <= 1 || hasExtraSquareRoot)
220  break;
221  hasExtraSquareRoot = true;
222  }
223  matrix_sqrt_triangular(T, sqrtT);
224  T = sqrtT.template triangularView<Upper>();
225  ++numberOfSquareRoots;
226  }
227  computePade(degree, IminusT, res);
228 
229  for (; numberOfSquareRoots; --numberOfSquareRoots) {
230  compute2x2(res, ldexp(m_p, -numberOfSquareRoots));
231  res = res.template triangularView<Upper>() * res;
232  }
233  compute2x2(res, m_p);
234 }
MatrixType
Matrix< float, 1, Dynamic > MatrixType
Eigen::MatrixPowerAtomic::computePade
void computePade(int degree, const MatrixType &IminusT, ResultType &res) const
Definition: MatrixPower.h:162
Eigen::MatrixPowerAtomic::getPadeDegree
static int getPadeDegree(float normIminusT)
Definition: MatrixPower.h:237
Eigen::MatrixPowerAtomic::RealScalar
MatrixType::RealScalar RealScalar
Definition: MatrixPower.h:96
Eigen::matrix_sqrt_triangular
void matrix_sqrt_triangular(const MatrixType &arg, ResultType &result)
Compute matrix square root of triangular matrix.
Definition: MatrixSquareRoot.h:206

◆ computePade()

template<typename MatrixType >
void Eigen::MatrixPowerAtomic< MatrixType >::computePade ( int  degree,
const MatrixType &  IminusT,
ResultType &  res 
) const
private

Definition at line 162 of file MatrixPower.h.

163 {
164  int i = 2*degree;
165  res = (m_p-RealScalar(degree)) / RealScalar(2*i-2) * IminusT;
166 
167  for (--i; i; --i) {
168  res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
169  .solve((i==1 ? -m_p : i&1 ? (-m_p-RealScalar(i/2))/RealScalar(2*i) : (m_p-RealScalar(i/2))/RealScalar(2*i-2)) * IminusT).eval();
170  }
171  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
172 }

◆ computeSuperDiag() [1/2]

template<typename MatrixType >
MatrixPowerAtomic< MatrixType >::ComplexScalar Eigen::MatrixPowerAtomic< MatrixType >::computeSuperDiag ( const ComplexScalar &  curr,
const ComplexScalar &  prev,
RealScalar  p 
)
inlinestaticprivate

Definition at line 293 of file MatrixPower.h.

294 {
295  using std::ceil;
296  using std::exp;
297  using std::log;
298  using std::sinh;
299 
300  ComplexScalar logCurr = log(curr);
301  ComplexScalar logPrev = log(prev);
302  RealScalar unwindingNumber = ceil((numext::imag(logCurr - logPrev) - RealScalar(EIGEN_PI)) / RealScalar(2*EIGEN_PI));
303  ComplexScalar w = numext::log1p((curr-prev)/prev)/RealScalar(2) + ComplexScalar(0, RealScalar(EIGEN_PI)*unwindingNumber);
304  return RealScalar(2) * exp(RealScalar(0.5) * p * (logCurr + logPrev)) * sinh(p * w) / (curr - prev);
305 }
w
RowVector3d w
Eigen::MatrixPowerAtomic::ComplexScalar
std::complex< RealScalar > ComplexScalar
Definition: MatrixPower.h:97
Eigen::sinh
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Eigen::internal::remove_all_t< DerType >, typename Eigen::internal::traits< Eigen::internal::remove_all_t< DerType >>::Scalar, product) > sinh(const Eigen::AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:695
Eigen::log
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Eigen::internal::remove_all_t< DerType >, typename Eigen::internal::traits< Eigen::internal::remove_all_t< DerType >>::Scalar, product) > log(const Eigen::AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:641
Eigen::exp
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Eigen::internal::remove_all_t< DerType >, typename Eigen::internal::traits< Eigen::internal::remove_all_t< DerType >>::Scalar, product) > exp(const Eigen::AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:637
Eigen::ceil
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_ceil_op< typename Derived::Scalar >, const Derived > ceil(const Eigen::ArrayBase< Derived > &x)
Eigen::exp
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_exp_op< typename Derived::Scalar >, const Derived > exp(const Eigen::ArrayBase< Derived > &x)
Eigen::log
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_log_op< typename Derived::Scalar >, const Derived > log(const Eigen::ArrayBase< Derived > &x)
Eigen::sinh
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sinh_op< typename Derived::Scalar >, const Derived > sinh(const Eigen::ArrayBase< Derived > &x)
adtl::imag
adouble imag(const adouble &)
Definition: AdolcForward:73

◆ computeSuperDiag() [2/2]

template<typename MatrixType >
MatrixPowerAtomic< MatrixType >::RealScalar Eigen::MatrixPowerAtomic< MatrixType >::computeSuperDiag ( RealScalar  curr,
RealScalar  prev,
RealScalar  p 
)
inlinestaticprivate

Definition at line 309 of file MatrixPower.h.

310 {
311  using std::exp;
312  using std::log;
313  using std::sinh;
314 
315  RealScalar w = numext::log1p((curr-prev)/prev)/RealScalar(2);
316  return 2 * exp(p * (log(curr) + log(prev)) / 2) * sinh(p * w) / (curr - prev);
317 }

◆ getPadeDegree() [1/3]

template<typename MatrixType >
int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( double  normIminusT)
inlinestaticprivate

Definition at line 248 of file MatrixPower.h.

249 {
250  const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
251  1.999045567181744e-1, 2.789358995219730e-1 };
252  int degree = 3;
253  for (; degree <= 7; ++degree)
254  if (normIminusT <= maxNormForPade[degree - 3])
255  break;
256  return degree;
257 }

◆ getPadeDegree() [2/3]

template<typename MatrixType >
int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( float  normIminusT)
inlinestaticprivate

Definition at line 237 of file MatrixPower.h.

238 {
239  const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
240  int degree = 3;
241  for (; degree <= 4; ++degree)
242  if (normIminusT <= maxNormForPade[degree - 3])
243  break;
244  return degree;
245 }

◆ getPadeDegree() [3/3]

template<typename MatrixType >
int Eigen::MatrixPowerAtomic< MatrixType >::getPadeDegree ( long double  normIminusT)
inlinestaticprivate

Definition at line 260 of file MatrixPower.h.

261 {
262 #if LDBL_MANT_DIG == 53
263  const int maxPadeDegree = 7;
264  const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
265  1.999045567181744e-1L, 2.789358995219730e-1L };
266 #elif LDBL_MANT_DIG <= 64
267  const int maxPadeDegree = 8;
268  const long double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
269  6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
270 #elif LDBL_MANT_DIG <= 106
271  const int maxPadeDegree = 10;
272  const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
273  1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
274  2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
275  1.1016843812851143391275867258512e-1L };
276 #else
277  const int maxPadeDegree = 10;
278  const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
279  6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
280  9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
281  3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
282  9.134603732914548552537150753385375e-2L };
283 #endif
284  int degree = 3;
285  for (; degree <= maxPadeDegree; ++degree)
286  if (normIminusT <= static_cast<long double>(maxNormForPade[degree - 3]))
287  break;
288  return degree;
289 }

Member Data Documentation

◆ m_A

template<typename MatrixType >
const MatrixType& Eigen::MatrixPowerAtomic< MatrixType >::m_A
private

Definition at line 100 of file MatrixPower.h.

◆ m_p

template<typename MatrixType >
RealScalar Eigen::MatrixPowerAtomic< MatrixType >::m_p
private

Definition at line 101 of file MatrixPower.h.

The documentation for this class was generated from the following file:
Generated on Tue Jul 18 2023 15:28:51 for Eigen-unsupported by Doxygen 1.9.1 and the doxyYoda theme. This is an unofficial themed variant maintained by Rohit Goswami, based on the official documentation.