CFD/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

bug #1380: for Map<> as input of matrix exponential

Browse Source 
(grafted from d8b1f6cebd
)
This commit is contained in:
Gael Guennebaud 2017-02-20 14:06:06 +01:00
parent e1385337ff
commit 17bbd82f7d
1 changed files with 25 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
View File

@ -61,10 +61,11 @@ struct MatrixExponentialScalingOp
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*/
template <typename MatrixType>
void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V)
{
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
const MatrixType A2 = A * A;
const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
@ -77,9 +78,10 @@ void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V)
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*/
template <typename MatrixType>
void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V)
{
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
const MatrixType A2 = A * A;
@ -94,9 +96,10 @@ void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V)
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*/
template <typename MatrixType>
void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V)
{
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
const MatrixType A2 = A * A;
@ -114,9 +117,10 @@ void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V)
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*/
template <typename MatrixType>
void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V)
{
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
2162160.L, 110880.L, 3960.L, 90.L, 1.L};
@ -135,9 +139,10 @@ void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V)
* After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
* approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
*/
template <typename MatrixType>
void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V)
{
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
@ -162,9 +167,10 @@ void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V)
* This function activates only if your long double is double-double or quadruple.
*/
#if LDBL_MANT_DIG > 64
template <typename MatrixType>
void matrix_exp_pade17(const MatrixType &A, MatrixType &U, MatrixType &V)
template <typename MatA, typename MatU, typename MatV>
void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V)
{
typedef typename MatA::PlainObject MatrixType;
typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
100610229646136770560000.L, 15720348382208870400000.L,
@ -342,9 +348,10 @@ struct matrix_exp_computeUV<MatrixType, long double>
* \param arg argument of matrix exponential (should be plain object)
* \param result variable in which result will be stored
*/
template <typename MatrixType, typename ResultType>
void matrix_exp_compute(const MatrixType& arg, ResultType &result)
template <typename ArgType, typename ResultType>
void matrix_exp_compute(const ArgType& arg, ResultType &result)
{
typedef typename ArgType::PlainObject MatrixType;
#if LDBL_MANT_DIG > 112 // rarely happens
typedef typename traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
@ -354,11 +361,11 @@ void matrix_exp_compute(const MatrixType& arg, ResultType &result)
return;
}
#endif
typename MatrixType::PlainObject U, V;
MatrixType U, V;
int squarings;
matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
typename MatrixType::PlainObject numer = U + V;
typename MatrixType::PlainObject denom = -U + V;
MatrixType numer = U + V;
MatrixType denom = -U + V;
result = denom.partialPivLu().solve(numer);
for (int i=0; i<squarings; i++)
result *= result; // undo scaling by repeated squaring