MatrixExponential.h
Go to the documentation of this file.
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 
00018 #if defined(_MSC_VER) || defined(__FreeBSD__)
00019   template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
00020 #endif
00021 
00022 
00028 template <typename MatrixType>
00029 class MatrixExponential {
00030 
00031   public:
00032 
00040     MatrixExponential(const MatrixType &M);
00041 
00046     template <typename ResultType> 
00047     void compute(ResultType &result);
00048 
00049   private:
00050 
00051     // Prevent copying
00052     MatrixExponential(const MatrixExponential&);
00053     MatrixExponential& operator=(const MatrixExponential&);
00054 
00062     void pade3(const MatrixType &A);
00063 
00071     void pade5(const MatrixType &A);
00072 
00080     void pade7(const MatrixType &A);
00081 
00089     void pade9(const MatrixType &A);
00090 
00098     void pade13(const MatrixType &A);
00099 
00109     void pade17(const MatrixType &A);
00110 
00124     void computeUV(double);
00125 
00130     void computeUV(float);
00131     
00136     void computeUV(long double);
00137 
00138     typedef typename internal::traits<MatrixType>::Scalar Scalar;
00139     typedef typename NumTraits<Scalar>::Real RealScalar;
00140     typedef typename std::complex<RealScalar> ComplexScalar;
00141 
00143     typename internal::nested<MatrixType>::type m_M;
00144 
00146     MatrixType m_U;
00147 
00149     MatrixType m_V;
00150 
00152     MatrixType m_tmp1;
00153 
00155     MatrixType m_tmp2;
00156 
00158     MatrixType m_Id;
00159 
00161     int m_squarings;
00162 
00164     RealScalar m_l1norm;
00165 };
00166 
00167 template <typename MatrixType>
00168 MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) :
00169   m_M(M),
00170   m_U(M.rows(),M.cols()),
00171   m_V(M.rows(),M.cols()),
00172   m_tmp1(M.rows(),M.cols()),
00173   m_tmp2(M.rows(),M.cols()),
00174   m_Id(MatrixType::Identity(M.rows(), M.cols())),
00175   m_squarings(0),
00176   m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff())
00177 {
00178   /* empty body */
00179 }
00180 
00181 template <typename MatrixType>
00182 template <typename ResultType> 
00183 void MatrixExponential<MatrixType>::compute(ResultType &result)
00184 {
00185 #if LDBL_MANT_DIG > 112 // rarely happens
00186   if(sizeof(RealScalar) > 14) {
00187     result = m_M.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
00188     return;
00189   }
00190 #endif
00191   computeUV(RealScalar());
00192   m_tmp1 = m_U + m_V;   // numerator of Pade approximant
00193   m_tmp2 = -m_U + m_V;  // denominator of Pade approximant
00194   result = m_tmp2.partialPivLu().solve(m_tmp1);
00195   for (int i=0; i<m_squarings; i++)
00196     result *= result;   // undo scaling by repeated squaring
00197 }
00198 
00199 template <typename MatrixType>
00200 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A)
00201 {
00202   const RealScalar b[] = {120., 60., 12., 1.};
00203   m_tmp1.noalias() = A * A;
00204   m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id;
00205   m_U.noalias() = A * m_tmp2;
00206   m_V = b[2]*m_tmp1 + b[0]*m_Id;
00207 }
00208 
00209 template <typename MatrixType>
00210 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A)
00211 {
00212   const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.};
00213   MatrixType A2 = A * A;
00214   m_tmp1.noalias() = A2 * A2;
00215   m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id;
00216   m_U.noalias() = A * m_tmp2;
00217   m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id;
00218 }
00219 
00220 template <typename MatrixType>
00221 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A)
00222 {
00223   const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
00224   MatrixType A2 = A * A;
00225   MatrixType A4 = A2 * A2;
00226   m_tmp1.noalias() = A4 * A2;
00227   m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00228   m_U.noalias() = A * m_tmp2;
00229   m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00230 }
00231 
00232 template <typename MatrixType>
00233 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A)
00234 {
00235   const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
00236                       2162160., 110880., 3960., 90., 1.};
00237   MatrixType A2 = A * A;
00238   MatrixType A4 = A2 * A2;
00239   MatrixType A6 = A4 * A2;
00240   m_tmp1.noalias() = A6 * A2;
00241   m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00242   m_U.noalias() = A * m_tmp2;
00243   m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00244 }
00245 
00246 template <typename MatrixType>
00247 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A)
00248 {
00249   const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
00250                       1187353796428800., 129060195264000., 10559470521600., 670442572800.,
00251                       33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
00252   MatrixType A2 = A * A;
00253   MatrixType A4 = A2 * A2;
00254   m_tmp1.noalias() = A4 * A2;
00255   m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage
00256   m_tmp2.noalias() = m_tmp1 * m_V;
00257   m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00258   m_U.noalias() = A * m_tmp2;
00259   m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2;
00260   m_V.noalias() = m_tmp1 * m_tmp2;
00261   m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00262 }
00263 
00264 #if LDBL_MANT_DIG > 64
00265 template <typename MatrixType>
00266 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade17(const MatrixType &A)
00267 {
00268   const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
00269             100610229646136770560000.L, 15720348382208870400000.L,
00270             1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
00271             595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
00272             33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
00273             46512.L, 306.L, 1.L};
00274   MatrixType A2 = A * A;
00275   MatrixType A4 = A2 * A2;
00276   MatrixType A6 = A4 * A2;
00277   m_tmp1.noalias() = A4 * A4;
00278   m_V = b[17]*m_tmp1 + b[15]*A6 + b[13]*A4 + b[11]*A2; // used for temporary storage
00279   m_tmp2.noalias() = m_tmp1 * m_V;
00280   m_tmp2 += b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00281   m_U.noalias() = A * m_tmp2;
00282   m_tmp2 = b[16]*m_tmp1 + b[14]*A6 + b[12]*A4 + b[10]*A2;
00283   m_V.noalias() = m_tmp1 * m_tmp2;
00284   m_V += b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00285 }
00286 #endif
00287 
00288 template <typename MatrixType>
00289 void MatrixExponential<MatrixType>::computeUV(float)
00290 {
00291   using std::max;
00292   using std::pow;
00293   using std::ceil;
00294   if (m_l1norm < 4.258730016922831e-001) {
00295     pade3(m_M);
00296   } else if (m_l1norm < 1.880152677804762e+000) {
00297     pade5(m_M);
00298   } else {
00299     const float maxnorm = 3.925724783138660f;
00300     m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
00301     MatrixType A = m_M / pow(Scalar(2), m_squarings);
00302     pade7(A);
00303   }
00304 }
00305 
00306 template <typename MatrixType>
00307 void MatrixExponential<MatrixType>::computeUV(double)
00308 {
00309   using std::max;
00310   using std::pow;
00311   using std::ceil;
00312   if (m_l1norm < 1.495585217958292e-002) {
00313     pade3(m_M);
00314   } else if (m_l1norm < 2.539398330063230e-001) {
00315     pade5(m_M);
00316   } else if (m_l1norm < 9.504178996162932e-001) {
00317     pade7(m_M);
00318   } else if (m_l1norm < 2.097847961257068e+000) {
00319     pade9(m_M);
00320   } else {
00321     const double maxnorm = 5.371920351148152;
00322     m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
00323     MatrixType A = m_M / pow(Scalar(2), m_squarings);
00324     pade13(A);
00325   }
00326 }
00327 
00328 template <typename MatrixType>
00329 void MatrixExponential<MatrixType>::computeUV(long double)
00330 {
00331   using std::max;
00332   using std::pow;
00333   using std::ceil;
00334 #if   LDBL_MANT_DIG == 53   // double precision
00335   computeUV(double());
00336 #elif LDBL_MANT_DIG <= 64   // extended precision
00337   if (m_l1norm < 4.1968497232266989671e-003L) {
00338     pade3(m_M);
00339   } else if (m_l1norm < 1.1848116734693823091e-001L) {
00340     pade5(m_M);
00341   } else if (m_l1norm < 5.5170388480686700274e-001L) {
00342     pade7(m_M);
00343   } else if (m_l1norm < 1.3759868875587845383e+000L) {
00344     pade9(m_M);
00345   } else {
00346     const long double maxnorm = 4.0246098906697353063L;
00347     m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
00348     MatrixType A = m_M / pow(Scalar(2), m_squarings);
00349     pade13(A);
00350   }
00351 #elif LDBL_MANT_DIG <= 106  // double-double
00352   if (m_l1norm < 3.2787892205607026992947488108213e-005L) {
00353     pade3(m_M);
00354   } else if (m_l1norm < 6.4467025060072760084130906076332e-003L) {
00355     pade5(m_M);
00356   } else if (m_l1norm < 6.8988028496595374751374122881143e-002L) {
00357     pade7(m_M);
00358   } else if (m_l1norm < 2.7339737518502231741495857201670e-001L) {
00359     pade9(m_M);
00360   } else if (m_l1norm < 1.3203382096514474905666448850278e+000L) {
00361     pade13(m_M);
00362   } else {
00363     const long double maxnorm = 3.2579440895405400856599663723517L;
00364     m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
00365     MatrixType A = m_M / pow(Scalar(2), m_squarings);
00366     pade17(A);
00367   }
00368 #elif LDBL_MANT_DIG <= 112  // quadruple precison
00369   if (m_l1norm < 1.639394610288918690547467954466970e-005L) {
00370     pade3(m_M);
00371   } else if (m_l1norm < 4.253237712165275566025884344433009e-003L) {
00372     pade5(m_M);
00373   } else if (m_l1norm < 5.125804063165764409885122032933142e-002L) {
00374     pade7(m_M);
00375   } else if (m_l1norm < 2.170000765161155195453205651889853e-001L) {
00376     pade9(m_M);
00377   } else if (m_l1norm < 1.125358383453143065081397882891878e+000L) {
00378     pade13(m_M);
00379   } else {
00380     const long double maxnorm = 2.884233277829519311757165057717815L;
00381     m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
00382     MatrixType A = m_M / pow(Scalar(2), m_squarings);
00383     pade17(A);
00384   }
00385 #else
00386   // this case should be handled in compute()
00387   eigen_assert(false && "Bug in MatrixExponential"); 
00388 #endif  // LDBL_MANT_DIG
00389 }
00390 
00403 template<typename Derived> struct MatrixExponentialReturnValue
00404 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
00405 {
00406     typedef typename Derived::Index Index;
00407   public:
00413     MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
00414 
00420     template <typename ResultType>
00421     inline void evalTo(ResultType& result) const
00422     {
00423       const typename Derived::PlainObject srcEvaluated = m_src.eval();
00424       MatrixExponential<typename Derived::PlainObject> me(srcEvaluated);
00425       me.compute(result);
00426     }
00427 
00428     Index rows() const { return m_src.rows(); }
00429     Index cols() const { return m_src.cols(); }
00430 
00431   protected:
00432     const Derived& m_src;
00433   private:
00434     MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&);
00435 };
00436 
00437 namespace internal {
00438 template<typename Derived>
00439 struct traits<MatrixExponentialReturnValue<Derived> >
00440 {
00441   typedef typename Derived::PlainObject ReturnType;
00442 };
00443 }
00444 
00445 template <typename Derived>
00446 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
00447 {
00448   eigen_assert(rows() == cols());
00449   return MatrixExponentialReturnValue<Derived>(derived());
00450 }
00451 
00452 } // end namespace Eigen
00453 
00454 #endif // EIGEN_MATRIX_EXPONENTIAL


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:18