Re: [eigen] Eigen containers cannot be vectorized

You might try one of the following pragmas to inform the compiler to ignore aliasing in the loop nest. I haven't tested your code, so I apologize if none of these work.

#pragma ivdep (or #pragma GCC ivdep, in the case of GCC) is supported by many compilers. I don't have know the full list (https://github.com/jeffhammond/nwchem-tce-triples-kernels/blob/master/src/pragma_vendor.h lists all the ones I knew at some point, and Google tells me MSVC supports it too). In any case, it is easy enough to write build system tests for this sort of thing.

#pragma simd is an Intel compiler keyword. See https://software.intel.com/en-us/node/524555 for details. The Intel compiler takes this pragma (and the OpenMP equivalent noted below) seriously and will generate vector code when you use these, even in cases where the heuristics say it isn't beneficial.

#pragma omp simd is part of OpenMP 4 but you can get it and only it from GCC and Intel via -fopenmp-simd and -qopenmp-simd, respectively. See https://gcc.gnu.org/onlinedocs/gcc-5.3.0/gcc/C-Dialect-Options.html for details on GCC.

Best,

Jeff

On Thu, Dec 8, 2016 at 7:05 AM, Francois Fayard <fayard@xxxxxxxxxxxxx> wrote:

(Side note - pretty sure you were the math teacher of my younger brother Pierre at Janson-de-Sailly ;-) )

Yes I was. The world is so small. I really remember him as he was both good in maths and a wonderful student to have. Donne-lui le bonjour.

For the scalar type, yes my intent was to have pointer aliasing. So I have used a matrix of Eigen::Index types. But it turns out that even with “double”, the following code does not vectorize:

=======================================================
#include <Eigen/Dense>

void f(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& B) {
Eigen::Index cols = B.cols();
Eigen::Index rows = B.rows();
for (Eigen::Index j = 0; j < cols; ++j) {
for (Eigen::Index i = 0; i < rows; ++i) {
B(i, j) = 1.0 / (1 + i + j);
}
}
}
=======================================================

I don’t know Eigen enough, but I fear that the code "const Derived& derived() const { return *static_cast<const Derived*>(this) }” does not help the compiler in analyzing possible pointer aliasing..

So the problem seems to come with any type.

François Fayard
Founder & Consultant - Inside Loop
Applied Mathematics & High Performance Computing
Tel: +33 (0)6 01 44 06 93
Web: www.insideloop.io

On Dec 8, 2016, at 3:45 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

2016-12-08 9:33 GMT-05:00 Francois Fayard <fayard@xxxxxxxxxxxxx>:

On Dec 8, 2016, at 3:21 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Interesting analysis!
I would draw a different conclusion from it though: compiler auto-vectorization is hopeless :-)
Which is why explicitly-vectorized code like Eigen exists :-)

Cheers,
Benoit

Hi Benoit.

Automatic vectorization is difficult, but not hopeless. I agree that the compiler needs help from time to time, especially in C/C++ because of pointer aliasing. But Intel compilers do a lot with automatic vectorization. Since you look french, here is a video I have done for Intel showing the optimization of a code that relies completely on automatic vectorization : https://www.youtube.com/watch?v=pJbo6Z4NrgM&feature=youtu.be . Sorry, I don’t have any english version yet.

Thanks for the links! I'll watch it.

I'm serious though about how your analysis here is a pretty solid piece of evidence /against/ compiler auto-vectorization. Having to change the basic representation of the number or rows/columns to be pointers instead of integers would be a very intrusive change to make; regardless of the compatibility issue raised by Christoph, it simply isn't the way that the human programmer /wants/ to represent a matrix.

Also a question about your original email:

void f(Eigen::Matrix<Eigen::Index, Eigen::Dynamic, Eigen::Dynamic>& A)

here you set the Scalar type to be the same as the Index type. That's a highly unusual usage pattern. The Index type is a pointer-sized integer, while in most applications the Scalar type is either floating-point or a narrower integer type. Do I understand correctly that the whole issue being discussed here is specific to this special case?

If one wanted to ensure that the Scalar type is a different C++ type from the Index type even when they both are integer types of the same size and signedness, one could wrap that integer type into a class providing the right arithmetic operators, and use that scalar type with Eigen::Matrix --- which supports arbitrary scalar types given the right interface:

https://eigen.tuxfamily.org/dox-devel/TopicCustomizing_CustomScalar.html

François Fayard
Founder & Consultant - Inside Loop
Applied Mathematics & High Performance Computing
Tel: +33 (0)6 01 44 06 93
Web: www.insideloop.io

Jeff Hammond
jeff.science@xxxxxxxxx
http://jeffhammond.github.io/