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 MatA,
typename MatU,
typename MatV>
65 void matrix_exp_pade3(
const MatA& A, MatU& U, MatV& V)
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());
81 template <
typename MatA,
typename MatU,
typename MatV>
82 void matrix_exp_pade5(
const MatA& A, MatU& U, MatV& V)
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());
99 template <
typename MatA,
typename MatU,
typename MatV>
100 void matrix_exp_pade7(
const MatA& A, MatU& U, MatV& V)
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());
120 template <
typename MatA,
typename MatU,
typename MatV>
121 void matrix_exp_pade9(
const MatA& A, MatU& U, MatV& V)
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());
142 template <
typename MatA,
typename MatU,
typename MatV>
143 void matrix_exp_pade13(
const MatA& A, MatU& U, MatV& V)
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;
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());
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)
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;
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());
197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
198 struct matrix_exp_computeUV
207 static void run(
const MatrixType&
arg, MatrixType& U, MatrixType& V,
int& squarings);
210 template <
typename MatrixType>
211 struct matrix_exp_computeUV<MatrixType, float>
213 template <
typename ArgType>
214 static void run(
const ArgType&
arg, MatrixType& U, MatrixType& V,
int& squarings)
218 const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
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);
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);
234 template <
typename MatrixType>
235 struct matrix_exp_computeUV<MatrixType, double>
237 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
238 template <
typename ArgType>
239 static void run(
const ArgType& arg, MatrixType& U, MatrixType& V,
int& squarings)
243 const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
245 if (l1norm < 1.495585217958292e-002) {
246 matrix_exp_pade3(arg, U, V);
247 }
else if (l1norm < 2.539398330063230e-001) {
248 matrix_exp_pade5(arg, U, V);
249 }
else if (l1norm < 9.504178996162932e-001) {
250 matrix_exp_pade7(arg, U, V);
251 }
else if (l1norm < 2.097847961257068e+000) {
252 matrix_exp_pade9(arg, U, V);
254 const RealScalar maxnorm = 5.371920351148152;
255 frexp(l1norm / maxnorm, &squarings);
256 if (squarings < 0) squarings = 0;
257 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<RealScalar>(squarings));
258 matrix_exp_pade13(A, U, V);
263 template <
typename MatrixType>
264 struct matrix_exp_computeUV<MatrixType, long double>
266 template <
typename ArgType>
267 static void run(
const ArgType& arg, MatrixType& U, MatrixType& V,
int& squarings)
269 #if LDBL_MANT_DIG == 53 // double precision 270 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
276 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
279 #if LDBL_MANT_DIG <= 64 // extended precision 281 if (l1norm < 4.1968497232266989671e-003L) {
282 matrix_exp_pade3(arg, U, V);
283 }
else if (l1norm < 1.1848116734693823091e-001L) {
284 matrix_exp_pade5(arg, U, V);
285 }
else if (l1norm < 5.5170388480686700274e-001L) {
286 matrix_exp_pade7(arg, U, V);
287 }
else if (l1norm < 1.3759868875587845383e+000L) {
288 matrix_exp_pade9(arg, U, V);
290 const long double maxnorm = 4.0246098906697353063L;
291 frexp(l1norm / maxnorm, &squarings);
292 if (squarings < 0) squarings = 0;
293 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
294 matrix_exp_pade13(A, U, V);
297 #elif LDBL_MANT_DIG <= 106 // double-double 299 if (l1norm < 3.2787892205607026992947488108213e-005L) {
300 matrix_exp_pade3(arg, U, V);
301 }
else if (l1norm < 6.4467025060072760084130906076332e-003L) {
302 matrix_exp_pade5(arg, U, V);
303 }
else if (l1norm < 6.8988028496595374751374122881143e-002L) {
304 matrix_exp_pade7(arg, U, V);
305 }
else if (l1norm < 2.7339737518502231741495857201670e-001L) {
306 matrix_exp_pade9(arg, U, V);
307 }
else if (l1norm < 1.3203382096514474905666448850278e+000L) {
308 matrix_exp_pade13(arg, U, V);
310 const long double maxnorm = 3.2579440895405400856599663723517L;
311 frexp(l1norm / maxnorm, &squarings);
312 if (squarings < 0) squarings = 0;
313 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
314 matrix_exp_pade17(A, U, V);
317 #elif LDBL_MANT_DIG <= 113 // quadruple precision 319 if (l1norm < 1.639394610288918690547467954466970e-005L) {
320 matrix_exp_pade3(arg, U, V);
321 }
else if (l1norm < 4.253237712165275566025884344433009e-003L) {
322 matrix_exp_pade5(arg, U, V);
323 }
else if (l1norm < 5.125804063165764409885122032933142e-002L) {
324 matrix_exp_pade7(arg, U, V);
325 }
else if (l1norm < 2.170000765161155195453205651889853e-001L) {
326 matrix_exp_pade9(arg, U, V);
327 }
else if (l1norm < 1.125358383453143065081397882891878e+000L) {
328 matrix_exp_pade13(arg, U, V);
330 const long double maxnorm = 2.884233277829519311757165057717815L;
331 frexp(l1norm / maxnorm, &squarings);
332 if (squarings < 0) squarings = 0;
333 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
334 matrix_exp_pade17(A, U, V);
340 eigen_assert(
false &&
"Bug in MatrixExponential");
343 #endif // LDBL_MANT_DIG 347 template<
typename T>
struct is_exp_known_type : false_type {};
348 template<>
struct is_exp_known_type<float> : true_type {};
349 template<>
struct is_exp_known_type<double> : true_type {};
350 #if LDBL_MANT_DIG <= 113 351 template<>
struct is_exp_known_type<long double> : true_type {};
354 template <
typename ArgType,
typename ResultType>
355 void matrix_exp_compute(
const ArgType& arg, ResultType &result, true_type)
357 typedef typename ArgType::PlainObject MatrixType;
360 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings);
361 MatrixType numer = U + V;
362 MatrixType denom = -U + V;
363 result = denom.partialPivLu().solve(numer);
364 for (
int i=0; i<squarings; i++)
374 template <
typename ArgType,
typename ResultType>
375 void matrix_exp_compute(
const ArgType& arg, ResultType &result, false_type)
377 typedef typename ArgType::PlainObject MatrixType;
378 typedef typename traits<MatrixType>::Scalar Scalar;
379 typedef typename NumTraits<Scalar>::Real RealScalar;
380 typedef typename std::complex<RealScalar> ComplexScalar;
381 result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
397 :
public ReturnByValue<MatrixExponentialReturnValue<Derived> >
410 template <
typename ResultType>
411 inline void evalTo(ResultType& result)
const 413 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
414 internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
417 Index rows()
const {
return m_src.rows(); }
418 Index cols()
const {
return m_src.cols(); }
421 const typename internal::ref_selector<Derived>::type m_src;
425 template<
typename Derived>
426 struct traits<MatrixExponentialReturnValue<Derived> >
428 typedef typename Derived::PlainObject ReturnType;
432 template <
typename Derived>
435 eigen_assert(rows() == cols());
441 #endif // EIGEN_MATRIX_EXPONENTIAL const Eigen::CwiseUnaryOp< Eigen::internal::scalar_arg_op< typename Derived::Scalar >, const Derived > arg(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:411
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
const MatrixExponentialReturnValue< Derived > exp() const
Definition: MatrixExponential.h:433
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:396
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:404