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 上使用带有 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

首先，Ubuntu16.04提供了Eigen 3.3-beta1，不太推荐使用。我建议升级到更新的版本。此外，要包含 Eigen，请写（例如）：

#include <Eigen/Eigenvalues>

并用-I /usr/include/eigen3编译（如果你使用OS提供的版本），或者更好-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 提供的测试版上运行（我仍然建议更新）。

一些不相关的笔记：

以下是错误的，因为您应该从索引 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);

这条线不适合（除非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() 第二个因素，而不是第一个因素。
你可以避免'big'矩阵local_p和centered，直接用[=31累加Eigen::Matrix3d sum2和Eigen::Vector3d sum =] 和 sum +=v 以及计算
```
Eigen::Vector3d mu = sum / diff;
Eigen::Matrix3d cov_e = (sum2 - mu*mu.adjoint()*diff)/(diff-1);
```

Eigen JacobiSVD cuda编译错误

Eigen JacobiSVD cuda compile error

nvcc

eigen