RE: [eigen] Intel (R) MKL IE SpBLAS support in Eigen

["Zhukova, Maria" <maria.zhukova@xxxxxxxxx>]

Hello Edward and thanks much for your feedback!

This is a good point. 
Actually, the correct place for creating IE SpBLAS handle is one of the key moments which I want to discuss.
The thing is -- all Intel (R) MKL IE SpBLAS functionality requires matrix to be stored in compressed row/column storage format. Necessary storage format can be obtained by calling SparseMatrix::makeCompressed() function. 
So, I wonder if the better decision will be to create sparse handle inside the makeCompressed() and destroy in SparseMatrix destructor. What do you think about this option?

Best regards,
Maria

-----Original Message-----
From: Edward Lam [mailto:edward@xxxxxxxxxx] 
Sent: Tuesday, April 3, 2018 4:28 PM
To: eigen@xxxxxxxxxxxxxxxxxxx
Subject: Re: [eigen] Intel (R) MKL IE SpBLAS support in Eigen

Hi Maria,

I'm only an Eigen user but I'm extremely happy at this effort to have more MKL wrapping in Eigen! As our application depends on Intel TBB, the only real way to use Eigen with multithreaded support is via MKL.

 > A.createSparseHandle(); /* *NEW*: is used to create handle required for all IE  > SpBLAS routines */  > ...
 > A.destroySparseHandle(); /* *NEW*: is used to delete created handle */

I think that this is very un-C++ like. I'm not sure why it was decided to be done this way. Why wouldn't the sparse handle be allocated within A itself and then deallocated during A's constructor? Note also that in C++, we need to worry about what happens when the matrix is copied so this state management needs to be within the matrix itself or else the semantics of these create/destroy handles are very troublesome.

Best Regards,
-Edward

On 4/3/2018 5:39 PM, Zhukova, Maria wrote:
> 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-e
> xecutor-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
>