![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk> 00005 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@ms63.hinet.net> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_MATRIX_EXPONENTIAL 00012 #define EIGEN_MATRIX_EXPONENTIAL 00013 00014 #include "StemFunction.h" 00015 00016 namespace Eigen { 00017 namespace internal { 00018 00023 template <typename RealScalar> 00024 struct MatrixExponentialScalingOp 00025 { 00030 MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { } 00031 00032 00037 inline const RealScalar operator() (const RealScalar& x) const 00038 { 00039 using std::ldexp; 00040 return ldexp(x, -m_squarings); 00041 } 00042 00043 typedef std::complex<RealScalar> ComplexScalar; 00044 00049 inline const ComplexScalar operator() (const ComplexScalar& x) const 00050 { 00051 using std::ldexp; 00052 return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings)); 00053 } 00054 00055 private: 00056 int m_squarings; 00057 }; 00058 00064 template <typename MatA, typename MatU, typename MatV> 00065 void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V) 00066 { 00067 typedef typename MatA::PlainObject MatrixType; 00068 typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar; 00069 const RealScalar b[] = {120.L, 60.L, 12.L, 1.L}; 00070 const MatrixType A2 = A * A; 00071 const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00072 U.noalias() = A * tmp; 00073 V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00074 } 00075 00081 template <typename MatA, typename MatU, typename MatV> 00082 void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V) 00083 { 00084 typedef typename MatA::PlainObject MatrixType; 00085 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; 00086 const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L}; 00087 const MatrixType A2 = A * A; 00088 const MatrixType A4 = A2 * A2; 00089 const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00090 U.noalias() = A * tmp; 00091 V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00092 } 00093 00099 template <typename MatA, typename MatU, typename MatV> 00100 void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V) 00101 { 00102 typedef typename MatA::PlainObject MatrixType; 00103 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; 00104 const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L}; 00105 const MatrixType A2 = A * A; 00106 const MatrixType A4 = A2 * A2; 00107 const MatrixType A6 = A4 * A2; 00108 const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 00109 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00110 U.noalias() = A * tmp; 00111 V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00112 00113 } 00114 00120 template <typename MatA, typename MatU, typename MatV> 00121 void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V) 00122 { 00123 typedef typename MatA::PlainObject MatrixType; 00124 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; 00125 const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L, 00126 2162160.L, 110880.L, 3960.L, 90.L, 1.L}; 00127 const MatrixType A2 = A * A; 00128 const MatrixType A4 = A2 * A2; 00129 const MatrixType A6 = A4 * A2; 00130 const MatrixType A8 = A6 * A2; 00131 const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 00132 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00133 U.noalias() = A * tmp; 00134 V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00135 } 00136 00142 template <typename MatA, typename MatU, typename MatV> 00143 void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V) 00144 { 00145 typedef typename MatA::PlainObject MatrixType; 00146 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; 00147 const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L, 00148 1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L, 00149 33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L}; 00150 const MatrixType A2 = A * A; 00151 const MatrixType A4 = A2 * A2; 00152 const MatrixType A6 = A4 * A2; 00153 V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage 00154 MatrixType tmp = A6 * V; 00155 tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00156 U.noalias() = A * tmp; 00157 tmp = b[12] * A6 + b[10] * A4 + b[8] * A2; 00158 V.noalias() = A6 * tmp; 00159 V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00160 } 00161 00169 #if LDBL_MANT_DIG > 64 00170 template <typename MatA, typename MatU, typename MatV> 00171 void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V) 00172 { 00173 typedef typename MatA::PlainObject MatrixType; 00174 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; 00175 const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L, 00176 100610229646136770560000.L, 15720348382208870400000.L, 00177 1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L, 00178 595373117923584000.L, 27563570274240000.L, 1060137318240000.L, 00179 33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L, 00180 46512.L, 306.L, 1.L}; 00181 const MatrixType A2 = A * A; 00182 const MatrixType A4 = A2 * A2; 00183 const MatrixType A6 = A4 * A2; 00184 const MatrixType A8 = A4 * A4; 00185 V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage 00186 MatrixType tmp = A8 * V; 00187 tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 00188 + b[1] * MatrixType::Identity(A.rows(), A.cols()); 00189 U.noalias() = A * tmp; 00190 tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2; 00191 V.noalias() = tmp * A8; 00192 V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 00193 + b[0] * MatrixType::Identity(A.rows(), A.cols()); 00194 } 00195 #endif 00196 00197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real> 00198 struct matrix_exp_computeUV 00199 { 00207 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings); 00208 }; 00209 00210 template <typename MatrixType> 00211 struct matrix_exp_computeUV<MatrixType, float> 00212 { 00213 template <typename ArgType> 00214 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) 00215 { 00216 using std::frexp; 00217 using std::pow; 00218 const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); 00219 squarings = 0; 00220 if (l1norm < 4.258730016922831e-001f) { 00221 matrix_exp_pade3(arg, U, V); 00222 } else if (l1norm < 1.880152677804762e+000f) { 00223 matrix_exp_pade5(arg, U, V); 00224 } else { 00225 const float maxnorm = 3.925724783138660f; 00226 frexp(l1norm / maxnorm, &squarings); 00227 if (squarings < 0) squarings = 0; 00228 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings)); 00229 matrix_exp_pade7(A, U, V); 00230 } 00231 } 00232 }; 00233 00234 template <typename MatrixType> 00235 struct matrix_exp_computeUV<MatrixType, double> 00236 { 00237 template <typename ArgType> 00238 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) 00239 { 00240 using std::frexp; 00241 using std::pow; 00242 const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); 00243 squarings = 0; 00244 if (l1norm < 1.495585217958292e-002) { 00245 matrix_exp_pade3(arg, U, V); 00246 } else if (l1norm < 2.539398330063230e-001) { 00247 matrix_exp_pade5(arg, U, V); 00248 } else if (l1norm < 9.504178996162932e-001) { 00249 matrix_exp_pade7(arg, U, V); 00250 } else if (l1norm < 2.097847961257068e+000) { 00251 matrix_exp_pade9(arg, U, V); 00252 } else { 00253 const double maxnorm = 5.371920351148152; 00254 frexp(l1norm / maxnorm, &squarings); 00255 if (squarings < 0) squarings = 0; 00256 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<double>(squarings)); 00257 matrix_exp_pade13(A, U, V); 00258 } 00259 } 00260 }; 00261 00262 template <typename MatrixType> 00263 struct matrix_exp_computeUV<MatrixType, long double> 00264 { 00265 template <typename ArgType> 00266 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) 00267 { 00268 #if LDBL_MANT_DIG == 53 // double precision 00269 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings); 00270 00271 #else 00272 00273 using std::frexp; 00274 using std::pow; 00275 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); 00276 squarings = 0; 00277 00278 #if LDBL_MANT_DIG <= 64 // extended precision 00279 00280 if (l1norm < 4.1968497232266989671e-003L) { 00281 matrix_exp_pade3(arg, U, V); 00282 } else if (l1norm < 1.1848116734693823091e-001L) { 00283 matrix_exp_pade5(arg, U, V); 00284 } else if (l1norm < 5.5170388480686700274e-001L) { 00285 matrix_exp_pade7(arg, U, V); 00286 } else if (l1norm < 1.3759868875587845383e+000L) { 00287 matrix_exp_pade9(arg, U, V); 00288 } else { 00289 const long double maxnorm = 4.0246098906697353063L; 00290 frexp(l1norm / maxnorm, &squarings); 00291 if (squarings < 0) squarings = 0; 00292 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings)); 00293 matrix_exp_pade13(A, U, V); 00294 } 00295 00296 #elif LDBL_MANT_DIG <= 106 // double-double 00297 00298 if (l1norm < 3.2787892205607026992947488108213e-005L) { 00299 matrix_exp_pade3(arg, U, V); 00300 } else if (l1norm < 6.4467025060072760084130906076332e-003L) { 00301 matrix_exp_pade5(arg, U, V); 00302 } else if (l1norm < 6.8988028496595374751374122881143e-002L) { 00303 matrix_exp_pade7(arg, U, V); 00304 } else if (l1norm < 2.7339737518502231741495857201670e-001L) { 00305 matrix_exp_pade9(arg, U, V); 00306 } else if (l1norm < 1.3203382096514474905666448850278e+000L) { 00307 matrix_exp_pade13(arg, U, V); 00308 } else { 00309 const long double maxnorm = 3.2579440895405400856599663723517L; 00310 frexp(l1norm / maxnorm, &squarings); 00311 if (squarings < 0) squarings = 0; 00312 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings)); 00313 matrix_exp_pade17(A, U, V); 00314 } 00315 00316 #elif LDBL_MANT_DIG <= 112 // quadruple precison 00317 00318 if (l1norm < 1.639394610288918690547467954466970e-005L) { 00319 matrix_exp_pade3(arg, U, V); 00320 } else if (l1norm < 4.253237712165275566025884344433009e-003L) { 00321 matrix_exp_pade5(arg, U, V); 00322 } else if (l1norm < 5.125804063165764409885122032933142e-002L) { 00323 matrix_exp_pade7(arg, U, V); 00324 } else if (l1norm < 2.170000765161155195453205651889853e-001L) { 00325 matrix_exp_pade9(arg, U, V); 00326 } else if (l1norm < 1.125358383453143065081397882891878e+000L) { 00327 matrix_exp_pade13(arg, U, V); 00328 } else { 00329 frexp(l1norm / maxnorm, &squarings); 00330 if (squarings < 0) squarings = 0; 00331 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings)); 00332 matrix_exp_pade17(A, U, V); 00333 } 00334 00335 #else 00336 00337 // this case should be handled in compute() 00338 eigen_assert(false && "Bug in MatrixExponential"); 00339 00340 #endif 00341 #endif // LDBL_MANT_DIG 00342 } 00343 }; 00344 00345 00346 /* Computes the matrix exponential 00347 * 00348 * \param arg argument of matrix exponential (should be plain object) 00349 * \param result variable in which result will be stored 00350 */ 00351 template <typename ArgType, typename ResultType> 00352 void matrix_exp_compute(const ArgType& arg, ResultType &result) 00353 { 00354 typedef typename ArgType::PlainObject MatrixType; 00355 #if LDBL_MANT_DIG > 112 // rarely happens 00356 typedef typename traits<MatrixType>::Scalar Scalar; 00357 typedef typename NumTraits<Scalar>::Real RealScalar; 00358 typedef typename std::complex<RealScalar> ComplexScalar; 00359 if (sizeof(RealScalar) > 14) { 00360 result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>); 00361 return; 00362 } 00363 #endif 00364 MatrixType U, V; 00365 int squarings; 00366 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V) 00367 MatrixType numer = U + V; 00368 MatrixType denom = -U + V; 00369 result = denom.partialPivLu().solve(numer); 00370 for (int i=0; i<squarings; i++) 00371 result *= result; // undo scaling by repeated squaring 00372 } 00373 00374 } // end namespace Eigen::internal 00375 00386 template<typename Derived> struct MatrixExponentialReturnValue 00387 : public ReturnByValue<MatrixExponentialReturnValue<Derived> > 00388 { 00389 typedef typename Derived::Index Index; 00390 public: 00395 MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } 00396 00401 template <typename ResultType> 00402 inline void evalTo(ResultType& result) const 00403 { 00404 const typename internal::nested_eval<Derived, 10>::type tmp(m_src); 00405 internal::matrix_exp_compute(tmp, result); 00406 } 00407 00408 Index rows() const { return m_src.rows(); } 00409 Index cols() const { return m_src.cols(); } 00410 00411 protected: 00412 const typename internal::ref_selector<Derived>::type m_src; 00413 }; 00414 00415 namespace internal { 00416 template<typename Derived> 00417 struct traits<MatrixExponentialReturnValue<Derived> > 00418 { 00419 typedef typename Derived::PlainObject ReturnType; 00420 }; 00421 } 00422 00423 template <typename Derived> 00424 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 00425 { 00426 eigen_assert(rows() == cols()); 00427 return MatrixExponentialReturnValue<Derived>(derived()); 00428 } 00429 00430 } // end namespace Eigen 00431 00432 #endif // EIGEN_MATRIX_EXPONENTIAL