Hi !

Working with SparseMatrix<double>.
Trying to compute A - B where B=M - diag(M)

I'am not even able to compile the very first step of this computation :

Code: Select all

SparseMatrix<double> diagM = M.diagonal( ).asDiagonal( ) ;


Receiving the compiler error :
/home/rousseau/MPI/proj/src/elemFinSeq.cpp|27|erreur: conversion from ‘const Eigen::DiagonalWrapper<const Eigen::Diagonal<const Eigen::SparseMatrix<double>, 0> >’ to non-scalar type ‘Eigen::SparseMatrix<double>’ requested|

May I ask you some hints please ?

All the best

To make it efficient one option would be:

Code: Select all

C = A - M;
C.diagonal() += M.diagonal();


In case you need B=M-diag(M) multiple times:

Code: Select all

B=M;
B.prune([] (int i, int j, double&) { return i!=j; });
C = A-B;


Many thanks Gaël for your precious support.

May I ask you what's wrong with my naive way to extract the sparse matrix diagonal in my snippet , I don't actually understand the meaning of the compiler complaint ?

Tried to use your hint to fetch a sparse matrix without its own diagonal :

Code: Select all

  SparseMatrix<double> M(3,3) ,C(3,3);
    C = M;
   C.diagonal() -= M.diagonal() ;


Does not compile neither.

All the best

Sylvain

Your version does not work because there is no conversion from a diagonal matrix to a sparse one. It should not be difficult to add, but this approach will still be slower than the ones I suggested because computing M-diagM is a O(nnz) operation with a somewhat complicated logic and dynamic memory management (as A-B). On the other hand, C.diagonal() -= M.diagonal() exploits the fact that the coefficients are sorted to quickly find the entries which as to be modified.

Regarding, "C.diagonal() -= M.diagonal()", indeed, I've forgotten that SparseMatrix::diagonal() was read-only. In the devel branch it is now read-write with the restriction that the diagonal entries must exist if the user writes to the diagonal.