MatrixExponential.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
5 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_MATRIX_EXPONENTIAL
12 #define EIGEN_MATRIX_EXPONENTIAL
13 
14 #include "StemFunction.h"
15 
16 namespace Eigen {
17 namespace internal {
18 
23 template <typename RealScalar>
24 struct MatrixExponentialScalingOp
25 {
30  MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
31 
32 
37  inline const RealScalar operator() (const RealScalar& x) const
38  {
39  using std::ldexp;
40  return ldexp(x, -m_squarings);
41  }
42 
43  typedef std::complex<RealScalar> ComplexScalar;
44 
49  inline const ComplexScalar operator() (const ComplexScalar& x) const
50  {
51  using std::ldexp;
52  return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
53  }
54 
55  private:
56  int m_squarings;
57 };
58 
64 template <typename MatA, typename MatU, typename MatV>
65 void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V)
66 {
67  typedef typename MatA::PlainObject MatrixType;
68  typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
69  const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
70  const MatrixType A2 = A * A;
71  const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
72  U.noalias() = A * tmp;
73  V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
74 }
75 
81 template <typename MatA, typename MatU, typename MatV>
82 void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V)
83 {
84  typedef typename MatA::PlainObject MatrixType;
85  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
86  const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
87  const MatrixType A2 = A * A;
88  const MatrixType A4 = A2 * A2;
89  const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
90  U.noalias() = A * tmp;
91  V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
92 }
93 
99 template <typename MatA, typename MatU, typename MatV>
100 void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V)
101 {
102  typedef typename MatA::PlainObject MatrixType;
103  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
104  const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
105  const MatrixType A2 = A * A;
106  const MatrixType A4 = A2 * A2;
107  const MatrixType A6 = A4 * A2;
108  const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2
109  + b[1] * MatrixType::Identity(A.rows(), A.cols());
110  U.noalias() = A * tmp;
111  V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
112 
113 }
114 
120 template <typename MatA, typename MatU, typename MatV>
121 void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V)
122 {
123  typedef typename MatA::PlainObject MatrixType;
124  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
125  const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
126  2162160.L, 110880.L, 3960.L, 90.L, 1.L};
127  const MatrixType A2 = A * A;
128  const MatrixType A4 = A2 * A2;
129  const MatrixType A6 = A4 * A2;
130  const MatrixType A8 = A6 * A2;
131  const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
132  + b[1] * MatrixType::Identity(A.rows(), A.cols());
133  U.noalias() = A * tmp;
134  V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
135 }
136 
142 template <typename MatA, typename MatU, typename MatV>
143 void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V)
144 {
145  typedef typename MatA::PlainObject MatrixType;
146  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
147  const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
148  1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
149  33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
150  const MatrixType A2 = A * A;
151  const MatrixType A4 = A2 * A2;
152  const MatrixType A6 = A4 * A2;
153  V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
154  MatrixType tmp = A6 * V;
155  tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
156  U.noalias() = A * tmp;
157  tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
158  V.noalias() = A6 * tmp;
159  V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
160 }
161 
169 #if LDBL_MANT_DIG > 64
170 template <typename MatA, typename MatU, typename MatV>
171 void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V)
172 {
173  typedef typename MatA::PlainObject MatrixType;
174  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
175  const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
176  100610229646136770560000.L, 15720348382208870400000.L,
177  1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
178  595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
179  33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
180  46512.L, 306.L, 1.L};
181  const MatrixType A2 = A * A;
182  const MatrixType A4 = A2 * A2;
183  const MatrixType A6 = A4 * A2;
184  const MatrixType A8 = A4 * A4;
185  V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
186  MatrixType tmp = A8 * V;
187  tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
188  + b[1] * MatrixType::Identity(A.rows(), A.cols());
189  U.noalias() = A * tmp;
190  tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
191  V.noalias() = tmp * A8;
192  V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2
193  + b[0] * MatrixType::Identity(A.rows(), A.cols());
194 }
195 #endif
196 
197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
198 struct matrix_exp_computeUV
199 {
207  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
208 };
209 
210 template <typename MatrixType>
211 struct matrix_exp_computeUV<MatrixType, float>
212 {
213  template <typename ArgType>
214  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
215  {
216  using std::frexp;
217  using std::pow;
218  const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
219  squarings = 0;
220  if (l1norm < 4.258730016922831e-001f) {
221  matrix_exp_pade3(arg, U, V);
222  } else if (l1norm < 1.880152677804762e+000f) {
223  matrix_exp_pade5(arg, U, V);
224  } else {
225  const float maxnorm = 3.925724783138660f;
226  frexp(l1norm / maxnorm, &squarings);
227  if (squarings < 0) squarings = 0;
228  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
229  matrix_exp_pade7(A, U, V);
230  }
231  }
232 };
233 
234 template <typename MatrixType>
235 struct matrix_exp_computeUV<MatrixType, double>
236 {
237  template <typename ArgType>
238  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
239  {
240  using std::frexp;
241  using std::pow;
242  const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
243  squarings = 0;
244  if (l1norm < 1.495585217958292e-002) {
245  matrix_exp_pade3(arg, U, V);
246  } else if (l1norm < 2.539398330063230e-001) {
247  matrix_exp_pade5(arg, U, V);
248  } else if (l1norm < 9.504178996162932e-001) {
249  matrix_exp_pade7(arg, U, V);
250  } else if (l1norm < 2.097847961257068e+000) {
251  matrix_exp_pade9(arg, U, V);
252  } else {
253  const double maxnorm = 5.371920351148152;
254  frexp(l1norm / maxnorm, &squarings);
255  if (squarings < 0) squarings = 0;
256  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<double>(squarings));
257  matrix_exp_pade13(A, U, V);
258  }
259  }
260 };
261 
262 template <typename MatrixType>
263 struct matrix_exp_computeUV<MatrixType, long double>
264 {
265  template <typename ArgType>
266  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
267  {
268 #if LDBL_MANT_DIG == 53 // double precision
269  matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
270 
271 #else
272 
273  using std::frexp;
274  using std::pow;
275  const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
276  squarings = 0;
277 
278 #if LDBL_MANT_DIG <= 64 // extended precision
279 
280  if (l1norm < 4.1968497232266989671e-003L) {
281  matrix_exp_pade3(arg, U, V);
282  } else if (l1norm < 1.1848116734693823091e-001L) {
283  matrix_exp_pade5(arg, U, V);
284  } else if (l1norm < 5.5170388480686700274e-001L) {
285  matrix_exp_pade7(arg, U, V);
286  } else if (l1norm < 1.3759868875587845383e+000L) {
287  matrix_exp_pade9(arg, U, V);
288  } else {
289  const long double maxnorm = 4.0246098906697353063L;
290  frexp(l1norm / maxnorm, &squarings);
291  if (squarings < 0) squarings = 0;
292  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
293  matrix_exp_pade13(A, U, V);
294  }
295 
296 #elif LDBL_MANT_DIG <= 106 // double-double
297 
298  if (l1norm < 3.2787892205607026992947488108213e-005L) {
299  matrix_exp_pade3(arg, U, V);
300  } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
301  matrix_exp_pade5(arg, U, V);
302  } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
303  matrix_exp_pade7(arg, U, V);
304  } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
305  matrix_exp_pade9(arg, U, V);
306  } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
307  matrix_exp_pade13(arg, U, V);
308  } else {
309  const long double maxnorm = 3.2579440895405400856599663723517L;
310  frexp(l1norm / maxnorm, &squarings);
311  if (squarings < 0) squarings = 0;
312  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
313  matrix_exp_pade17(A, U, V);
314  }
315 
316 #elif LDBL_MANT_DIG <= 112 // quadruple precison
317 
318  if (l1norm < 1.639394610288918690547467954466970e-005L) {
319  matrix_exp_pade3(arg, U, V);
320  } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
321  matrix_exp_pade5(arg, U, V);
322  } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
323  matrix_exp_pade7(arg, U, V);
324  } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
325  matrix_exp_pade9(arg, U, V);
326  } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
327  matrix_exp_pade13(arg, U, V);
328  } else {
329  frexp(l1norm / maxnorm, &squarings);
330  if (squarings < 0) squarings = 0;
331  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
332  matrix_exp_pade17(A, U, V);
333  }
334 
335 #else
336 
337  // this case should be handled in compute()
338  eigen_assert(false && "Bug in MatrixExponential");
339 
340 #endif
341 #endif // LDBL_MANT_DIG
342  }
343 };
344 
345 
346 /* Computes the matrix exponential
347  *
348  * \param arg argument of matrix exponential (should be plain object)
349  * \param result variable in which result will be stored
350  */
351 template <typename ArgType, typename ResultType>
352 void matrix_exp_compute(const ArgType& arg, ResultType &result)
353 {
354  typedef typename ArgType::PlainObject MatrixType;
355 #if LDBL_MANT_DIG > 112 // rarely happens
356  typedef typename traits<MatrixType>::Scalar Scalar;
357  typedef typename NumTraits<Scalar>::Real RealScalar;
358  typedef typename std::complex<RealScalar> ComplexScalar;
359  if (sizeof(RealScalar) > 14) {
360  result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
361  return;
362  }
363 #endif
364  MatrixType U, V;
365  int squarings;
366  matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
367  MatrixType numer = U + V;
368  MatrixType denom = -U + V;
369  result = denom.partialPivLu().solve(numer);
370  for (int i=0; i<squarings; i++)
371  result *= result; // undo scaling by repeated squaring
372 }
373 
374 } // end namespace Eigen::internal
375 
386 template<typename Derived> struct MatrixExponentialReturnValue
387 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
388 {
389  typedef typename Derived::Index Index;
390  public:
395  MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
396 
401  template <typename ResultType>
402  inline void evalTo(ResultType& result) const
403  {
404  const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
405  internal::matrix_exp_compute(tmp, result);
406  }
407 
408  Index rows() const { return m_src.rows(); }
409  Index cols() const { return m_src.cols(); }
410 
411  protected:
412  const typename internal::ref_selector<Derived>::type m_src;
413 };
414 
415 namespace internal {
416 template<typename Derived>
417 struct traits<MatrixExponentialReturnValue<Derived> >
418 {
419  typedef typename Derived::PlainObject ReturnType;
420 };
421 }
422 
423 template <typename Derived>
424 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
425 {
426  eigen_assert(rows() == cols());
427  return MatrixExponentialReturnValue<Derived>(derived());
428 }
429 
430 } // end namespace Eigen
431 
432 #endif // EIGEN_MATRIX_EXPONENTIAL
Eigen::half_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half pow(const half &a, const half &b)
Definition: Half.h:407
Eigen::MatrixExponentialReturnValue::evalTo
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:402
Eigen::internal::matrix_exp_compute
void matrix_exp_compute(const ArgType &arg, ResultType &result)
Definition: MatrixExponential.h:352
mpfr::ldexp
const mpreal ldexp(const mpreal &v, mp_exp_t exp)
Definition: mpreal.h:2020
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:238
Eigen
Definition: LDLT.h:16
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Eigen::internal::matrix_exp_computeUV
Compute the (17,17)-Padé approximant to the exponential.
Definition: MatrixExponential.h:198
Eigen::MatrixBase::exp
const MatrixExponentialReturnValue< Derived > exp() const
Definition: MatrixExponential.h:424
Eigen::internal::matrix_exp_pade13
void matrix_exp_pade13(const MatA &A, MatU &U, MatV &V)
Compute the (13,13)-Padé approximant to the exponential.
Definition: MatrixExponential.h:143
Eigen::internal::matrix_exp_computeUV::run
static void run(const MatrixType &arg, MatrixType &U, MatrixType &V, int &squarings)
Compute Padé approximant to the exponential.
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Eigen::internal::matrix_exp_computeUV< MatrixType, double >::run
static void run(const ArgType &arg, MatrixType &U, MatrixType &V, int &squarings)
Definition: MatrixExponential.h:238
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:577
Eigen::MatrixExponentialReturnValue::m_src
const internal::ref_selector< Derived >::type m_src
Definition: MatrixExponential.h:412
Eigen::internal::matrix_exp_computeUV< MatrixType, float >::run
static void run(const ArgType &arg, MatrixType &U, MatrixType &V, int &squarings)
Definition: MatrixExponential.h:214
Eigen::internal::matrix_exp_pade3
void matrix_exp_pade3(const MatA &A, MatU &U, MatV &V)
Compute the (3,3)-Padé approximant to the exponential.
Definition: MatrixExponential.h:65
mpfr::frexp
const mpreal frexp(const mpreal &x, mp_exp_t *exp, mp_rnd_t mode=mpreal::get_default_rnd())
Definition: mpreal.h:2008
Eigen::internal::matrix_exp_pade5
void matrix_exp_pade5(const MatA &A, MatU &U, MatV &V)
Compute the (5,5)-Padé approximant to the exponential.
Definition: MatrixExponential.h:82
Eigen::internal::matrix_exp_pade9
void matrix_exp_pade9(const MatA &A, MatU &U, MatV &V)
Compute the (9,9)-Padé approximant to the exponential.
Definition: MatrixExponential.h:121
arg
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArgReturnType arg() const
Definition: ArrayCwiseUnaryOps.h:57
Eigen::internal::MatrixExponentialScalingOp::MatrixExponentialScalingOp
MatrixExponentialScalingOp(int squarings)
Constructor.
Definition: MatrixExponential.h:30
Eigen::internal::matrix_exp_pade7
void matrix_exp_pade7(const MatA &A, MatU &U, MatV &V)
Compute the (7,7)-Padé approximant to the exponential.
Definition: MatrixExponential.h:100
Eigen::MatrixExponentialReturnValue
Proxy for the matrix exponential of some matrix (expression).
Definition: ForwardDeclarations.h:284
Eigen::internal::MatrixExponentialScalingOp::operator()
const RealScalar operator()(const RealScalar &x) const
Scale a matrix coefficient.
Definition: MatrixExponential.h:37
Eigen::internal::true_type
Definition: Meta.h:54
Eigen::internal::matrix_exp_computeUV< MatrixType, long double >::run
static void run(const ArgType &arg, MatrixType &U, MatrixType &V, int &squarings)
Definition: MatrixExponential.h:266
Eigen::MatrixExponentialReturnValue::MatrixExponentialReturnValue
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:395
Eigen::TensorSycl::run
void run(Expr &expr, Dev &dev)
Definition: TensorSyclRun.h:33
Eigen::numext::b
EIGEN_DEVICE_FUNC const Scalar & b
Definition: SpecialFunctionsImpl.h:1556

hebiros
Author(s): Xavier Artache , Matthew Tesch
autogenerated on Thu Sep 3 2020 04:08:27