使用可重定位设备代码构建 R 包

Question

我正在编写一个 R 包，它使用 Thrust 来处理内存分配并避免编写我自己的 CUDA 内核。

在某些情况下，我从设备代码而不是主机代码调用 cuBLAS 例程。这改变了编译要求。虽然代码使用下面的nvcc 命令编译，但可能需要显式调用主机链接器（g++）。我怎样才能修改当前的构建过程来完成这个？

我使用的步骤是：

1. 使用 -dc 开关编译包含设备可重定位代码的输出文件 (max.o)

2. 创建要链接的库 (libmax.a)

3. 使用 -c 开关编译包含包装函数 (somePackage.o) 的输出文件

4. 使用-shared 开关创建链接到libmax.a 的共享库(somePackage.so)

如下所示的工作示例：

iterator.h：它定义了一些类型，包括strideAccessor。

max.h: max.cu 中的函数声明

max.cu：定义一个函数，该函数在每个 n 维数为 d 的串联数组中查找最大元素的索引。

somePackage.cu：处理 R/C++ 接口的包装器

$ cat iterator.h
#ifndef ITER_H
#define ITER_H

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/tuple.h>
#include <thrust/iterator/zip_iterator.h>

typedef thrust::device_vector<int> ivec_d;
typedef thrust::device_vector<double> fvec_d;
typedef thrust::device_vector<int>::iterator intIter;
typedef thrust::device_vector<double>::iterator realIter;
typedef thrust::host_vector<int> ivec_h;
typedef thrust::host_vector<double> fvec_h;

typedef thrust::counting_iterator<int> countIter;

//Used for generating rep( (1:len)*incr, times=infinity)
struct stride: public thrust::unary_function<int, int>{

  int incr;

  __host__ __device__ stride(int incr=1): incr(incr){}

  __host__ __device__ int operator()(int x){

    return x*incr;
  }
};

typedef thrust::transform_iterator<stride, countIter> strideIter;
typedef thrust::permutation_iterator<realIter, strideIter> strideAccessor;


#endif

$ cat max.h
#include "iterator.h"

void cublas_max(fvec_d &x, ivec_d &result, int n, int d);

$ cat max.cu
#include "iterator.h"
#include <thrust/functional.h>
#include <thrust/transform.h>
#include <cublas_v2.h>
#include <iostream>

struct whichMax : thrust::unary_function<double, int>{
  int dim;

  __host__ __device__ whichMax(int dim): dim(dim){}

  __host__ __device__ int operator()(double &vec){

    cublasHandle_t handle;
    cublasCreate_v2(&handle);
    int incx=1, n = dim, result =0;
    double *vec_ptr = thrust::raw_pointer_cast(&vec);

    //find the first index of a maximal element
    cublasIdamax(handle, n, vec_ptr, incx, &result);
    cublasDestroy_v2(handle);
    return result;
  }
};

void cublas_max(fvec_d &x, ivec_d &result, int n, int d){

  stride f(d);
  strideIter siter = thrust::transform_iterator<stride, countIter>(thrust::make_counting_iterator<int>(0), f);
  strideAccessor stridex = thrust::permutation_iterator<realIter, strideIter>(x.begin(), siter);

  whichMax g(d);

  //find the index of maximum for each of n subvectors
  thrust::copy(result.begin(), result.end(), std::ostream_iterator<int>(std::cout, " "));
  std::cout << std::endl;
  thrust::transform(stridex, stridex + n, result.begin(),  g);
  thrust::copy(result.begin(), result.end(), std::ostream_iterator<int>(std::cout, " "));
  std::cout << std::endl;
}

$ cat somePackage.cu
#include "iterator.h"
#include "max.h"
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include <iostream>

extern "C" SEXP Rcublas_max(SEXP x, SEXP n, SEXP dim){

  double *xptr = REAL(x);
  int N = INTEGER(n)[0], D = INTEGER(n)[0];

  fvec_d dx(xptr, xptr+N*D);
  ivec_d dresult(N);

  cublas_max(dx, dresult, N, D);

  ivec_h hresult(N);
  thrust::copy(dresult.begin(), dresult.end(), hresult.begin());

  SEXP indices = PROTECT(allocVector(INTSXP, N));

  for(int i=0; i<N; ++i)
    INTEGER(indices)[i] = hresult[i];

  UNPROTECT(1);
  return indices;
}

$ make
nvcc -dc -arch=sm_35 -Xcompiler -fPIC -lcublas_device -lcublas_device max.cu -o max.o
nvcc -lib -arch=sm_35 -Xcompiler -fPIC -lcublas_device -lcublas_device max.o -o libmax.a
nvcc -c -arch=sm_35 -Xcompiler -fPIC -lcublas_device somePackage.cu -lmax -I/home/emittman/src/R-3.3.1/builddir/include -I. -o somePackage.o
nvcc -shared -arch=sm_35 -Xcompiler -fPIC -lcublas_device somePackage.o -I/home/emittman/src/R-3.3.1/builddir/include -I. -L. -lcublas_device -lmax -o somePackage.so
ptxas info    : 'device-function-maxrregcount' is a BETA feature

Answer 1

我用 Rcpp 创建了一个 R 包，从 C++ 共享库中调用一些外部函数，然后调用 CUDA 内核来执行所需的计算。

您在这里尝试做的是将您的 CUDA 代码编译成静态库，然后将其链接到您的 R 包（它本身将被编译成一个共享库）。我的方法与你的不同，我对我的方法进行描述只是为了给你一个不同的想法。

这是一个简化的例子。

包含CUDA代码的共享库的kernels.cu：

__global__
void my_cuda_kernel( ... ) {
    // ......
}

包含CUDA代码的共享库的main.cu：

extern "C" {
    void do_cuda_work( ... ) {
        thrust::copy( ... );
        my_cuda_kernel <<< ... >>> ( ... );
    }
}

R包中的package.cpp：

extern void do_cuda_work( ... );

// [[Rcpp::export]]
void call_cuda_code( ... ) {
    do_cuda_work( ... );
}

要将CUDA代码编译成共享库，需要使用：

nvcc -arch=sm_35 -dc ... kernels.cu -o kernels.o
nvcc -arch=sm_35 -dc ... main.cu -o main.o
nvcc --shared -arch=sm_35 ... kernels.o main.o ... libMyCUDALibrary.so

请注意，要使单独编译工作，您需要为编译器和链接器指定-arch=sm_35，为编译器指定-dc。成功创建共享库后，将 R 包链接到它是相当简单的。 （但是，您可能需要在 R 包的 src 文件夹下创建一个 Makevars 文件来指定包含和库路径，也可能是 RPATH）：

CXX_STD= CXX11
PKG_CPPFLAGS= -I../../../CPP/include
PKG_LIBS= -L../../../CPP/bin/Release -lMyCUDALibrary -Wl,-rpath=$$HOME/MyCUDALibrary/CPP/bin/Release `$(R_HOME)/bin/Rscript -e "Rcpp::LdFlags()"`