Re: [eigen] matrix::linspace

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/1/26 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.

Hm. "high = (size-1)*s" is true to machine precision which is good
enough. But I guess that what you mean is that

  s + s + s + ... + s    (size-1 times)

can be farther away from 'high'.

The only idea that I see is that every, say, 8 iterations, you
recompute the product with a multiplication (thus you get something
good to machine precision) and you forget about the 'accumulated'
value.

Benoit

>
> 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
>>
>>
>>
>
>
>