11 #ifndef EIGEN_MATRIX_EXPONENTIAL
12 #define EIGEN_MATRIX_EXPONENTIAL
14 #include "StemFunction.h"
23 template <
typename RealScalar>
24 struct MatrixExponentialScalingOp
30 MatrixExponentialScalingOp(
int squarings) : m_squarings(squarings) { }
37 inline const RealScalar operator() (
const RealScalar& x)
const
40 return ldexp(x, -m_squarings);
43 typedef std::complex<RealScalar> ComplexScalar;
49 inline const ComplexScalar operator() (
const ComplexScalar& x)
const
52 return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
64 template <
typename MatrixType>
65 void matrix_exp_pade3(
const MatrixType &A, MatrixType &U, MatrixType &V)
67 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
68 const RealScalar b[] = {120., 60., 12., 1.};
69 const MatrixType A2 = A * A;
70 const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
71 U.noalias() = A * tmp;
72 V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
80 template <
typename MatrixType>
81 void matrix_exp_pade5(
const MatrixType &A, MatrixType &U, MatrixType &V)
83 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
84 const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.};
85 const MatrixType A2 = A * A;
86 const MatrixType A4 = A2 * A2;
87 const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
88 U.noalias() = A * tmp;
89 V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
97 template <
typename MatrixType>
98 void matrix_exp_pade7(
const MatrixType &A, MatrixType &U, MatrixType &V)
100 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
101 const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
102 const MatrixType A2 = A * A;
103 const MatrixType A4 = A2 * A2;
104 const MatrixType A6 = A4 * A2;
105 const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2
106 + b[1] * MatrixType::Identity(A.rows(), A.cols());
107 U.noalias() = A * tmp;
108 V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
117 template <
typename MatrixType>
118 void matrix_exp_pade9(
const MatrixType &A, MatrixType &U, MatrixType &V)
120 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
121 const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
122 2162160., 110880., 3960., 90., 1.};
123 const MatrixType A2 = A * A;
124 const MatrixType A4 = A2 * A2;
125 const MatrixType A6 = A4 * A2;
126 const MatrixType A8 = A6 * A2;
127 const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
128 + b[1] * MatrixType::Identity(A.rows(), A.cols());
129 U.noalias() = A * tmp;
130 V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
138 template <
typename MatrixType>
139 void matrix_exp_pade13(
const MatrixType &A, MatrixType &U, MatrixType &V)
141 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
142 const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
143 1187353796428800., 129060195264000., 10559470521600., 670442572800.,
144 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
145 const MatrixType A2 = A * A;
146 const MatrixType A4 = A2 * A2;
147 const MatrixType A6 = A4 * A2;
148 V = b[13] * A6 + b[11] * A4 + b[9] * A2;
149 MatrixType tmp = A6 * V;
150 tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
151 U.noalias() = A * tmp;
152 tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
153 V.noalias() = A6 * tmp;
154 V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
164 #if LDBL_MANT_DIG > 64
165 template <
typename MatrixType>
166 void matrix_exp_pade17(
const MatrixType &A, MatrixType &U, MatrixType &V)
168 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
169 const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
170 100610229646136770560000.L, 15720348382208870400000.L,
171 1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
172 595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
173 33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
174 46512.L, 306.L, 1.L};
175 const MatrixType A2 = A * A;
176 const MatrixType A4 = A2 * A2;
177 const MatrixType A6 = A4 * A2;
178 const MatrixType A8 = A4 * A4;
179 V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2;
180 MatrixType tmp = A8 * V;
181 tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
182 + b[1] * MatrixType::Identity(A.rows(), A.cols());
183 U.noalias() = A * tmp;
184 tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
185 V.noalias() = tmp * A8;
186 V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2
187 + b[0] * MatrixType::Identity(A.rows(), A.cols());
191 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
192 struct matrix_exp_computeUV
201 static void run(
const MatrixType& arg, MatrixType& U, MatrixType& V,
int& squarings);
204 template <
typename MatrixType>
205 struct matrix_exp_computeUV<MatrixType, float>
207 static void run(
const MatrixType& arg, MatrixType& U, MatrixType& V,
int& squarings)
211 const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
213 if (l1norm < 4.258730016922831e-001) {
214 matrix_exp_pade3(arg, U, V);
215 }
else if (l1norm < 1.880152677804762e+000) {
216 matrix_exp_pade5(arg, U, V);
218 const float maxnorm = 3.925724783138660f;
219 frexp(l1norm / maxnorm, &squarings);
220 if (squarings < 0) squarings = 0;
221 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
222 matrix_exp_pade7(A, U, V);
227 template <
typename MatrixType>
228 struct matrix_exp_computeUV<MatrixType, double>
230 static void run(
const MatrixType& arg, MatrixType& U, MatrixType& V,
int& squarings)
234 const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
236 if (l1norm < 1.495585217958292e-002) {
237 matrix_exp_pade3(arg, U, V);
238 }
else if (l1norm < 2.539398330063230e-001) {
239 matrix_exp_pade5(arg, U, V);
240 }
else if (l1norm < 9.504178996162932e-001) {
241 matrix_exp_pade7(arg, U, V);
242 }
else if (l1norm < 2.097847961257068e+000) {
243 matrix_exp_pade9(arg, U, V);
245 const double maxnorm = 5.371920351148152;
246 frexp(l1norm / maxnorm, &squarings);
247 if (squarings < 0) squarings = 0;
248 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<double>(squarings));
249 matrix_exp_pade13(A, U, V);
254 template <
typename MatrixType>
255 struct matrix_exp_computeUV<MatrixType, long double>
257 static void run(
const MatrixType& arg, MatrixType& U, MatrixType& V,
int& squarings)
259 #if LDBL_MANT_DIG == 53 // double precision
261 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
267 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
270 #if LDBL_MANT_DIG <= 64 // extended precision
272 if (l1norm < 4.1968497232266989671e-003L) {
273 matrix_exp_pade3(arg, U, V);
274 }
else if (l1norm < 1.1848116734693823091e-001L) {
275 matrix_exp_pade5(arg, U, V);
276 }
else if (l1norm < 5.5170388480686700274e-001L) {
277 matrix_exp_pade7(arg, U, V);
278 }
else if (l1norm < 1.3759868875587845383e+000L) {
279 matrix_exp_pade9(arg, U, V);
281 const long double maxnorm = 4.0246098906697353063L;
282 frexp(l1norm / maxnorm, &squarings);
283 if (squarings < 0) squarings = 0;
284 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
285 matrix_exp_pade13(A, U, V);
288 #elif LDBL_MANT_DIG <= 106 // double-double
290 if (l1norm < 3.2787892205607026992947488108213e-005L) {
291 matrix_exp_pade3(arg, U, V);
292 }
else if (l1norm < 6.4467025060072760084130906076332e-003L) {
293 matrix_exp_pade5(arg, U, V);
294 }
else if (l1norm < 6.8988028496595374751374122881143e-002L) {
295 matrix_exp_pade7(arg, U, V);
296 }
else if (l1norm < 2.7339737518502231741495857201670e-001L) {
297 matrix_exp_pade9(arg, U, V);
298 }
else if (l1norm < 1.3203382096514474905666448850278e+000L) {
299 matrix_exp_pade13(arg, U, V);
301 const long double maxnorm = 3.2579440895405400856599663723517L;
302 frexp(l1norm / maxnorm, &squarings);
303 if (squarings < 0) squarings = 0;
304 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
305 matrix_exp_pade17(A, U, V);
308 #elif LDBL_MANT_DIG <= 112 // quadruple precison
310 if (l1norm < 1.639394610288918690547467954466970e-005L) {
311 matrix_exp_pade3(arg, U, V);
312 }
else if (l1norm < 4.253237712165275566025884344433009e-003L) {
313 matrix_exp_pade5(arg, U, V);
314 }
else if (l1norm < 5.125804063165764409885122032933142e-002L) {
315 matrix_exp_pade7(arg, U, V);
316 }
else if (l1norm < 2.170000765161155195453205651889853e-001L) {
317 matrix_exp_pade9(arg, U, V);
318 }
else if (l1norm < 1.125358383453143065081397882891878e+000L) {
319 matrix_exp_pade13(arg, U, V);
321 frexp(l1norm / maxnorm, &squarings);
322 if (squarings < 0) squarings = 0;
323 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
324 matrix_exp_pade17(A, U, V);
330 eigen_assert(
false &&
"Bug in MatrixExponential");
333 #endif // LDBL_MANT_DIG
343 template <
typename MatrixType,
typename ResultType>
344 void matrix_exp_compute(
const MatrixType& arg, ResultType &result)
346 #if LDBL_MANT_DIG > 112 // rarely happens
347 typedef typename traits<MatrixType>::Scalar Scalar;
348 typedef typename NumTraits<Scalar>::Real RealScalar;
349 typedef typename std::complex<RealScalar> ComplexScalar;
350 if (
sizeof(RealScalar) > 14) {
351 result = arg.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
357 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings);
358 MatrixType numer = U + V;
359 MatrixType denom = -U + V;
360 result = denom.partialPivLu().solve(numer);
361 for (
int i=0; i<squarings; i++)
378 :
public ReturnByValue<MatrixExponentialReturnValue<Derived> >
380 typedef typename Derived::Index Index;
392 template <
typename ResultType>
393 inline void evalTo(ResultType& result)
const
395 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
396 internal::matrix_exp_compute(tmp, result);
399 Index rows()
const {
return m_src.rows(); }
400 Index cols()
const {
return m_src.cols(); }
403 const typename internal::ref_selector<Derived>::type m_src;
407 template<
typename Derived>
408 struct traits<MatrixExponentialReturnValue<Derived> >
410 typedef typename Derived::PlainObject ReturnType;
414 template <
typename Derived>
415 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp()
const
417 eigen_assert(rows() == cols());
418 return MatrixExponentialReturnValue<Derived>(derived());
423 #endif // EIGEN_MATRIX_EXPONENTIAL
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:393
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:377
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:386