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>
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>
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>
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>
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>
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) {
222 }
else if (l1norm < 1.880152677804762e+000f) {
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));
234 template <
typename MatrixType>
235 struct matrix_exp_computeUV<MatrixType, double>
237 template <
typename ArgType>
238 static void run(
const ArgType& arg, MatrixType& U, MatrixType& V,
int& squarings)
242 const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
244 if (l1norm < 1.495585217958292e-002) {
246 }
else if (l1norm < 2.539398330063230e-001) {
248 }
else if (l1norm < 9.504178996162932e-001) {
250 }
else if (l1norm < 2.097847961257068e+000) {
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));
262 template <
typename MatrixType>
263 struct matrix_exp_computeUV<MatrixType, long double>
265 template <
typename ArgType>
266 static void run(
const ArgType& arg, MatrixType& U, MatrixType& V,
int& squarings)
268 #if LDBL_MANT_DIG == 53 // double precision 269 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
275 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
278 #if LDBL_MANT_DIG <= 64 // extended precision 280 if (l1norm < 4.1968497232266989671e-003L) {
282 }
else if (l1norm < 1.1848116734693823091e-001L) {
284 }
else if (l1norm < 5.5170388480686700274e-001L) {
286 }
else if (l1norm < 1.3759868875587845383e+000L) {
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));
296 #elif LDBL_MANT_DIG <= 106 // double-double 298 if (l1norm < 3.2787892205607026992947488108213e-005L) {
300 }
else if (l1norm < 6.4467025060072760084130906076332e-003L) {
302 }
else if (l1norm < 6.8988028496595374751374122881143e-002L) {
304 }
else if (l1norm < 2.7339737518502231741495857201670e-001L) {
306 }
else if (l1norm < 1.3203382096514474905666448850278e+000L) {
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);
316 #elif LDBL_MANT_DIG <= 112 // quadruple precison 318 if (l1norm < 1.639394610288918690547467954466970e-005L) {
320 }
else if (l1norm < 4.253237712165275566025884344433009e-003L) {
322 }
else if (l1norm < 5.125804063165764409885122032933142e-002L) {
324 }
else if (l1norm < 2.170000765161155195453205651889853e-001L) {
326 }
else if (l1norm < 1.125358383453143065081397882891878e+000L) {
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);
338 eigen_assert(
false &&
"Bug in MatrixExponential");
341 #endif // LDBL_MANT_DIG 351 template <
typename ArgType,
typename ResultType>
352 void matrix_exp_compute(
const ArgType& arg, ResultType &result)
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>);
366 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings);
367 MatrixType numer = U + V;
368 MatrixType denom = -U + V;
369 result = denom.partialPivLu().solve(numer);
370 for (
int i=0; i<squarings; i++)
387 :
public ReturnByValue<MatrixExponentialReturnValue<Derived> >
389 typedef typename Derived::Index Index;
401 template <
typename ResultType>
402 inline void evalTo(ResultType& result)
const 404 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
405 internal::matrix_exp_compute(tmp, result);
408 Index rows()
const {
return m_src.rows(); }
409 Index cols()
const {
return m_src.cols(); }
412 const typename internal::ref_selector<Derived>::type m_src;
416 template<
typename Derived>
417 struct traits<MatrixExponentialReturnValue<Derived> >
419 typedef typename Derived::PlainObject ReturnType;
423 template <
typename Derived>
426 eigen_assert(rows() == cols());
432 #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:402
void matrix_exp_pade13(const MatA &A, MatU &U, MatV &V)
Compute the (13,13)-Padé approximant to the exponential.
Definition: MatrixExponential.h:143
void matrix_exp_pade3(const MatA &A, MatU &U, MatV &V)
Compute the (3,3)-Padé approximant to the exponential.
Definition: MatrixExponential.h:65
void matrix_exp_pade5(const MatA &A, MatU &U, MatV &V)
Compute the (5,5)-Padé approximant to the exponential.
Definition: MatrixExponential.h:82
void matrix_exp_pade9(const MatA &A, MatU &U, MatV &V)
Compute the (9,9)-Padé approximant to the exponential.
Definition: MatrixExponential.h:121
void matrix_exp_pade7(const MatA &A, MatU &U, MatV &V)
Compute the (7,7)-Padé approximant to the exponential.
Definition: MatrixExponential.h:100
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:386
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:395