[PROJ] Vector/SIMD acceleration

[Andrew Bell]

I'll add that Eigen, or something like that, has the advantages that 1) no
special expertise is needed to maintain the code and 2) as CPUs change, the
optimization comes without a code change on the part of the user -- someone
who's a relative expert in optimizing numeric computations is doing the
work.

On Fri, Apr 17, 2020 at 4:49 PM Charles Karney <charles at karney.com> wrote:

> Eigen is more capable than you imply (and this is the point of the
> link Andrew sent).  An Eigen "matrix" is for linear algebra.  But
> an Eigen "array" is for element-wise operations on arrays of numbers.
> and "scalar * array" and "sin(array)" do the "right things".  I'll
> have to check with comparisons.  I assume the ">" operator will
> return an array of bools and then you'll have to arrange that different
> portions of your array are handled differently (like MATLAB/Octave?
> see also the old Cray vector merge intrinsics).
>
> I'm not sure that Eigen will be the right approach.  But I think it's
> close.  (And probably my assertion that the code can remain essentially
> the same, except for the declarations, is a pipe dream.)
>
> By the way, Eigen is an awesome package.  It implements everything with
> header code and so the end result is often very efficient.  The downside
> is that because of all the flexibility that Eigen provides with how
> arrays are laid out and allocated, you can sometimes end up in template
> hell.
>
> On 4/17/20 3:36 PM, Even Rouault wrote:
> > On vendredi 17 avril 2020 13:29:17 CEST Charles Karney wrote:
> >
> >  > I wonder whether it's possible to get the benefits of vectorization
> >
> >  > without massive of changes to the code.
> >
> > That's exactly what I had in mind.
> >
> >  > Perhaps the basic projection
> >
> >  > functions could be templated to allow Eigen arrays of PJ_LP's to be
> >
> >  > passed. Eigen already has overloads for componentwise arithmetic
> >
> >  > operations and handles SIMD vectorization automatically.
> >
> > I didn't know Eigen, but from a quick look & try, it seems I've done
> > something similar.
> >
> > At least for types like Vector2d, Vector4d. But Eigen is quite limited
> > in the operations it supports. Just +, -, *, /. No comparisons, sqrt,
> > sin, cos, etc.
> >
> > Which my prototype supports (through sleef for sin, cos, etc).
> >
> > Adding support for functions that take Eigen types, like sqrt(Vector4d)
> > should be possible, but comparison operators require to be able to
> > modify the class definition.
> >
> > Another pain point I found with Eigen, is that it refuses to do things
> like
> >
> > a * b on vectors because of potential ambiguity of the operator in the
> > context of linear algebra (matrix multiplication, cross product,
> > member-to-member multiplication ?). But in a PROJ world where we'd port
> > from double to multiple-double, we just want member-to-member
> > multipication, and with Eigen you have to do it explictly with
> > a.cwiseProduct(b). Would require more changes and would clutter the
> > readability of the code.
> >
> >  > A major advantage of this approach is that PROJ doesn't need to get
> into
> >
> >  > the weeds with SIMD instructions. So when a new instruction set comes
> >
> >  > along we can (I hope) rely on the maintains of Eigen to do the hard
> >
> >  > work. (I see that there's already some interoperability with Eigen and
> >
> >  > CUDA.)
> >
> > Admitedly I've only prototyped SSE and AVX. But supporting other archs
> > would only require to add a specialization in the generic template stuff
> > (not a big deal. It is of the order of 10 intrinsincs per arch), not
> > projection code. And for non-supported archs, the templated operations
> > would just expand to a 1-coordinate at a time implementation, like the
> > current serial code. But I guess if Intel is covered, we have already >
> > 90% use cases that matter (nobody ever paid me to optimize something for
> > non-Intel archs).
> >
> > One potential downside (in the context of PROJ licencing) of Eigen is
> > its license. It is MPL 2, weak copyleft, which would require binary
> > distributions of PROJ to provide a way of accessing Eigen source code.
> > Not a huge deal, but this would put additional requirements for
> > proprietary shops using PROJ.
> >
> >  > If this approach is followed, I would also recommend that the basic
> >
> >  > floating point type, double, be either templated or typedef'ed.
> >
> > Yep, exactly my idea of template<typename T, int N> struct VF
> >
> > --
> >
> > Spatialys - Geospatial professional services
> >
> > http://www.spatialys.com
> >
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
>


-- 
Andrew Bell
andrew.bell.ia at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200417/3fb1ea79/attachment-0001.html>