Re: [eigen] Recursion and block matrices

[ Thread Index | Date Index | More Archives ]

Attached are two short plugins that I wrote to implement a sub-matrix method, subm(), on a PlainObject and a MapBase. The interface is

// (startRow, startCol) is the upper left corner of the sub-matrix.
// blockRows is the number of rows in the sub-matrix
// blockCols is the number of columns in the sub-matrix
View subm( Index startRow,
                    Index startCol,
                    Index blockRows,
                    Index blockCols );

and View is the respective return type in the two plugins. The test routine, trec.cpp, shows that this works.
The correct output is

 10  20  30  40
 50  60  70  80
 90 100 110 120
 11  20  30  40
 50  61  70  80
 90 100 111 120

which shows that the identity has been added to the main diagonal of the "Before" matrix.

valgrind shows no errors and no leaks:
==32578== HEAP SUMMARY:
==32578==     in use at exit: 0 bytes in 0 blocks
==32578==   total heap usage: 56 allocs, 56 frees, 13,022 bytes allocated
==32578== All heap blocks were freed -- no leaks are possible
==32578== For counts of detected and suppressed errors, rerun with: -v
==32578== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

However, there is a software engineering issue that I have not been able to resolve:

There are two places in trec.cpp, right before calling the "incr" function, where I need to explicitly reference a Map object as a MapBase object,
      or else the code will not compile / run correctly.

For example, when changing the "incr(sub)" call in main to incr(sub0), the compilation error is:
error: no matching function for call to

incr( Eigen::Map< Eigen::Matrix<double, -0x000000001, -0x000000001>,
Eigen::Stride<-0x000000001,-0x000000001> >& )

note: candidate is:

template< class Derived > void incr( Eigen::MapBase< Derived >& )

Or, when changing the reference to be a copy constructor, the code compiles,
but produces a core dump:


Thank you,

PS The subm notation is from DLIB.

On 09/18/2012 10:03 PM, Norman Goldstein wrote:
Thanks for explaining Stride<0,0>.

A few weeks, ago, when I read Gael's email, I was already sufficiently overwhelmed with Eigen, and did not want to get involved in the Devel branch. Today I took a look, and the Ref certainly seems related, but I can't say that I understand it well enough to use it. It seems that the inner stride must be 1?

What I am aiming for is conceptually simple, and it is allowing me to learn the Eigen class hierarchy well enough to implement sub-matrices that will be less efficient as expressions, but will not have infinite recursion, and will share the same memory as the original matrix. In a nutshell, a sub-matrix is a Map of a rectangular subset of a PlainObject.

Will let you guys know when I am done. I'm sure I will not have done this in the best way, and any suggestions you have for improvement will be taken to heart. I hope will not have too many question to bother you along the way. If you have any other ideas for defining submatrices with the qualities I described above, please shout. Although I enjoy the challenge, I would not classify this as fun, and am looking forward to using this functionality.

Thanks, again,

On 09/18/2012 08:59 AM, Christoph Hertzberg wrote:
On 18.09.2012 17:51, Norman Goldstein wrote:
I figured this one out by decomposing the C++ statements until the
compiler message was more specific.
It turns out that the compiler was being asked to match Stride<0,0> with
Stride<-1,-1>.  The former is
defaulted in some declarations, while the latter is for dynamic stride.
When would one want to use

Stride<0,0> means "natural" stride, so colstride is number of rows, rowstride is 1 (assuming colmajor matrices). Have you yet tried to use the Ref<MatrixXd> class, as Gael suggested some mails ago? Or are you implementing this more for the fun of it?


On 08/28/2012 01:01 AM, Gael Guennebaud wrote:

With the devel branch you might try:

void incr(Ref<MatrixXd> mat) {

// no cast needed



typedef Derived View;

View subm( Index startRow2, 
	   Index startCol2, 
	   Index blockRows2, 
	   Index blockCols2 )
    Derived( &this->operator()(startRow2, startCol2),
	     blockRows2, blockCols2, 
  	     Stride< Dynamic, Dynamic >( this->outerStride(), 
  					 this->innerStride() ) );

}// subm

typedef typename StridedMapType< Stride<Dynamic,Dynamic> >::type View;

View subm( Index startRow2, 
	   Index startCol2, 
	   Index blockRows2, 
	   Index blockCols2 )
    Map( &this->operator()(startRow2, startCol2),
	 blockRows2, blockCols2, 
	 Stride< Dynamic, Dynamic >( this->outerStride(), 1 ) );
}// subm
#include <Eigen/Core>
using namespace Eigen;

#include <iostream>
using namespace std;

template< class Derived >
void incr( MapBase< Derived >& matV )

  if( ( 1 == matV.rows() ) ||
      ( 1 == matV.cols() ) )

  auto sub0 = matV.subm( 1,
			 matV.rows() - 1,
			 matV.cols() - 1 );
  MapBase< Derived >& sub = sub0;
  incr( sub );
}// incr

int main( int argc, char* argv[] )
  MatrixXd mat( 3, 4 );
  mat << 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120;
  cout << "Before:\n" << mat << endl;

  Map< MatrixXd, Unaligned, Stride<-1,-1>> sub0 
    = mat.subm( 0, 0, mat.rows(), mat.cols() );

  // Compiles, but causes a core dump
  // MapBase< Map< MatrixXd, Unaligned, Stride<-1,-1>> > sub = sub0;
  MapBase< Map< MatrixXd, Unaligned, Stride<-1,-1>> >& sub = sub0;

  // Does not compile
  // incr( sub0 );  
  incr( sub );

  cout << "After:\n" << mat << endl;
  return 0;
}// main

Mail converted by MHonArc 2.6.19+