Hello Eigen community,
My name is Maria Zhukova and I’m a software development engineer at Intel ® MKL Sparse team.
My team is interested in contributing into Eigen, so I’ve investigated our possibilities and so far this is what I have:
Eigen support different operations for sparse matrices stored in CSR and CSC format which can be implemented on a basis of IE SpBLAS kernels (please, refer to
https://software.intel.com/en-us/mkl-developer-reference-c-inspector-executor-sparse-blas-routines for the general idea of interfaces)
, basically we want to implement calls to our IE SpBLAS into next operations:
SparseMatrix + SparseMatrix (mkl_sparse_?_add)
SparseMatrix * DenseVector (mkl_sparse_?_mv)
SparseMatrix * DenseMatrix (mkl_sparse_?_mm)
SparseMatrix * SparseMatrix (mkl_sparse_spmm),
and Triangular solve (mkl_sparse_?_trsv).
I’ve already started with implementation of sparse_time_dense_impl_mkl kernel which is based on mkl_sparse_?_mv (included in patch).
This is how it will look like for user:
#include <Eigen/SpBLASSupport>
<--
NEW: IE SpBLAS include module
void main () {
SparseMatrix<double, RowMajor> A;
Matrix<double, Dynamic, 1> x, y;
A.makeCompressed(); /* Convert matrix A into CSR/CSC format */
A.createSparseHandle();
/*
NEW: is used to create handle required for all IE SpBLAS routines */
// support of IE SpBLAS is here
y = beta*y + alpha*A*x;
/* call to mkl_sparse_?_mv with operation = SPARSE_OPERATION_NON_TRANSPOSE */
y = beta*y + alpha*A.transpose()*x;
/* call to mkl_sparse_?_mv with operation = SPARSE_OPERATION_TRANSPOSE */
y = beta*y + alpha*A.adjoint()*x;
/* call to mkl_sparse_?_mv with operation = SPARSE_OPERATION_CONJUGATE_TRANSPOSE */
A.destroySparseHandle(); /*
NEW: is used to delete created handle */
}
I’ve attached a draft patch including all necessary changes and would like to hear your feedback.
Please, let me know if you have any questions and comments.
Best regards,
Maria