How to work with Eigen in CUDA kernels

Question

Eigen is a c++ linear algebra library http://eigen.tuxfamily.org.

It's easy to work with basic data types, like basic float arrays, and just copy it to device memory and pass the pointer to cuda kernels. But Eigen matrix are complex type so how to copy it to device memory and let cuda kernels read/write with it?

Answer 1

Since November 2016 (Release of Eigen 3.3), a new option exists: Using Eigen directly inside CUDA kernels - see this question.

Example from the linked question:

__global__ void cu_dot(Eigen::Vector3f *v1, Eigen::Vector3f *v2, double *out, size_t N)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx < N)
    {
        out[idx] = v1[idx].dot(v2[idx]);
    }
    return;
}

Copying an array of Eigen::Vector3f to device:

Eigen::Vector3f *host_vectors = new Eigen::Vector3f[N];
Eigen::Vector3f *dev_vectors;
cudaMalloc((void **)&dev_vectors, sizeof(Eigen::Vector3f)*N)
cudaMemcpy(dev_vectors, host_vectors, sizeof(Eigen::Vector3f)*N, cudaMemcpyHostToDevice)

Answer 2

If all you want is to access the data of an Eigen::Matrix via a raw C pointer, then you can use the .data() function. Coefficient are stored sequentially in memory in a column major order by default, or row major if you asked for:

MatrixXd A(10,10);
double *A_data = A.data();

Answer 3

Apart from rewriting and refitting the code, there is an Eigen-compatible library written as a byproduct of a research project that performs matrix calculations on the GPU, and you can use multiple backends: https://github.com/rudaoshi/gpumatrix

I cannot vouch for it, but if it works it is probably exactly what you are looking for.

If you'd like a more general-purpose solution, this thread seems to contain very useful info

Answer 4

There are two ways.

Make eigen work on GPU, which probably is hard and won't perform good. At least if work on GPU means only get it to compile and produce results. Eigen is practically hand optimized for modern CPUs. Internally Eigen uses its own allocators and memory layouts which are most probably not going to work well on CUDA.

2nd way is easier to do and should not break legacy Eigen code, and probaly is the only suitable in your case. Switch your underlying matrices to plain matrices (i.e. double**) use Eigen::Map. This way, you will have Eigen interface to plain data-type so the codes should not break, and you can send the matrix to the GPU as normal c-array, like its usually done. The drawback is that probably you won't utilize Eigen to the full potential, however if you offload most of work to GPU it's ok.

It's actually reversing things a bit. Instead of making Eigen arrays work on CUDA, you can make Eigen work on normal arrays.