$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009, 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 00005 // Copyright (C) 2011 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 00023 template <typename MatrixType> 00024 class MatrixExponential { 00025 00026 public: 00027 00035 MatrixExponential(const MatrixType &M); 00036 00041 template <typename ResultType> 00042 void compute(ResultType &result); 00043 00044 private: 00045 00046 // Prevent copying 00047 MatrixExponential(const MatrixExponential&); 00048 MatrixExponential& operator=(const MatrixExponential&); 00049 00057 void pade3(const MatrixType &A); 00058 00066 void pade5(const MatrixType &A); 00067 00075 void pade7(const MatrixType &A); 00076 00084 void pade9(const MatrixType &A); 00085 00093 void pade13(const MatrixType &A); 00094 00104 void pade17(const MatrixType &A); 00105 00119 void computeUV(double); 00120 00125 void computeUV(float); 00126 00131 void computeUV(long double); 00132 00133 typedef typename internal::traits<MatrixType>::Scalar Scalar; 00134 typedef typename NumTraits<Scalar>::Real RealScalar; 00135 typedef typename std::complex<RealScalar> ComplexScalar; 00136 00138 typename internal::nested<MatrixType>::type m_M; 00139 00141 MatrixType m_U; 00142 00144 MatrixType m_V; 00145 00147 MatrixType m_tmp1; 00148 00150 MatrixType m_tmp2; 00151 00153 MatrixType m_Id; 00154 00156 int m_squarings; 00157 00159 RealScalar m_l1norm; 00160 }; 00161 00162 template <typename MatrixType> 00163 MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) : 00164 m_M(M), 00165 m_U(M.rows(),M.cols()), 00166 m_V(M.rows(),M.cols()), 00167 m_tmp1(M.rows(),M.cols()), 00168 m_tmp2(M.rows(),M.cols()), 00169 m_Id(MatrixType::Identity(M.rows(), M.cols())), 00170 m_squarings(0), 00171 m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff()) 00172 { 00173 /* empty body */ 00174 } 00175 00176 template <typename MatrixType> 00177 template <typename ResultType> 00178 void MatrixExponential<MatrixType>::compute(ResultType &result) 00179 { 00180 #if LDBL_MANT_DIG > 112 // rarely happens 00181 if(sizeof(RealScalar) > 14) { 00182 result = m_M.matrixFunction(StdStemFunctions<ComplexScalar>::exp); 00183 return; 00184 } 00185 #endif 00186 computeUV(RealScalar()); 00187 m_tmp1 = m_U + m_V; // numerator of Pade approximant 00188 m_tmp2 = -m_U + m_V; // denominator of Pade approximant 00189 result = m_tmp2.partialPivLu().solve(m_tmp1); 00190 for (int i=0; i<m_squarings; i++) 00191 result *= result; // undo scaling by repeated squaring 00192 } 00193 00194 template <typename MatrixType> 00195 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A) 00196 { 00197 const RealScalar b[] = {120., 60., 12., 1.}; 00198 m_tmp1.noalias() = A * A; 00199 m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id; 00200 m_U.noalias() = A * m_tmp2; 00201 m_V = b[2]*m_tmp1 + b[0]*m_Id; 00202 } 00203 00204 template <typename MatrixType> 00205 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A) 00206 { 00207 const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.}; 00208 MatrixType A2 = A * A; 00209 m_tmp1.noalias() = A2 * A2; 00210 m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id; 00211 m_U.noalias() = A * m_tmp2; 00212 m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id; 00213 } 00214 00215 template <typename MatrixType> 00216 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A) 00217 { 00218 const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; 00219 MatrixType A2 = A * A; 00220 MatrixType A4 = A2 * A2; 00221 m_tmp1.noalias() = A4 * A2; 00222 m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 00223 m_U.noalias() = A * m_tmp2; 00224 m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 00225 } 00226 00227 template <typename MatrixType> 00228 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A) 00229 { 00230 const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., 00231 2162160., 110880., 3960., 90., 1.}; 00232 MatrixType A2 = A * A; 00233 MatrixType A4 = A2 * A2; 00234 MatrixType A6 = A4 * A2; 00235 m_tmp1.noalias() = A6 * A2; 00236 m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 00237 m_U.noalias() = A * m_tmp2; 00238 m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 00239 } 00240 00241 template <typename MatrixType> 00242 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A) 00243 { 00244 const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., 00245 1187353796428800., 129060195264000., 10559470521600., 670442572800., 00246 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; 00247 MatrixType A2 = A * A; 00248 MatrixType A4 = A2 * A2; 00249 m_tmp1.noalias() = A4 * A2; 00250 m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage 00251 m_tmp2.noalias() = m_tmp1 * m_V; 00252 m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 00253 m_U.noalias() = A * m_tmp2; 00254 m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2; 00255 m_V.noalias() = m_tmp1 * m_tmp2; 00256 m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 00257 } 00258 00259 #if LDBL_MANT_DIG > 64 00260 template <typename MatrixType> 00261 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade17(const MatrixType &A) 00262 { 00263 const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L, 00264 100610229646136770560000.L, 15720348382208870400000.L, 00265 1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L, 00266 595373117923584000.L, 27563570274240000.L, 1060137318240000.L, 00267 33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L, 00268 46512.L, 306.L, 1.L}; 00269 MatrixType A2 = A * A; 00270 MatrixType A4 = A2 * A2; 00271 MatrixType A6 = A4 * A2; 00272 m_tmp1.noalias() = A4 * A4; 00273 m_V = b[17]*m_tmp1 + b[15]*A6 + b[13]*A4 + b[11]*A2; // used for temporary storage 00274 m_tmp2.noalias() = m_tmp1 * m_V; 00275 m_tmp2 += b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 00276 m_U.noalias() = A * m_tmp2; 00277 m_tmp2 = b[16]*m_tmp1 + b[14]*A6 + b[12]*A4 + b[10]*A2; 00278 m_V.noalias() = m_tmp1 * m_tmp2; 00279 m_V += b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 00280 } 00281 #endif 00282 00283 template <typename MatrixType> 00284 void MatrixExponential<MatrixType>::computeUV(float) 00285 { 00286 using std::frexp; 00287 using std::pow; 00288 if (m_l1norm < 4.258730016922831e-001) { 00289 pade3(m_M); 00290 } else if (m_l1norm < 1.880152677804762e+000) { 00291 pade5(m_M); 00292 } else { 00293 const float maxnorm = 3.925724783138660f; 00294 frexp(m_l1norm / maxnorm, &m_squarings); 00295 if (m_squarings < 0) m_squarings = 0; 00296 MatrixType A = m_M / pow(Scalar(2), m_squarings); 00297 pade7(A); 00298 } 00299 } 00300 00301 template <typename MatrixType> 00302 void MatrixExponential<MatrixType>::computeUV(double) 00303 { 00304 using std::frexp; 00305 using std::pow; 00306 if (m_l1norm < 1.495585217958292e-002) { 00307 pade3(m_M); 00308 } else if (m_l1norm < 2.539398330063230e-001) { 00309 pade5(m_M); 00310 } else if (m_l1norm < 9.504178996162932e-001) { 00311 pade7(m_M); 00312 } else if (m_l1norm < 2.097847961257068e+000) { 00313 pade9(m_M); 00314 } else { 00315 const double maxnorm = 5.371920351148152; 00316 frexp(m_l1norm / maxnorm, &m_squarings); 00317 if (m_squarings < 0) m_squarings = 0; 00318 MatrixType A = m_M / pow(Scalar(2), m_squarings); 00319 pade13(A); 00320 } 00321 } 00322 00323 template <typename MatrixType> 00324 void MatrixExponential<MatrixType>::computeUV(long double) 00325 { 00326 using std::frexp; 00327 using std::pow; 00328 #if LDBL_MANT_DIG == 53 // double precision 00329 computeUV(double()); 00330 #elif LDBL_MANT_DIG <= 64 // extended precision 00331 if (m_l1norm < 4.1968497232266989671e-003L) { 00332 pade3(m_M); 00333 } else if (m_l1norm < 1.1848116734693823091e-001L) { 00334 pade5(m_M); 00335 } else if (m_l1norm < 5.5170388480686700274e-001L) { 00336 pade7(m_M); 00337 } else if (m_l1norm < 1.3759868875587845383e+000L) { 00338 pade9(m_M); 00339 } else { 00340 const long double maxnorm = 4.0246098906697353063L; 00341 frexp(m_l1norm / maxnorm, &m_squarings); 00342 if (m_squarings < 0) m_squarings = 0; 00343 MatrixType A = m_M / pow(Scalar(2), m_squarings); 00344 pade13(A); 00345 } 00346 #elif LDBL_MANT_DIG <= 106 // double-double 00347 if (m_l1norm < 3.2787892205607026992947488108213e-005L) { 00348 pade3(m_M); 00349 } else if (m_l1norm < 6.4467025060072760084130906076332e-003L) { 00350 pade5(m_M); 00351 } else if (m_l1norm < 6.8988028496595374751374122881143e-002L) { 00352 pade7(m_M); 00353 } else if (m_l1norm < 2.7339737518502231741495857201670e-001L) { 00354 pade9(m_M); 00355 } else if (m_l1norm < 1.3203382096514474905666448850278e+000L) { 00356 pade13(m_M); 00357 } else { 00358 const long double maxnorm = 3.2579440895405400856599663723517L; 00359 frexp(m_l1norm / maxnorm, &m_squarings); 00360 if (m_squarings < 0) m_squarings = 0; 00361 MatrixType A = m_M / pow(Scalar(2), m_squarings); 00362 pade17(A); 00363 } 00364 #elif LDBL_MANT_DIG <= 112 // quadruple precison 00365 if (m_l1norm < 1.639394610288918690547467954466970e-005L) { 00366 pade3(m_M); 00367 } else if (m_l1norm < 4.253237712165275566025884344433009e-003L) { 00368 pade5(m_M); 00369 } else if (m_l1norm < 5.125804063165764409885122032933142e-002L) { 00370 pade7(m_M); 00371 } else if (m_l1norm < 2.170000765161155195453205651889853e-001L) { 00372 pade9(m_M); 00373 } else if (m_l1norm < 1.125358383453143065081397882891878e+000L) { 00374 pade13(m_M); 00375 } else { 00376 const long double maxnorm = 2.884233277829519311757165057717815L; 00377 frexp(m_l1norm / maxnorm, &m_squarings); 00378 if (m_squarings < 0) m_squarings = 0; 00379 MatrixType A = m_M / pow(Scalar(2), m_squarings); 00380 pade17(A); 00381 } 00382 #else 00383 // this case should be handled in compute() 00384 eigen_assert(false && "Bug in MatrixExponential"); 00385 #endif // LDBL_MANT_DIG 00386 } 00387 00400 template<typename Derived> struct MatrixExponentialReturnValue 00401 : public ReturnByValue<MatrixExponentialReturnValue<Derived> > 00402 { 00403 typedef typename Derived::Index Index; 00404 public: 00410 MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } 00411 00417 template <typename ResultType> 00418 inline void evalTo(ResultType& result) const 00419 { 00420 const typename Derived::PlainObject srcEvaluated = m_src.eval(); 00421 MatrixExponential<typename Derived::PlainObject> me(srcEvaluated); 00422 me.compute(result); 00423 } 00424 00425 Index rows() const { return m_src.rows(); } 00426 Index cols() const { return m_src.cols(); } 00427 00428 protected: 00429 const Derived& m_src; 00430 private: 00431 MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&); 00432 }; 00433 00434 namespace internal { 00435 template<typename Derived> 00436 struct traits<MatrixExponentialReturnValue<Derived> > 00437 { 00438 typedef typename Derived::PlainObject ReturnType; 00439 }; 00440 } 00441 00442 template <typename Derived> 00443 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 00444 { 00445 eigen_assert(rows() == cols()); 00446 return MatrixExponentialReturnValue<Derived>(derived()); 00447 } 00448 00449 } // end namespace Eigen 00450 00451 #endif // EIGEN_MATRIX_EXPONENTIAL