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 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_MATRIX_EXPONENTIAL
00026 #define EIGEN_MATRIX_EXPONENTIAL
00027 
00028 #ifdef _MSC_VER
00029   template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
00030 #endif
00031 
00032 
00038 template <typename MatrixType>
00039 class MatrixExponential {
00040 
00041   public:
00042 
00050     MatrixExponential(const MatrixType &M);
00051 
00056     template <typename ResultType> 
00057     void compute(ResultType &result);
00058 
00059   private:
00060 
00061     // Prevent copying
00062     MatrixExponential(const MatrixExponential&);
00063     MatrixExponential& operator=(const MatrixExponential&);
00064 
00072     void pade3(const MatrixType &A);
00073 
00081     void pade5(const MatrixType &A);
00082 
00090     void pade7(const MatrixType &A);
00091 
00099     void pade9(const MatrixType &A);
00100 
00108     void pade13(const MatrixType &A);
00109 
00123     void computeUV(double);
00124 
00129     void computeUV(float);
00130 
00131     typedef typename internal::traits<MatrixType>::Scalar Scalar;
00132     typedef typename NumTraits<Scalar>::Real RealScalar;
00133 
00135     typename internal::nested<MatrixType>::type m_M;
00136 
00138     MatrixType m_U;
00139 
00141     MatrixType m_V;
00142 
00144     MatrixType m_tmp1;
00145 
00147     MatrixType m_tmp2;
00148 
00150     MatrixType m_Id;
00151 
00153     int m_squarings;
00154 
00156     float m_l1norm;
00157 };
00158 
00159 template <typename MatrixType>
00160 MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) :
00161   m_M(M),
00162   m_U(M.rows(),M.cols()),
00163   m_V(M.rows(),M.cols()),
00164   m_tmp1(M.rows(),M.cols()),
00165   m_tmp2(M.rows(),M.cols()),
00166   m_Id(MatrixType::Identity(M.rows(), M.cols())),
00167   m_squarings(0),
00168   m_l1norm(static_cast<float>(M.cwiseAbs().colwise().sum().maxCoeff()))
00169 {
00170   /* empty body */
00171 }
00172 
00173 template <typename MatrixType>
00174 template <typename ResultType> 
00175 void MatrixExponential<MatrixType>::compute(ResultType &result)
00176 {
00177   computeUV(RealScalar());
00178   m_tmp1 = m_U + m_V;   // numerator of Pade approximant
00179   m_tmp2 = -m_U + m_V;  // denominator of Pade approximant
00180   result = m_tmp2.partialPivLu().solve(m_tmp1);
00181   for (int i=0; i<m_squarings; i++)
00182     result *= result;           // undo scaling by repeated squaring
00183 }
00184 
00185 template <typename MatrixType>
00186 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A)
00187 {
00188   const Scalar b[] = {120., 60., 12., 1.};
00189   m_tmp1.noalias() = A * A;
00190   m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id;
00191   m_U.noalias() = A * m_tmp2;
00192   m_V = b[2]*m_tmp1 + b[0]*m_Id;
00193 }
00194 
00195 template <typename MatrixType>
00196 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A)
00197 {
00198   const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
00199   MatrixType A2 = A * A;
00200   m_tmp1.noalias() = A2 * A2;
00201   m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id;
00202   m_U.noalias() = A * m_tmp2;
00203   m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id;
00204 }
00205 
00206 template <typename MatrixType>
00207 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A)
00208 {
00209   const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
00210   MatrixType A2 = A * A;
00211   MatrixType A4 = A2 * A2;
00212   m_tmp1.noalias() = A4 * A2;
00213   m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00214   m_U.noalias() = A * m_tmp2;
00215   m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00216 }
00217 
00218 template <typename MatrixType>
00219 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A)
00220 {
00221   const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
00222                       2162160., 110880., 3960., 90., 1.};
00223   MatrixType A2 = A * A;
00224   MatrixType A4 = A2 * A2;
00225   MatrixType A6 = A4 * A2;
00226   m_tmp1.noalias() = A6 * A2;
00227   m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00228   m_U.noalias() = A * m_tmp2;
00229   m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00230 }
00231 
00232 template <typename MatrixType>
00233 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A)
00234 {
00235   const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
00236                       1187353796428800., 129060195264000., 10559470521600., 670442572800.,
00237                       33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
00238   MatrixType A2 = A * A;
00239   MatrixType A4 = A2 * A2;
00240   m_tmp1.noalias() = A4 * A2;
00241   m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage
00242   m_tmp2.noalias() = m_tmp1 * m_V;
00243   m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
00244   m_U.noalias() = A * m_tmp2;
00245   m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2;
00246   m_V.noalias() = m_tmp1 * m_tmp2;
00247   m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
00248 }
00249 
00250 template <typename MatrixType>
00251 void MatrixExponential<MatrixType>::computeUV(float)
00252 {
00253   using std::max;
00254   using std::pow;
00255   using std::ceil;
00256   if (m_l1norm < 4.258730016922831e-001) {
00257     pade3(m_M);
00258   } else if (m_l1norm < 1.880152677804762e+000) {
00259     pade5(m_M);
00260   } else {
00261     const float maxnorm = 3.925724783138660f;
00262     m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
00263     MatrixType A = m_M / pow(Scalar(2), Scalar(static_cast<RealScalar>(m_squarings)));
00264     pade7(A);
00265   }
00266 }
00267 
00268 template <typename MatrixType>
00269 void MatrixExponential<MatrixType>::computeUV(double)
00270 {
00271   using std::max;
00272   using std::pow;
00273   using std::ceil;
00274   if (m_l1norm < 1.495585217958292e-002) {
00275     pade3(m_M);
00276   } else if (m_l1norm < 2.539398330063230e-001) {
00277     pade5(m_M);
00278   } else if (m_l1norm < 9.504178996162932e-001) {
00279     pade7(m_M);
00280   } else if (m_l1norm < 2.097847961257068e+000) {
00281     pade9(m_M);
00282   } else {
00283     const double maxnorm = 5.371920351148152;
00284     m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
00285     MatrixType A = m_M / pow(Scalar(2), Scalar(m_squarings));
00286     pade13(A);
00287   }
00288 }
00289 
00302 template<typename Derived> struct MatrixExponentialReturnValue
00303 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
00304 {
00305     typedef typename Derived::Index Index;
00306   public:
00312     MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
00313 
00319     template <typename ResultType>
00320     inline void evalTo(ResultType& result) const
00321     {
00322       const typename Derived::PlainObject srcEvaluated = m_src.eval();
00323       MatrixExponential<typename Derived::PlainObject> me(srcEvaluated);
00324       me.compute(result);
00325     }
00326 
00327     Index rows() const { return m_src.rows(); }
00328     Index cols() const { return m_src.cols(); }
00329 
00330   protected:
00331     const Derived& m_src;
00332   private:
00333     MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&);
00334 };
00335 
00336 namespace internal {
00337 template<typename Derived>
00338 struct traits<MatrixExponentialReturnValue<Derived> >
00339 {
00340   typedef typename Derived::PlainObject ReturnType;
00341 };
00342 }
00343 
00344 template <typename Derived>
00345 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
00346 {
00347   eigen_assert(rows() == cols());
00348   return MatrixExponentialReturnValue<Derived>(derived());
00349 }
00350 
00351 #endif // EIGEN_MATRIX_EXPONENTIAL


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:06