Потоково-параллельное умножение параллельных разреженных матриц с Eigen и OpenMP

Question

Согласно документам , Eigen поддерживает многопоточные плотные v. Плотные или мажорные строки, разреженные v. Плотные, но не разреженные v. Разреженные умножения матриц. В моем приложении у меня есть 2 разреженных матрицы, и я хочу умножить их параллельно, то есть получить C = A * B, где C и B - основной столбец.

[Редактировать:] Для моего приложения я также заранее знаю шаблон разреженности всех матриц, и мне нужно многократно выполнять такое умножение матриц (каждый раз при обновлении ненулевого значения). значения C с использованием обновленных A и B). Это заставляет меня думать, что в этом случае должен существовать поточно-безопасный способ выполнения этой операции, поскольку фиксированное выделение памяти C [end edit]

Учитывая, что Eigenдает быстрое чтение / запись для блоков столбцов для главных разреженных матриц столбцов, я могу выполнить операцию путем подстановки B в блоки столбцов (предположим, что A, B и C все были инициализированы и имеют правильные значенияsize):

  int nchunks = 7; 
  int chunksize = B.cols() / nchunks;

#pragma omp parallel for
  for(int i=0; i<nchunks; i++){
    if(i == nchunks-1){
      C.rightCols(B.cols() - (nchunks-1)*chunksize) = A * B.rightCols(B.cols() - (nchunks-1)*chunksize);
    } else {
      C.middleCols(i*chunksize, chunksize) = A * B.middleCols(i*chunksize, chunksize);
    }
  }

В моем приложении я должен выполнить A.transpose() * S * A тысячи раз, где S симметрично положительно определен, и я могу разбить эту операцию на 2 шага, каждый из которых выполняется параллельно. Результат должен быть симметричным. Тем самым я получаю значительное ускорение.

На мой взгляд, это не должно быть проблемой, потому что все потоки имеют общий C, но каждый из них работает со своим набором столбцов. На самом деле я ожидаю, что мое приложение будет чрезвычайно чувствительным даже к крошечным проблемам, потому что оно должно поддерживать симметрию после выполнения 2 таких разреженных умножений матриц.

Однако этот код работает в Linux (Ubuntu 18.04, R 3.6.1 связан с Intel MKL), но не в Windows (я пробовал оба с vanilla R 3.6.1и MRO 3.5.3). В некоторых случаях в Windows происходит сбой моего приложения из-за сбоя Cholesky, что, в свою очередь, вызвано неправильным созданием C, что я проследил до указанного выше фрагмента кода. Простое удаление OMP решает проблему, конечно, за счет производительности. Я не смог воспроизвести это последовательно в Windows: в отдельности продукт работает правильно. Но это определенно проблема OMP, потому что удаление параллельного предложения решает проблему.

Так что же происходит?

Вот рабочий фрагмент кода

//[[Rcpp::export]]
Rcpp::List spmat_mult(const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
                      const Eigen::SparseMatrix<double>& B){

  int nchunks = 7;
  int chunksize = B.cols() / nchunks;

  std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
  Eigen::SparseMatrix<double> Cs(A.rows(), B.cols());
  Cs = A * B;

  std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
  clog << "Standard "
       << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
       << "us.\n";

  Eigen::SparseMatrix<double> C(A.rows(), B.cols());
  C = Cs;

  start = std::chrono::steady_clock::now();
#pragma omp parallel for
  for(int i=0; i<nchunks; i++){
    if(i == nchunks-1){
      C.rightCols(B.cols() - (nchunks-1)*chunksize) = A * B.rightCols(B.cols() - (nchunks-1)*chunksize);
    } else {
      C.middleCols(i*chunksize, chunksize) = A * B.middleCols(i*chunksize, chunksize);
    }
  }

  end = std::chrono::steady_clock::now();
  clog << "OMP "
       << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
       << "us.\n";

  return Rcpp::List::create(
    Rcpp::Named("Cs") = Cs,
    Rcpp::Named("C") = C
  );
}

и R

  library(Matrix)

  set.seed(1)

  n <- 5000 # choose size
  A <- B <- matrix(0, ncol=n, nrow=n)

  nnz <- 250
  A_rows <- sample(1:n, nnz, replace=F)
  A_cols <- sample(1:n, nnz, replace=F)
  B_rows <- sample(1:n, nnz, replace=F)
  B_cols <- sample(1:n, nnz, replace=F)

  A[A_rows, A_cols] <- apply(A[A_rows, A_cols], 1:2, function(x) runif(1))
  B[B_rows, B_cols] <- apply(B[B_rows, B_cols], 1:2, function(x) runif(1))

  A <- as(A, "dgRMatrix")
  B <- as(B, "dgCMatrix")

  C <- spmat_mult(A, B)