Eigen JacobiSVD cuda compile error

I've got an error, regarding calling JacobiSVD in my cuda function.

This is the part of the code that causing the error.

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

And this is the error message.

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

I've used the following command to compile it.

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

I'm using code 8.0 with eigen3 on ubuntu 16.04. It seems like other functions such as eigen value decomposition also gives the same error.

Anyone knows a solution? I'm enclosing my code below.

//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;
                }
            }
        }




    }
}

Solution

First of all, Ubuntu 16.04 provides Eigen 3.3-beta1, which is not really recommended to be used. I would suggest upgrading to a more recent version. Furthermore, to include Eigen, write (e.g.):

#include <Eigen/Eigenvalues>

and compile with -I /usr/include/eigen3 (if you use the version provided by the OS), or better -I /path/to/local/eigen-version.

Then, as talonmies noted, you can't call host-functions from kernels, (I'm not sure at the moment, why JacobiSVD is not marked as device function), but in your case it would make much more sense to use Eigen::SelfAdjointEigenSolver, anyway. Since the matrix you are decomposing is fixed-size 3x3 you should actually use the optimized computeDirect method:

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

It seems the computeDirect even works on the beta version provided by Ubuntu (I'd still recommend to update).

Some unrelated notes:

The following is wrong, since you should start with index 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];

Also, you can write this in one line:

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

This line will not fit (unless diff==3):

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

What you probably mean is (local_p is actually 3xn not nx3)

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

And when computing cov_e you need to .adjoint() the second factor, not the first.

You can avoid both 'big' matrices local_p and centered, by directly accumulating Eigen::Matrix3d sum2 and Eigen::Vector3d sum with sum2 += v*v.adjoint() and sum +=v and computing
```
Eigen::Vector3d mu = sum / diff;
Eigen::Matrix3d cov_e = (sum2 - mu*mu.adjoint()*diff)/(diff-1);
```