Re: [eigen] matrix::linspace |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] matrix::linspace*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Tue, 26 Jan 2010 13:56:27 -0500*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=L63K7tmNwjZuYQdvQH2QusQpM0vs9p4e8HU+NnsoDAs=; b=g+UvJfwNDk+E0SmGw/BQcRWNhwKONrh6t+pVKy+rckmIA0pRlpktb6JPgEsOF2HJTU Apk/+cYGV8Ufr3f+llijdOyQZg2tLoFwT9ADPR3teyai10hD+savcc9imG9LyOpLeus9 mbjVzPPWQEic5uqResLabS2KLPLvaA37B6drM=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=ZXmj1tcPgIrlNLk/cICFjP27sB1g2PUp1MB7Czdk6ye0rCeAlCu7+Ilkhv5ZWnPjqA ijVJdOJ9hufdLNfULTgDTpr1ELhg6ufuGchMo22N99XKGN7lAI7kK2yA+y/Q2c1kO7v7 Y3iukamyxG830+pOVBzJkNsR48nPjEIvPob6k=

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

**Follow-Ups**:**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**References**:**[eigen] matrix::linspace***From:*Hauke Heibel

**Re: [eigen] matrix::linspace***From:*Benoit Jacob

**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**Re: [eigen] matrix::linspace***From:*Benoit Jacob

**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**Re: [eigen] matrix::linspace***From:*Benoit Jacob

**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**Re: [eigen] matrix::linspace***From:*Benoit Jacob

**Re: [eigen] matrix::linspace***From:*Hauke Heibel

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] matrix::linspace** - Next by Date:
**Re: [eigen] std::vector specialization** - Previous by thread:
**Re: [eigen] matrix::linspace** - Next by thread:
**Re: [eigen] matrix::linspace**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |