aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
authorGravatar Chen-Pang He <jdh8@ms63.hinet.net>2011-08-24 18:26:38 +0800
committerGravatar Chen-Pang He <jdh8@ms63.hinet.net>2011-08-24 18:26:38 +0800
commit8ddd1e390b43b9d86c02fbdacdf7b70322eecfff (patch)
tree67866f37004cb9c2db1dd0e2f4ca5d4fe7a0559a
parent8414be739b78b2d750be03d41cb5017325f4ae49 (diff)
fix: <ctime> is necessary for srand(time(NULL))
-rw-r--r--unsupported/Eigen/MPRealSupport3
-rw-r--r--unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h29
2 files changed, 19 insertions, 13 deletions
diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport
index 30e8d66e4..8f6132946 100644
--- a/unsupported/Eigen/MPRealSupport
+++ b/unsupported/Eigen/MPRealSupport
@@ -30,6 +30,7 @@
#ifndef EIGEN_MPREALSUPPORT_MODULE_H
#define EIGEN_MPREALSUPPORT_MODULE_H
+#include <ctime>
#include <mpreal.h>
#include <Eigen/Core>
@@ -51,7 +52,7 @@ namespace Eigen {
*
\code
#include <iostream>
-#include <Eigen/Mpfrc++Support>
+#include <Eigen/MPRealSupport>
#include <Eigen/LU>
using namespace mpfr;
using namespace Eigen;
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 2b014682d..29ec7e393 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -26,6 +26,8 @@
#ifndef EIGEN_MATRIX_EXPONENTIAL
#define EIGEN_MATRIX_EXPONENTIAL
+#include "StemFunction.h"
+
#ifdef _MSC_VER
template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
#endif
@@ -148,6 +150,7 @@ class MatrixExponential {
typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef typename std::complex<RealScalar> ComplexScalar;
/** \brief Reference to matrix whose exponential is to be computed. */
typename internal::nested<MatrixType>::type m_M;
@@ -183,7 +186,7 @@ MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) :
m_tmp2(M.rows(),M.cols()),
m_Id(MatrixType::Identity(M.rows(), M.cols())),
m_squarings(0),
- m_l1norm(static_cast<float>(M.cwiseAbs().colwise().sum().maxCoeff()))
+ m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff())
{
/* empty body */
}
@@ -192,12 +195,16 @@ template <typename MatrixType>
template <typename ResultType>
void MatrixExponential<MatrixType>::compute(ResultType &result)
{
- computeUV(RealScalar());
- m_tmp1 = m_U + m_V; // numerator of Pade approximant
- m_tmp2 = -m_U + m_V; // denominator of Pade approximant
- result = m_tmp2.partialPivLu().solve(m_tmp1);
- for (int i=0; i<m_squarings; i++)
- result *= result; // undo scaling by repeated squaring
+ try {
+ computeUV(RealScalar());
+ m_tmp1 = m_U + m_V; // numerator of Pade approximant
+ m_tmp2 = -m_U + m_V; // denominator of Pade approximant
+ result = m_tmp2.partialPivLu().solve(m_tmp1);
+ for (int i=0; i<m_squarings; i++)
+ result *= result; // undo scaling by repeated squaring
+ } catch(int) {
+ result = m_M.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
+ }
}
template <typename MatrixType>
@@ -386,11 +393,9 @@ void MatrixExponential<MatrixType>::computeUV(long double)
MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade17(A);
}
-#else // should never happen
- MatrixType A = m_M / Scalar(2);
- m_U = m_M.sinh();
- m_V = m_M.cosh();
-#endif
+#else // rarely happens
+ throw 0; // will be caught
+#endif // LDBL_MANT_DIG
}
/** \ingroup MatrixFunctions_Module