并行计算的二次形式在RcppParallel

Parallel computation of a quadratic form in RcppParallel

本文关键字：RcppParallel 二次 并行计算      更新时间：2023-10-16

我想写一个代码来计算v^T Av使用RcppParallel。这里v是一个大小为n的向量a是一个n × n的矩阵。我的想法是以并行方式计算Av，然后计算这个向量与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。

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)
相关文章：