1

I'm trying to port a function that transforms all non-zero entries in a sparsematrix using RcppEigen to run in parallel with RcppParallel but I can't get it to work.

I.e. with this source:

#include <RcppEigen.h>
#include <RcppParallel.h>

using namespace Rcpp;

typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Map<SpMat> MSpMat;

// [[Rcpp::export]]
S4 serial_test( S4 m ) {
    SpMat src = as<MSpMat>( m );
    SpMat tgt = as<MSpMat>( clone( m ) );
    for( int i = 0; i < src.rows(); ++i ) {
        for( SpMat::InnerIterator srcIt( src, i ), tgtIt( tgt, i ); srcIt; ++srcIt, ++tgtIt ) {
            tgtIt.valueRef() = -1;
        }
    }
    S4 out = wrap( tgt );
    return out;
}

struct MyWorker : public RcppParallel::Worker {
    SpMat src;
    SpMat tgt;

    MyWorker( const SpMat& src, SpMat& tgt )
        : src( src ), tgt( tgt ) {}

    void operator()( std::size_t begin, std::size_t end ) {
        for( int i = begin; i < end; ++i ) {
            for( SpMat::InnerIterator srcIt( src, i ), tgtIt( tgt, i ); srcIt; ++srcIt, ++tgtIt ) {
                tgtIt.valueRef() = -1;
            }
        }
    }
};

// [[Rcpp::export]]
S4 para_test( S4 m ) {
    SpMat src = as<MSpMat>( m );
    SpMat tgt = as<MSpMat>( clone( m ) );
    MyWorker wrkr( src, tgt );
    RcppParallel::parallelFor( 0, src.outerSize(), wrkr );
    S4 out = wrap( tgt );
    return out;
}

And this input:

> m <- Matrix::rsparsematrix( 1000, 1000, 0.1 )
> m[ 1:10, 1:10 ]
10 x 10 sparse Matrix of class "dgCMatrix"

 [1,] . . . . 0.52 . . .     . .
 [2,] . . . . 0.59 . . .     . .
 [3,] . . . . .    . . 0.47  . .
 [4,] . . . . .    . . 0.25  . .
 [5,] . . . . .    . . .     . .
 [6,] . . . . .    . . .     . .
 [7,] . . . . .    . . .     . .
 [8,] . . . . .    . . .     . .
 [9,] . . . . .    . . .    -1 .
[10,] . . . . .    . . .     . .

Calling serial_test works:

> m1 <- serial_test( m )
> m1[ 1:10, 1:10 ]
10 x 10 sparse Matrix of class "dgCMatrix"

 [1,] . . . . -1 . .  .  . .
 [2,] . . . . -1 . .  .  . .
 [3,] . . . .  . . . -1  . .
 [4,] . . . .  . . . -1  . .
 [5,] . . . .  . . .  .  . .
 [6,] . . . .  . . .  .  . .
 [7,] . . . .  . . .  .  . .
 [8,] . . . .  . . .  .  . .
 [9,] . . . .  . . .  . -1 .
[10,] . . . .  . . .  .  . .

But para_test does nothing, even though it runs fine:

> m2 <- para_test( m )
> m2[ 1:10, 1:10 ]
10 x 10 sparse Matrix of class "dgCMatrix"

 [1,] . . . . 0.52 . . .     . .
 [2,] . . . . 0.59 . . .     . .
 [3,] . . . . .    . . 0.47  . .
 [4,] . . . . .    . . 0.25  . .
 [5,] . . . . .    . . .     . .
 [6,] . . . . .    . . .     . .
 [7,] . . . . .    . . .     . .
 [8,] . . . . .    . . .     . .
 [9,] . . . . .    . . .    -1 .
[10,] . . . . .    . . .     . .

Is it possible to use RcppParallel with Eigen's SparseMatrix class?

jtatria
  • 527
  • 3
  • 12
  • 1
    I just poked around a little. There is something odd happening but I can't quite put my finger on it. Print `tgt` in the `worker` and you see it is visited and changed; I just don't quite see yet why the changes do not persist. – Dirk Eddelbuettel Aug 08 '17 at 12:48
  • As a hunch, can you try with dense matrices to see if maybe just maybe the non-contiguous layout of sparse matrices gets in the way? – Dirk Eddelbuettel Aug 08 '17 at 15:53
  • The same using Eigen's dense matrices works flawlessly, it doesn't even need RcppParallel::RMatrix. I am not sure what you mean by non-contigous, though... I thought the InnerIterator would take care of that and that the actual memory region would be contigous as it is basically mapping the values vector of a Matrix::dgCMatrix, but this is beyond my c++ knowledge... Thanks! – jtatria Aug 08 '17 at 18:22
  • I am not entirely sure either but there must be "something" else going on with the sparse representation. If I were you, I'd try a simple direct C++ example with OpenMP or Intel TBB over such a sparse matrix and see if that works. If not, back to Eigen experts. If it does, we can dig into what Rcpp may do here. – Dirk Eddelbuettel Aug 08 '17 at 18:27

0 Answers0