Re: [eigen] matrix::linspace

[Hauke Heibel <hauke.heibel@xxxxxxxxxxxxxx>]

The initial commit is there, including documentation and new snippets.
There is one issue left...

VectorXd a = VectorXd::LinSpaced(low,high,size);

does not necessarily result in

a(0) = low

nor in

a(a.size()-1) = high

I have no real idea how to fix this - I mean none that does not impact
performance. The problem arises due to numerical instabilities. The
sequential version accumulates values and thus 'size' times the
floating point error inherent in adding values. The random access
method suffers from the fact that for s = (high-low)/(size-1)

high = (size-1)*s

is not necessarily true.

Any ideas are welcome.

- Hauke

On Tue, Jan 26, 2010 at 2:02 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> Hi,
>
> 2010/1/25 Hauke Heibel <hauke.heibel@xxxxxxxxxxxxxx>:
>> Ok, I am done. I have the whole patch ready for review. Since I have
>> no idea what else to do for reviewing larger commits, I am posting
>> here.
>>
>> The main new functionality offered is
>>
>>  DenseBase::setLinSpace(Scalar l,Scalar h,int s)
>>
>
> as discussed on IRC, LinSpace ---> LinSpaced.
>
>> which computes for each vector index
>>
>>  a(i) = l + i*(h-l)/(s-1)
>>
>> Now to the changes this involves.
>>
>> a) I added a new packet op ei_plset(Scalar a). It creates packages [a,
>> a+1, ..., a+packet_size]. The nice part as opposed to ei_pset(...) is
>> that for all types (float, double, int) there is only a single
>> argument and not a varying number of arguments. This also allows to
>> offer a GenericPacketMath implementation where this operator acts as
>> identity.
>
> sounds good,
>
>>
>> b) I added new constraints for functors that are meant to be used with
>> CwiseNullaryOp. They are now required to offer a vector indexing
>> operator()(int). For minimal invasive modifications, I simply changed
>> existing ones by defining operator()(int,int) as operator()(int,int =
>> 0) - where it was appropriate. That was already partially done.
>
> OK,
>
>>
>> c) As discussed, ei_linspace_op_impl has two specializations one for
>> random access and one for sequential or linear access. The vectorized
>> sequential versions are much faster than non-vectorized ones and the
>> random access versions seem to perform equally fast and faster (at
>> least not worse).
>
> I understand that the API for the sequential case is setLinSpaced.
> What is the API for the random case? Are you going to provide a static
> LinSpaced() ? Meanwhile, it's untested?
>
>>
>> d) I introduced a new proxy ojbect ei_unaligned_assign_impl. It
>> performs a simple loop assignment when required and it does not
>> inline. Preventing the inlining was crucial in order to allow MSVC to
>> perform proper optimization of packet ops. Maybe this should be
>> guarded by a pre-processor define and activated only for MSVC.
>
> Good stuff, it even decrements the number of for loops in Eigen.
>
>> If I am not wrong, this is all. Any comments, or am I good to synch?
>
> You can push, I'm just curious about the random access version,
> especially API wise and testing wise.
>
> Benoit
>
>
>