aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Core6
-rw-r--r--Eigen/src/Core/PacketMath.h205
-rw-r--r--Eigen/src/Core/PacketMath_Altivec.h113
-rw-r--r--Eigen/src/Core/PacketMath_SSE.h154
-rw-r--r--Eigen/src/Core/ProductWIP.h260
-rw-r--r--Eigen/src/Core/util/Macros.h6
-rw-r--r--Eigen/src/Core/util/Meta.h3
7 files changed, 527 insertions, 220 deletions
diff --git a/Eigen/Core b/Eigen/Core
index cfe7a1059..999bb96f2 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -39,6 +39,12 @@ namespace Eigen {
#include "src/Core/NumTraits.h"
#include "src/Core/MathFunctions.h"
#include "src/Core/PacketMath.h"
+#if defined EIGEN_VECTORIZE_SSE
+#include "src/Core/PacketMath_SSE.h"
+#elif defined EIGEN_VECTORIZE_ALTIVEC
+#include "src/Core/PacketMath_Altivec.h"
+#endif
+
#include "src/Core/Functors.h"
#include "src/Core/MatrixBase.h"
#include "src/Core/Coeffs.h"
diff --git a/Eigen/src/Core/PacketMath.h b/Eigen/src/Core/PacketMath.h
index cfa19eb6a..59c143ede 100644
--- a/Eigen/src/Core/PacketMath.h
+++ b/Eigen/src/Core/PacketMath.h
@@ -26,7 +26,7 @@
#define EIGEN_PACKET_MATH_H
#ifndef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16
+#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16
#endif
// Default implementation for types not supported by the vectorization.
@@ -35,211 +35,46 @@
// called, TODO so sould we raise an assertion or not ?
/** \internal \returns a + b (coeff-wise) */
template <typename Scalar> inline Scalar ei_padd(const Scalar& a, const Scalar& b) { return a + b; }
+
/** \internal \returns a - b (coeff-wise) */
template <typename Scalar> inline Scalar ei_psub(const Scalar& a, const Scalar& b) { return a - b; }
+
/** \internal \returns a * b (coeff-wise) */
template <typename Scalar> inline Scalar ei_pmul(const Scalar& a, const Scalar& b) { return a * b; }
+
/** \internal \returns a * b - c (coeff-wise) */
template <typename Scalar> inline Scalar ei_pmadd(const Scalar& a, const Scalar& b, const Scalar& c)
{ return ei_padd(ei_pmul(a, b),c); }
+
/** \internal \returns the min of \a a and \a b (coeff-wise) */
template <typename Scalar> inline Scalar ei_pmin(const Scalar& a, const Scalar& b) { return std::min(a,b); }
+
/** \internal \returns the max of \a a and \a b (coeff-wise) */
template <typename Scalar> inline Scalar ei_pmax(const Scalar& a, const Scalar& b) { return std::max(a,b); }
+
/** \internal \returns a packet version of \a *from, from must be 16 bytes aligned */
template <typename Scalar> inline Scalar ei_pload(const Scalar* from) { return *from; }
+
+/** \internal \returns a packet version of \a *from, (un-aligned load) */
+template <typename Scalar> inline Scalar ei_ploadu(const Scalar* from) { return *from; }
+
/** \internal \returns a packet with constant coefficients \a a, e.g.: (a,a,a,a) */
template <typename Scalar> inline Scalar ei_pset1(const Scalar& a) { return a; }
+
/** \internal copy the packet \a from to \a *to, \a to must be 16 bytes aligned */
template <typename Scalar> inline void ei_pstore(Scalar* to, const Scalar& from) { (*to) = from; }
-/** \internal \returns the first element of a packet */
-template <typename Scalar> inline Scalar ei_pfirst(const Scalar& a) { return a; }
-/** \internal \returns a packet where the element i contains the sum of the packet of \a vec[i] */
-template <typename Scalar> inline Scalar ei_predux(const Scalar vecs[1]) { return vecs[0]; }
-/** \internal \returns the sum of the elements of \a a*/
-template <typename Scalar> inline Scalar ei_predux(const Scalar& a) { return a; }
-
-#ifdef EIGEN_VECTORIZE_SSE
-
-#ifdef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-#undef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16
-#endif
-template<> struct ei_packet_traits<float> { typedef __m128 type; enum {size=4}; };
-template<> struct ei_packet_traits<double> { typedef __m128d type; enum {size=2}; };
-template<> struct ei_packet_traits<int> { typedef __m128i type; enum {size=4}; };
-
-template<> inline __m128 ei_padd(const __m128& a, const __m128& b) { return _mm_add_ps(a,b); }
-template<> inline __m128d ei_padd(const __m128d& a, const __m128d& b) { return _mm_add_pd(a,b); }
-template<> inline __m128i ei_padd(const __m128i& a, const __m128i& b) { return _mm_add_epi32(a,b); }
-
-template<> inline __m128 ei_psub(const __m128& a, const __m128& b) { return _mm_sub_ps(a,b); }
-template<> inline __m128d ei_psub(const __m128d& a, const __m128d& b) { return _mm_sub_pd(a,b); }
-template<> inline __m128i ei_psub(const __m128i& a, const __m128i& b) { return _mm_sub_epi32(a,b); }
-
-template<> inline __m128 ei_pmul(const __m128& a, const __m128& b) { return _mm_mul_ps(a,b); }
-template<> inline __m128d ei_pmul(const __m128d& a, const __m128d& b) { return _mm_mul_pd(a,b); }
-template<> inline __m128i ei_pmul(const __m128i& a, const __m128i& b)
-{
- return _mm_or_si128(
- _mm_and_si128(
- _mm_mul_epu32(a,b),
- _mm_setr_epi32(0xffffffff,0,0xffffffff,0)),
- _mm_slli_si128(
- _mm_and_si128(
- _mm_mul_epu32(_mm_srli_si128(a,4),_mm_srli_si128(b,4)),
- _mm_setr_epi32(0xffffffff,0,0xffffffff,0)), 4));
-}
-
-// for some weird raisons, it has to be overloaded for packet integer
-template<> inline __m128i ei_pmadd(const __m128i& a, const __m128i& b, const __m128i& c) { return ei_padd(ei_pmul(a,b), c); }
-
-template<> inline __m128 ei_pmin(const __m128& a, const __m128& b) { return _mm_min_ps(a,b); }
-template<> inline __m128d ei_pmin(const __m128d& a, const __m128d& b) { return _mm_min_pd(a,b); }
-// FIXME this vectorized min operator is likely to be slower than the standard one
-template<> inline __m128i ei_pmin(const __m128i& a, const __m128i& b)
-{
- __m128i mask = _mm_cmplt_epi32(a,b);
- return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
-}
-
-template<> inline __m128 ei_pmax(const __m128& a, const __m128& b) { return _mm_max_ps(a,b); }
-template<> inline __m128d ei_pmax(const __m128d& a, const __m128d& b) { return _mm_max_pd(a,b); }
-// FIXME this vectorized max operator is likely to be slower than the standard one
-template<> inline __m128i ei_pmax(const __m128i& a, const __m128i& b)
-{
- __m128i mask = _mm_cmpgt_epi32(a,b);
- return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
-}
-
-inline __m128 ei_pload(const float* from) { return _mm_load_ps(from); }
-inline __m128d ei_pload(const double* from) { return _mm_load_pd(from); }
-inline __m128i ei_pload(const int* from) { return _mm_load_si128(reinterpret_cast<const __m128i*>(from)); }
-
-inline __m128 ei_ploadu(const float* from) { return _mm_loadu_ps(from); }
-inline __m128d ei_ploadu(const double* from) { return _mm_loadu_pd(from); }
-inline __m128i ei_ploadu(const int* from) { return _mm_loadu_si128(reinterpret_cast<const __m128i*>(from)); }
-
-inline __m128 ei_pset1(const float& from) { return _mm_set1_ps(from); }
-inline __m128d ei_pset1(const double& from) { return _mm_set1_pd(from); }
-inline __m128i ei_pset1(const int& from) { return _mm_set1_epi32(from); }
-
-inline void ei_pstore(float* to, const __m128& from) { _mm_store_ps(to, from); }
-inline void ei_pstore(double* to, const __m128d& from) { _mm_store_pd(to, from); }
-inline void ei_pstore(int* to, const __m128i& from) { _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
-
-inline void ei_pstoreu(float* to, const __m128& from) { _mm_storeu_ps(to, from); }
-inline void ei_pstoreu(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
-inline void ei_pstoreu(int* to, const __m128i& from) { _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
-
-inline float ei_pfirst(const __m128& a) { return _mm_cvtss_f32(a); }
-inline double ei_pfirst(const __m128d& a) { return _mm_cvtsd_f64(a); }
-inline int ei_pfirst(const __m128i& a) { return _mm_cvtsi128_si32(a); }
-
-#ifdef __SSE3__
-// TODO implement SSE2 versions as well as integer versions
-inline __m128 ei_predux(const __m128* vecs)
-{
- return _mm_hadd_ps(_mm_hadd_ps(vecs[0], vecs[1]),_mm_hadd_ps(vecs[2], vecs[3]));
-}
-inline __m128d ei_predux(const __m128d* vecs)
-{
- return _mm_hadd_pd(vecs[0], vecs[1]);
-}
-
-inline float ei_predux(const __m128& a)
-{
- __m128 tmp0 = _mm_hadd_ps(a,a);
- return ei_pfirst(_mm_hadd_ps(tmp0, tmp0));
-}
-
-inline double ei_predux(const __m128d& a) { return ei_pfirst(_mm_hadd_pd(a, a)); }
-#endif
+/** \internal copy the packet \a from to \a *to, (un-aligned store) */
+template <typename Scalar> inline void ei_pstoreu(Scalar* to, const Scalar& from) { (*to) = from; }
-#elif defined(EIGEN_VECTORIZE_ALTIVEC)
+/** \internal \returns the first element of a packet */
+template <typename Scalar> inline Scalar ei_pfirst(const Scalar& a) { return a; }
-#ifdef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-#undef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 4
-#endif
+/** \internal \returns a packet where the element i contains the sum of the packet of \a vec[i] */
+// template <typename Scalar> inline Scalar ei_predux(const Scalar* vecs) { return vecs[0]; }
-static const vector int v0i = vec_splat_u32(0);
-static const vector int v16i_ = vec_splat_u32(-16);
-static const vector float v0f = (vector float) v0i;
-
-template<> struct ei_packet_traits<float> { typedef vector float type; enum {size=4}; };
-template<> struct ei_packet_traits<int> { typedef vector int type; enum {size=4}; };
-
-inline vector float ei_padd(const vector float a, const vector float b) { return vec_add(a,b); }
-inline vector int ei_padd(const vector int a, const vector int b) { return vec_add(a,b); }
-
-inline vector float ei_psub(const vector float a, const vector float b) { return vec_sub(a,b); }
-inline vector int ei_psub(const vector int a, const vector int b) { return vec_sub(a,b); }
-
-inline vector float ei_pmul(const vector float a, const vector float b) { return vec_madd(a,b, v0f); }
-inline vector int ei_pmul(const vector int a, const vector int b)
-{
- // Taken from http://
-
- //Set up constants
- vector int bswap, lowProduct, highProduct;
-
- //Do real work
- bswap = vec_rl( (vector unsigned int)b, (vector unsigned int)v16i_ );
- lowProduct = vec_mulo( (vector short)a,(vector short)b );
- highProduct = vec_msum((vector short)a,(vector short)bswap, v0i);
- highProduct = vec_sl( (vector unsigned int)highProduct, (vector unsigned int)v16i_ );
- return vec_add( lowProduct, highProduct );
-}
-
-inline vector float ei_pmadd(const vector float a, const vector float b, const vector float c) { return vec_madd(a, b, c); }
-
-inline vector float ei_pmin(const vector float a, const vector float b) { return vec_min(a,b); }
-inline vector int ei_pmin(const vector int a, const vector int b) { return vec_min(a,b); }
-
-inline vector float ei_pmax(const vector float a, const vector float b) { return vec_max(a,b); }
-inline vector int ei_pmax(const vector int a, const vector int b) { return vec_max(a,b); }
-
-inline vector float ei_pload(const float* from) { return vec_ld(0, from); }
-inline vector int ei_pload(const int* from) { return vec_ld(0, from); }
-
-inline vector float ei_pset1(const float& from)
-{
- static float __attribute__(aligned(16)) af[4];
- af[0] = from;
- vector float vc = vec_ld(0, af);
- vc = vec_splat(vc, 0);
- return vc;
-}
-
-inline vector int ei_pset1(const int& from)
-{
- static int __attribute__(aligned(16)) ai[4];
- ai[0] = from;
- vector int vc = vec_ld(0, ai);
- vc = vec_splat(vc, 0);
- return vc;
-}
-
-inline void ei_pstore(float* to, const vector float from) { vec_st(from, 0, to); }
-inline void ei_pstore(int* to, const vector int from) { vec_st(from, 0, to); }
-
-inline float ei_pfirst(const vector float a)
-{
- static float __attribute__(aligned(16)) af[4];
- vec_st(a, 0, af);
- return af[0];
-}
-
-inline int ei_pfirst(const vector int a)
-{
- static int __attribute__(aligned(16)) ai[4];
- vec_st(a, 0, ai);
- return ai[0];
-}
-
-#endif // EIGEN_VECTORIZE_ALTIVEC & SSE
+/** \internal \returns the sum of the elements of \a a*/
+// template <typename Scalar> inline Scalar ei_predux(const Scalar& a) { return a; }
#endif // EIGEN_PACKET_MATH_H
diff --git a/Eigen/src/Core/PacketMath_Altivec.h b/Eigen/src/Core/PacketMath_Altivec.h
new file mode 100644
index 000000000..eb702af8c
--- /dev/null
+++ b/Eigen/src/Core/PacketMath_Altivec.h
@@ -0,0 +1,113 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_PACKET_MATH_ALTIVEC_H
+#define EIGEN_PACKET_MATH_ALTIVEC_H
+
+#ifndef EIGEN_VECTORIZE_ALTIVEC
+#error include PacketMath_Altivec without EIGEN_VECTORIZE_ALTIVEC
+#endif
+
+#ifdef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+#undef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 4
+#endif
+
+static const vector int v0i = vec_splat_u32(0);
+static const vector int v16i_ = vec_splat_u32(-16);
+static const vector float v0f = (vector float) v0i;
+
+template<> struct ei_packet_traits<float> { typedef vector float type; enum {size=4}; };
+template<> struct ei_packet_traits<int> { typedef vector int type; enum {size=4}; };
+
+inline vector float ei_padd(const vector float a, const vector float b) { return vec_add(a,b); }
+inline vector int ei_padd(const vector int a, const vector int b) { return vec_add(a,b); }
+
+inline vector float ei_psub(const vector float a, const vector float b) { return vec_sub(a,b); }
+inline vector int ei_psub(const vector int a, const vector int b) { return vec_sub(a,b); }
+
+inline vector float ei_pmul(const vector float a, const vector float b) { return vec_madd(a,b, v0f); }
+inline vector int ei_pmul(const vector int a, const vector int b)
+{
+ // Taken from http://
+
+ //Set up constants
+ vector int bswap, lowProduct, highProduct;
+
+ //Do real work
+ bswap = vec_rl( (vector unsigned int)b, (vector unsigned int)v16i_ );
+ lowProduct = vec_mulo( (vector short)a,(vector short)b );
+ highProduct = vec_msum((vector short)a,(vector short)bswap, v0i);
+ highProduct = vec_sl( (vector unsigned int)highProduct, (vector unsigned int)v16i_ );
+ return vec_add( lowProduct, highProduct );
+}
+
+inline vector float ei_pmadd(const vector float a, const vector float b, const vector float c) { return vec_madd(a, b, c); }
+
+inline vector float ei_pmin(const vector float a, const vector float b) { return vec_min(a,b); }
+inline vector int ei_pmin(const vector int a, const vector int b) { return vec_min(a,b); }
+
+inline vector float ei_pmax(const vector float a, const vector float b) { return vec_max(a,b); }
+inline vector int ei_pmax(const vector int a, const vector int b) { return vec_max(a,b); }
+
+inline vector float ei_pload(const float* from) { return vec_ld(0, from); }
+inline vector int ei_pload(const int* from) { return vec_ld(0, from); }
+
+inline vector float ei_pset1(const float& from)
+{
+ static float __attribute__(aligned(16)) af[4];
+ af[0] = from;
+ vector float vc = vec_ld(0, af);
+ vc = vec_splat(vc, 0);
+ return vc;
+}
+
+inline vector int ei_pset1(const int& from)
+{
+ static int __attribute__(aligned(16)) ai[4];
+ ai[0] = from;
+ vector int vc = vec_ld(0, ai);
+ vc = vec_splat(vc, 0);
+ return vc;
+}
+
+inline void ei_pstore(float* to, const vector float from) { vec_st(from, 0, to); }
+inline void ei_pstore(int* to, const vector int from) { vec_st(from, 0, to); }
+
+inline float ei_pfirst(const vector float a)
+{
+ static float __attribute__(aligned(16)) af[4];
+ vec_st(a, 0, af);
+ return af[0];
+}
+
+inline int ei_pfirst(const vector int a)
+{
+ static int __attribute__(aligned(16)) ai[4];
+ vec_st(a, 0, ai);
+ return ai[0];
+}
+
+#endif // EIGEN_PACKET_MATH_ALTIVEC_H
+
diff --git a/Eigen/src/Core/PacketMath_SSE.h b/Eigen/src/Core/PacketMath_SSE.h
new file mode 100644
index 000000000..1872affea
--- /dev/null
+++ b/Eigen/src/Core/PacketMath_SSE.h
@@ -0,0 +1,154 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_PACKET_MATH_SSE_H
+#define EIGEN_PACKET_MATH_SSE_H
+
+#ifndef EIGEN_VECTORIZE_SSE
+#error include PacketMath_SSE without EIGEN_VECTORIZE_SSE
+#endif
+
+template<> struct ei_packet_traits<float> { typedef __m128 type; enum {size=4}; };
+template<> struct ei_packet_traits<double> { typedef __m128d type; enum {size=2}; };
+template<> struct ei_packet_traits<int> { typedef __m128i type; enum {size=4}; };
+
+template<> inline __m128 ei_padd(const __m128& a, const __m128& b) { return _mm_add_ps(a,b); }
+template<> inline __m128d ei_padd(const __m128d& a, const __m128d& b) { return _mm_add_pd(a,b); }
+template<> inline __m128i ei_padd(const __m128i& a, const __m128i& b) { return _mm_add_epi32(a,b); }
+
+template<> inline __m128 ei_psub(const __m128& a, const __m128& b) { return _mm_sub_ps(a,b); }
+template<> inline __m128d ei_psub(const __m128d& a, const __m128d& b) { return _mm_sub_pd(a,b); }
+template<> inline __m128i ei_psub(const __m128i& a, const __m128i& b) { return _mm_sub_epi32(a,b); }
+
+template<> inline __m128 ei_pmul(const __m128& a, const __m128& b) { return _mm_mul_ps(a,b); }
+template<> inline __m128d ei_pmul(const __m128d& a, const __m128d& b) { return _mm_mul_pd(a,b); }
+template<> inline __m128i ei_pmul(const __m128i& a, const __m128i& b)
+{
+ return _mm_or_si128(
+ _mm_and_si128(
+ _mm_mul_epu32(a,b),
+ _mm_setr_epi32(0xffffffff,0,0xffffffff,0)),
+ _mm_slli_si128(
+ _mm_and_si128(
+ _mm_mul_epu32(_mm_srli_si128(a,4),_mm_srli_si128(b,4)),
+ _mm_setr_epi32(0xffffffff,0,0xffffffff,0)), 4));
+}
+
+// for some weird raisons, it has to be overloaded for packet integer
+template<> inline __m128i ei_pmadd(const __m128i& a, const __m128i& b, const __m128i& c) { return ei_padd(ei_pmul(a,b), c); }
+
+template<> inline __m128 ei_pmin(const __m128& a, const __m128& b) { return _mm_min_ps(a,b); }
+template<> inline __m128d ei_pmin(const __m128d& a, const __m128d& b) { return _mm_min_pd(a,b); }
+// FIXME this vectorized min operator is likely to be slower than the standard one
+template<> inline __m128i ei_pmin(const __m128i& a, const __m128i& b)
+{
+ __m128i mask = _mm_cmplt_epi32(a,b);
+ return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
+}
+
+template<> inline __m128 ei_pmax(const __m128& a, const __m128& b) { return _mm_max_ps(a,b); }
+template<> inline __m128d ei_pmax(const __m128d& a, const __m128d& b) { return _mm_max_pd(a,b); }
+// FIXME this vectorized max operator is likely to be slower than the standard one
+template<> inline __m128i ei_pmax(const __m128i& a, const __m128i& b)
+{
+ __m128i mask = _mm_cmpgt_epi32(a,b);
+ return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
+}
+
+inline __m128 ei_pload(const float* from) { return _mm_load_ps(from); }
+inline __m128d ei_pload(const double* from) { return _mm_load_pd(from); }
+inline __m128i ei_pload(const int* from) { return _mm_load_si128(reinterpret_cast<const __m128i*>(from)); }
+
+inline __m128 ei_ploadu(const float* from) { return _mm_loadu_ps(from); }
+inline __m128d ei_ploadu(const double* from) { return _mm_loadu_pd(from); }
+inline __m128i ei_ploadu(const int* from) { return _mm_loadu_si128(reinterpret_cast<const __m128i*>(from)); }
+
+inline __m128 ei_pset1(const float& from) { return _mm_set1_ps(from); }
+inline __m128d ei_pset1(const double& from) { return _mm_set1_pd(from); }
+inline __m128i ei_pset1(const int& from) { return _mm_set1_epi32(from); }
+
+inline void ei_pstore(float* to, const __m128& from) { _mm_store_ps(to, from); }
+inline void ei_pstore(double* to, const __m128d& from) { _mm_store_pd(to, from); }
+inline void ei_pstore(int* to, const __m128i& from) { _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
+
+inline void ei_pstoreu(float* to, const __m128& from) { _mm_storeu_ps(to, from); }
+inline void ei_pstoreu(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
+inline void ei_pstoreu(int* to, const __m128i& from) { _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
+
+inline float ei_pfirst(const __m128& a) { return _mm_cvtss_f32(a); }
+inline double ei_pfirst(const __m128d& a) { return _mm_cvtsd_f64(a); }
+inline int ei_pfirst(const __m128i& a) { return _mm_cvtsi128_si32(a); }
+
+#ifdef __SSE3__
+// TODO implement SSE2 versions as well as integer versions
+inline __m128 ei_predux(const __m128* vecs)
+{
+ return _mm_hadd_ps(_mm_hadd_ps(vecs[0], vecs[1]),_mm_hadd_ps(vecs[2], vecs[3]));
+}
+inline __m128d ei_predux(const __m128d* vecs)
+{
+ return _mm_hadd_pd(vecs[0], vecs[1]);
+}
+
+inline float ei_predux(const __m128& a)
+{
+ __m128 tmp0 = _mm_hadd_ps(a,a);
+ return ei_pfirst(_mm_hadd_ps(tmp0, tmp0));
+}
+
+inline double ei_predux(const __m128d& a) { return ei_pfirst(_mm_hadd_pd(a, a)); }
+#else
+// SSE2 versions
+inline float ei_predux(const __m128& a)
+{
+ __m128 tmp = _mm_add_ps(a, _mm_movehl_ps(a,a));
+ return ei_pfirst(_mm_add_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1)));
+}
+inline double ei_predux(const __m128d& a)
+{
+ return ei_pfirst(_mm_add_sd(a, _mm_unpackhi_pd(a,a)));
+}
+
+inline __m128 ei_predux(const __m128* vecs)
+{
+ __m128 tmp0, tmp1, tmp2;
+ tmp0 = _mm_unpacklo_ps(vecs[0], vecs[1]);
+ tmp1 = _mm_unpackhi_ps(vecs[0], vecs[1]);
+ tmp2 = _mm_unpackhi_ps(vecs[2], vecs[3]);
+ tmp0 = _mm_add_ps(tmp0, tmp1);
+ tmp1 = _mm_unpacklo_ps(vecs[2], vecs[3]);
+ tmp1 = _mm_add_ps(tmp1, tmp2);
+ tmp2 = _mm_movehl_ps(tmp1, tmp0);
+ tmp0 = _mm_movelh_ps(tmp0, tmp1);
+ return _mm_add_ps(tmp0, tmp2);
+}
+
+inline __m128d ei_predux(const __m128d* vecs)
+{
+ return _mm_add_pd(_mm_unpacklo_pd(vecs[0], vecs[1]), _mm_unpackhi_pd(vecs[0], vecs[1]));
+}
+#endif // SSE3
+
+#endif // EIGEN_PACKET_MATH_SSE_H
+
diff --git a/Eigen/src/Core/ProductWIP.h b/Eigen/src/Core/ProductWIP.h
index a1c10d5d8..57bb899d6 100644
--- a/Eigen/src/Core/ProductWIP.h
+++ b/Eigen/src/Core/ProductWIP.h
@@ -193,6 +193,17 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
typedef typename ei_traits<Product>::_LhsNested _LhsNested;
typedef typename ei_traits<Product>::_RhsNested _RhsNested;
+ enum {
+ PacketSize = ei_packet_traits<Scalar>::size,
+ #if (defined __i386__)
+ // i386 architectures provides only 8 xmmm register,
+ // so let's reduce the max number of rows processed at once
+ MaxBlockRows = 4,
+ #else
+ MaxBlockRows = 8,
+ #endif
+ };
+
Product(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{
@@ -200,7 +211,18 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
}
/** \internal */
- template<typename DestDerived> void _cacheFriendlyEval(DestDerived& res) const;
+ template<typename DestDerived>
+ void _cacheFriendlyEval(DestDerived& res) const;
+
+ /** \internal */
+ template<typename DestDerived, int RhsAlignment, int ResAlignment>
+ void _cacheFriendlyEvalImpl(DestDerived& res) const __attribute__ ((noinline));
+
+ /** \internal */
+ template<typename DestDerived, int RhsAlignment, int ResAlignment, int BlockRows>
+ void _cacheFriendlyEvalKernel(DestDerived& res,
+ int l2i, int l2j, int l2k, int l1i,
+ int l2blockRowEnd, int l2blockColEnd, int l2blockSizeEnd, const Scalar* block) const EIGEN_DONT_INLINE;
private:
@@ -299,7 +321,7 @@ template<typename Derived>
template<typename Lhs, typename Rhs>
Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheOptimalProduct>& product)
{
- product._cacheFriendlyEval(*this);
+ product._cacheFriendlyEval(derived());
return derived();
}
@@ -307,6 +329,127 @@ template<typename Lhs, typename Rhs, int EvalMode>
template<typename DestDerived>
void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
{
+ const bool rhsIsAligned = (m_lhs.cols()%PacketSize == 0);
+ const bool resIsAligned = ((_rows()%PacketSize) == 0);
+
+ if (rhsIsAligned && resIsAligned)
+ _cacheFriendlyEvalImpl<DestDerived, Aligned, Aligned>(res);
+ else if (rhsIsAligned && (!resIsAligned))
+ _cacheFriendlyEvalImpl<DestDerived, Aligned, UnAligned>(res);
+ else if ((!rhsIsAligned) && resIsAligned)
+ _cacheFriendlyEvalImpl<DestDerived, UnAligned, Aligned>(res);
+ else
+ _cacheFriendlyEvalImpl<DestDerived, UnAligned, UnAligned>(res);
+
+}
+
+template<typename Lhs, typename Rhs, int EvalMode>
+template<typename DestDerived, int RhsAlignment, int ResAlignment, int BlockRows>
+void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEvalKernel(DestDerived& res,
+ int l2i, int l2j, int l2k, int l1i,
+ int l2blockRowEnd, int l2blockColEnd, int l2blockSizeEnd, const Scalar* block) const
+{
+ asm("#eigen begin kernel");
+
+ ei_internal_assert(BlockRows<=8);
+
+ // NOTE: sounds like we cannot rely on meta-unrolling to access dst[I] without enforcing GCC
+ // to create the dst's elements in memory, hence killing the performance.
+
+ for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
+ {
+ int offsetblock = l2k * (l2blockRowEnd-l2i) + (l1i-l2i)*(l2blockSizeEnd-l2k) - l2k*BlockRows;
+ const Scalar* localB = &block[offsetblock];
+
+ int l1jsize = l1j * m_lhs.cols(); //TODO find a better way to optimize address computation ?
+
+ // don't worry, dst is a set of registers
+ PacketScalar dst[BlockRows];
+ dst[0] = ei_pset1(Scalar(0.));
+ switch(BlockRows)
+ {
+ case 8: dst[7] = dst[0];
+ case 7: dst[6] = dst[0];
+ case 6: dst[5] = dst[0];
+ case 5: dst[4] = dst[0];
+ case 4: dst[3] = dst[0];
+ case 3: dst[2] = dst[0];
+ case 2: dst[1] = dst[0];
+ default: break;
+ }
+
+ // let's declare a few other temporary registers
+ PacketScalar tmp, tmp1;
+
+ // unaligned loads are expensive, therefore let's preload the next element in advance
+ if (RhsAlignment==UnAligned)
+ tmp1 = ei_ploadu(&m_rhs.derived().data()[l1jsize+l2k]);
+
+ for(int k=l2k; k<l2blockSizeEnd; k+=PacketSize)
+ {
+ // FIXME if we don't cache l1j*m_lhs.cols() then the performance are poor,
+ // let's directly access to the data
+ //PacketScalar tmp = m_rhs.template packetCoeff<Aligned>(k, l1j);
+ if (RhsAlignment==Aligned)
+ {
+ tmp = ei_pload(&m_rhs.data()[l1jsize + k]);
+ }
+ else
+ {
+ tmp = tmp1;
+ if (k+PacketSize<l2blockSizeEnd)
+ tmp1 = ei_ploadu(&m_rhs.data()[l1jsize + k+PacketSize]);
+ }
+
+ dst[0] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows ])), dst[0]);
+ if (BlockRows>=2) dst[1] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+ PacketSize])), dst[1]);
+ if (BlockRows>=3) dst[2] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+2*PacketSize])), dst[2]);
+ if (BlockRows>=4) dst[3] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+3*PacketSize])), dst[3]);
+ if (BlockRows>=5) dst[4] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+4*PacketSize])), dst[4]);
+ if (BlockRows>=6) dst[5] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+5*PacketSize])), dst[5]);
+ if (BlockRows>=7) dst[6] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+6*PacketSize])), dst[6]);
+ if (BlockRows>=8) dst[7] = ei_pmadd(tmp, ei_pload(&(localB[k*BlockRows+7*PacketSize])), dst[7]);
+ }
+
+ enum {
+ // Number of rows we can reduce per packet
+ PacketRows = (ResAlignment==Aligned && PacketSize>1) ? (BlockRows / PacketSize) : 0,
+ // First row index from which we have to to do redux once at a time
+ RemainingStart = PacketSize * PacketRows
+ };
+
+ // we have up to 4 packets (for doubles: 8 rows / 2)
+ if (PacketRows>=1)
+ res.template writePacketCoeff<Aligned>(l1i, l1j,
+ ei_padd(res.template packetCoeff<Aligned>(l1i, l1j), ei_predux(&(dst[0]))));
+ if (PacketRows>=2)
+ res.template writePacketCoeff<Aligned>(l1i+PacketSize, l1j,
+ ei_padd(res.template packetCoeff<Aligned>(l1i+PacketSize, l1j), ei_predux(&(dst[PacketSize]))));
+ if (PacketRows>=3)
+ res.template writePacketCoeff<Aligned>(l1i+2*PacketSize, l1j,
+ ei_padd(res.template packetCoeff<Aligned>(l1i+2*PacketSize, l1j), ei_predux(&(dst[2*PacketSize]))));
+ if (PacketRows>=4)
+ res.template writePacketCoeff<Aligned>(l1i+3*PacketSize, l1j,
+ ei_padd(res.template packetCoeff<Aligned>(l1i+3*PacketSize, l1j), ei_predux(&(dst[3*PacketSize]))));
+
+ // process the remaining rows one at a time
+ if (RemainingStart<=0 && BlockRows>=1) res.coeffRef(l1i+0, l1j) += ei_predux(dst[0]);
+ if (RemainingStart<=1 && BlockRows>=2) res.coeffRef(l1i+1, l1j) += ei_predux(dst[1]);
+ if (RemainingStart<=2 && BlockRows>=3) res.coeffRef(l1i+2, l1j) += ei_predux(dst[2]);
+ if (RemainingStart<=3 && BlockRows>=4) res.coeffRef(l1i+3, l1j) += ei_predux(dst[3]);
+ if (RemainingStart<=4 && BlockRows>=5) res.coeffRef(l1i+4, l1j) += ei_predux(dst[4]);
+ if (RemainingStart<=5 && BlockRows>=6) res.coeffRef(l1i+5, l1j) += ei_predux(dst[5]);
+ if (RemainingStart<=6 && BlockRows>=7) res.coeffRef(l1i+6, l1j) += ei_predux(dst[6]);
+ if (RemainingStart<=7 && BlockRows>=8) res.coeffRef(l1i+7, l1j) += ei_predux(dst[7]);
+
+ asm("#eigen end kernel");
+ }
+}
+
+template<typename Lhs, typename Rhs, int EvalMode>
+template<typename DestDerived, int RhsAlignment, int ResAlignment>
+void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEvalImpl(DestDerived& res) const
+{
// allow direct access to data for benchmark purpose
const Scalar* __restrict__ a = m_lhs.derived().data();
const Scalar* __restrict__ b = m_rhs.derived().data();
@@ -316,21 +459,13 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
// then we don't need to clear res and avoid and additional mat-mat sum
// res.setZero();
- const int ps = ei_packet_traits<Scalar>::size; // size of a packet
- #if (defined __i386__)
- // i386 architectures provides only 8 xmmm register,
- // so let's reduce the max number of rows processed at once
- const int bw = 4; // number of rows treated at once
- #else
- const int bw = 8; // number of rows treated at once
- #endif
- const int bs = ps * bw; // total number of elements treated at once
+ const int bs = PacketSize * MaxBlockRows; // total number of elements treated at once
const int rows = _rows();
const int cols = _cols();
- const int size = m_lhs.cols(); // third dimension of the product
+ const int remainingSize = m_lhs.cols()%PacketSize;
+ const int size = m_lhs.cols() - remainingSize; // third dimension of the product clamped to packet boundaries
const int l2blocksize = 256 > _cols() ? _cols() : 256;
- const bool rhsIsAligned = ((size%ps) == 0);
- const bool resIsAligned = ((cols%ps) == 0);
+ // FIXME use calloca ?? (allocation on the stack)
Scalar* __restrict__ block = new Scalar[l2blocksize*size];
// loops on each L2 cache friendly blocks of the result
@@ -348,23 +483,23 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
{
const int l2blockSizeEnd = std::min(l2k+l2blocksize, size);
- for (int i = l2i; i<l2blockRowEndBW; i+=bw)
+ for (int i = l2i; i<l2blockRowEndBW; i+=MaxBlockRows)
{
- for (int k=l2k; k<l2blockSizeEnd; k+=ps)
+ for (int k=l2k; k<l2blockSizeEnd; k+=PacketSize)
{
// TODO write these two loops using meta unrolling
// negligible for large matrices but useful for small ones
- for (int w=0; w<bw; ++w)
- for (int s=0; s<ps; ++s)
+ for (int w=0; w<MaxBlockRows; ++w)
+ for (int s=0; s<PacketSize; ++s)
block[count++] = m_lhs.coeff(i+w,k+s);
}
}
if (l2blockRowRemaining>0)
{
- for (int k=l2k; k<l2blockSizeEnd; k+=ps)
+ for (int k=l2k; k<l2blockSizeEnd; k+=PacketSize)
{
for (int w=0; w<l2blockRowRemaining; ++w)
- for (int s=0; s<ps; ++s)
+ for (int s=0; s<PacketSize; ++s)
block[count++] = m_lhs.coeff(l2blockRowEndBW+w,k+s);
}
}
@@ -376,19 +511,21 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
for(int l2k=0; l2k<size; l2k+=l2blocksize)
{
- // acumulate a full row of current a block time 4 cols of current a block
- // to a 1x4 c block
+ // acumulate a bw rows of lhs time a single column of rhs to a bw x 1 block of res
int l2blockSizeEnd = std::min(l2k+l2blocksize, size);
- // for each 4x1 result's block sub blocks...
- for(int l1i=l2i; l1i<l2blockRowEndBW; l1i+=bw)
+ // for each bw x 1 result's block
+ for(int l1i=l2i; l1i<l2blockRowEndBW; l1i+=MaxBlockRows)
{
+ _cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, MaxBlockRows>(
+ res, l2i, l2j, l2k, l1i, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block);
+#if 0
for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
{
int offsetblock = l2k * (l2blockRowEnd-l2i) + (l1i-l2i)*(l2blockSizeEnd-l2k) - l2k*bw/*bs*/;
const Scalar* localB = &block[offsetblock];
- int l1jsize = l1j * size; //TODO find a better way to optimize address computation ?
+ int l1jsize = l1j * m_lhs.cols(); //TODO find a better way to optimize address computation ?
PacketScalar dst[bw];
dst[0] = ei_pset1(Scalar(0.));
@@ -408,7 +545,8 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
asm("#eigen begincore");
for(int k=l2k; k<l2blockSizeEnd; k+=ps)
{
- //PacketScalar tmp = m_rhs.packetCoeff(k, l1j);
+// PacketScalar tmp = m_rhs.template packetCoeff<Aligned>(k, l1j);
+ // TODO make this branching compile time (costly for doubles)
if (rhsIsAligned)
tmp = ei_pload(&m_rhs.derived().data()[l1jsize + k]);
else
@@ -436,21 +574,61 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
}
}
- res.template writePacketCoeff<Aligned>(l1i, l1j, ei_padd(res.template packetCoeff<Aligned>(l1i, l1j), ei_predux(dst)));
- if (ps==2)
- res.template writePacketCoeff<Aligned>(l1i+2,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+2,l1j), ei_predux(&(dst[2]))));
- if (bw==8)
+// if (resIsAligned)
{
- res.template writePacketCoeff<Aligned>(l1i+4,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+4,l1j), ei_predux(&(dst[4]))));
+ res.template writePacketCoeff<Aligned>(l1i, l1j, ei_padd(res.template packetCoeff<Aligned>(l1i, l1j), ei_predux(dst)));
if (ps==2)
- res.template writePacketCoeff<Aligned>(l1i+6,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+6,l1j), ei_predux(&(dst[6]))));
+ res.template writePacketCoeff<Aligned>(l1i+2,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+2,l1j), ei_predux(&(dst[2]))));
+ if (bw==8)
+ {
+ res.template writePacketCoeff<Aligned>(l1i+4,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+4,l1j), ei_predux(&(dst[4]))));
+ if (ps==2)
+ res.template writePacketCoeff<Aligned>(l1i+6,l1j, ei_padd(res.template packetCoeff<Aligned>(l1i+6,l1j), ei_predux(&(dst[6]))));
+ }
}
+// else
+// {
+// // TODO uncommenting this code kill the perf, even though it is never called !!
+// // TODO optimize this loop
+// // TODO is it better to do one redux at once or packet reduxes + unaligned store ?
+// for (int w = 0; w<bw; ++w)
+// res.coeffRef(l1i+w, l1j) += ei_predux(dst[w]);
+// std::cout << "!\n";
+// }
asm("#eigen endcore");
}
+#endif
}
if (l2blockRowRemaining>0)
{
+ // this is an attempt to build an array of kernels, but I did not manage to get it compiles
+// typedef void (*Kernel)(DestDerived& , int, int, int, int, int, int, int, const Scalar*);
+// Kernel kernels[8];
+// kernels[0] = (Kernel)(&Product<Lhs,Rhs,EvalMode>::template _cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 1>);
+// kernels[l2blockRowRemaining](res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block);
+
+ switch(l2blockRowRemaining)
+ {
+ case 1:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 1>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 2:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 2>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 3:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 3>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 4:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 4>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 5:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 5>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 6:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 6>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ case 7:_cacheFriendlyEvalKernel<DestDerived, RhsAlignment, ResAlignment, 7>(
+ res, l2i, l2j, l2k, l2blockRowEndBW, l2blockRowEnd, l2blockColEnd, l2blockSizeEnd, block); break;
+ default:
+ ei_internal_assert(false && "internal error"); break;
+ }
+
+#if 0
// TODO optimize this part using a generic templated function that processes N rows
// here we process the remaining l2blockRowRemaining rows
for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
@@ -460,13 +638,13 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
int l1jsize = l1j * size;
- PacketScalar dst[bw];
+ PacketScalar dst[MaxBlockRows];
dst[0] = ei_pset1(Scalar(0.));
for (int w = 1; w<l2blockRowRemaining; ++w)
dst[w] = dst[0];
PacketScalar b0, b1, tmp;
asm("#eigen begincore dynamic");
- for(int k=l2k; k<l2blockSizeEnd; k+=ps)
+ for(int k=l2k; k<l2blockSizeEnd; k+=PacketSize)
{
//PacketScalar tmp = m_rhs.packetCoeff(k, l1j);
if (rhsIsAligned)
@@ -476,7 +654,7 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
// TODO optimize this loop
for (int w = 0; w<l2blockRowRemaining; ++w)
- dst[w] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRowRemaining+w*ps])), dst[w]);
+ dst[w] = ei_pmadd(tmp, ei_pload(&(localB[k*l2blockRowRemaining+w*PacketSize])), dst[w]);
}
// TODO optimize this loop
@@ -485,11 +663,23 @@ void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
asm("#eigen endcore dynamic");
}
+#endif
}
}
}
}
+ // handle the part which cannot be processed by the vectorized path
+ if (remainingSize)
+ {
+ res += Product<
+ Block<typename ei_unconst<_LhsNested>::type,Dynamic,Dynamic>,
+ Block<typename ei_unconst<_RhsNested>::type,Dynamic,Dynamic>,
+ NormalProduct>(
+ m_lhs.block(0,size, _rows(), remainingSize),
+ m_rhs.block(size,0, remainingSize, _cols())).lazy();
+ }
+
delete[] block;
}
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index be5e7bba5..613ab9617 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -85,6 +85,12 @@ using Eigen::MatrixBase;
#endif
#if (defined __GNUC__)
+#define EIGEN_DONT_INLINE __attribute__((noinline))
+#else
+#define EIGEN_DONT_INLINE
+#endif
+
+#if (defined __GNUC__)
#define EIGEN_ALIGN_128 __attribute__ ((aligned(16)))
#else
#define EIGEN_ALIGN_128
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index 19768c1ca..264802b9b 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -189,6 +189,9 @@ template<typename T> class ei_eval
template<typename T> struct ei_unref { typedef T type; };
template<typename T> struct ei_unref<T&> { typedef T type; };
+template<typename T> struct ei_unconst { typedef T type; };
+template<typename T> struct ei_unconst<const T> { typedef T type; };
+
template<typename T> struct ei_is_temporary
{
enum { ret = 0 };
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-05 00:34:42 +0000