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.L, 60.L, 12.L, 1.L};
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.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
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.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
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.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
122 2162160.L, 110880.L, 3960.L, 90.L, 1.L};
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.L, 32382376266240000.L, 7771770303897600.L,
143 1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
144 33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
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-001f) {
214 matrix_exp_pade3(arg, U, V);
215 }
else if (l1norm < 1.880152677804762e+000f) {
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 260 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
266 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
269 #if LDBL_MANT_DIG <= 64 // extended precision 271 if (l1norm < 4.1968497232266989671e-003L) {
272 matrix_exp_pade3(arg, U, V);
273 }
else if (l1norm < 1.1848116734693823091e-001L) {
274 matrix_exp_pade5(arg, U, V);
275 }
else if (l1norm < 5.5170388480686700274e-001L) {
276 matrix_exp_pade7(arg, U, V);
277 }
else if (l1norm < 1.3759868875587845383e+000L) {
278 matrix_exp_pade9(arg, U, V);
280 const long double maxnorm = 4.0246098906697353063L;
281 frexp(l1norm / maxnorm, &squarings);
282 if (squarings < 0) squarings = 0;
283 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
284 matrix_exp_pade13(A, U, V);
287 #elif LDBL_MANT_DIG <= 106 // double-double 289 if (l1norm < 3.2787892205607026992947488108213e-005L) {
290 matrix_exp_pade3(arg, U, V);
291 }
else if (l1norm < 6.4467025060072760084130906076332e-003L) {
292 matrix_exp_pade5(arg, U, V);
293 }
else if (l1norm < 6.8988028496595374751374122881143e-002L) {
294 matrix_exp_pade7(arg, U, V);
295 }
else if (l1norm < 2.7339737518502231741495857201670e-001L) {
296 matrix_exp_pade9(arg, U, V);
297 }
else if (l1norm < 1.3203382096514474905666448850278e+000L) {
298 matrix_exp_pade13(arg, U, V);
300 const long double maxnorm = 3.2579440895405400856599663723517L;
301 frexp(l1norm / maxnorm, &squarings);
302 if (squarings < 0) squarings = 0;
303 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
304 matrix_exp_pade17(A, U, V);
307 #elif LDBL_MANT_DIG <= 112 // quadruple precison 309 if (l1norm < 1.639394610288918690547467954466970e-005L) {
310 matrix_exp_pade3(arg, U, V);
311 }
else if (l1norm < 4.253237712165275566025884344433009e-003L) {
312 matrix_exp_pade5(arg, U, V);
313 }
else if (l1norm < 5.125804063165764409885122032933142e-002L) {
314 matrix_exp_pade7(arg, U, V);
315 }
else if (l1norm < 2.170000765161155195453205651889853e-001L) {
316 matrix_exp_pade9(arg, U, V);
317 }
else if (l1norm < 1.125358383453143065081397882891878e+000L) {
318 matrix_exp_pade13(arg, U, V);
320 frexp(l1norm / maxnorm, &squarings);
321 if (squarings < 0) squarings = 0;
322 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
323 matrix_exp_pade17(A, U, V);
329 eigen_assert(
false &&
"Bug in MatrixExponential");
332 #endif // LDBL_MANT_DIG 342 template <
typename MatrixType,
typename ResultType>
343 void matrix_exp_compute(
const MatrixType& arg, ResultType &result)
345 #if LDBL_MANT_DIG > 112 // rarely happens 346 typedef typename traits<MatrixType>::Scalar Scalar;
347 typedef typename NumTraits<Scalar>::Real RealScalar;
348 typedef typename std::complex<RealScalar> ComplexScalar;
349 if (
sizeof(RealScalar) > 14) {
350 result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
356 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings);
357 MatrixType numer = U + V;
358 MatrixType denom = -U + V;
359 result = denom.partialPivLu().solve(numer);
360 for (
int i=0; i<squarings; i++)
377 :
public ReturnByValue<MatrixExponentialReturnValue<Derived> >
379 typedef typename Derived::Index Index;
391 template <
typename ResultType>
392 inline void evalTo(ResultType& result)
const 394 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
395 internal::matrix_exp_compute(tmp, result);
398 Index rows()
const {
return m_src.rows(); }
399 Index cols()
const {
return m_src.cols(); }
402 const typename internal::ref_selector<Derived>::type m_src;
406 template<
typename Derived>
407 struct traits<MatrixExponentialReturnValue<Derived> >
409 typedef typename Derived::PlainObject ReturnType;
413 template <
typename Derived>
416 eigen_assert(rows() == cols());
422 #endif // EIGEN_MATRIX_EXPONENTIAL Namespace containing all symbols from the Eigen library.
Definition: AdolcForward:45
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:392
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:376
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:385