[PROJ] Vector/SIMD acceleration

[Charles Karney]

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
>