From d60cca32e54ea02757713935e852cd53343d52c1 Mon Sep 17 00:00:00 2001 From: ermak Date: Sat, 17 Dec 2016 00:45:13 +0700 Subject: [PATCH 01/19] Transformation methods added to ParametrizedLine class. --- Eigen/src/Geometry/ParametrizedLine.h | 39 ++++++++++++++++++++++++++- test/geo_parametrizedline.cpp | 27 +++++++++++++++++++ 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Geometry/ParametrizedLine.h b/Eigen/src/Geometry/ParametrizedLine.h index 1e985d8cd..3929ca87f 100644 --- a/Eigen/src/Geometry/ParametrizedLine.h +++ b/Eigen/src/Geometry/ParametrizedLine.h @@ -104,7 +104,44 @@ public: template EIGEN_DEVICE_FUNC VectorType intersectionPoint(const Hyperplane<_Scalar, _AmbientDim, OtherOptions>& hyperplane) const; - /** \returns \c *this with scalar type casted to \a NewScalarType + /** Applies the transformation matrix \a mat to \c *this and returns a reference to \c *this. + * + * \param mat the Dim x Dim transformation matrix + * \param traits specifies whether the matrix \a mat represents an #Isometry + * or a more generic #Affine transformation. The default is #Affine. + */ + template + EIGEN_DEVICE_FUNC inline ParametrizedLine& transform(const MatrixBase& mat, TransformTraits traits = Affine) + { + if (traits==Affine) + direction() = (mat * direction()).normalized(); + else if (traits==Isometry) + direction() = mat * direction(); + else + { + eigen_assert(0 && "invalid traits value in ParametrizedLine::transform()"); + } + origin() = mat * origin(); + return *this; + } + + /** Applies the transformation \a t to \c *this and returns a reference to \c *this. + * + * \param t the transformation of dimension Dim + * \param traits specifies whether the transformation \a t represents an #Isometry + * or a more generic #Affine transformation. The default is #Affine. + * Other kind of transformations are not supported. + */ + template + EIGEN_DEVICE_FUNC inline ParametrizedLine& transform(const Transform& t, + TransformTraits traits = Affine) + { + transform(t.linear(), traits); + origin() += t.translation(); + return *this; + } + +/** \returns \c *this with scalar type casted to \a NewScalarType * * Note that if \a NewScalarType is equal to the current scalar type of \c *this * then this function smartly returns a const reference to \c *this. diff --git a/test/geo_parametrizedline.cpp b/test/geo_parametrizedline.cpp index 9bf5f3c1d..29c1b105c 100644 --- a/test/geo_parametrizedline.cpp +++ b/test/geo_parametrizedline.cpp @@ -25,6 +25,8 @@ template void parametrizedline(const LineType& _line) typedef typename NumTraits::Real RealScalar; typedef Matrix VectorType; typedef Hyperplane HyperplaneType; + typedef Matrix MatrixType; VectorType p0 = VectorType::Random(dim); VectorType p1 = VectorType::Random(dim); @@ -59,6 +61,31 @@ template void parametrizedline(const LineType& _line) VERIFY_IS_MUCH_SMALLER_THAN(hp.signedDistance(pi), RealScalar(1)); VERIFY_IS_MUCH_SMALLER_THAN(l0.distance(pi), RealScalar(1)); VERIFY_IS_APPROX(l0.intersectionPoint(hp), pi); + + // transform + if (!NumTraits::IsComplex) + { + MatrixType rot = MatrixType::Random(dim,dim).householderQr().householderQ(); + DiagonalMatrix scaling(VectorType::Random()); + Translation translation(VectorType::Random()); + + while(scaling.diagonal().cwiseAbs().minCoeff() void parametrizedline_alignment() From fc94258e77ca8a69dca81ab29aa011b173ed15cd Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 18 Dec 2016 22:11:48 +0000 Subject: [PATCH 02/19] Fix unused warning --- Eigen/src/Core/arch/AltiVec/MathFunctions.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Eigen/src/Core/arch/AltiVec/MathFunctions.h b/Eigen/src/Core/arch/AltiVec/MathFunctions.h index 5511245dd..c5e4bede7 100644 --- a/Eigen/src/Core/arch/AltiVec/MathFunctions.h +++ b/Eigen/src/Core/arch/AltiVec/MathFunctions.h @@ -84,8 +84,10 @@ static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_q3, 2.00000000000000000009e0); static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C1, 0.693145751953125); static _EIGEN_DECLARE_CONST_Packet2d(cephes_exp_C2, 1.42860682030941723212e-6); +#ifdef __POWER8_VECTOR__ static Packet2l p2l_1023 = { 1023, 1023 }; static Packet2ul p2ul_52 = { 52, 52 }; +#endif #endif From 8c0e70150433e8fe50c980ff629a9f80162eaf92 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 18 Dec 2016 22:13:19 +0000 Subject: [PATCH 03/19] bug #1360: fix sign issue with pmull on altivec --- Eigen/src/Core/arch/AltiVec/Complex.h | 10 +++++----- Eigen/src/Core/arch/AltiVec/PacketMath.h | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 45213f791..59367ba29 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -15,14 +15,14 @@ namespace Eigen { namespace internal { -static Packet4ui p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_ZERO_);//{ 0x00000000, 0x80000000, 0x00000000, 0x80000000 }; +static Packet4ui p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_MZERO);//{ 0x00000000, 0x80000000, 0x00000000, 0x80000000 }; #ifdef __VSX__ #if defined(_BIG_ENDIAN) -static Packet2ul p2ul_CONJ_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2d_ZERO_, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; -static Packet2ul p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO, (Packet4ui) p2d_ZERO_, 8);//{ 0x8000000000000000, 0x0000000000000000 }; +static Packet2ul p2ul_CONJ_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2d_MZERO, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; +static Packet2ul p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO, (Packet4ui) p2d_MZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; #else -static Packet2ul p2ul_CONJ_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO, (Packet4ui) p2d_ZERO_, 8);//{ 0x8000000000000000, 0x0000000000000000 }; -static Packet2ul p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2d_ZERO_, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; +static Packet2ul p2ul_CONJ_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO, (Packet4ui) p2d_MZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; +static Packet2ul p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2d_MZERO, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; #endif #endif diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index cbfef3503..e7d4f4d8d 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -72,7 +72,7 @@ static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0); //{ 0, 0, 0, 0,} static _EIGEN_DECLARE_CONST_FAST_Packet4i(ONE,1); //{ 1, 1, 1, 1} static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS16,-16); //{ -16, -16, -16, -16} static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS1,-1); //{ -1, -1, -1, -1} -static Packet4f p4f_ZERO_ = (Packet4f) vec_sl((Packet4ui)p4i_MINUS1, (Packet4ui)p4i_MINUS1); //{ 0x80000000, 0x80000000, 0x80000000, 0x80000000} +static Packet4f p4f_MZERO = (Packet4f) vec_sl((Packet4ui)p4i_MINUS1, (Packet4ui)p4i_MINUS1); //{ 0x80000000, 0x80000000, 0x80000000, 0x80000000} #ifndef __VSX__ static Packet4f p4f_ONE = vec_ctf(p4i_ONE, 0); //{ 1.0, 1.0, 1.0, 1.0} #endif @@ -358,7 +358,7 @@ template<> EIGEN_STRONG_INLINE Packet4i pnegate(const Packet4i& a) { return p4i_ template<> EIGEN_STRONG_INLINE Packet4f pconj(const Packet4f& a) { return a; } template<> EIGEN_STRONG_INLINE Packet4i pconj(const Packet4i& a) { return a; } -template<> EIGEN_STRONG_INLINE Packet4f pmul(const Packet4f& a, const Packet4f& b) { return vec_madd(a,b, p4f_ZERO); } +template<> EIGEN_STRONG_INLINE Packet4f pmul(const Packet4f& a, const Packet4f& b) { return vec_madd(a,b, p4f_MZERO); } template<> EIGEN_STRONG_INLINE Packet4i pmul(const Packet4i& a, const Packet4i& b) { return a * b; } template<> EIGEN_STRONG_INLINE Packet4f pdiv(const Packet4f& a, const Packet4f& b) @@ -373,7 +373,7 @@ template<> EIGEN_STRONG_INLINE Packet4f pdiv(const Packet4f& a, const t = vec_nmsub(y_0, b, p4f_ONE); y_1 = vec_madd(y_0, t, y_0); - return vec_madd(a, y_1, p4f_ZERO); + return vec_madd(a, y_1, p4f_MZERO); #else return vec_div(a, b); #endif @@ -766,7 +766,7 @@ static Packet2l p2l_ONE = { 1, 1 }; static Packet2l p2l_ZERO = reinterpret_cast(p4i_ZERO); static Packet2d p2d_ONE = { 1.0, 1.0 }; static Packet2d p2d_ZERO = reinterpret_cast(p4f_ZERO); -static Packet2d p2d_ZERO_ = { -0.0, -0.0 }; +static Packet2d p2d_MZERO = { -0.0, -0.0 }; #ifdef _BIG_ENDIAN static Packet2d p2d_COUNTDOWN = reinterpret_cast(vec_sld(reinterpret_cast(p2d_ZERO), reinterpret_cast(p2d_ONE), 8)); @@ -904,7 +904,7 @@ template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a) { return p2d_ template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; } -template<> EIGEN_STRONG_INLINE Packet2d pmul(const Packet2d& a, const Packet2d& b) { return vec_madd(a,b,p2d_ZERO); } +template<> EIGEN_STRONG_INLINE Packet2d pmul(const Packet2d& a, const Packet2d& b) { return vec_madd(a,b,p2d_MZERO); } template<> EIGEN_STRONG_INLINE Packet2d pdiv(const Packet2d& a, const Packet2d& b) { return vec_div(a,b); } // for some weird raisons, it has to be overloaded for packet of integers From fb1d0138ece3206cd67c04893e2a9d6a71594816 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Mon, 19 Dec 2016 07:32:48 -0800 Subject: [PATCH 04/19] Include SSE packet instructions when compiling with avx512 enabled. --- Eigen/Core | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Eigen/Core b/Eigen/Core index 444c1c8d7..16be82ac2 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -150,6 +150,10 @@ #define EIGEN_VECTORIZE_AVX2 #define EIGEN_VECTORIZE_AVX #define EIGEN_VECTORIZE_FMA + #define EIGEN_VECTORIZE_SSE3 + #define EIGEN_VECTORIZE_SSSE3 + #define EIGEN_VECTORIZE_SSE4_1 + #define EIGEN_VECTORIZE_SSE4_2 #ifdef __AVX512DQ__ #define EIGEN_VECTORIZE_AVX512DQ #endif @@ -360,6 +364,8 @@ using std::ptrdiff_t; #include "src/Core/arch/SSE/PacketMath.h" #include "src/Core/arch/AVX/PacketMath.h" #include "src/Core/arch/AVX512/PacketMath.h" + #include "src/Core/arch/SSE/MathFunctions.h" + #include "src/Core/arch/AVX/MathFunctions.h" #include "src/Core/arch/AVX512/MathFunctions.h" #elif defined EIGEN_VECTORIZE_AVX // Use AVX for floats and doubles, SSE for integers From 751e097c57e84d368d782c4a18b960ed2350c2f0 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Mon, 19 Dec 2016 13:44:46 -0500 Subject: [PATCH 05/19] Use 32 registers on ARM64 --- Eigen/src/Core/arch/NEON/PacketMath.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h index 2a8f58d74..d392bf3ff 100644 --- a/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/Eigen/src/Core/arch/NEON/PacketMath.h @@ -28,11 +28,13 @@ namespace internal { #define EIGEN_HAS_SINGLE_INSTRUCTION_CJMADD #endif -// FIXME NEON has 16 quad registers, but since the current register allocator -// is so bad, it is much better to reduce it to 8 #ifndef EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS +#if EIGEN_ARCH_ARM64 +#define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS 32 +#else #define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS 16 #endif +#endif typedef float32x2_t Packet2f; typedef float32x4_t Packet4f; From 923acadfacef98ef234ed108cc6c3de877c0fe89 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Mon, 19 Dec 2016 13:02:27 -0800 Subject: [PATCH 06/19] Fixed compilation errors with gcc6 when compiling the AVX512 intrinsics --- Eigen/src/Core/arch/AVX512/PacketMath.h | 62 ++++++------------------- 1 file changed, 15 insertions(+), 47 deletions(-) diff --git a/Eigen/src/Core/arch/AVX512/PacketMath.h b/Eigen/src/Core/arch/AVX512/PacketMath.h index f6500a16e..0580b80f8 100644 --- a/Eigen/src/Core/arch/AVX512/PacketMath.h +++ b/Eigen/src/Core/arch/AVX512/PacketMath.h @@ -461,53 +461,21 @@ EIGEN_STRONG_INLINE Packet16i ploadu(const int* from) { // {a0, a0 a1, a1, a2, a2, a3, a3, a4, a4, a5, a5, a6, a6, a7, a7} template <> EIGEN_STRONG_INLINE Packet16f ploaddup(const float* from) { - Packet8f lane0 = _mm256_broadcast_ps((const __m128*)(const void*)from); - // mimic an "inplace" permutation of the lower 128bits using a blend - lane0 = _mm256_blend_ps( - lane0, _mm256_castps128_ps256(_mm_permute_ps( - _mm256_castps256_ps128(lane0), _MM_SHUFFLE(1, 0, 1, 0))), - 15); - // then we can perform a consistent permutation on the global register to get - // everything in shape: - lane0 = _mm256_permute_ps(lane0, _MM_SHUFFLE(3, 3, 2, 2)); - - Packet8f lane1 = _mm256_broadcast_ps((const __m128*)(const void*)(from + 4)); - // mimic an "inplace" permutation of the lower 128bits using a blend - lane1 = _mm256_blend_ps( - lane1, _mm256_castps128_ps256(_mm_permute_ps( - _mm256_castps256_ps128(lane1), _MM_SHUFFLE(1, 0, 1, 0))), - 15); - // then we can perform a consistent permutation on the global register to get - // everything in shape: - lane1 = _mm256_permute_ps(lane1, _MM_SHUFFLE(3, 3, 2, 2)); - -#ifdef EIGEN_VECTORIZE_AVX512DQ - Packet16f res = _mm512_undefined_ps(); - return _mm512_insertf32x8(res, lane0, 0); - return _mm512_insertf32x8(res, lane1, 1); - return res; -#else - Packet16f res = _mm512_undefined_ps(); - res = _mm512_insertf32x4(res, _mm256_extractf128_ps(lane0, 0), 0); - res = _mm512_insertf32x4(res, _mm256_extractf128_ps(lane0, 1), 1); - res = _mm512_insertf32x4(res, _mm256_extractf128_ps(lane1, 0), 2); - res = _mm512_insertf32x4(res, _mm256_extractf128_ps(lane1, 1), 3); - return res; -#endif + __m256i low_half = _mm256_load_si256(reinterpret_cast(from)); + __m512 even_elements = _mm512_castsi512_ps(_mm512_cvtepu32_epi64(low_half)); + __m512 pairs = _mm512_permute_ps(even_elements, _MM_SHUFFLE(2, 2, 0, 0)); + return pairs; } // Loads 4 doubles from memory a returns the packet {a0, a0 a1, a1, a2, a2, a3, // a3} template <> EIGEN_STRONG_INLINE Packet8d ploaddup(const double* from) { - Packet4d lane0 = _mm256_broadcast_pd((const __m128d*)(const void*)from); - lane0 = _mm256_permute_pd(lane0, 3 << 2); - - Packet4d lane1 = _mm256_broadcast_pd((const __m128d*)(const void*)(from + 2)); - lane1 = _mm256_permute_pd(lane1, 3 << 2); - - Packet8d res = _mm512_undefined_pd(); - res = _mm512_insertf64x4(res, lane0, 0); - return _mm512_insertf64x4(res, lane1, 1); + __m512d x = _mm512_setzero_pd(); + x = _mm512_insertf64x2(x, _mm_loaddup_pd(&from[0]), 0); + x = _mm512_insertf64x2(x, _mm_loaddup_pd(&from[1]), 1); + x = _mm512_insertf64x2(x, _mm_loaddup_pd(&from[2]), 2); + x = _mm512_insertf64x2(x, _mm_loaddup_pd(&from[3]), 3); + return x; } // Loads 4 floats from memory a returns the packet @@ -525,11 +493,11 @@ EIGEN_STRONG_INLINE Packet16f ploadquad(const float* from) { // {a0, a0 a0, a0, a1, a1, a1, a1} template <> EIGEN_STRONG_INLINE Packet8d ploadquad(const double* from) { - Packet8d tmp = _mm512_undefined_pd(); - Packet2d tmp0 = _mm_load_pd1(from); - Packet2d tmp1 = _mm_load_pd1(from + 1); - Packet4d lane0 = _mm256_broadcastsd_pd(tmp0); - Packet4d lane1 = _mm256_broadcastsd_pd(tmp1); + __m128d tmp0 = _mm_load_pd1(from); + __m256d lane0 = _mm256_broadcastsd_pd(tmp0); + __m128d tmp1 = _mm_load_pd1(from + 1); + __m256d lane1 = _mm256_broadcastsd_pd(tmp1); + __m512d tmp = _mm512_undefined_pd(); tmp = _mm512_insertf64x4(tmp, lane0, 0); return _mm512_insertf64x4(tmp, lane1, 1); } From 27ceb43bf6b06dda898e5d027097f33a970f7355 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Mon, 19 Dec 2016 15:34:42 -0800 Subject: [PATCH 07/19] Fixed race condition in the tensor_shuffling_sycl test --- unsupported/test/cxx11_tensor_shuffling_sycl.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/unsupported/test/cxx11_tensor_shuffling_sycl.cpp b/unsupported/test/cxx11_tensor_shuffling_sycl.cpp index b2b75cbde..c4521aac8 100644 --- a/unsupported/test/cxx11_tensor_shuffling_sycl.cpp +++ b/unsupported/test/cxx11_tensor_shuffling_sycl.cpp @@ -57,6 +57,7 @@ static void test_simple_shuffling_sycl(const Eigen::SyclDevice& sycl_device) gpu2.device(sycl_device)=gpu1.shuffle(shuffles); sycl_device.memcpyDeviceToHost(no_shuffle.data(), gpu_data2, buffSize); + sycl_device.synchronize(); VERIFY_IS_EQUAL(no_shuffle.dimension(0), sizeDim1); VERIFY_IS_EQUAL(no_shuffle.dimension(1), sizeDim2); @@ -82,8 +83,9 @@ static void test_simple_shuffling_sycl(const Eigen::SyclDevice& sycl_device) DataType* gpu_data3 = static_cast(sycl_device.allocate(buffSize)); TensorMap> gpu3(gpu_data3, tensorrangeShuffle); - gpu3.device(sycl_device)=gpu1.shuffle(shuffles); - sycl_device.memcpyDeviceToHost(shuffle.data(), gpu_data3, buffSize); + gpu3.device(sycl_device)=gpu1.shuffle(shuffles); + sycl_device.memcpyDeviceToHost(shuffle.data(), gpu_data3, buffSize); + sycl_device.synchronize(); VERIFY_IS_EQUAL(shuffle.dimension(0), sizeDim3); VERIFY_IS_EQUAL(shuffle.dimension(1), sizeDim4); From f5d644b4155e2cbf2a03988c2ef592a3d2857031 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 09:35:00 +0100 Subject: [PATCH 08/19] Make sure that HyperPlane::transform manitains a unit normal vector in the Affine case. --- Eigen/src/Geometry/Hyperplane.h | 3 +++ test/geo_hyperplane.cpp | 3 +++ 2 files changed, 6 insertions(+) diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h index 07f2659b2..05929b299 100644 --- a/Eigen/src/Geometry/Hyperplane.h +++ b/Eigen/src/Geometry/Hyperplane.h @@ -217,7 +217,10 @@ public: EIGEN_DEVICE_FUNC inline Hyperplane& transform(const MatrixBase& mat, TransformTraits traits = Affine) { if (traits==Affine) + { normal() = mat.inverse().transpose() * normal(); + m_coeffs /= normal().norm(); + } else if (traits==Isometry) normal() = mat * normal(); else diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index e77702bc7..27892850d 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -66,12 +66,15 @@ template void hyperplane(const HyperplaneType& _plane) VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot,Isometry).absDistance(rot * p1), Scalar(1) ); pl2 = pl1; VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling).absDistance((rot*scaling) * p1), Scalar(1) ); + VERIFY_IS_APPROX( pl2.normal().norm(), RealScalar(1) ); pl2 = pl1; VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling*translation) .absDistance((rot*scaling*translation) * p1), Scalar(1) ); + VERIFY_IS_APPROX( pl2.normal().norm(), RealScalar(1) ); pl2 = pl1; VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*translation,Isometry) .absDistance((rot*translation) * p1), Scalar(1) ); + VERIFY_IS_APPROX( pl2.normal().norm(), RealScalar(1) ); } // casting From 10c6bcdc2eac68d0b48ecfe4d391079610eb6124 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Mon, 19 Dec 2016 14:07:42 +0100 Subject: [PATCH 09/19] Add support for long indexes and for (real-valued) row-major matrices to CholmodSupport module --- Eigen/src/CholmodSupport/CholmodSupport.h | 71 ++++++++++++++++++----- test/cholmod_support.cpp | 42 +++++++++----- 2 files changed, 84 insertions(+), 29 deletions(-) diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index 571972023..a07efb1d5 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -32,7 +32,7 @@ template<> struct cholmod_configure_matrix > { } }; -// Other scalar types are not yet suppotred by Cholmod +// Other scalar types are not yet supported by Cholmod // template<> struct cholmod_configure_matrix { // template // static void run(CholmodType& mat) { @@ -124,6 +124,9 @@ cholmod_sparse viewAsCholmod(const SparseSelfAdjointView::IsComplex == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); + if(_Options & RowMajorBit) res.stype *=-1; return res; } @@ -159,6 +162,44 @@ MappedSparseMatrix viewAsEigen(cholmod_sparse& cm) static_cast(cm.p), static_cast(cm.i),static_cast(cm.x) ); } +namespace internal { + +// template specializations for int and long that call the correct cholmod method + +#define EIGEN_CHOLMOD_SPECIALIZE0(ret, name) \ + template ret cm_ ## name (cholmod_common &Common) { return cholmod_ ## name (&Common); } \ + template<> ret cm_ ## name (cholmod_common &Common) { return cholmod_l_ ## name (&Common); } + +#define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1) \ + template ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_ ## name (&a1, &Common); } \ + template<> ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_l_ ## name (&a1, &Common); } + +EIGEN_CHOLMOD_SPECIALIZE0(int, start) +EIGEN_CHOLMOD_SPECIALIZE0(int, finish) + +EIGEN_CHOLMOD_SPECIALIZE1(int, free_factor, cholmod_factor*, L) +EIGEN_CHOLMOD_SPECIALIZE1(int, free_dense, cholmod_dense*, X) +EIGEN_CHOLMOD_SPECIALIZE1(int, free_sparse, cholmod_sparse*, A) + +EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A) + +template cholmod_dense* cm_solve (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_solve (sys, &L, &B, &Common); } +template<> cholmod_dense* cm_solve (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_l_solve (sys, &L, &B, &Common); } + +template cholmod_sparse* cm_spsolve (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_spsolve (sys, &L, &B, &Common); } +template<> cholmod_sparse* cm_spsolve (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_l_spsolve (sys, &L, &B, &Common); } + +template +int cm_factorize_p (cholmod_sparse* A, double beta[2], _StorageIndex* fset, size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_factorize_p (A, beta, fset, fsize, L, &Common); } +template<> +int cm_factorize_p (cholmod_sparse* A, double beta[2], long* fset, size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_l_factorize_p (A, beta, fset, fsize, L, &Common); } + +#undef EIGEN_CHOLMOD_SPECIALIZE0 +#undef EIGEN_CHOLMOD_SPECIALIZE1 + +} // namespace internal + + enum CholmodMode { CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt }; @@ -195,7 +236,7 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + internal::cm_start(m_cholmod); } explicit CholmodBase(const MatrixType& matrix) @@ -203,15 +244,15 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + internal::cm_start(m_cholmod); compute(matrix); } ~CholmodBase() { if(m_cholmodFactor) - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); - cholmod_finish(&m_cholmod); + internal::cm_free_factor(m_cholmodFactor, m_cholmod); + internal::cm_finish(m_cholmod); } inline StorageIndex cols() const { return internal::convert_index(m_cholmodFactor->n); } @@ -219,7 +260,7 @@ class CholmodBase : public SparseSolverBase /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, + * \returns \c Success if computation was successful, * \c NumericalIssue if the matrix.appears to be negative. */ ComputationInfo info() const @@ -246,11 +287,11 @@ class CholmodBase : public SparseSolverBase { if(m_cholmodFactor) { - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); + internal::cm_free_factor(m_cholmodFactor, m_cholmod); m_cholmodFactor = 0; } cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - m_cholmodFactor = cholmod_analyze(&A, &m_cholmod); + m_cholmodFactor = internal::cm_analyze(A, m_cholmod); this->m_isInitialized = true; this->m_info = Success; @@ -268,7 +309,7 @@ class CholmodBase : public SparseSolverBase { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod); + internal::cm_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod); // If the factorization failed, minor is the column at which it did. On success minor == n. this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue); @@ -289,19 +330,20 @@ class CholmodBase : public SparseSolverBase EIGEN_UNUSED_VARIABLE(size); eigen_assert(size==b.rows()); - // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref. + // Cholmod needs column-major storage without inner-stride, which corresponds to the default behavior of Ref. Ref > b_ref(b.derived()); cholmod_dense b_cd = viewAsCholmod(b_ref); - cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod); + cholmod_dense* x_cd = internal::cm_solve(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod); if(!x_cd) { this->m_info = NumericalIssue; return; } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) + // NOTE Actually, the copy can be avoided by calling cholmod_solve2 instead of cholmod_solve dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows(),b.cols()); - cholmod_free_dense(&x_cd, &m_cholmod); + internal::cm_free_dense(x_cd, m_cholmod); } /** \internal */ @@ -316,15 +358,16 @@ class CholmodBase : public SparseSolverBase // note: cs stands for Cholmod Sparse Ref > b_ref(b.const_cast_derived()); cholmod_sparse b_cs = viewAsCholmod(b_ref); - cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); + cholmod_sparse* x_cs = internal::cm_spsolve(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod); if(!x_cs) { this->m_info = NumericalIssue; return; } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) + // NOTE cholmod_spsolve in fact just calls the dense solver for blocks of 4 columns at a time (similar to Eigen's sparse solver) dest.derived() = viewAsEigen(*x_cs); - cholmod_free_sparse(&x_cs, &m_cholmod); + internal::cm_free_sparse(x_cs, m_cholmod); } #endif // EIGEN_PARSED_BY_DOXYGEN diff --git a/test/cholmod_support.cpp b/test/cholmod_support.cpp index a7eda28f7..931207334 100644 --- a/test/cholmod_support.cpp +++ b/test/cholmod_support.cpp @@ -12,21 +12,21 @@ #include -template void test_cholmod_T() +template void test_cholmod_ST() { - CholmodDecomposition, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); - CholmodDecomposition, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); - CholmodDecomposition, Lower> g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); - CholmodDecomposition, Upper> g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); - CholmodDecomposition, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); - CholmodDecomposition, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); + CholmodDecomposition g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); + CholmodDecomposition g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); + CholmodDecomposition g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); + CholmodDecomposition g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); + CholmodDecomposition g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); + CholmodDecomposition g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); - CholmodSupernodalLLT, Lower> chol_colmajor_lower; - CholmodSupernodalLLT, Upper> chol_colmajor_upper; - CholmodSimplicialLLT, Lower> llt_colmajor_lower; - CholmodSimplicialLLT, Upper> llt_colmajor_upper; - CholmodSimplicialLDLT, Lower> ldlt_colmajor_lower; - CholmodSimplicialLDLT, Upper> ldlt_colmajor_upper; + CholmodSupernodalLLT chol_colmajor_lower; + CholmodSupernodalLLT chol_colmajor_upper; + CholmodSimplicialLLT llt_colmajor_lower; + CholmodSimplicialLLT llt_colmajor_upper; + CholmodSimplicialLDLT ldlt_colmajor_lower; + CholmodSimplicialLDLT ldlt_colmajor_upper; check_sparse_spd_solving(g_chol_colmajor_lower); check_sparse_spd_solving(g_chol_colmajor_upper); @@ -50,8 +50,20 @@ template void test_cholmod_T() check_sparse_spd_determinant(ldlt_colmajor_upper); } +template void test_cholmod_T() +{ + test_cholmod_ST >(); +} + void test_cholmod_support() { - CALL_SUBTEST_1(test_cholmod_T()); - CALL_SUBTEST_2(test_cholmod_T >()); + CALL_SUBTEST_11( (test_cholmod_T()) ); + CALL_SUBTEST_12( (test_cholmod_T()) ); + CALL_SUBTEST_13( (test_cholmod_T()) ); + CALL_SUBTEST_14( (test_cholmod_T()) ); + CALL_SUBTEST_21( (test_cholmod_T, ColMajor, int >()) ); + CALL_SUBTEST_22( (test_cholmod_T, ColMajor, long>()) ); + // TODO complex row-major matrices do not work at the moment: + // CALL_SUBTEST_23( (test_cholmod_T, RowMajor, int >()) ); + // CALL_SUBTEST_24( (test_cholmod_T, RowMajor, long>()) ); } From 316673bbdeb27688d454480405ca4b7eb0b10a4c Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 14:38:05 +0100 Subject: [PATCH 10/19] Clean-up usage of ExpressionTraits in all/any implementation. --- Eigen/src/Core/BooleanRedux.h | 38 +++++++++++++++++------------------ 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/Eigen/src/Core/BooleanRedux.h b/Eigen/src/Core/BooleanRedux.h index 8409d8749..ed607d5d8 100644 --- a/Eigen/src/Core/BooleanRedux.h +++ b/Eigen/src/Core/BooleanRedux.h @@ -14,56 +14,54 @@ namespace Eigen { namespace internal { -template +template struct all_unroller { - typedef typename Derived::ExpressionTraits Traits; enum { - col = (UnrollCount-1) / Traits::RowsAtCompileTime, - row = (UnrollCount-1) % Traits::RowsAtCompileTime + col = (UnrollCount-1) / Rows, + row = (UnrollCount-1) % Rows }; static inline bool run(const Derived &mat) { - return all_unroller::run(mat) && mat.coeff(row, col); + return all_unroller::run(mat) && mat.coeff(row, col); } }; -template -struct all_unroller +template +struct all_unroller { static inline bool run(const Derived &/*mat*/) { return true; } }; -template -struct all_unroller +template +struct all_unroller { static inline bool run(const Derived &) { return false; } }; -template +template struct any_unroller { - typedef typename Derived::ExpressionTraits Traits; enum { - col = (UnrollCount-1) / Traits::RowsAtCompileTime, - row = (UnrollCount-1) % Traits::RowsAtCompileTime + col = (UnrollCount-1) / Rows, + row = (UnrollCount-1) % Rows }; static inline bool run(const Derived &mat) { - return any_unroller::run(mat) || mat.coeff(row, col); + return any_unroller::run(mat) || mat.coeff(row, col); } }; -template -struct any_unroller +template +struct any_unroller { static inline bool run(const Derived & /*mat*/) { return false; } }; -template -struct any_unroller +template +struct any_unroller { static inline bool run(const Derived &) { return false; } }; @@ -87,7 +85,7 @@ inline bool DenseBase::all() const }; Evaluator evaluator(derived()); if(unroll) - return internal::all_unroller::run(evaluator); + return internal::all_unroller::RowsAtCompileTime>::run(evaluator); else { for(Index j = 0; j < cols(); ++j) @@ -111,7 +109,7 @@ inline bool DenseBase::any() const }; Evaluator evaluator(derived()); if(unroll) - return internal::any_unroller::run(evaluator); + return internal::any_unroller::RowsAtCompileTime>::run(evaluator); else { for(Index j = 0; j < cols(); ++j) From 1c024e5585ba78e79b19d04081ab1d0eef5ab8f7 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Tue, 20 Dec 2016 14:45:44 +0100 Subject: [PATCH 11/19] Added some possible temporaries to .hgignore --- .hgignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.hgignore b/.hgignore index 769a47f1f..dcd9f4431 100644 --- a/.hgignore +++ b/.hgignore @@ -28,7 +28,11 @@ activity.png *.rej log patch +*.patch a a.* lapack/testing lapack/reference +.*project +.settings +Makefile From 5271474b1512370eb09bd8021aaba04e08586310 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 15:51:30 +0100 Subject: [PATCH 12/19] Remove common "noncopyable" base class from evaluator_base to get a chance to get EBO (Empty Base Optimization) Note: we should probbaly get rid of this class and define a macro instead. --- Eigen/src/Core/CoreEvaluators.h | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index 1d14af652..0f05ea76e 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -106,7 +106,7 @@ struct evaluator // ---------- base class for all evaluators ---------- template -struct evaluator_base : public noncopyable +struct evaluator_base { // TODO that's not very nice to have to propagate all these traits. They are currently only needed to handle outer,inner indices. typedef traits ExpressionTraits; @@ -114,6 +114,14 @@ struct evaluator_base : public noncopyable enum { Alignment = 0 }; + // noncopyable: + // Don't make this class inherit noncopyable as this kills EBO (Empty Base Optimization) + // and make complex evaluator much larger than then should do. + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE evaluator_base() {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ~evaluator_base() {} +private: + EIGEN_DEVICE_FUNC evaluator_base(const evaluator_base&); + EIGEN_DEVICE_FUNC const evaluator_base& operator=(const evaluator_base&); }; // -------------------- Matrix and Array -------------------- From 11f55b2979a0044d663151643f8e59d9285cb5b4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 15:55:40 +0100 Subject: [PATCH 13/19] Optimize storage layout of Cwise* and PlainObjectBase evaluator to remove the functor or outer-stride if they are empty. For instance, sizeof("(A-B).cwiseAbs2()") with A,B Vector4f is now 16 bytes, instead of 48 before this optimization. In theory, evaluators should be completely optimized away by the compiler, but this might help in some cases. --- Eigen/src/Core/CoreEvaluators.h | 193 ++++++++++++++++++++------------ 1 file changed, 121 insertions(+), 72 deletions(-) diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index 0f05ea76e..24fc7835b 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -131,6 +131,27 @@ private: // Here we directly specialize evaluator. This is not really a unary expression, and it is, by definition, dense, // so no need for more sophisticated dispatching. +// this helper permits to completely eliminate m_outerStride if it is known at compiletime. +template class plainobjectbase_evaluator_data { +public: + plainobjectbase_evaluator_data(const Scalar* ptr, Index outerStride) : data(ptr) + { + EIGEN_ONLY_USED_FOR_DEBUG(outerStride); + eigen_internal_assert(outerStride==OuterStride); + } + Index outerStride() const { return OuterStride; } + const Scalar *data; +}; + +template class plainobjectbase_evaluator_data { +public: + plainobjectbase_evaluator_data(const Scalar* ptr, Index outerStride) : data(ptr), m_outerStride(outerStride) {} + Index outerStride() const { return m_outerStride; } + const Scalar *data; +protected: + Index m_outerStride; +}; + template struct evaluator > : evaluator_base @@ -149,18 +170,21 @@ struct evaluator > Flags = traits::EvaluatorFlags, Alignment = traits::Alignment }; - + enum { + // We do not need to know the outer stride for vectors + OuterStrideAtCompileTime = IsVectorAtCompileTime ? 0 + : int(IsRowMajor) ? ColsAtCompileTime + : RowsAtCompileTime + }; + EIGEN_DEVICE_FUNC evaluator() - : m_data(0), - m_outerStride(IsVectorAtCompileTime ? 0 - : int(IsRowMajor) ? ColsAtCompileTime - : RowsAtCompileTime) + : m_d(0,OuterStrideAtCompileTime) { EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); } - + EIGEN_DEVICE_FUNC explicit evaluator(const PlainObjectType& m) - : m_data(m.data()), m_outerStride(IsVectorAtCompileTime ? 0 : m.outerStride()) + : m_d(m.data(),IsVectorAtCompileTime ? 0 : m.outerStride()) { EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); } @@ -169,30 +193,30 @@ struct evaluator > CoeffReturnType coeff(Index row, Index col) const { if (IsRowMajor) - return m_data[row * m_outerStride.value() + col]; + return m_d.data[row * m_d.outerStride() + col]; else - return m_data[row + col * m_outerStride.value()]; + return m_d.data[row + col * m_d.outerStride()]; } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { - return m_data[index]; + return m_d.data[index]; } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar& coeffRef(Index row, Index col) { if (IsRowMajor) - return const_cast(m_data)[row * m_outerStride.value() + col]; + return const_cast(m_d.data)[row * m_d.outerStride() + col]; else - return const_cast(m_data)[row + col * m_outerStride.value()]; + return const_cast(m_d.data)[row + col * m_d.outerStride()]; } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar& coeffRef(Index index) { - return const_cast(m_data)[index]; + return const_cast(m_d.data)[index]; } template @@ -200,16 +224,16 @@ struct evaluator > PacketType packet(Index row, Index col) const { if (IsRowMajor) - return ploadt(m_data + row * m_outerStride.value() + col); + return ploadt(m_d.data + row * m_d.outerStride() + col); else - return ploadt(m_data + row + col * m_outerStride.value()); + return ploadt(m_d.data + row + col * m_d.outerStride()); } template EIGEN_STRONG_INLINE PacketType packet(Index index) const { - return ploadt(m_data + index); + return ploadt(m_d.data + index); } template @@ -218,26 +242,22 @@ struct evaluator > { if (IsRowMajor) return pstoret - (const_cast(m_data) + row * m_outerStride.value() + col, x); + (const_cast(m_d.data) + row * m_d.outerStride() + col, x); else return pstoret - (const_cast(m_data) + row + col * m_outerStride.value(), x); + (const_cast(m_d.data) + row + col * m_d.outerStride(), x); } template EIGEN_STRONG_INLINE void writePacket(Index index, const PacketType& x) { - return pstoret(const_cast(m_data) + index, x); + return pstoret(const_cast(m_d.data) + index, x); } protected: - const Scalar *m_data; - // We do not need to know the outer stride for vectors - variable_if_dynamic m_outerStride; + plainobjectbase_evaluator_data m_d; }; template @@ -535,9 +555,7 @@ struct unary_evaluator, IndexBased > }; EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE - explicit unary_evaluator(const XprType& op) - : m_functor(op.functor()), - m_argImpl(op.nestedExpression()) + explicit unary_evaluator(const XprType& op) : m_d(op) { EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); @@ -548,32 +566,43 @@ struct unary_evaluator, IndexBased > EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const { - return m_functor(m_argImpl.coeff(row, col)); + return m_d.func()(m_d.argImpl.coeff(row, col)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { - return m_functor(m_argImpl.coeff(index)); + return m_d.func()(m_d.argImpl.coeff(index)); } template EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const { - return m_functor.packetOp(m_argImpl.template packet(row, col)); + return m_d.func().packetOp(m_d.argImpl.template packet(row, col)); } template EIGEN_STRONG_INLINE PacketType packet(Index index) const { - return m_functor.packetOp(m_argImpl.template packet(index)); + return m_d.func().packetOp(m_d.argImpl.template packet(index)); } protected: - const UnaryOp m_functor; - evaluator m_argImpl; + + // this helper permits to completely eliminate the functor if it is empty + class Data : private UnaryOp + { + public: + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Data(const XprType& xpr) : UnaryOp(xpr.functor()), argImpl(xpr.nestedExpression()) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const UnaryOp& func() const { return static_cast(*this); } + evaluator argImpl; + }; + + Data m_d; }; // -------------------- CwiseTernaryOp -------------------- @@ -617,11 +646,7 @@ struct ternary_evaluator, IndexBased evaluator::Alignment) }; - EIGEN_DEVICE_FUNC explicit ternary_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_arg1Impl(xpr.arg1()), - m_arg2Impl(xpr.arg2()), - m_arg3Impl(xpr.arg3()) + EIGEN_DEVICE_FUNC explicit ternary_evaluator(const XprType& xpr) : m_d(xpr) { EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); @@ -632,38 +657,47 @@ struct ternary_evaluator, IndexBased EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const { - return m_functor(m_arg1Impl.coeff(row, col), m_arg2Impl.coeff(row, col), m_arg3Impl.coeff(row, col)); + return m_d.func()(m_d.arg1Impl.coeff(row, col), m_d.arg2Impl.coeff(row, col), m_d.arg3Impl.coeff(row, col)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { - return m_functor(m_arg1Impl.coeff(index), m_arg2Impl.coeff(index), m_arg3Impl.coeff(index)); + return m_d.func()(m_d.arg1Impl.coeff(index), m_d.arg2Impl.coeff(index), m_d.arg3Impl.coeff(index)); } template EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const { - return m_functor.packetOp(m_arg1Impl.template packet(row, col), - m_arg2Impl.template packet(row, col), - m_arg3Impl.template packet(row, col)); + return m_d.func().packetOp(m_d.arg1Impl.template packet(row, col), + m_d.arg2Impl.template packet(row, col), + m_d.arg3Impl.template packet(row, col)); } template EIGEN_STRONG_INLINE PacketType packet(Index index) const { - return m_functor.packetOp(m_arg1Impl.template packet(index), - m_arg2Impl.template packet(index), - m_arg3Impl.template packet(index)); + return m_d.func().packetOp(m_d.arg1Impl.template packet(index), + m_d.arg2Impl.template packet(index), + m_d.arg3Impl.template packet(index)); } protected: - const TernaryOp m_functor; - evaluator m_arg1Impl; - evaluator m_arg2Impl; - evaluator m_arg3Impl; + // this helper permits to completely eliminate the functor if it is empty + struct Data : private TernaryOp + { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Data(const XprType& xpr) : TernaryOp(xpr.functor()), arg1Impl(xpr.arg1()), arg2Impl(xpr.arg2()), arg3Impl(xpr.arg3()) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const TernaryOp& func() const { return static_cast(*this); } + evaluator arg1Impl; + evaluator arg2Impl; + evaluator arg3Impl; + }; + + Data m_d; }; // -------------------- CwiseBinaryOp -------------------- @@ -704,10 +738,7 @@ struct binary_evaluator, IndexBased, IndexBase Alignment = EIGEN_PLAIN_ENUM_MIN(evaluator::Alignment,evaluator::Alignment) }; - EIGEN_DEVICE_FUNC explicit binary_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()) + EIGEN_DEVICE_FUNC explicit binary_evaluator(const XprType& xpr) : m_d(xpr) { EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); @@ -718,35 +749,45 @@ struct binary_evaluator, IndexBased, IndexBase EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const { - return m_functor(m_lhsImpl.coeff(row, col), m_rhsImpl.coeff(row, col)); + return m_d.func()(m_d.lhsImpl.coeff(row, col), m_d.rhsImpl.coeff(row, col)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { - return m_functor(m_lhsImpl.coeff(index), m_rhsImpl.coeff(index)); + return m_d.func()(m_d.lhsImpl.coeff(index), m_d.rhsImpl.coeff(index)); } template EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const { - return m_functor.packetOp(m_lhsImpl.template packet(row, col), - m_rhsImpl.template packet(row, col)); + return m_d.func().packetOp(m_d.lhsImpl.template packet(row, col), + m_d.rhsImpl.template packet(row, col)); } template EIGEN_STRONG_INLINE PacketType packet(Index index) const { - return m_functor.packetOp(m_lhsImpl.template packet(index), - m_rhsImpl.template packet(index)); + return m_d.func().packetOp(m_d.lhsImpl.template packet(index), + m_d.rhsImpl.template packet(index)); } protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; + + // this helper permits to completely eliminate the functor if it is empty + struct Data : private BinaryOp + { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Data(const XprType& xpr) : BinaryOp(xpr.functor()), lhsImpl(xpr.lhs()), rhsImpl(xpr.rhs()) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const BinaryOp& func() const { return static_cast(*this); } + evaluator lhsImpl; + evaluator rhsImpl; + }; + + Data m_d; }; // -------------------- CwiseUnaryView -------------------- @@ -765,9 +806,7 @@ struct unary_evaluator, IndexBased> Alignment = 0 // FIXME it is not very clear why alignment is necessarily lost... }; - EIGEN_DEVICE_FUNC explicit unary_evaluator(const XprType& op) - : m_unaryOp(op.functor()), - m_argImpl(op.nestedExpression()) + EIGEN_DEVICE_FUNC explicit unary_evaluator(const XprType& op) : m_d(op) { EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); @@ -779,30 +818,40 @@ struct unary_evaluator, IndexBased> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const { - return m_unaryOp(m_argImpl.coeff(row, col)); + return m_d.func()(m_d.argImpl.coeff(row, col)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { - return m_unaryOp(m_argImpl.coeff(index)); + return m_d.func()(m_d.argImpl.coeff(index)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar& coeffRef(Index row, Index col) { - return m_unaryOp(m_argImpl.coeffRef(row, col)); + return m_d.func()(m_d.argImpl.coeffRef(row, col)); } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar& coeffRef(Index index) { - return m_unaryOp(m_argImpl.coeffRef(index)); + return m_d.func()(m_d.argImpl.coeffRef(index)); } protected: - const UnaryOp m_unaryOp; - evaluator m_argImpl; + + // this helper permits to completely eliminate the functor if it is empty + struct Data : private UnaryOp + { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Data(const XprType& xpr) : UnaryOp(xpr.functor()), argImpl(xpr.nestedExpression()) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const UnaryOp& func() const { return static_cast(*this); } + evaluator argImpl; + }; + + Data m_d; }; // -------------------- Map -------------------- From 684cfc762d70e8ab766bc94968d8d5e462c44074 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 16:33:53 +0100 Subject: [PATCH 14/19] Add transpose, adjoint, conjugate methods to SelfAdjointView (useful to write generic code) --- Eigen/src/Core/SelfAdjointView.h | 36 ++++++++++++++++++++++++++++++-- test/product_symm.cpp | 18 ++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 62d4180da..06484ab30 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -45,7 +45,7 @@ struct traits > : traits }; } -// FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ?? + template class SelfAdjointView : public TriangularBase > { @@ -60,10 +60,12 @@ template class SelfAdjointView /** \brief The type of coefficients in this matrix */ typedef typename internal::traits::Scalar Scalar; typedef typename MatrixType::StorageIndex StorageIndex; + typedef typename internal::remove_all::type MatrixConjugateReturnType; enum { Mode = internal::traits::Mode, - Flags = internal::traits::Flags + Flags = internal::traits::Flags, + TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0) }; typedef typename MatrixType::PlainObject PlainObject; @@ -187,6 +189,36 @@ template class SelfAdjointView TriangularView >::type(tmp2); } + typedef SelfAdjointView ConjugateReturnType; + /** \sa MatrixBase::conjugate() const */ + EIGEN_DEVICE_FUNC + inline const ConjugateReturnType conjugate() const + { return ConjugateReturnType(m_matrix.conjugate()); } + + typedef SelfAdjointView AdjointReturnType; + /** \sa MatrixBase::adjoint() const */ + EIGEN_DEVICE_FUNC + inline const AdjointReturnType adjoint() const + { return AdjointReturnType(m_matrix.adjoint()); } + + typedef SelfAdjointView TransposeReturnType; + /** \sa MatrixBase::transpose() */ + EIGEN_DEVICE_FUNC + inline TransposeReturnType transpose() + { + EIGEN_STATIC_ASSERT_LVALUE(MatrixType) + typename MatrixType::TransposeReturnType tmp(m_matrix); + return TransposeReturnType(tmp); + } + + typedef SelfAdjointView ConstTransposeReturnType; + /** \sa MatrixBase::transpose() const */ + EIGEN_DEVICE_FUNC + inline const ConstTransposeReturnType transpose() const + { + return ConstTransposeReturnType(m_matrix.transpose()); + } + /** \returns a const expression of the main diagonal of the matrix \c *this * * This method simply returns the diagonal of the nested expression, thus by-passing the SelfAdjointView decorator. diff --git a/test/product_symm.cpp b/test/product_symm.cpp index 74d7329b1..8c44383f9 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -39,6 +39,24 @@ template void symm(int size = Size, in VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView() * (s2*rhs1), rhs13 = (s1*m1) * (s2*rhs1)); + VERIFY_IS_APPROX(rhs12 = (s1*m2).transpose().template selfadjointView() * (s2*rhs1), + rhs13 = (s1*m1.transpose()) * (s2*rhs1)); + + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView().transpose() * (s2*rhs1), + rhs13 = (s1*m1.transpose()) * (s2*rhs1)); + + VERIFY_IS_APPROX(rhs12 = (s1*m2).conjugate().template selfadjointView() * (s2*rhs1), + rhs13 = (s1*m1).conjugate() * (s2*rhs1)); + + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView().conjugate() * (s2*rhs1), + rhs13 = (s1*m1).conjugate() * (s2*rhs1)); + + VERIFY_IS_APPROX(rhs12 = (s1*m2).adjoint().template selfadjointView() * (s2*rhs1), + rhs13 = (s1*m1).adjoint() * (s2*rhs1)); + + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView().adjoint() * (s2*rhs1), + rhs13 = (s1*m1).adjoint() * (s2*rhs1)); + m2 = m1.template triangularView(); rhs12.setRandom(); rhs13 = rhs12; m3 = m2.template selfadjointView(); VERIFY_IS_EQUAL(m1, m3); From e2f4ee1c2b0933bd60d169dc328b321f2db40605 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 21:56:21 +0100 Subject: [PATCH 15/19] Speed up parsing of sparse Market file. --- unsupported/Eigen/src/SparseExtra/MarketIO.h | 80 +++++++++++--------- unsupported/test/sparse_extra.cpp | 23 ++++++ 2 files changed, 66 insertions(+), 37 deletions(-) diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index cdc14f86e..9dbd0488d 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -12,38 +12,38 @@ #define EIGEN_SPARSE_MARKET_IO_H #include +#include namespace Eigen { namespace internal { - template - inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, Scalar& value) + template + inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, Scalar& value) { - line >> i >> j >> value; - i--; - j--; - if(i>=0 && j>=0 && i> i >> j >> value; } - template - inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, std::complex& value) + + template<> inline void GetMarketLine (const char* line, int& i, int& j, float& value) + { std::sscanf(line, "%d %d %g", &i, &j, &value); } + + template<> inline void GetMarketLine (const char* line, int& i, int& j, double& value) + { std::sscanf(line, "%d %d %lg", &i, &j, &value); } + + template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex& value) + { std::sscanf(line, "%d %d %g %g", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); } + + template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex& value) + { std::sscanf(line, "%d %d %lg %lg", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); } + + template + inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, std::complex& value) { + std::stringstream sline(line); Scalar valR, valI; - line >> i >> j >> valR >> valI; - i--; - j--; - if(i>=0 && j>=0 && i(valR, valI); - return true; - } - else - return false; + sline >> i >> j >> valR >> valI; + value = std::complex(valR,valI); } template @@ -81,13 +81,13 @@ namespace internal } } - template - inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out) + template + inline void PutMatrixElt(Scalar value, StorageIndex row, StorageIndex col, std::ofstream& out) { out << row << " "<< col << " " << value << "\n"; } - template - inline void PutMatrixElt(std::complex value, int row, int col, std::ofstream& out) + template + inline void PutMatrixElt(std::complex value, StorageIndex row, StorageIndex col, std::ofstream& out) { out << row << " " << col << " " << value.real() << " " << value.imag() << "\n"; } @@ -133,17 +133,20 @@ template bool loadMarket(SparseMatrixType& mat, const std::string& filename) { typedef typename SparseMatrixType::Scalar Scalar; - typedef typename SparseMatrixType::Index Index; + typedef typename SparseMatrixType::StorageIndex StorageIndex; std::ifstream input(filename.c_str(),std::ios::in); if(!input) return false; + + char rdbuffer[4096]; + input.rdbuf()->pubsetbuf(rdbuffer, 4096); const int maxBuffersize = 2048; char buffer[maxBuffersize]; bool readsizes = false; - typedef Triplet T; + typedef Triplet T; std::vector elements; Index M(-1), N(-1), NNZ(-1); @@ -154,33 +157,36 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) //NOTE An appropriate test should be done on the header to get the symmetry if(buffer[0]=='%') continue; - - std::stringstream line(buffer); - + if(!readsizes) { + std::stringstream line(buffer); line >> M >> N >> NNZ; if(M > 0 && N > 0 && NNZ > 0) { readsizes = true; - //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n"; mat.resize(M,N); mat.reserve(NNZ); } } else { - Index i(-1), j(-1); + StorageIndex i(-1), j(-1); Scalar value; - if( internal::GetMarketLine(line, M, N, i, j, value) ) + internal::GetMarketLine(buffer, i, j, value); + + i--; + j--; + if(i>=0 && j>=0 && i void sparse_extra(const SparseMatrixType& re } +template +void check_marketio() +{ + typedef Matrix DenseMatrix; + Index rows = internal::random(1,100); + Index cols = internal::random(1,100); + SparseMatrixType m1, m2; + m1 = DenseMatrix::Random(rows, cols).sparseView(); + saveMarket(m1, "sparse_extra.mtx"); + loadMarket(m2, "sparse_extra.mtx"); + VERIFY_IS_APPROX(m1,m2); +} + void test_sparse_extra() { for(int i = 0; i < g_repeat; i++) { @@ -143,5 +156,15 @@ void test_sparse_extra() CALL_SUBTEST_3( (sparse_product >()) ); CALL_SUBTEST_3( (sparse_product >()) ); + + CALL_SUBTEST_4( (check_marketio >()) ); + CALL_SUBTEST_4( (check_marketio >()) ); + CALL_SUBTEST_4( (check_marketio,ColMajor,int> >()) ); + CALL_SUBTEST_4( (check_marketio,ColMajor,int> >()) ); + CALL_SUBTEST_4( (check_marketio >()) ); + CALL_SUBTEST_4( (check_marketio >()) ); + CALL_SUBTEST_4( (check_marketio,ColMajor,long int> >()) ); + CALL_SUBTEST_4( (check_marketio,ColMajor,long int> >()) ); + TEST_SET_BUT_UNUSED_VARIABLE(s); } } From e8d6862f14445c28a299ce087bbbffdf5af0ebf6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 22:10:33 +0100 Subject: [PATCH 16/19] Properly adjust precision when saving to Market format. --- unsupported/Eigen/src/SparseExtra/MarketIO.h | 9 +++++---- unsupported/test/sparse_extra.cpp | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index 9dbd0488d..fc70a24d8 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -231,12 +231,13 @@ template bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0) { typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::RealScalar RealScalar; std::ofstream out(filename.c_str(),std::ios::out); if(!out) return false; out.flags(std::ios_base::scientific); - out.precision(64); + out.precision(std::numeric_limits::digits10 + 2); std::string header; internal::putMarketHeader(header, sym); out << header << std::endl; @@ -247,7 +248,6 @@ bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sy { ++ count; internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); - // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n"; } out.close(); return true; @@ -256,13 +256,14 @@ bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sy template bool saveMarketVector (const VectorType& vec, const std::string& filename) { - typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::RealScalar RealScalar; std::ofstream out(filename.c_str(),std::ios::out); if(!out) return false; out.flags(std::ios_base::scientific); - out.precision(64); + out.precision(std::numeric_limits::digits10 + 2); if(internal::is_same >::value || internal::is_same >::value) out << "%%MatrixMarket matrix array complex general\n"; else diff --git a/unsupported/test/sparse_extra.cpp b/unsupported/test/sparse_extra.cpp index ceec3d8a8..4f6723d6d 100644 --- a/unsupported/test/sparse_extra.cpp +++ b/unsupported/test/sparse_extra.cpp @@ -139,7 +139,7 @@ void check_marketio() m1 = DenseMatrix::Random(rows, cols).sparseView(); saveMarket(m1, "sparse_extra.mtx"); loadMarket(m2, "sparse_extra.mtx"); - VERIFY_IS_APPROX(m1,m2); + VERIFY_IS_EQUAL(DenseMatrix(m1),DenseMatrix(m2)); } void test_sparse_extra() From 94e8d8902f882058ff9912238e5934ab5632217d Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 22:17:01 +0100 Subject: [PATCH 17/19] Fix bug #1367: compilation fix for gcc 4.1! --- Eigen/src/Core/util/Macros.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 4ce427478..56b7f0dbb 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -501,10 +501,11 @@ // attribute to maximize inlining. This should only be used when really necessary: in particular, // it uses __attribute__((always_inline)) on GCC, which most of the time is useless and can severely harm compile times. // FIXME with the always_inline attribute, -// gcc 3.4.x reports the following compilation error: +// gcc 3.4.x and 4.1 reports the following compilation error: // Eval.h:91: sorry, unimplemented: inlining failed in call to 'const Eigen::Eval Eigen::MatrixBase::eval() const' // : function body not available -#if EIGEN_GNUC_AT_LEAST(4,0) +// See also bug 1367 +#if EIGEN_GNUC_AT_LEAST(4,2) #define EIGEN_ALWAYS_INLINE __attribute__((always_inline)) inline #else #define EIGEN_ALWAYS_INLINE EIGEN_STRONG_INLINE From 2b3fc981b8603820982429520d4362d6b9607bda Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 22:52:27 +0100 Subject: [PATCH 18/19] bug #1362: workaround constant conditional warning produced by MSVC --- Eigen/src/Core/util/Macros.h | 8 ++++++++ Eigen/src/Geometry/AlignedBox.h | 2 +- Eigen/src/Geometry/Transform.h | 14 +++++++------- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 56b7f0dbb..95960b448 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -628,6 +628,14 @@ namespace Eigen { #endif +#if EIGEN_COMP_MSVC + // NOTE MSVC often gives C4127 warnings with compiletime if statements. See bug 1362. + // This workaround is ugly, but it does the job. +# define EIGEN_CONST_CONDITIONAL(cond) (void)0, cond +#else +# define EIGEN_CONST_CONDITIONAL(cond) cond +#endif + //------------------------------------------------------------------------------------------ // Static and dynamic alignment control // diff --git a/Eigen/src/Geometry/AlignedBox.h b/Eigen/src/Geometry/AlignedBox.h index 066eae4f9..c902d8f0a 100644 --- a/Eigen/src/Geometry/AlignedBox.h +++ b/Eigen/src/Geometry/AlignedBox.h @@ -63,7 +63,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) /** Default constructor initializing a null box. */ EIGEN_DEVICE_FUNC inline AlignedBox() - { if (AmbientDimAtCompileTime!=Dynamic) setEmpty(); } + { if (EIGEN_CONST_CONDITIONAL(AmbientDimAtCompileTime!=Dynamic)) setEmpty(); } /** Constructs a null box with \a _dim the dimension of the ambient space. */ EIGEN_DEVICE_FUNC inline explicit AlignedBox(Index _dim) : m_min(_dim), m_max(_dim) diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 3f31ee45d..2d36dfadf 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -335,7 +335,7 @@ public: OtherModeIsAffineCompact = OtherMode == int(AffineCompact) }; - if(ModeIsAffineCompact == OtherModeIsAffineCompact) + if(EIGEN_CONST_CONDITIONAL(ModeIsAffineCompact == OtherModeIsAffineCompact)) { // We need the block expression because the code is compiled for all // combinations of transformations and will trigger a compile time error @@ -343,7 +343,7 @@ public: m_matrix.template block(0,0) = other.matrix().template block(0,0); makeAffine(); } - else if(OtherModeIsAffineCompact) + else if(EIGEN_CONST_CONDITIONAL(OtherModeIsAffineCompact)) { typedef typename Transform::MatrixType OtherMatrixType; internal::transform_construct_from_matrix::run(this, other.matrix()); @@ -481,7 +481,7 @@ public: TransformTimeDiagonalReturnType res; res.linear().noalias() = a*b.linear(); res.translation().noalias() = a*b.translation(); - if (Mode!=int(AffineCompact)) + if (EIGEN_CONST_CONDITIONAL(Mode!=int(AffineCompact))) res.matrix().row(Dim) = b.matrix().row(Dim); return res; } @@ -755,7 +755,7 @@ template Transform& Transform::operator=(const QMatrix& other) { EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) - if (Mode == int(AffineCompact)) + if (EIGEN_CONST_CONDITIONAL(Mode == int(AffineCompact))) m_matrix << other.m11(), other.m21(), other.dx(), other.m12(), other.m22(), other.dy(); else @@ -801,7 +801,7 @@ Transform& Transform::operator { check_template_params(); EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) - if (Mode == int(AffineCompact)) + if (EIGEN_CONST_CONDITIONAL(Mode == int(AffineCompact))) m_matrix << other.m11(), other.m21(), other.dx(), other.m12(), other.m22(), other.dy(); else @@ -819,7 +819,7 @@ template QTransform Transform::toQTransform(void) const { EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) - if (Mode == int(AffineCompact)) + if (EIGEN_CONST_CONDITIONAL(Mode == int(AffineCompact))) return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(0,2), m_matrix.coeff(1,2)); @@ -912,7 +912,7 @@ EIGEN_DEVICE_FUNC Transform& Transform::pretranslate(const MatrixBase &other) { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,int(Dim)) - if(int(Mode)==int(Projective)) + if(EIGEN_CONST_CONDITIONAL(int(Mode)==int(Projective))) affine() += other * m_matrix.row(Dim); else translation() += other; From f2f9df8aa57d7a303eb113c251245e315f2ad2b7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 20 Dec 2016 22:53:19 +0100 Subject: [PATCH 19/19] Remove MSVC warning 4127 - conditional expression is constant from the disabled list as we now have a local workaround. --- Eigen/src/Core/util/DisableStupidWarnings.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Eigen/src/Core/util/DisableStupidWarnings.h b/Eigen/src/Core/util/DisableStupidWarnings.h index 21f80d86b..4431f2fc4 100755 --- a/Eigen/src/Core/util/DisableStupidWarnings.h +++ b/Eigen/src/Core/util/DisableStupidWarnings.h @@ -4,7 +4,6 @@ #ifdef _MSC_VER // 4100 - unreferenced formal parameter (occurred e.g. in aligned_allocator::destroy(pointer p)) // 4101 - unreferenced local variable - // 4127 - conditional expression is constant // 4181 - qualifier applied to reference type ignored // 4211 - nonstandard extension used : redefined extern to static // 4244 - 'argument' : conversion from 'type1' to 'type2', possible loss of data @@ -20,7 +19,7 @@ #ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS #pragma warning( push ) #endif - #pragma warning( disable : 4100 4101 4127 4181 4211 4244 4273 4324 4503 4512 4522 4700 4714 4717 4800) + #pragma warning( disable : 4100 4101 4181 4211 4244 4273 4324 4503 4512 4522 4700 4714 4717 4800) #elif defined __INTEL_COMPILER // 2196 - routine is both "inline" and "noinline" ("noinline" assumed)