CODEWITHSUNDEEP
c++11sparse-matrixeigenintel-mkl
DrStS
DrStS

Reputation: 21

Use Intel MKL Sparse BLAS extensions through C++ Eigen library

I would like to use Intel MKL Sparse BLAS through the Eigen C++ library. I checked the documentation of Eigen https://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html. It seems that this is not applicable for Intel MKL Sparse BLAS. The only MKL Sparse example which I found is the Intel MKL PARDISO interface.

Therefore, I would like to ask the question if this is supported?

Thank a lot!

Upvotes: 1

Views: 860

Answers (1)

Charlie S
Charlie S

Reputation: 315

MKL Sparse BLAS routines are not currently supported in Eigen. However, someone started integrating the newer MKL Inspector-Executor interface, though it hasn't been active since 2018. I haven't tried it myself.

https://bitbucket.org/eigen/eigen/pull-requests/473/add-support-of-intel-r-mkl-inspector/diff

Alternatively, you can create an MKL sparse matrix handle for your Eigen matrix and do some MKL stuff on it. My solution isn't elegant, but its good enough for tests.

static void EigenToMKLSparseCSR(const SparseMatrix<double> &A_Eigen, sparse_matrix_t &A_MKL)
{
    int m = (int)A_Eigen.rows();
    int n = (int)A_Eigen.cols();
    double *values = (double*)A_Eigen.valuePtr();

    if (A_Eigen.IsRowMajor)
    {
        int *rows_start = (int*)A_Eigen.outerIndexPtr();
        int *rows_end = (int*)A_Eigen.outerIndexPtr() + 1;
        int *col_indx = (int*)A_Eigen.innerIndexPtr();
        mkl_sparse_d_create_csr(&A_MKL, SPARSE_INDEX_BASE_ZERO, m, n, rows_start, rows_end, col_indx, values);
        mkl_sparse
    }
    else
    {
        int *cols_start = (int*)A_Eigen.outerIndexPtr();
        int *cols_end = (int*)A_Eigen.outerIndexPtr() + 1;
        int *row_indx = (int*)A_Eigen.innerIndexPtr();
        mkl_sparse_d_create_csc(&A_MKL, SPARSE_INDEX_BASE_ZERO, m, n, cols_start, cols_end, row_indx, values);
        mkl_sparse_convert_csr(A_MKL, SPARSE_OPERATION_NON_TRANSPOSE, &A_MKL);
    }
    mkl_sparse_set_memory_hint(A_MKL, SPARSE_MEMORY_AGGRESSIVE);
    mkl_sparse_optimize(A_MKL);
}

Note that I am assuming that your scalar type is a double. You could probably make it a template that accepts float, double, MKL_Complex8, and MKL_Complex16. I honestly have no idea what SPARSE_MEMORY_AGGRESSIVE does, but I threw it in there for good measure.

Note that Eigen defaults to column-major storage, whereas MKL (mostly) requires you to use row-major. Ironically, this routine does NOT work with the intel compiler! It does work with Visual Studio 2017.

Upvotes: 2

Related Questions