a***@fitzgibbon.ie
2018-04-05 14:54:29 UTC
It's great to see discussions like the MKL sparse one below, but may I make a meta-point that it would be very valuable to have it in an "issues" format, as is available on github e.g. https://github.com/JuliaLang/julia/issues/9147
Is it an option to do this for Eigen? Would it be [another] prod towards moving to github?
-----Original Message-----
From: Edward Lam <***@sidefx.com>
Sent: 05 April 2018 14:32
To: ***@lists.tuxfamily.org
Subject: Re: [eigen] Intel (R) MKL IE SpBLAS support in Eigen
Would it be useful to incorporate lambda's into the interface to avoid begin/end pairs? So from the user side, we would write code like this instead:
1) Analyze and run
SimplicialLDLT<MklSparseMatrix<double>> llt(A);
int it = 0;
for (int it = 0; ...; ++it) {
if (it == 0)
llt.matrixL().sparseAnalyzeAndRun(100, [&] { llt.solve(b); });
else
llt.solve(b);
// ...
}
2) Analyze only
SimplicialLDLT<MklSparseMatrix<double>> llt(A);
llt.matrixL().sparseAnalyze(100, [&] { llt.solve(b); });
// and solve as usual
For more complicated algorithms, one can always outline the lambda and pass it into the analysis.
Cheers,
-Edward
Is it an option to do this for Eigen? Would it be [another] prod towards moving to github?
-----Original Message-----
From: Edward Lam <***@sidefx.com>
Sent: 05 April 2018 14:32
To: ***@lists.tuxfamily.org
Subject: Re: [eigen] Intel (R) MKL IE SpBLAS support in Eigen
Would it be useful to incorporate lambda's into the interface to avoid begin/end pairs? So from the user side, we would write code like this instead:
1) Analyze and run
SimplicialLDLT<MklSparseMatrix<double>> llt(A);
int it = 0;
for (int it = 0; ...; ++it) {
if (it == 0)
llt.matrixL().sparseAnalyzeAndRun(100, [&] { llt.solve(b); });
else
llt.solve(b);
// ...
}
2) Analyze only
SimplicialLDLT<MklSparseMatrix<double>> llt(A);
llt.matrixL().sparseAnalyze(100, [&] { llt.solve(b); });
// and solve as usual
For more complicated algorithms, one can always outline the lambda and pass it into the analysis.
Cheers,
-Edward
Thank you for opening this discussion on the public mailing list.
So let's discuss about the public API, which currently is not very
(i1) - Storing MKL's handle in SparseMatrix breaks ABI and does not
sounds very generic.
(i2) - which operations are going to be analyzed/optimized,
(i3) - and specify the 'expected_calls' parameter.
(e.g., non-linear optimization, eigenvalues, ...)
SimplicialLDLT<SparseMatrix<double> > llt(A);
while(...) {
...
x = llt.solve(b);
...
}
Here the triangular L factor is going to be used for triangular and
transposed-triangular solves dozens to hundreds of time but only the
user of SimplicialLDLT knowns that, not SimplicialLDLT, nor
SparseMatrix. Moreover, the user does not own the SparseMatrix that we
want to analyze/optimize for. Other patterns are likely easier to handle, so let's focus on it for now.
Regarding (i1), I would suggest to introduce a new type, say
MklSparseMatrix<> that would enhance SparseMatrix<> through
MklSparseMatrix::beginAnalysis(Index expected_calls) const { // turn
*this to compressed mode // create handle // store expected_calls //
enable recording mode }
MklSparseMatrix::endAnalysis() const { // disable recording mode //
[optional] call mkl_sparse_optimize }
All states in MklSparseMatrix would be mutable.
Between a pair of beginAnalysis/endAnalysis each call to a supported
operation would trigger calls to mkl_sparse_set_*_hint()/mkl_sparse_optimize.
Optionally, we could even add a "dryrun" mode for which no operation
would be performed, only calls to mkl_sparse_set_*_hint() and then
mkl_sparse_optimize would be called in endAnalysis(). This way
mkl_sparse_optimize() would be called only once.
SimplicialLDLT<MklSparseMatrix<double> > llt(A); int it=0;
while(...) {
...
if(it==0) llt.matrixL().beginAnalysis(100); x = llt.solve(b);
if(it==0) llt.matrixL().endAnalysis(); ...
++it;
}
SimplicialLDLT<MklSparseMatrix<double> > llt(A);
llt.matrixL().beginAnalysis(100, DryRun); x = llt.solve(b); //
permutation and division by the diagonal matrix D would still be
performed, but calls to actual triangular solves would be by-passed
llt.matrixL().endAnalysis();
while(...) {
...
x = llt.solve(b);
...
}
If someone directly deal with the factor L, then we could follow the
SimplicialLLT<SparseMatrix<double> > llt(A);
MklSparseMatrix L(llt.matrixL());
L.beginAnalysis(100,DryRun);
y = L.triangularView<Lower>() * x;
L.endAnalysis();
while(...) {
...
y = L.triangularView<Lower>() * x;
...
}
This design in quite general and expendable to any sparse-optimizers,
even built-in ones in the future.
In contrast to the current proposal, only selected operations would be
passed to MKL (need to use a MklSparseMatrix + begin/end recording phase).
What do you think?
gael
On Tue, Apr 3, 2018 at 11:39 PM, Zhukova, Maria
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
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
<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).____
*#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____
__ __
__ __
__ __
So let's discuss about the public API, which currently is not very
(i1) - Storing MKL's handle in SparseMatrix breaks ABI and does not
sounds very generic.
(i2) - which operations are going to be analyzed/optimized,
(i3) - and specify the 'expected_calls' parameter.
(e.g., non-linear optimization, eigenvalues, ...)
SimplicialLDLT<SparseMatrix<double> > llt(A);
while(...) {
...
x = llt.solve(b);
...
}
Here the triangular L factor is going to be used for triangular and
transposed-triangular solves dozens to hundreds of time but only the
user of SimplicialLDLT knowns that, not SimplicialLDLT, nor
SparseMatrix. Moreover, the user does not own the SparseMatrix that we
want to analyze/optimize for. Other patterns are likely easier to handle, so let's focus on it for now.
Regarding (i1), I would suggest to introduce a new type, say
MklSparseMatrix<> that would enhance SparseMatrix<> through
MklSparseMatrix::beginAnalysis(Index expected_calls) const { // turn
*this to compressed mode // create handle // store expected_calls //
enable recording mode }
MklSparseMatrix::endAnalysis() const { // disable recording mode //
[optional] call mkl_sparse_optimize }
All states in MklSparseMatrix would be mutable.
Between a pair of beginAnalysis/endAnalysis each call to a supported
operation would trigger calls to mkl_sparse_set_*_hint()/mkl_sparse_optimize.
Optionally, we could even add a "dryrun" mode for which no operation
would be performed, only calls to mkl_sparse_set_*_hint() and then
mkl_sparse_optimize would be called in endAnalysis(). This way
mkl_sparse_optimize() would be called only once.
SimplicialLDLT<MklSparseMatrix<double> > llt(A); int it=0;
while(...) {
...
if(it==0) llt.matrixL().beginAnalysis(100); x = llt.solve(b);
if(it==0) llt.matrixL().endAnalysis(); ...
++it;
}
SimplicialLDLT<MklSparseMatrix<double> > llt(A);
llt.matrixL().beginAnalysis(100, DryRun); x = llt.solve(b); //
permutation and division by the diagonal matrix D would still be
performed, but calls to actual triangular solves would be by-passed
llt.matrixL().endAnalysis();
while(...) {
...
x = llt.solve(b);
...
}
If someone directly deal with the factor L, then we could follow the
SimplicialLLT<SparseMatrix<double> > llt(A);
MklSparseMatrix L(llt.matrixL());
L.beginAnalysis(100,DryRun);
y = L.triangularView<Lower>() * x;
L.endAnalysis();
while(...) {
...
y = L.triangularView<Lower>() * x;
...
}
This design in quite general and expendable to any sparse-optimizers,
even built-in ones in the future.
In contrast to the current proposal, only selected operations would be
passed to MKL (need to use a MklSparseMatrix + begin/end recording phase).
What do you think?
gael
On Tue, Apr 3, 2018 at 11:39 PM, Zhukova, Maria
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
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
<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).____
*#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____
__ __
__ __
__ __