Eigen JacobiSVD cuda编译错误

Question

强调文本我有一个错误，关于在我的cuda函数中调用JacobiSVD。

这是导致错误的代码的一部分。

Eigen::JacobiSVD<Eigen::Matrix3d> svd(cov_e, Eigen::ComputeThinU | Eigen::ComputeThinV); 

并且这是错误消息。

CUDA_voxel_building.cu(43): error: calling a __host__ function("Eigen::JacobiSVD , (int)2> ::JacobiSVD") from a __global__ function("kernel") is not allowed

我用下面的命令来编译它。

nvcc -std=c++11 -D_MWAITXINTRIN_H_INCLUDED -D__STRICT_ANSI__ -ptx CUDA_voxel_building.cu 

我在ubuntu 16.04上使用代码8.0和eigen3。 似乎像其他功能，如特征值分解也给出了相同的错误。

任何人都知道解决方案吗？我在下面附上我的代码。

//nvcc -ptx CUDA_voxel_building.cu 
#include </usr/include/eigen3/Eigen/Core> 
#include </usr/include/eigen3/Eigen/SVD> 
/* 
#include </usr/include/eigen3/Eigen/Sparse> 

#include </usr/include/eigen3/Eigen/Dense> 
#include </usr/include/eigen3/Eigen/Eigenvalues> 

*/ 





__global__ void kernel(double *p, double *breaks,double *ind, double *mu, double *cov, double *e,double *v, int *n, char *isgood, int minpts, int maxgpu){ 
    bool debuginfo = false; 
    int idx = threadIdx.x + blockIdx.x * blockDim.x; 
    if(debuginfo)printf("Thread %d got pointer\n",idx); 
    if(idx < maxgpu){ 


     int s_ind = breaks[idx]; 
     int e_ind = breaks[idx+1]; 
     int diff = e_ind-s_ind; 

     if(diff >minpts){ 
      int cnt = 0; 
      Eigen::MatrixXd local_p(3,diff) ; 
      for(int k = s_ind;k<e_ind;k++){ 
       int temp_ind=ind[k]; 

       //Eigen::Matrix<double, 3, diff> local_p; 
       local_p(1,cnt) = p[temp_ind*3]; 
       local_p(2,cnt) = p[temp_ind*3+1]; 
       local_p(3,cnt) = p[temp_ind*3+2]; 
       cnt++; 
      } 

      Eigen::Matrix3d centered = local_p.rowwise() - local_p.colwise().mean(); 
      Eigen::Matrix3d cov_e = (centered.adjoint() * centered)/double(local_p.rows() - 1); 

      Eigen::JacobiSVD<Eigen::Matrix3d> svd(cov_e, Eigen::ComputeThinU | Eigen::ComputeThinV); 
    /*   Eigen::Matrix3d Cp = svd.matrixU() * svd.singularValues().asDiagonal() * svd.matrixV().transpose(); 


      mu[idx]=p[ind[s_ind]*3]; 
      mu[idx+1]=p[ind[s_ind+1]*3]; 
      mu[idx+2]=p[ind[s_ind+2]*3]; 

      e[idx]=svd.singularValues()(0); 
      e[idx+1]=svd.singularValues()(1); 
      e[idx+2]=svd.singularValues()(2); 

      n[idx] = diff; 
      isgood[idx] = 1; 

      for(int x = 0; x < 3; x++) 
      { 
       for(int y = 0; y < 3; y++) 
       { 
        v[x+ 3*y +idx*9]=svd.matrixV()(x, y); 
        cov[x+ 3*y +idx*9]=cov_e(x, y); 
        //if(debuginfo)printf("%f ",R[x+ 3*y +i*9]); 
        if(debuginfo)printf("%f ",Rm(x, y)); 
       } 
      } 
*/ 

     } else { 
      mu[idx]=0; 
      mu[idx+1]=0; 
      mu[idx+2]=0; 

      e[idx]=0; 
      e[idx+1]=0; 
      e[idx+2]=0; 

      n[idx] = 0; 
      isgood[idx] = 0; 

      for(int x = 0; x < 3; x++) 
      { 
       for(int y = 0; y < 3; y++) 
       { 
        v[x+ 3*y +idx*9]=0; 
        cov[x+ 3*y +idx*9]=0; 
       } 
      } 
     } 




    } 
} 

Answer 1

首先，Ubuntu 16.04提供了Eigen 3.3-beta1，这并不真正推荐使用。我会建议升级到更新的版本。此外，包括艾根，写（例如）：

#include <Eigen/Eigenvalues> 

与-I /usr/include/eigen3编译（如果使用操作系统提供的版本），或更好的-I /path/to/local/eigen-version。

然后，正如talonmies所说，你不能从内核调用主机函数，（我不确定此刻，为什么JacobiSVD没有标记为设备函数），但在你的情况下，它会做得更多无论如何，感觉要使用Eigen::SelfAdjointEigenSolver。既然你是分解矩阵是固定大小的3×3，你应该实际使用的优化computeDirect方法：

Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig; // default constructor 
eig.computeDirect(cov_e); // works for 2x2 and 3x3 matrices, does not require loops 

看来computeDirect甚至适用于被Ubuntu提供的测试版（我还是建议更新）。

一些不相关的注意事项：

1. 下列哪项是错误的，因为你应该与指数0开始：

   local_p(1,cnt) = p[temp_ind*3]; 
   local_p(2,cnt) = p[temp_ind*3+1]; 
   local_p(3,cnt) = p[temp_ind*3+2]; 
   

   此外，您还可以在一行写：

   local_p.col(cnt) = Eigen::Vector3d::Map(p+temp_ind*3); 
   

2. 此行不适合（除非diff==3）：

   Eigen::Matrix3d centered = local_p.rowwise() - local_p.colwise().mean(); 
   

   什么你大概的意思是（local_p实际上是3XN不NX3）

   Eigen::Matrix<double, 3, Eigen::Dynamic> centered = local_p.colwise() - local_p.rowwise().mean(); 
   

   和计算cov_e当你需要.adjoint()第二个因素，不是第一次。

3. 你可以同时避免 '大' 矩阵local_p和centered，直接累加Eigen::Matrix3d sum2和Eigen::Vector3d sum与sum2 += v*v.adjoint()和sum +=v和计算

   Eigen::Vector3d mu = sum/diff; 
   Eigen::Matrix3d cov_e = (sum2 - mu*mu.adjoint()*diff)/(diff-1);