c++ - RcppParallel 中二次型的并行计算

Question

我想编写代码使用 RcppParallel 计算 v^T Av。这里 v 是一个大小为 n 的 vector ，A 是一个 n×n 矩阵。我的想法是以并行方式评估 Av，然后计算该 vector 与 v 的内积。这是我的代码:

#include <Rcpp.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include<numeric>

using namespace std;
using namespace Rcpp;
using namespace RcppParallel;

struct Par_MatVec_Mult: public Worker {

  const RMatrix<double> Mat;
  const vector<double> Vec;
  vector<double> output;

  Par_MatVec_Mult(RMatrix<double> A, vector<double> v, vector<double> Av): \ 
  vector<double> Av): Mat(A), Vec(v), output(Av) { }
  void operator() ( size_t begin, size_t end) {

    for( size_t i = begin; i < end; i++ ) {
      RMatrix<double>::Row rowi = Mat.row(i);
       output.at(i) = inner_product( rowi.begin(), rowi.end(),Vec.begin(), \
       0.0 );
    }
  }
};

// [[Rcpp::export]]
double Parallel_MatVec_Mult( NumericMatrix A, vector<double> v ){

  vector<double> b( A.nrow(), 0.0 );
  Par_MatVec_Mult Av(A, v, b);
  parallelFor(0, A.nrow(), Av);
  return inner_product( Av.output.begin(), Av.output.end(), v.begin(), 0.0   );
}

编译时出现以下错误。

no known conversion for argument 1 from 'Rcpp::NumericMatrix {akaRcpp::Matrix<14>}' 
to 'RcppParallel::RMatrix<double>'

我很清楚 NumericMatrix 和 RMatrix 是两个不同的对象。但是我不知道它们到底有什么区别，以及如何更改我的代码以消除此错误。

顺便说一句，我在 Windows 10 上使用 RStudio 0.99.903。

Answer 1

single-argument constructor RMatrix 被标记为 explicit:

inline explicit RMatrix(const Source& source) 
   : data_(const_cast<Source&>(source).begin()),
     nrow_(source.nrow()),
     ncol_(source.ncol())
{}

您已为 Par_MatVec_Mult 构造函数提供此签名

Par_MatVec_Mult(RMatrix<double> A, vector<double> v, vector<double> Av)

并试图稍后将 NumericMatrix 传递给它。这将需要在 RMatrix 的构造中进行隐式转换，但由于原本合适的构造函数被标记为 explicit，因此这是不允许的，并且会出现错误.

这是一个简单的例子来证明这一点:

#include <Rcpp.h>

class Wrapper {
private:
    std::size_t nr;
    std::size_t nc;

public:
    template <typename T>
    explicit Wrapper(const T& t)
        : nr(t.nrow()),
          nc(t.ncol())
    {}

    std::size_t size() const { return nr * nc; }
};

// [[Rcpp::export]]
int test(Rcpp::NumericMatrix x) {
    // Wrapper w = x;   // Error: conversion from ‘Rcpp::NumericMatrix
                        // {aka Rcpp::Matrix<14>}’ to non-scalar
                        // type ‘Wrapper’ requested

    Wrapper w(x);       // Ok

    return w.size();
}

---

在您的情况下，修复很简单:只需将 Par_MatVec_Mult 的构造函数的签名更改为:

Par_MatVec_Mult(NumericMatrix A, vector<double> v, vector<double> Av)