CODEWITHSUNDEEP
c++matrixsparse-matrixeigen
Alberto
Alberto

Reputation: 2982

Extract a block from a sparse matrix as another sparse matric

How to extract a block from a Eigen::SparseMatrix<double>. It seems there aren't the methods I used for the dense ones.

‘class Eigen::SparseMatrix<double>’ has no member named ‘topLeftCorner’
‘class Eigen::SparseMatrix<double>’ has no member named ‘block’

There is a way to extract a block as a Eigen::SparseMatrix<double> ?

Upvotes: 5

Views: 3395

Answers (3)

Akavall
Akavall

Reputation: 86168

This is now supported in Eigen 3.2.2 Docs (though maybe earlier versions support it too).

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>

using namespace Eigen;

int main()
{
  MatrixXd silly(6, 3);

  silly << 0, 1, 2,
           0, 3, 0,
           2, 0, 0,
           3, 2, 1,
           0, 1, 0,
           2, 0, 0;



  SparseMatrix<double, RowMajor> sparse_silly = silly.sparseView();

  std::cout <<"Whole Matrix" << std::endl;
  std::cout << sparse_silly << std::endl;

  std::cout << "block of matrix" << std::endl;
  std::cout << sparse_silly.block(1,1,3,2) << std::endl;

  return 0;
}

Upvotes: 4

Jakob
Jakob

Reputation: 2380

There is very sparse support (sorry, no pun intended) for submatrices in sparse matrices. Effectively you can only access a continuous set of rows for row-major, and columns for column major. The reason for that is not that the matrices could be empty, but rather that the indexing scheme is somewhat more complicated than with dense matrices. With dense matrices you only need an additional stride number in order to support sub-matrix support.

Upvotes: 1

Alberto
Alberto

Reputation: 2982

I made this function to extract blocks from a Eigen::SparseMatrix<double,ColMaior>

typedef Triplet<double> Tri;
SparseMatrix<double> sparseBlock(SparseMatrix<double,ColMajor> M,
        int ibegin, int jbegin, int icount, int jcount){
        //only for ColMajor Sparse Matrix
    assert(ibegin+icount <= M.rows());
    assert(jbegin+jcount <= M.cols());
    int Mj,Mi,i,j,currOuterIndex,nextOuterIndex;
    vector<Tri> tripletList;
    tripletList.reserve(M.nonZeros());

    for(j=0; j<jcount; j++){
        Mj=j+jbegin;
        currOuterIndex = M.outerIndexPtr()[Mj];
        nextOuterIndex = M.outerIndexPtr()[Mj+1];

        for(int a = currOuterIndex; a<nextOuterIndex; a++){
            Mi=M.innerIndexPtr()[a];

            if(Mi < ibegin) continue;
            if(Mi >= ibegin + icount) break;

            i=Mi-ibegin;    
            tripletList.push_back(Tri(i,j,M.valuePtr()[a]));
        }
    }
    SparseMatrix<double> matS(icount,jcount);
    matS.setFromTriplets(tripletList.begin(), tripletList.end());
    return matS;
}

And these if the sub-matrix is in one of the four corners:

SparseMatrix<double> sparseTopLeftBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,0,0,icount,jcount);
}
SparseMatrix<double> sparseTopRightBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,0,M.cols()-jcount,icount,jcount);
}
SparseMatrix<double> sparseBottomLeftBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,M.rows()-icount,0,icount,jcount);
}
SparseMatrix<double> sparseBottomRightBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,M.rows()-icount,M.cols()-jcount,icount,jcount);
}

Upvotes: 5

Related Questions