GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-15 02:43:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Optimize getting exponent from IEEE floating points.

Browse Source
This commit is contained in:
jdh8 2012-08-08 01:27:11 +08:00
commit 8cddcaf839
3 changed files with 24 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Core/Transpose.h
View File

@ -353,7 +353,7 @@ struct check_transpose_aliasing_run_time_selector
{ {
static bool run(const Scalar* dest, const OtherDerived& src) static bool run(const Scalar* dest, const OtherDerived& src)
{ {
return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src)); return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src));
} }
}; };
@ -362,8 +362,8 @@ struct check_transpose_aliasing_run_time_selector<Scalar,DestIsTransposed,CwiseB
{ {
static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src) static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src)
{ {
return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.lhs()))) return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.lhs())))
|| ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.rhs()))); || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.rhs())));
} }
}; };

11
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -743,7 +743,16 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
// RealScalar e2 = abs2(subdiag[end-1]); // RealScalar e2 = abs2(subdiag[end-1]);
// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
// This explain the following, somewhat more complicated, version: // This explain the following, somewhat more complicated, version:
RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e)); RealScalar mu = diag[end];
if(td==0)
mu -= abs(e);
else
{
RealScalar e2 = abs2(subdiag[end-1]);
RealScalar h = hypot(td,e);
if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
else mu -= e2 / (td + (td>0 ? h : -h));
}
RealScalar x = diag[start] - mu; RealScalar x = diag[start] - mu;
RealScalar z = subdiag[start]; RealScalar z = subdiag[start];

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

@ -15,11 +15,6 @@
namespace Eigen { namespace Eigen {
#if defined(_MSC_VER) || defined(__FreeBSD__)
template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
#endif
/** \ingroup MatrixFunctions_Module /** \ingroup MatrixFunctions_Module
* \brief Class for computing the matrix exponential. * \brief Class for computing the matrix exponential.
* \tparam MatrixType type of the argument of the exponential, * \tparam MatrixType type of the argument of the exponential,
@ -297,7 +292,8 @@ void MatrixExponential<MatrixType>::computeUV(float)
pade5(m_M); pade5(m_M);
} else { } else {
const float maxnorm = 3.925724783138660f; const float maxnorm = 3.925724783138660f;
m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); frexp(m_l1norm / maxnorm, &m_squarings);
if (m_squarings < 0) m_squarings = 0;
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade7(A); pade7(A);
} }
@ -319,7 +315,8 @@ void MatrixExponential<MatrixType>::computeUV(double)
pade9(m_M); pade9(m_M);
} else { } else {
const double maxnorm = 5.371920351148152; const double maxnorm = 5.371920351148152;
m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); frexp(m_l1norm / maxnorm, &m_squarings);
if (m_squarings < 0) m_squarings = 0;
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade13(A); pade13(A);
} }
@ -344,7 +341,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
pade9(m_M); pade9(m_M);
} else { } else {
const long double maxnorm = 4.0246098906697353063L; const long double maxnorm = 4.0246098906697353063L;
m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); frexp(m_l1norm / maxnorm, &m_squarings);
if (m_squarings < 0) m_squarings = 0;
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade13(A); pade13(A);
} }
@ -361,7 +359,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
pade13(m_M); pade13(m_M);
} else { } else {
const long double maxnorm = 3.2579440895405400856599663723517L; const long double maxnorm = 3.2579440895405400856599663723517L;
m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); frexp(m_l1norm / maxnorm, &m_squarings);
if (m_squarings < 0) m_squarings = 0;
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade17(A); pade17(A);
} }
@ -378,7 +377,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
pade13(m_M); pade13(m_M);
} else { } else {
const long double maxnorm = 2.884233277829519311757165057717815L; const long double maxnorm = 2.884233277829519311757165057717815L;
m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); frexp(m_l1norm / maxnorm, &m_squarings);
if (m_squarings < 0) m_squarings = 0;
MatrixType A = m_M / pow(Scalar(2), m_squarings); MatrixType A = m_M / pow(Scalar(2), m_squarings);
pade17(A); pade17(A);
} }