|  | // This file is part of Eigen, a lightweight C++ template library | 
|  | // for linear algebra. | 
|  | // | 
|  | // Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk> | 
|  | // | 
|  | // Eigen is free software; you can redistribute it and/or | 
|  | // modify it under the terms of the GNU Lesser General Public | 
|  | // License as published by the Free Software Foundation; either | 
|  | // version 3 of the License, or (at your option) any later version. | 
|  | // | 
|  | // Alternatively, you can redistribute it and/or | 
|  | // modify it under the terms of the GNU General Public License as | 
|  | // published by the Free Software Foundation; either version 2 of | 
|  | // the License, or (at your option) any later version. | 
|  | // | 
|  | // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY | 
|  | // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS | 
|  | // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the | 
|  | // GNU General Public License for more details. | 
|  | // | 
|  | // You should have received a copy of the GNU Lesser General Public | 
|  | // License and a copy of the GNU General Public License along with | 
|  | // Eigen. If not, see <http://www.gnu.org/licenses/>. | 
|  |  | 
|  | #ifndef EIGEN_MATRIX_FUNCTIONS | 
|  | #define EIGEN_MATRIX_FUNCTIONS | 
|  |  | 
|  | #include <list> | 
|  | #include <functional> | 
|  | #include <iterator> | 
|  |  | 
|  | #include <Eigen/Core> | 
|  | #include <Eigen/LU> | 
|  | #include <Eigen/Eigenvalues> | 
|  |  | 
|  | namespace Eigen { | 
|  |  | 
|  | /** \ingroup Unsupported_modules | 
|  | * \defgroup MatrixFunctions_Module Matrix functions module | 
|  | * \brief This module aims to provide various methods for the computation of | 
|  | * matrix functions. | 
|  | * | 
|  | * To use this module, add | 
|  | * \code | 
|  | * #include <unsupported/Eigen/MatrixFunctions> | 
|  | * \endcode | 
|  | * at the start of your source file. | 
|  | * | 
|  | * This module defines the following MatrixBase methods. | 
|  | *  - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine | 
|  | *  - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine | 
|  | *  - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential | 
|  | *  - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions | 
|  | *  - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine | 
|  | *  - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine | 
|  | * | 
|  | * These methods are the main entry points to this module. | 
|  | * | 
|  | * %Matrix functions are defined as follows.  Suppose that \f$ f \f$ | 
|  | * is an entire function (that is, a function on the complex plane | 
|  | * that is everywhere complex differentiable).  Then its Taylor | 
|  | * series | 
|  | * \f[ f(0) + f'(0) x + \frac{f''(0)}{2} x^2 + \frac{f'''(0)}{3!} x^3 + \cdots \f] | 
|  | * converges to \f$ f(x) \f$. In this case, we can define the matrix | 
|  | * function by the same series: | 
|  | * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] | 
|  | * | 
|  | */ | 
|  |  | 
|  | #include "src/MatrixFunctions/MatrixExponential.h" | 
|  | #include "src/MatrixFunctions/MatrixFunction.h" | 
|  |  | 
|  |  | 
|  |  | 
|  | /** | 
|  | \page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module | 
|  | \ingroup MatrixFunctions_Module | 
|  |  | 
|  | The remainder of the page documents the following MatrixBase methods | 
|  | which are defined in the MatrixFunctions module. | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_cos MatrixBase::cos() | 
|  |  | 
|  | Compute the matrix cosine. | 
|  |  | 
|  | \code | 
|  | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  a square matrix. | 
|  | \returns  expression representing \f$ \cos(M) \f$. | 
|  |  | 
|  | This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). | 
|  |  | 
|  | \sa \ref matrixbase_sin "sin()" for an example. | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_cosh MatrixBase::cosh() | 
|  |  | 
|  | Compute the matrix hyberbolic cosine. | 
|  |  | 
|  | \code | 
|  | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  a square matrix. | 
|  | \returns  expression representing \f$ \cosh(M) \f$ | 
|  |  | 
|  | This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). | 
|  |  | 
|  | \sa \ref matrixbase_sinh "sinh()" for an example. | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_exp MatrixBase::exp() | 
|  |  | 
|  | Compute the matrix exponential. | 
|  |  | 
|  | \code | 
|  | const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  matrix whose exponential is to be computed. | 
|  | \returns    expression representing the matrix exponential of \p M. | 
|  |  | 
|  | The matrix exponential of \f$ M \f$ is defined by | 
|  | \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] | 
|  | The matrix exponential can be used to solve linear ordinary | 
|  | differential equations: the solution of \f$ y' = My \f$ with the | 
|  | initial condition \f$ y(0) = y_0 \f$ is given by | 
|  | \f$ y(t) = \exp(M) y_0 \f$. | 
|  |  | 
|  | The cost of the computation is approximately \f$ 20 n^3 \f$ for | 
|  | matrices of size \f$ n \f$. The number 20 depends weakly on the | 
|  | norm of the matrix. | 
|  |  | 
|  | The matrix exponential is computed using the scaling-and-squaring | 
|  | method combined with Padé approximation. The matrix is first | 
|  | rescaled, then the exponential of the reduced matrix is computed | 
|  | approximant, and then the rescaling is undone by repeated | 
|  | squaring. The degree of the Padé approximant is chosen such | 
|  | that the approximation error is less than the round-off | 
|  | error. However, errors may accumulate during the squaring phase. | 
|  |  | 
|  | Details of the algorithm can be found in: Nicholas J. Higham, "The | 
|  | scaling and squaring method for the matrix exponential revisited," | 
|  | <em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193, | 
|  | 2005. | 
|  |  | 
|  | Example: The following program checks that | 
|  | \f[ \exp \left[ \begin{array}{ccc} | 
|  | 0 & \frac14\pi & 0 \\ | 
|  | -\frac14\pi & 0 & 0 \\ | 
|  | 0 & 0 & 0 | 
|  | \end{array} \right] = \left[ \begin{array}{ccc} | 
|  | \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ | 
|  | \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ | 
|  | 0 & 0 & 1 | 
|  | \end{array} \right]. \f] | 
|  | This corresponds to a rotation of \f$ \frac14\pi \f$ radians around | 
|  | the z-axis. | 
|  |  | 
|  | \include MatrixExponential.cpp | 
|  | Output: \verbinclude MatrixExponential.out | 
|  |  | 
|  | \note \p M has to be a matrix of \c float, \c double, | 
|  | \c complex<float> or \c complex<double> . | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_matrixfunction MatrixBase::matrixFunction() | 
|  |  | 
|  | Compute a matrix function. | 
|  |  | 
|  | \code | 
|  | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  argument of matrix function, should be a square matrix. | 
|  | \param[in]  f  an entire function; \c f(x,n) should compute the n-th | 
|  | derivative of f at x. | 
|  | \returns  expression representing \p f applied to \p M. | 
|  |  | 
|  | Suppose that \p M is a matrix whose entries have type \c Scalar. | 
|  | Then, the second argument, \p f, should be a function with prototype | 
|  | \code | 
|  | ComplexScalar f(ComplexScalar, int) | 
|  | \endcode | 
|  | where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is | 
|  | real (e.g., \c float or \c double) and \c ComplexScalar = | 
|  | \c Scalar if \c Scalar is complex. The return value of \c f(x,n) | 
|  | should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. | 
|  |  | 
|  | This routine uses the algorithm described in: | 
|  | Philip Davies and Nicholas J. Higham, | 
|  | "A Schur-Parlett algorithm for computing matrix functions", | 
|  | <em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003. | 
|  |  | 
|  | The actual work is done by the MatrixFunction class. | 
|  |  | 
|  | Example: The following program checks that | 
|  | \f[ \exp \left[ \begin{array}{ccc} | 
|  | 0 & \frac14\pi & 0 \\ | 
|  | -\frac14\pi & 0 & 0 \\ | 
|  | 0 & 0 & 0 | 
|  | \end{array} \right] = \left[ \begin{array}{ccc} | 
|  | \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ | 
|  | \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ | 
|  | 0 & 0 & 1 | 
|  | \end{array} \right]. \f] | 
|  | This corresponds to a rotation of \f$ \frac14\pi \f$ radians around | 
|  | the z-axis. This is the same example as used in the documentation | 
|  | of \ref matrixbase_exp "exp()". | 
|  |  | 
|  | \include MatrixFunction.cpp | 
|  | Output: \verbinclude MatrixFunction.out | 
|  |  | 
|  | Note that the function \c expfn is defined for complex numbers | 
|  | \c x, even though the matrix \c A is over the reals. Instead of | 
|  | \c expfn, we could also have used StdStemFunctions::exp: | 
|  | \code | 
|  | A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); | 
|  | \endcode | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_sin MatrixBase::sin() | 
|  |  | 
|  | Compute the matrix sine. | 
|  |  | 
|  | \code | 
|  | const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  a square matrix. | 
|  | \returns  expression representing \f$ \sin(M) \f$. | 
|  |  | 
|  | This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). | 
|  |  | 
|  | Example: \include MatrixSine.cpp | 
|  | Output: \verbinclude MatrixSine.out | 
|  |  | 
|  |  | 
|  |  | 
|  | \section matrixbase_sinh const MatrixBase::sinh() | 
|  |  | 
|  | Compute the matrix hyperbolic sine. | 
|  |  | 
|  | \code | 
|  | MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const | 
|  | \endcode | 
|  |  | 
|  | \param[in]  M  a square matrix. | 
|  | \returns  expression representing \f$ \sinh(M) \f$ | 
|  |  | 
|  | This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). | 
|  |  | 
|  | Example: \include MatrixSinh.cpp | 
|  | Output: \verbinclude MatrixSinh.out | 
|  |  | 
|  | */ | 
|  |  | 
|  | } | 
|  |  | 
|  | #endif // EIGEN_MATRIX_FUNCTIONS | 
|  |  |