Eigen Iterate over inner iterator in sparse matrix with SIMD

Question

Apply simd over eigen inner iterator of sparse matrix:

for(auto i = 0; i < smat.outerSize(); i++){
    #pragma omp simd
    for(SMat::InnerIterator iter(smat,i); it; ++it){
        it.valueRef() = value;
    }
}

This does not work due to an error with parenthetical initialization in for loop is incompatible with simd. Next I try:

    SMat::InnerIterator iter(smat,i);
    #pragma omp simd
    for(;it;++it){ // error, declaration or initialization expected

    for(it;it;++it){ // error, declaration or initialization expected

Then I google and search the documentation, only to encounter a phrase mentioning that simd is implicit when adding sparse matrices (so I know this is possible, and that somewhere in the templated bowels of eigen, there is a simd loop over the inner vector; but I don't know how to do it).

Next, I check and discover that Eigen has only three calls to omp in the entire code. Does this mean that Eigen is reliant only on compiler flags for simd activation?

Finally, I try changing the loop to canonical form (per comment below), and get a different error:

for(auto it = typename SMat::InnerIterator(smat,i); it; ++it)

// error: '#pragma omp simd' used with class iteration variable 'it'

What is the expected way to trigger or iterate over the inner vector in an Eigen::SparseMatrix<double> with simd?

What version of OpenMP does your compiler support? Anyway, in OpenMP, loops have to be in so-called canonical form, which should in your case look likefor (auto it = Eigen::SparseMatrix<double>::InnerIterator(smat, i); it != false; ++it). — Daniel Langr
SIMD in Eigen is almost entirely done using compiler intrinsics. Where did you find "that simd is implicit when adding sparse matrices"? (Depends on what you are actually adding, but especially adding two sparse matrices with different non-zero patterns, this sounds difficult. — chtz
InnerIterator is indeed not documented at the moment. It could be made compatible to a standard random-access iterator (it lacks some operators like it-diff, --it, it-=diff). For reference, the source is here: bitbucket.org/eigen/eigen/src/0fdf0e56237f92e5/Eigen/src/… — chtz
@chtz: InnerIterator is mentioned in eigen.tuxfamily.org/dox/group__TutorialSparse.html, but that's a tutorial and doesn't detail the type properties. — Peter Cordes
@DanielLangr the canonical loop form does not allow for != as an operator. — Zulan

Zulan Zulan · Accepted Answer · 2018-06-06T10:32:51

It is impossible to apply #pragma omp simd in this context. As per the OpenMP specifications (2.6 Canonical Loop Form), "in the simd construct the only random access iterator types that are allowed [...] are pointer types.". The involved iterators are clearly not pointer types. It might be possible to change them to such to allow an OpenMP simd loop, but that would require insight into the implementations and data layout of the involved types.

Eigen Iterate over inner iterator in sparse matrix with SIMD

2 Answers