Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:57:04 

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
0005 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010 
0011 #ifndef EIGEN_MATRIX_EXPONENTIAL
0012 #define EIGEN_MATRIX_EXPONENTIAL
0013 
0014 #include "StemFunction.h"
0015 
0016 namespace Eigen {
0017 namespace internal {
0018 
0019 /** \brief Scaling operator.
0020  *
0021  * This struct is used by CwiseUnaryOp to scale a matrix by \f$ 2^{-s} \f$.
0022  */
0023 template <typename RealScalar>
0024 struct MatrixExponentialScalingOp
0025 {
0026   /** \brief Constructor.
0027    *
0028    * \param[in] squarings  The integer \f$ s \f$ in this document.
0029    */
0030   MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
0031 
0032 
0033   /** \brief Scale a matrix coefficient.
0034    *
0035    * \param[in,out] x  The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
0036    */
0037   inline const RealScalar operator() (const RealScalar& x) const
0038   {
0039     using std::ldexp;
0040     return ldexp(x, -m_squarings);
0041   }
0042 
0043   typedef std::complex<RealScalar> ComplexScalar;
0044 
0045   /** \brief Scale a matrix coefficient.
0046    *
0047    * \param[in,out] x  The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
0048    */
0049   inline const ComplexScalar operator() (const ComplexScalar& x) const
0050   {
0051     using std::ldexp;
0052     return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
0053   }
0054 
0055   private:
0056     int m_squarings;
0057 };
0058 
0059 /** \brief Compute the (3,3)-Pad&eacute; approximant to the exponential.
0060  *
0061  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0062  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0063  */
0064 template <typename MatA, typename MatU, typename MatV>
0065 void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V)
0066 {
0067   typedef typename MatA::PlainObject MatrixType;
0068   typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
0069   const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
0070   const MatrixType A2 = A * A;
0071   const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
0072   U.noalias() = A * tmp;
0073   V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
0074 }
0075 
0076 /** \brief Compute the (5,5)-Pad&eacute; approximant to the exponential.
0077  *
0078  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0079  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0080  */
0081 template <typename MatA, typename MatU, typename MatV>
0082 void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V)
0083 {
0084   typedef typename MatA::PlainObject MatrixType;
0085   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0086   const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
0087   const MatrixType A2 = A * A;
0088   const MatrixType A4 = A2 * A2;
0089   const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
0090   U.noalias() = A * tmp;
0091   V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
0092 }
0093 
0094 /** \brief Compute the (7,7)-Pad&eacute; approximant to the exponential.
0095  *
0096  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0097  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0098  */
0099 template <typename MatA, typename MatU, typename MatV>
0100 void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V)
0101 {
0102   typedef typename MatA::PlainObject MatrixType;
0103   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0104   const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
0105   const MatrixType A2 = A * A;
0106   const MatrixType A4 = A2 * A2;
0107   const MatrixType A6 = A4 * A2;
0108   const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 
0109     + b[1] * MatrixType::Identity(A.rows(), A.cols());
0110   U.noalias() = A * tmp;
0111   V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
0112 
0113 }
0114 
0115 /** \brief Compute the (9,9)-Pad&eacute; approximant to the exponential.
0116  *
0117  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0118  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0119  */
0120 template <typename MatA, typename MatU, typename MatV>
0121 void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V)
0122 {
0123   typedef typename MatA::PlainObject MatrixType;
0124   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0125   const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
0126                           2162160.L, 110880.L, 3960.L, 90.L, 1.L};
0127   const MatrixType A2 = A * A;
0128   const MatrixType A4 = A2 * A2;
0129   const MatrixType A6 = A4 * A2;
0130   const MatrixType A8 = A6 * A2;
0131   const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 
0132     + b[1] * MatrixType::Identity(A.rows(), A.cols());
0133   U.noalias() = A * tmp;
0134   V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
0135 }
0136 
0137 /** \brief Compute the (13,13)-Pad&eacute; approximant to the exponential.
0138  *
0139  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0140  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0141  */
0142 template <typename MatA, typename MatU, typename MatV>
0143 void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V)
0144 {
0145   typedef typename MatA::PlainObject MatrixType;
0146   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0147   const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
0148                           1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
0149                           33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
0150   const MatrixType A2 = A * A;
0151   const MatrixType A4 = A2 * A2;
0152   const MatrixType A6 = A4 * A2;
0153   V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
0154   MatrixType tmp = A6 * V;
0155   tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
0156   U.noalias() = A * tmp;
0157   tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
0158   V.noalias() = A6 * tmp;
0159   V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
0160 }
0161 
0162 /** \brief Compute the (17,17)-Pad&eacute; approximant to the exponential.
0163  *
0164  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
0165  *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
0166  *
0167  *  This function activates only if your long double is double-double or quadruple.
0168  */
0169 #if LDBL_MANT_DIG > 64
0170 template <typename MatA, typename MatU, typename MatV>
0171 void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V)
0172 {
0173   typedef typename MatA::PlainObject MatrixType;
0174   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0175   const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
0176                           100610229646136770560000.L, 15720348382208870400000.L,
0177                           1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
0178                           595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
0179                           33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
0180                           46512.L, 306.L, 1.L};
0181   const MatrixType A2 = A * A;
0182   const MatrixType A4 = A2 * A2;
0183   const MatrixType A6 = A4 * A2;
0184   const MatrixType A8 = A4 * A4;
0185   V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
0186   MatrixType tmp = A8 * V;
0187   tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 
0188     + b[1] * MatrixType::Identity(A.rows(), A.cols());
0189   U.noalias() = A * tmp;
0190   tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
0191   V.noalias() = tmp * A8;
0192   V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 
0193     + b[0] * MatrixType::Identity(A.rows(), A.cols());
0194 }
0195 #endif
0196 
0197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
0198 struct matrix_exp_computeUV
0199 {
0200   /** \brief Compute Pad&eacute; approximant to the exponential.
0201     *
0202     * Computes \c U, \c V and \c squarings such that \f$ (V+U)(V-U)^{-1} \f$ is a Pad&eacute;
0203     * approximant of \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$, where \f$ M \f$
0204     * denotes the matrix \c arg. The degree of the Pad&eacute; approximant and the value of squarings
0205     * are chosen such that the approximation error is no more than the round-off error.
0206     */
0207   static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
0208 };
0209 
0210 template <typename MatrixType>
0211 struct matrix_exp_computeUV<MatrixType, float>
0212 {
0213   template <typename ArgType>
0214   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
0215   {
0216     using std::frexp;
0217     using std::pow;
0218     const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
0219     squarings = 0;
0220     if (l1norm < 4.258730016922831e-001f) {
0221       matrix_exp_pade3(arg, U, V);
0222     } else if (l1norm < 1.880152677804762e+000f) {
0223       matrix_exp_pade5(arg, U, V);
0224     } else {
0225       const float maxnorm = 3.925724783138660f;
0226       frexp(l1norm / maxnorm, &squarings);
0227       if (squarings < 0) squarings = 0;
0228       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
0229       matrix_exp_pade7(A, U, V);
0230     }
0231   }
0232 };
0233 
0234 template <typename MatrixType>
0235 struct matrix_exp_computeUV<MatrixType, double>
0236 {
0237   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
0238   template <typename ArgType>
0239   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
0240   {
0241     using std::frexp;
0242     using std::pow;
0243     const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
0244     squarings = 0;
0245     if (l1norm < 1.495585217958292e-002) {
0246       matrix_exp_pade3(arg, U, V);
0247     } else if (l1norm < 2.539398330063230e-001) {
0248       matrix_exp_pade5(arg, U, V);
0249     } else if (l1norm < 9.504178996162932e-001) {
0250       matrix_exp_pade7(arg, U, V);
0251     } else if (l1norm < 2.097847961257068e+000) {
0252       matrix_exp_pade9(arg, U, V);
0253     } else {
0254       const RealScalar maxnorm = 5.371920351148152;
0255       frexp(l1norm / maxnorm, &squarings);
0256       if (squarings < 0) squarings = 0;
0257       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<RealScalar>(squarings));
0258       matrix_exp_pade13(A, U, V);
0259     }
0260   }
0261 };
0262   
0263 template <typename MatrixType>
0264 struct matrix_exp_computeUV<MatrixType, long double>
0265 {
0266   template <typename ArgType>
0267   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
0268   {
0269 #if   LDBL_MANT_DIG == 53   // double precision
0270     matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
0271   
0272 #else
0273   
0274     using std::frexp;
0275     using std::pow;
0276     const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
0277     squarings = 0;
0278   
0279 #if LDBL_MANT_DIG <= 64   // extended precision
0280   
0281     if (l1norm < 4.1968497232266989671e-003L) {
0282       matrix_exp_pade3(arg, U, V);
0283     } else if (l1norm < 1.1848116734693823091e-001L) {
0284       matrix_exp_pade5(arg, U, V);
0285     } else if (l1norm < 5.5170388480686700274e-001L) {
0286       matrix_exp_pade7(arg, U, V);
0287     } else if (l1norm < 1.3759868875587845383e+000L) {
0288       matrix_exp_pade9(arg, U, V);
0289     } else {
0290       const long double maxnorm = 4.0246098906697353063L;
0291       frexp(l1norm / maxnorm, &squarings);
0292       if (squarings < 0) squarings = 0;
0293       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
0294       matrix_exp_pade13(A, U, V);
0295     }
0296   
0297 #elif LDBL_MANT_DIG <= 106  // double-double
0298   
0299     if (l1norm < 3.2787892205607026992947488108213e-005L) {
0300       matrix_exp_pade3(arg, U, V);
0301     } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
0302       matrix_exp_pade5(arg, U, V);
0303     } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
0304       matrix_exp_pade7(arg, U, V);
0305     } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
0306       matrix_exp_pade9(arg, U, V);
0307     } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
0308       matrix_exp_pade13(arg, U, V);
0309     } else {
0310       const long double maxnorm = 3.2579440895405400856599663723517L;
0311       frexp(l1norm / maxnorm, &squarings);
0312       if (squarings < 0) squarings = 0;
0313       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
0314       matrix_exp_pade17(A, U, V);
0315     }
0316   
0317 #elif LDBL_MANT_DIG <= 113  // quadruple precision
0318   
0319     if (l1norm < 1.639394610288918690547467954466970e-005L) {
0320       matrix_exp_pade3(arg, U, V);
0321     } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
0322       matrix_exp_pade5(arg, U, V);
0323     } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
0324       matrix_exp_pade7(arg, U, V);
0325     } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
0326       matrix_exp_pade9(arg, U, V);
0327     } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
0328       matrix_exp_pade13(arg, U, V);
0329     } else {
0330       const long double maxnorm = 2.884233277829519311757165057717815L;
0331       frexp(l1norm / maxnorm, &squarings);
0332       if (squarings < 0) squarings = 0;
0333       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
0334       matrix_exp_pade17(A, U, V);
0335     }
0336   
0337 #else
0338   
0339     // this case should be handled in compute()
0340     eigen_assert(false && "Bug in MatrixExponential"); 
0341   
0342 #endif
0343 #endif  // LDBL_MANT_DIG
0344   }
0345 };
0346 
0347 template<typename T> struct is_exp_known_type : false_type {};
0348 template<> struct is_exp_known_type<float> : true_type {};
0349 template<> struct is_exp_known_type<double> : true_type {};
0350 #if LDBL_MANT_DIG <= 113
0351 template<> struct is_exp_known_type<long double> : true_type {};
0352 #endif
0353 
0354 template <typename ArgType, typename ResultType>
0355 void matrix_exp_compute(const ArgType& arg, ResultType &result, true_type) // natively supported scalar type
0356 {
0357   typedef typename ArgType::PlainObject MatrixType;
0358   MatrixType U, V;
0359   int squarings;
0360   matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
0361   MatrixType numer = U + V;
0362   MatrixType denom = -U + V;
0363   result = denom.partialPivLu().solve(numer);
0364   for (int i=0; i<squarings; i++)
0365     result *= result;   // undo scaling by repeated squaring
0366 }
0367 
0368 
0369 /* Computes the matrix exponential
0370  *
0371  * \param arg    argument of matrix exponential (should be plain object)
0372  * \param result variable in which result will be stored
0373  */
0374 template <typename ArgType, typename ResultType>
0375 void matrix_exp_compute(const ArgType& arg, ResultType &result, false_type) // default
0376 {
0377   typedef typename ArgType::PlainObject MatrixType;
0378   typedef typename traits<MatrixType>::Scalar Scalar;
0379   typedef typename NumTraits<Scalar>::Real RealScalar;
0380   typedef typename std::complex<RealScalar> ComplexScalar;
0381   result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
0382 }
0383 
0384 } // end namespace Eigen::internal
0385 
0386 /** \ingroup MatrixFunctions_Module
0387   *
0388   * \brief Proxy for the matrix exponential of some matrix (expression).
0389   *
0390   * \tparam Derived  Type of the argument to the matrix exponential.
0391   *
0392   * This class holds the argument to the matrix exponential until it is assigned or evaluated for
0393   * some other reason (so the argument should not be changed in the meantime). It is the return type
0394   * of MatrixBase::exp() and most of the time this is the only way it is used.
0395   */
0396 template<typename Derived> struct MatrixExponentialReturnValue
0397 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
0398 {
0399   public:
0400     /** \brief Constructor.
0401       *
0402       * \param src %Matrix (expression) forming the argument of the matrix exponential.
0403       */
0404     MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
0405 
0406     /** \brief Compute the matrix exponential.
0407       *
0408       * \param result the matrix exponential of \p src in the constructor.
0409       */
0410     template <typename ResultType>
0411     inline void evalTo(ResultType& result) const
0412     {
0413       const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
0414       internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
0415     }
0416 
0417     Index rows() const { return m_src.rows(); }
0418     Index cols() const { return m_src.cols(); }
0419 
0420   protected:
0421     const typename internal::ref_selector<Derived>::type m_src;
0422 };
0423 
0424 namespace internal {
0425 template<typename Derived>
0426 struct traits<MatrixExponentialReturnValue<Derived> >
0427 {
0428   typedef typename Derived::PlainObject ReturnType;
0429 };
0430 }
0431 
0432 template <typename Derived>
0433 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
0434 {
0435   eigen_assert(rows() == cols());
0436   return MatrixExponentialReturnValue<Derived>(derived());
0437 }
0438 
0439 } // end namespace Eigen
0440 
0441 #endif // EIGEN_MATRIX_EXPONENTIAL
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!