mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-31 19:00:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

* add ei_pdiv intrinsic, make quotient functor vectorizable

Browse Source 
* add vdw benchmark from Tim's real-world use case
This commit is contained in:
Benoit Jacob 2008-06-23 22:00:18 +00:00
parent ac9aa47bbc
commit c9560df4a0
6 changed files with 105 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Core/DummyPacketMath.h
View File

@ -38,6 +38,9 @@ template <typename Scalar> inline Scalar ei_psub(const Scalar& a, const Scalar&
/** \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 (coeff-wise) */
template <typename Scalar> inline Scalar ei_pdiv(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); }

18
Eigen/src/Core/Functors.h
View File

@ -131,12 +131,18 @@ struct ei_functor_traits<ei_scalar_difference_op<Scalar> > {
* \sa class CwiseBinaryOp, MatrixBase::cwiseQuotient()
*/
template<typename Scalar> struct ei_scalar_quotient_op EIGEN_EMPTY_STRUCT {
inline const Scalar operator() (const Scalar& a, const Scalar& b) const { return a / b; }
inline const Scalar operator() (const Scalar& a, const Scalar& b) const { return a / b; }
template<typename PacketScalar>
inline const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const
{ return ei_pdiv(a,b); }
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_quotient_op<Scalar> >
{ enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = false }; };
struct ei_functor_traits<ei_scalar_quotient_op<Scalar> > {
enum {
Cost = 2 * NumTraits<Scalar>::MulCost,
PacketAccess = ei_packet_traits<Scalar>::size>1
};
};
// unary functors:
@ -179,7 +185,7 @@ template<typename Scalar> struct ei_scalar_abs2_op EIGEN_EMPTY_STRUCT {
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_abs2_op<Scalar> >
{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = NumTraits<Scalar>::IsComplex==false && int(ei_packet_traits<Scalar>::size)>1 }; };
{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = int(ei_packet_traits<Scalar>::size)>1 }; };
/** \internal
* \brief Template functor to compute the conjugate of a complex value
@ -272,7 +278,7 @@ struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> >
* \brief Template functor to divide a scalar by a fixed other one
*
* This functor is used to implement the quotient of a matrix by
* a scalar where the scalar type is not a floating point type.
* a scalar where the scalar type is not necessarily a floating point type.
*
* \sa class CwiseUnaryOp, MatrixBase::operator/
*/

2
Eigen/src/Core/Sum.h
View File

@ -216,7 +216,7 @@ struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling>
if(alignedSize == size) return res;
}
else // too small to vectorize anything.
// since this is dynamic-size hence inefficient anyway, don't try to optimize.
// since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
{
res = Scalar(0);
}

2
Eigen/src/Core/arch/AltiVec/PacketMath.h
View File

@ -59,6 +59,8 @@ inline vector int ei_pmul(const vector int a, const vector int b)
return vec_add( lowProduct, highProduct );
}
inline vector float ei_pdiv(const vector float a, const vector float b) { return vec_div(a,b); }
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); }

3
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@ -55,6 +55,9 @@ template<> inline __m128i ei_pmul(const __m128i& a, const __m128i& b)
_mm_setr_epi32(0xffffffff,0,0xffffffff,0)), 4));
}
template<> inline __m128 ei_pdiv(const __m128& a, const __m128& b) { return _mm_div_ps(a,b); }
template<> inline __m128d ei_pdiv(const __m128d& a, const __m128d& b) { return _mm_div_pd(a,b); }
// 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); }

84
bench/vdw_new.cpp Normal file
View File

@ -0,0 +1,84 @@
#include <Eigen/Core>
namespace Eigen {
template<typename Scalar> struct pow12_op EIGEN_EMPTY_STRUCT {
inline const Scalar operator() (const Scalar& a) const
{
Scalar b = a*a*a;
Scalar c = b*b;
return c*c;
}
template<typename PacketScalar>
inline const PacketScalar packetOp(const PacketScalar& a) const
{
PacketScalar b = ei_pmul(a, ei_pmul(a, a));
PacketScalar c = ei_pmul(b, b);
return ei_pmul(c, c);
}
};
template<typename Scalar>
struct ei_functor_traits<pow12_op<Scalar> >
{
enum {
Cost = 4 * NumTraits<Scalar>::MulCost,
PacketAccess = int(ei_packet_traits<Scalar>::size) > 1
};
};
} // namespace Eigen
using Eigen::pow12_op;
USING_PART_OF_NAMESPACE_EIGEN
#ifndef SCALAR
#define SCALAR float
#endif
#ifndef SIZE
#define SIZE 10000
#endif
#ifndef REPEAT
#define REPEAT 10000
#endif
typedef Matrix<SCALAR, Eigen::Dynamic, 1> Vec;
using namespace std;
SCALAR E_VDW(const Vec &interactions1, const Vec &interactions2)
{
return interactions2
.cwiseQuotient(interactions1)
.cwise(pow12_op<SCALAR>())
.sum();
}
int main()
{
//
// 1 2 3 4 ... (interactions)
// ka . . . . ...
// rab . . . . ...
// energy . . . . ...
// ... ... ... ... ... ...
// (variables
// for
// interaction)
//
Vec interactions1(SIZE), interactions2(SIZE); // SIZE is the number of vdw interactions in our system
// SetupCalculations()
SCALAR rab = 1.0;
interactions1.setConstant(2.4);
interactions2.setConstant(rab);
// Energy()
SCALAR energy = 0.0;
for (unsigned int i = 0; i<REPEAT; ++i) {
energy += E_VDW(interactions1, interactions2);
energy *= 1 + 1e-20 * i; // prevent compiler from optimizing the loop
}
cout << "energy = " << energy << endl;
}