diff --git a/Eigen/src/Core/functors/NullaryFunctors.h b/Eigen/src/Core/functors/NullaryFunctors.h index cd9fbf267..71629af4c 100644 --- a/Eigen/src/Core/functors/NullaryFunctors.h +++ b/Eigen/src/Core/functors/NullaryFunctors.h @@ -37,7 +37,7 @@ template struct functor_traits > { enum { Cost = NumTraits::AddCost, PacketAccess = false, IsRepeatable = true }; }; -template struct linspaced_op_impl; +template struct linspaced_op_impl; // linear access for packet ops: // 1) initialization @@ -48,12 +48,12 @@ template 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 -struct linspaced_op_impl +struct linspaced_op_impl { - linspaced_op_impl(const Scalar& low, const Scalar& step) : - m_low(low), m_step(step), - m_packetStep(pset1(unpacket_traits::size*step)), - m_base(padd(pset1(low), pmul(pset1(step),plset(-unpacket_traits::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(unpacket_traits::size*m_step)), + m_base(padd(pset1(low), pmul(pset1(m_step),plset(-unpacket_traits::size)))) {} template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (Index i) const @@ -75,11 +75,11 @@ struct linspaced_op_impl // 1) each step // [low, ..., low] + ( [step, ..., step] * ( [i, ..., i] + [0, ..., size] ) ) template -struct linspaced_op_impl +struct linspaced_op_impl { - linspaced_op_impl(const Scalar& low, const Scalar& step) : - m_low(low), m_step(step), - m_lowPacket(pset1(m_low)), m_stepPacket(pset1(m_step)), m_interPacket(plset(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(m_low)), m_stepPacket(pset1(m_step)), m_interPacket(plset(0)) {} template 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 const Packet m_interPacket; }; +template +struct linspaced_op_impl +{ + 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(0)) + {} + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const Scalar operator() (Index i) const { + return m_low + (m_length*Scalar(i))/(m_numSteps-1); + } + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + const Packet packetOp(Index i) const { + return internal::padd(pset1(m_low), pdiv(pmul(pset1(m_length), padd(pset1(Scalar(i)),m_interPacket)), + pset1(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 // nested expressions). template struct linspaced_op; template struct functor_traits< linspaced_op > -{ enum { Cost = 1, PacketAccess = packet_traits::HasSetLinear, IsRepeatable = true }; }; +{ + enum + { + Cost = 1, + PacketAccess = packet_traits::HasSetLinear + && ((!NumTraits::IsInteger) || packet_traits::HasDiv), + IsRepeatable = true + }; +}; template 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 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); } @@ -134,7 +169,9 @@ template 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 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::IsInteger?true:RandomAccess),NumTraits::IsInteger> impl; }; // all functors allow linear access, except scalar_identity_op. So we fix here a quick meta diff --git a/test/nullary.cpp b/test/nullary.cpp index 4844f2952..8d65910eb 100644 --- a/test/nullary.cpp +++ b/test/nullary.cpp @@ -48,30 +48,32 @@ void testVectorType(const VectorType& base) VectorType m(base); m.setLinSpaced(size,low,high); + if(!NumTraits::IsInteger) + { + VectorType n(size); + for (int i=0; i::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 row_vector(size); @@ -126,5 +128,12 @@ void test_nullary() CALL_SUBTEST_8( testVectorType(Vector4f()) ); CALL_SUBTEST_8( testVectorType(Matrix()) ); CALL_SUBTEST_8( testVectorType(Matrix()) ); + + CALL_SUBTEST_9( testVectorType(VectorXi(internal::random(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::epsilon() ); +#endif }