tst-arma-spmat.R

library(Rcpp)
library(Matrix)

n = 100000
m = 10000
nnz = 0.001 * n * m
set.seed(1)

x = sparseMatrix(i = sample(n, nnz, T), j = sample(m, nnz, T), x = 1, dims = c(n, m))
i = sample(n, nnz * 10, T)
j = sample(m, nnz * 10, T) 

install.packages("~/Downloads/RcppArmadillo_0.7.960.1.2.tar.gz", repos = NULL, type = "source")
sourceCpp("~/Downloads/tst-arma.cpp", rebuild = T)

system.time(temp <- test_spmat(x, i, j, 1))
# user  system elapsed 
# 0.568   0.003   0.572
temp
# 9830
system.time(temp <- test_spmat(x, i, j, 4))
# user  system elapsed 
# 0.636   0.004   0.164 
temp
# 9830

install.packages("~/Downloads/RcppArmadillo_0.8.100.1.0.tar.gz", repos = NULL, type = "source")
sourceCpp("~/Downloads/tst-arma.cpp", rebuild = T)

system.time(temp <- test_spmat(x, i, j, 1))
# user  system elapsed 
# 6.199   0.037   6.253 
temp
# 9830

# this one crash R session
system.time(temp <- test_spmat(x, i, j, 4))

tst-arma.cpp

#include <RcppArmadillo.h>
#include <queue>
#include <iostream>
#include <vector>
#ifdef _OPENMP
#include <omp.h>
#endif

#define GRAIN_SIZE 10

using namespace Rcpp;
using namespace RcppArmadillo;
using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
double test_spmat(const arma::sp_mat &x, IntegerVector I, IntegerVector J, int n_threads) {
  int *i_ptr = I.begin();
  int *j_ptr = J.begin();

  double sum = 0;
  #ifdef _OPENMP
  #pragma omp parallel for num_threads(n_threads) schedule(dynamic, GRAIN_SIZE) reduction(+:sum)
  #endif
  for(int k = 0; k < J.size(); k++) {
    //adjust to 0-based indexes
    int i = i_ptr[k] - 1;
    int j = j_ptr[k] - 1;
    sum += x(i, j);
  }
  return(sum);
}