fix: <ctime> is necessary for srand(time(NULL))

This commit is contained in:
Chen-Pang He 2011-08-24 18:26:38 +08:00
parent 8414be739b
commit 8ddd1e390b
2 changed files with 19 additions and 13 deletions

View File

@ -30,6 +30,7 @@
#ifndef EIGEN_MPREALSUPPORT_MODULE_H #ifndef EIGEN_MPREALSUPPORT_MODULE_H
#define EIGEN_MPREALSUPPORT_MODULE_H #define EIGEN_MPREALSUPPORT_MODULE_H
#include <ctime>
#include <mpreal.h> #include <mpreal.h>
#include <Eigen/Core> #include <Eigen/Core>
@ -51,7 +52,7 @@ namespace Eigen {
* *
\code \code
#include <iostream> #include <iostream>
#include <Eigen/Mpfrc++Support> #include <Eigen/MPRealSupport>
#include <Eigen/LU> #include <Eigen/LU>
using namespace mpfr; using namespace mpfr;
using namespace Eigen; using namespace Eigen;

View File

@ -26,6 +26,8 @@
#ifndef EIGEN_MATRIX_EXPONENTIAL #ifndef EIGEN_MATRIX_EXPONENTIAL
#define EIGEN_MATRIX_EXPONENTIAL #define EIGEN_MATRIX_EXPONENTIAL
#include "StemFunction.h"
#ifdef _MSC_VER #ifdef _MSC_VER
template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); } template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
#endif #endif
@ -148,6 +150,7 @@ class MatrixExponential {
typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename std::complex<RealScalar> ComplexScalar;
/** \brief Reference to matrix whose exponential is to be computed. */ /** \brief Reference to matrix whose exponential is to be computed. */
typename internal::nested<MatrixType>::type m_M; typename internal::nested<MatrixType>::type m_M;
@ -183,7 +186,7 @@ MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) :
m_tmp2(M.rows(),M.cols()), m_tmp2(M.rows(),M.cols()),
m_Id(MatrixType::Identity(M.rows(), M.cols())), m_Id(MatrixType::Identity(M.rows(), M.cols())),
m_squarings(0), m_squarings(0),
m_l1norm(static_cast<float>(M.cwiseAbs().colwise().sum().maxCoeff())) m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff())
{ {
/* empty body */ /* empty body */
} }
@ -192,12 +195,16 @@ template <typename MatrixType>
template <typename ResultType> template <typename ResultType>
void MatrixExponential<MatrixType>::compute(ResultType &result) void MatrixExponential<MatrixType>::compute(ResultType &result)
{ {
computeUV(RealScalar()); try {
m_tmp1 = m_U + m_V; // numerator of Pade approximant computeUV(RealScalar());
m_tmp2 = -m_U + m_V; // denominator of Pade approximant m_tmp1 = m_U + m_V; // numerator of Pade approximant
result = m_tmp2.partialPivLu().solve(m_tmp1); m_tmp2 = -m_U + m_V; // denominator of Pade approximant
for (int i=0; i<m_squarings; i++) result = m_tmp2.partialPivLu().solve(m_tmp1);
result *= result; // undo scaling by repeated squaring 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> template <typename MatrixType>
@ -386,11 +393,9 @@ void MatrixExponential<MatrixType>::computeUV(long double)
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade17(A); pade17(A);
} }
#else // should never happen #else // rarely happens
MatrixType A = m_M / Scalar(2); throw 0; // will be caught
m_U = m_M.sinh(); #endif // LDBL_MANT_DIG
m_V = m_M.cosh();
#endif
} }
/** \ingroup MatrixFunctions_Module /** \ingroup MatrixFunctions_Module