diff options
| -rw-r--r-- | unsupported/Eigen/MatrixFunctions | 57 | ||||
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h | 2 | ||||
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 22 | ||||
| -rw-r--r-- | unsupported/doc/examples/MatrixPower.cpp | 16 |
4 files changed, 85 insertions, 12 deletions
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 9917efeed..5133afb83 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -216,6 +216,63 @@ Output: \verbinclude MatrixLogarithm.out class MatrixLogarithmAtomic, MatrixBase::sqrt(). +\section matrixbase_pow MatrixBase::pow() + +Compute the matrix raised to arbitrary real power. + +\code +template <typename ExponentType> +const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const +\endcode + +\param[in] M base of the matrix power, should be a square matrix. +\param[in] p exponent of the matrix power, should be an integer or +the same type as the real scalar in \p M. + +The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$, +where exp denotes the matrix exponential, and log denotes the matrix +logarithm. + +The matrix \f$ M \f$ should meet the conditions to be an argument of +matrix logarithm. + +This function computes the matrix logarithm using the +Schur-Padé algorithm as implemented by MatrixBase::pow(). +The exponent is split into integral part and fractional part, where +the fractional part is in the interval \f$ (-1, 1) \f$. The main +diagonal and the first super-diagonal is directly computed. + +The actual work is done by the MatrixPower class, which can compute +\f$ M^p v \f$, where \p v is another matrix with the same rows as +\p M. The matrix \p v is set to be the identity matrix by default. + +Details of the algorithm can be found in: Nicholas J. Higham and +Lijing Lin, "A Schur-Padé algorithm for fractional powers of a +matrix," <em>SIAM J. %Matrix Anal. Applic.</em>, +<b>32(3)</b>:1056–1078, 2011. + +Example: The following program checks that +\f[ \left[ \begin{array}{ccc} + \cos1 & -\sin1 & 0 \\ + \sin1 & \cos1 & 0 \\ + 0 & 0 & 1 + \end{array} \right]^{\frac14\pi} = \left[ \begin{array}{ccc} + \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ + \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ + 0 & 0 & 1 + \end{array} \right]. \f] +This corresponds to \f$ \frac14\pi \f$ rotations of 1 radian around +the z-axis. + +\include MatrixPower.cpp +Output: \verbinclude MatrixPower.out + +\note \p M has to be a matrix of \c float, \c double, \c long double +\c complex<float>, \c complex<double>, or \c complex<long double> . + +\sa MatrixBase::exp(), MatrixBase::log(), class MatrixPower. + + \section matrixbase_matrixfunction MatrixBase::matrixFunction() Compute a matrix function. diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index f237b9cdc..781d7bf93 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -217,7 +217,7 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(long double normTminusI) 3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L, 8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L }; #endif - int degree = 3 + int degree = 3; for (; degree <= maxPadeDegree; ++degree) if (normTminusI <= maxNormForPade[degree - minPadeDegree]) break; diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 918db63b4..86ef24eac 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -71,8 +71,8 @@ class MatrixPower /** * \brief Compute the matrix power. * - * If \c b is \em fatter than \c A, it computes \f$ A^{p_{\textrm int}} - * \f$ first, and then multiplies it with \c b. Otherwise, + * If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}} + * \f$ first, and then multiplies it with \p b. Otherwise, * #computeChainProduct optimizes the expression. * * \sa computeChainProduct(ResultType&); @@ -124,13 +124,13 @@ class MatrixPower */ void computeBig(); - /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c double) */ + /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */ inline int getPadeDegree(double); - /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c float) */ + /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = float) */ inline int getPadeDegree(float); - /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c long double) */ + /** \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */ inline int getPadeDegree(long double); /** \brief Compute Padé approximation to matrix fractional power. */ @@ -196,8 +196,8 @@ class MatrixPower<MatrixType, IntExponent, PlainObject, 1> /** * \brief Compute the matrix power. * - * If \c b is \em fatter than \c A, it computes \f$ A^p \f$ first, and - * then multiplies it with \c b. Otherwise, #computeChainProduct + * If \p b is \em fatter than \p A, it computes \f$ A^p \f$ first, and + * then multiplies it with \p b. Otherwise, #computeChainProduct * optimizes the expression. * * \param[out] result \f$ A^p b \f$, as specified in the constructor. @@ -646,7 +646,7 @@ template<typename MatrixType, typename ExponentType, typename Derived> class Mat /** * \brief Compute the matrix exponential. * - * \param[out] result \f$ A^p b \f$ where \c A ,\c p and \c b are as in + * \param[out] result \f$ A^p b \f$ where \p A ,\p p and \p b are as in * the constructor. */ template <typename ResultType> @@ -700,12 +700,12 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue : m_A(A), m_p(p) { } /** - * \brief Return the matrix power multiplied by %Matrix \c b. + * \brief Return the matrix power multiplied by %Matrix \p b. * * The %MatrixPower class can optimize \f$ A^p b \f$ computing, and this * method provides an elegant way to call it: * - * \param[in] b %Matrix (exporession), the multiplier. + * \param[in] b %Matrix (expression), the multiplier. */ template <typename OtherDerived> const MatrixPowerMultiplied<Derived, ExponentType, OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const @@ -714,7 +714,7 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue /** * \brief Compute the matrix power. * - * \param[out] result \f$ A^p \f$ where \c A and \c p are as in the + * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the * constructor. */ template <typename ResultType> diff --git a/unsupported/doc/examples/MatrixPower.cpp b/unsupported/doc/examples/MatrixPower.cpp new file mode 100644 index 000000000..6ade0b8af --- /dev/null +++ b/unsupported/doc/examples/MatrixPower.cpp @@ -0,0 +1,16 @@ +#include <unsupported/Eigen/MatrixFunctions> +#include <iostream> + +using namespace Eigen; + +int main() +{ + const double pi = std::acos(-1.0); + Matrix3d A; + A << cos(1), -sin(1), 0, + sin(1), cos(1), 0, + 0 , 0 , 1; + std::cout << "The matrix A is:\n" << A << "\n\n" + << "The matrix power A^(pi/4) is:\n" << A.pow(pi/4) << std::endl; + return 0; +} |