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

1 files changed, 25 insertions, 18 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 9ad2b9cc8..bb6d9e1fe 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -61,10 +61,11 @@ struct MatrixExponentialScalingOp
  *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé
  *  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