Eigen JacobiSVD cuda 编译错误

Question

关于在我的 cuda 函数中调用 JacobiSVD 时出现错误。

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

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

这是错误信息。

CUDA_voxel_building.cu(43)：错误：调用 __host__ 来自 __global__ 的函数（“Eigen::JacobiSVD, (int)2> ::JacobiSVD”） 不允许使用函数（“内核”）

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

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

我在 ubuntu 16.04 上使用带有 eigen3 的代码 8.0。 似乎特征值分解等其他函数也给出了相同的错误。

有人知道解决方案吗？我在下面附上了我的代码。

//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，其实不推荐使用。我建议升级到更新的版本。此外，要包含 Eigen，请编写（例如）：

#include <Eigen/Eigenvalues>

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

然后，正如 talonmies 所指出的，你不能从内核调用主机函数，（我现在不确定，为什么 JacobiSVD 没有被标记为设备函数），但在你的情况下它会做很多无论如何，使用Eigen::SelfAdjointEigenSolver 更有意义。由于您要分解的矩阵是固定大小的 3x3，因此您实际上应该使用优化的 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. 您可以通过直接将Eigen::Matrix3d sum2 和Eigen::Vector3d sum 与sum2 += v*v.adjoint() 和sum +=v 累加并计算来避免'大'矩阵local_p 和centered

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