[PROJ] Vector/SIMD acceleration

[Even Rouault]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200417/1952cf49/attachment.html>