Eigen JacobiSVD cuda编译错误

Question

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

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

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

这是错误信息。

CUDA_voxel_building.cu（43）：错误：不允许从__global__函数（“内核”）调用__host__函数（“Eigen :: JacobiSVD，（int）2> :: JacobiSVD”）

我使用以下命令编译它。

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，这并不是真正推荐使用的。我建议升级到更新的版本。此外，要包括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. 你可以通过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);