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

bug #698: fix linspaced for integer types.

Browse Source
This commit is contained in:
Gael Guennebaud 2016-02-01 14:25:34 +01:00
parent 2c3224924b
commit e1d219e5c9
2 changed files with 70 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
Eigen/src/Core/functors/NullaryFunctors.h
View File

@ -37,7 +37,7 @@ template<typename Scalar>
struct functor_traits<scalar_identity_op<Scalar> >
{ enum { Cost = NumTraits<Scalar>::AddCost, PacketAccess = false, IsRepeatable = true }; };
template <typename Scalar, typename Packet, bool RandomAccess> struct linspaced_op_impl;
template <typename Scalar, typename Packet, bool RandomAccess, bool IsInteger> struct linspaced_op_impl;
// linear access for packet ops:
// 1) initialization
@ -48,12 +48,12 @@ template <typename Scalar, typename Packet, bool RandomAccess> struct linspaced_
// TODO: Perhaps it's better to initialize lazily (so not in the constructor but in packetOp)
// in order to avoid the padd() in operator() ?
template <typename Scalar, typename Packet>
struct linspaced_op_impl<Scalar,Packet,false>
struct linspaced_op_impl<Scalar,Packet,/*RandomAccess*/false,/*IsInteger*/false>
{
linspaced_op_impl(const Scalar& low, const Scalar& step) :
m_low(low), m_step(step),
m_packetStep(pset1<Packet>(unpacket_traits<Packet>::size*step)),
m_base(padd(pset1<Packet>(low), pmul(pset1<Packet>(step),plset<Packet>(-unpacket_traits<Packet>::size)))) {}
linspaced_op_impl(const Scalar& low, const Scalar& high, Index num_steps) :
m_low(low), m_step(num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1)),
m_packetStep(pset1<Packet>(unpacket_traits<Packet>::size*m_step)),
m_base(padd(pset1<Packet>(low), pmul(pset1<Packet>(m_step),plset<Packet>(-unpacket_traits<Packet>::size)))) {}
template<typename Index>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (Index i) const
@ -75,11 +75,11 @@ struct linspaced_op_impl<Scalar,Packet,false>
// 1) each step
// [low, ..., low] + ( [step, ..., step] * ( [i, ..., i] + [0, ..., size] ) )
template <typename Scalar, typename Packet>
struct linspaced_op_impl<Scalar,Packet,true>
struct linspaced_op_impl<Scalar,Packet,/*RandomAccess*/true,/*IsInteger*/false>
{
linspaced_op_impl(const Scalar& low, const Scalar& step) :
m_low(low), m_step(step),
m_lowPacket(pset1<Packet>(m_low)), m_stepPacket(pset1<Packet>(m_step)), m_interPacket(plset<Packet>(0)) {}
linspaced_op_impl(const Scalar& low, const Scalar& high, Index num_steps) :
m_low(low), m_step(num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1)),
m_lowPacket(pset1<Packet>(m_low)), m_stepPacket(pset1<Packet>(m_step)), m_interPacket(plset<Packet>(0)) {}
template<typename Index>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return m_low+i*m_step; }
@ -95,6 +95,31 @@ struct linspaced_op_impl<Scalar,Packet,true>
const Packet m_interPacket;
};
template <typename Scalar, typename Packet>
struct linspaced_op_impl<Scalar,Packet,/*RandomAccess*/true,/*IsInteger*/true>
{
linspaced_op_impl(const Scalar& low, const Scalar& high, Index num_steps) :
m_low(low), m_length(high-low), m_numSteps(num_steps), m_interPacket(plset<Packet>(0))
{}
template<typename Index>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
const Scalar operator() (Index i) const {
return m_low + (m_length*Scalar(i))/(m_numSteps-1);
}
template<typename Index>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
const Packet packetOp(Index i) const {
return internal::padd(pset1<Packet>(m_low), pdiv(pmul(pset1<Packet>(m_length), padd(pset1<Packet>(Scalar(i)),m_interPacket)),
pset1<Packet>(m_numSteps-1))); }
const Scalar m_low;
const Scalar m_length;
const Index m_numSteps;
const Packet m_interPacket;
};
// ----- Linspace functor ----------------------------------------------------------------
// Forward declaration (we default to random access which does not really give
@ -102,10 +127,20 @@ struct linspaced_op_impl<Scalar,Packet,true>
// nested expressions).
template <typename Scalar, typename PacketType, bool RandomAccess = true> struct linspaced_op;
template <typename Scalar, typename PacketType, bool RandomAccess> struct functor_traits< linspaced_op<Scalar,PacketType,RandomAccess> >
{ enum { Cost = 1, PacketAccess = packet_traits<Scalar>::HasSetLinear, IsRepeatable = true }; };
{
enum
{
Cost = 1,
PacketAccess = packet_traits<Scalar>::HasSetLinear
&& ((!NumTraits<Scalar>::IsInteger) || packet_traits<Scalar>::HasDiv),
IsRepeatable = true
};
};
template <typename Scalar, typename PacketType, bool RandomAccess> struct linspaced_op
{
linspaced_op(const Scalar& low, const Scalar& high, Index num_steps) : impl((num_steps==1 ? high : low), (num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1))) {}
linspaced_op(const Scalar& low, const Scalar& high, Index num_steps)
: impl((num_steps==1 ? high : low),high,num_steps)
{}
template<typename Index>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); }
@ -134,7 +169,9 @@ template <typename Scalar, typename PacketType, bool RandomAccess> struct linspa
// This proxy object handles the actual required temporaries, the different
// implementations (random vs. sequential access) as well as the
// correct piping to size 2/4 packet operations.
const linspaced_op_impl<Scalar,PacketType,RandomAccess> impl;
// As long as we don't have a Bresenham-like implementation for linear-access and integer types,
// we have to by-pass RandomAccess for integer types. See bug 698.
const linspaced_op_impl<Scalar,PacketType,(NumTraits<Scalar>::IsInteger?true:RandomAccess),NumTraits<Scalar>::IsInteger> impl;
};
// all functors allow linear access, except scalar_identity_op. So we fix here a quick meta

31
test/nullary.cpp
View File

@ -48,30 +48,32 @@ void testVectorType(const VectorType& base)
VectorType m(base);
m.setLinSpaced(size,low,high);
if(!NumTraits<Scalar>::IsInteger)
{
VectorType n(size);
for (int i=0; i<size; ++i)
n(i) = low+i*step;
VERIFY_IS_APPROX(m,n);
}
VectorType n(size);
for (int i=0; i<size; ++i)
n(i) = low+i*step;
n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/(size-1));
VERIFY_IS_APPROX(m,n);
// random access version
m = VectorType::LinSpaced(size,low,high);
VERIFY_IS_APPROX(m,n);
// Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<Scalar>::epsilon() );
// These guys sometimes fail! This is not good. Any ideas how to fix them!?
//VERIFY( m(m.size()-1) == high );
//VERIFY( m(0) == low );
VERIFY( internal::isApprox(m(m.size()-1),high) );
VERIFY( size==1 || internal::isApprox(m(0),low) );
// sequential access version
m = VectorType::LinSpaced(Sequential,size,low,high);
VERIFY_IS_APPROX(m,n);
// These guys sometimes fail! This is not good. Any ideas how to fix them!?
//VERIFY( m(m.size()-1) == high );
//VERIFY( m(0) == low );
VERIFY( internal::isApprox(m(m.size()-1),high) );
VERIFY( size==1 || internal::isApprox(m(0),low) );
// check whether everything works with row and col major vectors
Matrix<Scalar,Dynamic,1> row_vector(size);
@ -126,5 +128,12 @@ void test_nullary()
CALL_SUBTEST_8( testVectorType(Vector4f()) );
CALL_SUBTEST_8( testVectorType(Matrix<float,8,1>()) );
CALL_SUBTEST_8( testVectorType(Matrix<float,1,1>()) );
CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(1,300))) );
}
#ifdef EIGEN_TEST_PART_6
// Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<double>::epsilon() );
#endif
}