Passing argument as `Ref<const _MatrixType>` instead of `const _MatrixType&`

[dariomangoni]

I'm a first-time user of Eigen and Spectra, so what I'm experiencing may not be critical, even not correct, anyway...

I need to solve a Generalized Eigenvalue Problem on CSR-format matrices. I already have CSR arrays so I would like to wrap them in a Eigen::SparseMatrix in order to feed the Spectra::SymGEigsSolver by mean of Spectra::SparseSymMatProd and Spectra::SparseCholesky.

To wrap CSR C-style arrays I would use Map<SparseMatrix<double>>(...CSRmatrixdata...), but then I came to an empasse. How can I call the SparseSymMatProd constructor? There are few ways...

1. pass the Map<> itself
2. pass a Ref<> to the Map<>
3. pass an Eigen::SparseMatrix& referring to Map<>
4. pass an Eigen::SparseMatrix& referring to Ref<> that refers to the Map<>

What happens? Fact 1) inside the SparseSymMatProd constructor, a const Eigen::SparseMatrix& is initialized with the argument of the constructor. i.e. const Eigen::SparseMatrix& m_mat = _argument_

Fact 2a) these passages do not call copy* Map<>(CSRarrays) Ref<> = Map<> and obviously Eigen::SparseMatrix& = Eigen::SparseMatrix&

Fact 2b) these passages do call copy* Eigen::SparseMatrix& = Ref<> Eigen::SparseMatrix& = Map<>

*in order to prove that, I compared the valuePtr() address with the original CSR value address

Now: in [1] and [2]: the reference-initialization inside the constructor forces a copy, but suddenly the destructors of SparseMatrix first and CompressedStorage after are called and they destroy the arrays of m_mat that have just been created! The matrix is invalidated an no product can be actually performed.

in [3] and [4]: everything works fine, the reference-initialization called inside the constructor is of the last case of Fact 2a, so there is no copy and avoid this strange behaviour above, but there do is a copy in the reference initialization that is made before calling the constructor.

So: I would recommend to change the first lines of Spectra::SparseSymMatProd from

private:
  const SparseMatrix& m_mat;
public:
    SparseSymMatProd(const SparseMatrix& mat_) :
        m_mat(mat_) { }

to

private:
  Eigen::Ref<const SparseMatrix> m_mat;
public:
    SparseSymMatProd(Eigen::Ref<const SparseMatrix> mat_) :
        m_mat(mat_) { }

This is retro-compatible: it works with Eigen::SparseMatrix and it should works also with Refs and it doesn't force a copy.

I'm not aware of all the implications and if this is actually correct. If you, instead, have a better way to pass some arrays to the eigen solver without this mess and that doesn't force copies, please tell me!

[yixuan]

Hi @dariomangoni , thanks for the comments! This is indeed a good point. Map<SparseMatrix<double>> and Ref<SparseMatrix<double>> are new features introduced in the recent Eigen 3.3, and they didn't exist before. I need to think of a way to make it compatible with older versions of Eigen, but at least for your own use, it's perfectly a good idea to make that change.

[yixuan]

This has been solved in https://github.com/yixuan/spectra/commit/dec0d1fc1ad40f860c3d38e4835fa94135e7f65a. Thank you for the suggestion!