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