[eigen] Extract sparse sub-matrices from a sparse matrix

[Alberto Minetti <bangirasu@xxxxxxxxx>]

I need to extract a block from a Eigen::SparseMatrix<double> but 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’

So I made this function to extract blocks from the
Eigen::SparseMatrix<double> and other four functions if the sub-matrix
is in one of the four corners. The code below has a main with tests.

#include <iostream>
#include <vector>
#include <stdio.h>
#include <assert.h>
#include <Eigen/Sparse>

using namespace std;
using namespace Eigen;

#define N 30

typedef Triplet<double> Tri;

SparseMatrix<double> sparseBlock(SparseMatrix<double> M,
int ibegin, int jbegin, int icount, int jcount){
	//tested using 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;
}

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);
}

int main(int argc, char *argv[]){
	
	std::vector<Tri> tripletList;
	tripletList.reserve(N*N);

	for(int i=1; i<N; i=i+2)
		for(int j=1; j<N; j=j+2)
			tripletList.push_back(Tri(i,j,i*N+j));
	
	SparseMatrix<double> M(N,N);
	M.setFromTriplets(tripletList.begin(), tripletList.end());
	// M is ready to go!
	
	cout << "Here is the matrix M("
	<<M.rows()<<"x"<<M.cols()<<"):\n" << M << endl;

	SparseMatrix<double> B=sparseBlock(M,5,3,4,9);
	cout << "Here is the matrix B("
	<<B.rows()<<"x"<<B.cols()<<"):\n" << B << endl;

	B=sparseBottomLeftBlock(M,4,2);
	cout << "Here is the Bottom Left matrix B("
	<<B.rows()<<"x"<<B.cols()<<"):\n" << B << endl;
	
	B=sparseBottomRightBlock(M,3,2);
	cout << "Here is the Bottom Right matrix B("
	<<B.rows()<<"x"<<B.cols()<<"):\n" << B << endl;
	
	B=sparseTopLeftBlock(M,6,2);
	cout << "Here is the Top Left matrix B("
	<<B.rows()<<"x"<<B.cols()<<"):\n" << B << endl;
	
	B=sparseTopRightBlock(M,2,4);
	cout << "Here is the Top Right matrix B("
	<<B.rows()<<"x"<<B.cols()<<"):\n" << B << endl;
	
	return 0;
}


--
Alberto Minetti