GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-30 07:44:10 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge upstream updates.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2016-04-14 13:59:50 -07:00
commit d2e95492e7
7 changed files with 89 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/src/Core/arch/CUDA/Half.h
View File

@ -510,6 +510,11 @@ static EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC int (isfinite)(const Eigen::half& a
namespace std {
EIGEN_ALWAYS_INLINE ostream& operator << (ostream& os, const Eigen::half& v) {
os << static_cast<float>(v);
return os;
}
#if __cplusplus > 199711L
template <>
struct hash<Eigen::half> {

40
Eigen/src/SVD/JacobiSVD.h
View File

@ -350,7 +350,8 @@ template<typename MatrixType, int QRPreconditioner>
struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
{
typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, const typename MatrixType::RealScalar&) { return true; }
typedef typename MatrixType::RealScalar RealScalar;
static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
};
template<typename MatrixType, int QRPreconditioner>
@ -359,25 +360,30 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, const typename MatrixType::RealScalar& precision)
static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
{
using std::sqrt;
using std::abs;
Scalar z;
JacobiRotation<Scalar> rot;
RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
const RealScalar precision = NumTraits<Scalar>::epsilon();
if(n==0)
{
// make sure first column is zero
work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
if(work_matrix.coeff(p,q)!=Scalar(0))
if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
{
// work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when computing n
z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.row(p) *= z;
if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
}
if(work_matrix.coeff(q,q)!=Scalar(0))
if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
{
z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
work_matrix.row(q) *= z;
@ -391,13 +397,13 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
rot.s() = work_matrix.coeff(q,p) / n;
work_matrix.applyOnTheLeft(p,q,rot);
if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
if(work_matrix.coeff(p,q) != Scalar(0))
if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
{
z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.col(q) *= z;
if(svd.computeV()) svd.m_matrixV.col(q) *= z;
}
if(work_matrix.coeff(q,q) != Scalar(0))
if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
{
z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
work_matrix.row(q) *= z;
@ -405,11 +411,11 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
}
}
const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
RealScalar threshold = numext::maxi<RealScalar>(considerAsZero,
precision * numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
// return true if we still have some work to do
return numext::abs(work_matrix(p,q)) > threshold || numext::abs(work_matrix(q,p)) > threshold;
// update largest diagonal entry
maxDiagEntry = numext::maxi(maxDiagEntry,numext::maxi(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
// and check whether the 2x2 block is already diagonal
RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
}
};
@ -426,7 +432,6 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> rot1;
RealScalar t = m.coeff(0,0) + m.coeff(1,1);
RealScalar d = m.coeff(1,0) - m.coeff(0,1);
if(d == RealScalar(0))
{
rot1.s() = RealScalar(0);
@ -719,6 +724,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
}
/*** step 2. The main Jacobi SVD iteration. ***/
RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
bool finished = false;
while(!finished)
@ -734,16 +740,13 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
// if this 2x2 sub-matrix is not diagonal already...
// notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
// keep us iterating forever. Similarly, small denormal numbers are considered zero.
RealScalar threshold = numext::maxi<RealScalar>(considerAsZero,
precision * numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)),
abs(m_workMatrix.coeff(q,q))));
// We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
{
finished = false;
// perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
// the complex to real operation returns true is the updated 2x2 block is not already diagonal
if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, precision))
if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
{
JacobiRotation<RealScalar> j_left, j_right;
internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
@ -754,6 +757,9 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
m_workMatrix.applyOnTheRight(p,q,j_right);
if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
// keep track of the largest diagonal coefficient
maxDiagEntry = numext::maxi(maxDiagEntry,numext::maxi(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
}
}
}

12
test/main.h
View File

@ -316,9 +316,9 @@ inline bool test_isMuchSmallerThan(const float& a, const float& b)
{ return internal::isMuchSmallerThan(a, b, test_precision<float>()); }
inline bool test_isApproxOrLessThan(const float& a, const float& b)
{ return internal::isApproxOrLessThan(a, b, test_precision<float>()); }
inline bool test_isApprox(const double& a, const double& b)
{ return internal::isApprox(a, b, test_precision<double>()); }
inline bool test_isMuchSmallerThan(const double& a, const double& b)
{ return internal::isMuchSmallerThan(a, b, test_precision<double>()); }
inline bool test_isApproxOrLessThan(const double& a, const double& b)
@ -359,6 +359,12 @@ inline bool test_isApproxOrLessThan(const long double& a, const long double& b)
{ return internal::isApproxOrLessThan(a, b, test_precision<long double>()); }
#endif // EIGEN_TEST_NO_LONGDOUBLE
inline bool test_isApprox(const half& a, const half& b)
{ return internal::isApprox(a, b, test_precision<half>()); }
inline bool test_isMuchSmallerThan(const half& a, const half& b)
{ return internal::isMuchSmallerThan(a, b, test_precision<half>()); }
inline bool test_isApproxOrLessThan(const half& a, const half& b)
{ return internal::isApproxOrLessThan(a, b, test_precision<half>()); }
// test_relative_error returns the relative difference between a and b as a real scalar as used in isApprox.
template<typename T1,typename T2>
@ -426,9 +432,7 @@ template<typename T1,typename T2>
typename NumTraits<T1>::Real test_relative_error(const T1 &a, const T2 &b, typename internal::enable_if<internal::is_arithmetic<typename NumTraits<T1>::Real>::value, T1>::type* = 0)
{
typedef typename NumTraits<T1>::Real RealScalar;
using std::min;
using std::sqrt;
return sqrt(RealScalar(numext::abs2(a-b))/RealScalar((min)(numext::abs2(a),numext::abs2(b))));
return numext::sqrt(RealScalar(numext::abs2(a-b))/RealScalar((numext::mini)(numext::abs2(a),numext::abs2(b))));
}
template<typename T>

12
test/svd_fill.h
View File

@ -80,6 +80,8 @@ void svd_fill_random(MatrixType &m, int Option = 0)
Index i = internal::random<Index>(0,m.rows()-1);
Index j = internal::random<Index>(0,m.cols()-1);
m(j,i) = m(i,j) = samples(internal::random<Index>(0,samples.size()-1));
if(NumTraits<Scalar>::IsComplex)
*(&numext::real_ref(m(j,i))+1) = *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1));
}
}
}
@ -91,8 +93,14 @@ void svd_fill_random(MatrixType &m, int Option = 0)
if(!(dup && unit_uv))
{
Index n = internal::random<Index>(0,m.size()-1);
for(Index i=0; i<n; ++i)
m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = samples(internal::random<Index>(0,samples.size()-1));
for(Index k=0; k<n; ++k)
{
Index i = internal::random<Index>(0,m.rows()-1);
Index j = internal::random<Index>(0,m.cols()-1);
m(i,j) = samples(internal::random<Index>(0,samples.size()-1));
if(NumTraits<Scalar>::IsComplex)
*(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1));
}
}
}
}

4
unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h
View File

@ -64,7 +64,7 @@ struct scalar_sigmoid_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_sigmoid_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const {
const T one = T(1);
return one / (one + std::exp(-x));
return one / (one + numext::exp(-x));
}
template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
@ -803,7 +803,7 @@ class GaussianGenerator {
T offset = coordinates[i] - m_means[i];
tmp += offset * offset / m_two_sigmas[i];
}
return std::exp(-tmp);
return numext::exp(-tmp);
}
private:

4
unsupported/Eigen/CXX11/src/Tensor/TensorUInt128.h
View File

@ -53,9 +53,7 @@ struct TensorUInt128
template<typename T>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
explicit TensorUInt128(const T& x) : high(0), low(x) {
typedef typename conditional<sizeof(T) == 8, uint64_t, uint32_t>::type UnsignedT;
typedef typename conditional<sizeof(LOW) == 8, uint64_t, uint32_t>::type UnsignedLow;
eigen_assert(static_cast<UnsignedT>(x) <= static_cast<UnsignedLow>(NumTraits<LOW>::highest()));
eigen_assert((static_cast<typename conditional<sizeof(T) == 8, uint64_t, uint32_t>::type>(x) <= static_cast<typename conditional<sizeof(LOW) == 8, uint64_t, uint32_t>::type>(NumTraits<LOW>::highest())));
eigen_assert(x >= 0);
}

41
unsupported/test/cxx11_float16.cpp
View File

@ -122,6 +122,8 @@ void test_comparison()
VERIFY(half(1.0f) != half(2.0f));
// Comparisons with NaNs and infinities.
#if !EIGEN_COMP_MSVC
// Visual Studio errors out on divisions by 0
VERIFY(!(half(0.0 / 0.0) == half(0.0 / 0.0)));
VERIFY(half(0.0 / 0.0) != half(0.0 / 0.0));
@ -132,13 +134,26 @@ void test_comparison()
VERIFY(half(1.0) < half(1.0 / 0.0));
VERIFY(half(1.0) > half(-1.0 / 0.0));
#endif
}
void test_functions()
void test_basic_functions()
{
VERIFY_IS_EQUAL(float(numext::abs(half(3.5f))), 3.5f);
VERIFY_IS_EQUAL(float(numext::abs(half(-3.5f))), 3.5f);
VERIFY_IS_EQUAL(float(numext::floor(half(3.5f))), 3.0f);
VERIFY_IS_EQUAL(float(numext::floor(half(-3.5f))), -4.0f);
VERIFY_IS_EQUAL(float(numext::ceil(half(3.5f))), 4.0f);
VERIFY_IS_EQUAL(float(numext::ceil(half(-3.5f))), -3.0f);
VERIFY_IS_APPROX(float(numext::sqrt(half(0.0f))), 0.0f);
VERIFY_IS_APPROX(float(numext::sqrt(half(4.0f))), 2.0f);
VERIFY_IS_APPROX(float(numext::pow(half(0.0f), half(1.0f))), 0.0f);
VERIFY_IS_APPROX(float(numext::pow(half(2.0f), half(2.0f))), 4.0f);
VERIFY_IS_EQUAL(float(numext::exp(half(0.0f))), 1.0f);
VERIFY_IS_APPROX(float(numext::exp(half(EIGEN_PI))), float(20.0 + EIGEN_PI));
@ -146,10 +161,32 @@ void test_functions()
VERIFY_IS_APPROX(float(numext::log(half(10.0f))), 2.30273f);
}
void test_trigonometric_functions()
{
VERIFY_IS_APPROX(numext::cos(half(0.0f)), half(cosf(0.0f)));
VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI)), half(cosf(EIGEN_PI)));
//VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI/2)), half(cosf(EIGEN_PI/2)));
//VERIFY_IS_APPROX(numext::cos(half(3*EIGEN_PI/2)), half(cosf(3*EIGEN_PI/2)));
VERIFY_IS_APPROX(numext::cos(half(3.5f)), half(cosf(3.5f)));
VERIFY_IS_APPROX(numext::sin(half(0.0f)), half(sinf(0.0f)));
// VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI)), half(sinf(EIGEN_PI)));
VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI/2)), half(sinf(EIGEN_PI/2)));
VERIFY_IS_APPROX(numext::sin(half(3*EIGEN_PI/2)), half(sinf(3*EIGEN_PI/2)));
VERIFY_IS_APPROX(numext::sin(half(3.5f)), half(sinf(3.5f)));
VERIFY_IS_APPROX(numext::tan(half(0.0f)), half(tanf(0.0f)));
// VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI)), half(tanf(EIGEN_PI)));
// VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI/2)), half(tanf(EIGEN_PI/2)));
//VERIFY_IS_APPROX(numext::tan(half(3*EIGEN_PI/2)), half(tanf(3*EIGEN_PI/2)));
VERIFY_IS_APPROX(numext::tan(half(3.5f)), half(tanf(3.5f)));
}
void test_cxx11_float16()
{
CALL_SUBTEST(test_conversion());
CALL_SUBTEST(test_arithmetic());
CALL_SUBTEST(test_comparison());
CALL_SUBTEST(test_functions());
CALL_SUBTEST(test_basic_functions());
CALL_SUBTEST(test_trigonometric_functions());
}