Re: [eigen] Recursion and block matrices

[Norman Goldstein <normvcr@xxxxxxxxx>]

OK, I added


#   ifdef EIGEN_BLOCK_TRUE_PLUGIN
#     include EIGEN_BLOCK_TRUE_PLUGIN
#   endif

just before the protected part of  HasDirectAccess==true   in Block.h,
and have this in Block_true_plugin.h:

////////////////////// Block_true_plugin code ////////////////
inline
Block<XprType, BlockRows, BlockCols, InnerPanel, true>
block( Index startRow2,
       Index startCol2,
       Index blockRows2,
       Index blockCols2 )
{
  return
Block<XprType, BlockRows, BlockCols, InnerPanel, true>
    ( derived(),
      startRow2 + startRow(),   // <-- Accessor does not exist
      startCol2 + startCol(),      // <-- Accessor does not exist
      blockRows2,
      blockCols2 );
}
///////////////////////////////////////////////////////////////////

but, as I discovered, this specialization inherits from MapBase<Block>, and
unlike the generic Block definition: there is no startRow() and startCol()
accessors.  Any suggestions, please?



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

there is a specialization for Block (right after the generic
implementation) for HasDirectAccess==true, so you need to add your
specialization there too.

With the devel branch you might try:

void incr(Ref<MatrixXd> mat) {

// no cast needed
....

}

gael

On Tue, Aug 28, 2012 at 9:07 AM, Norman Goldstein <normvcr@xxxxxxxxx> wrote:

Thank you for the pointers.  It is reassuring to know that it was a bug that
stumped me.
I coded up an override Block::block() method to avoid the infinite
recursion:

////////////////////// Plugin Code //////////////////////////////////////
// EIGEN_BLOCK_PLUGIN

inline
Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess>
block( Index startRow2,
       Index startCol2,
       Index blockRows2,
       Index blockCols2 )
{
  return
Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess>
    ( derived(),
      startRow2 + startRow(),
      startCol2 + startCol(),
      blockRows2,
      blockCols2 );
}
/////////////////////////////////////////////////////////////////////////////

There is no Block plugin hook, so I put one in right before the protected
section of Block.h .
Also, I slightly modified the test driver:

//////////////////////////////// Code ////////////////////////////////////

#include <Eigen/Core>
using namespace Eigen;

#include <iostream>
using namespace std;

template< class Derived >
void incr( const Block< Derived >& mat )     // <-- Using Block< Derived >
{
  Block< Derived >& matV =
    const_cast< Block< Derived >& >( mat );

  ++matV(0,0);


  if( ( 1 == mat.rows() ) ||
      ( 1 == mat.cols() ) )
  {
    return;
  }
  else
  {
    incr( matV.block( 1,                                      // <-- Using
dynamic blocks, only
              1,

              mat.rows() - 1,
              mat.cols() - 1 ) );
  }
}// incr


int main( int argc, char* argv[] )
{
  MatrixXd mat( 2, 3 );
  mat << 1, 2, 3, 4, 5, 6;
  cout << "Before: " << mat << endl;

  incr( mat.block( 0, 0, mat.rows(), mat.cols() ) );   // <-- Start off with
a block


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

But there is still infinite recursion, even though the Block::block
method returns the same type of Block.  What am I missing?




On 08/27/2012 12:13 PM, Christoph Hertzberg wrote:

On 27.08.2012 19:33, Norman Goldstein wrote:

I haven't been able to figure out how to write a routine that calls
itself recursively
with successively smaller sub-blocks of the original matrix. Here is a
simple example,

With the current stable branch this is not possible, because this would
possibly generate infinitely many template instanciations. See also:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=408
There is an experimental solution here:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=481.

[...]

template< class Derived >
void incr( const DenseBase< Derived >& mat )
{
   DenseBase< Derived >& matV =
     const_cast< DenseBase< Derived >& >( mat );
[...]
   {
     incr( mat.bottomRightCorner( mat.rows() - 1,
                                                      mat.cols() - 1 ) );

Your problem here, however, is that mat is const, so
mat.bottomRightCorner(...) is a Block<const ...> and that inner const is not
removed by a const_cast. If you do matV.bottomRightCorner(...) instead, you
should get to the problem above ^^.

Btw:
If you just do tail recursion, you don't need to do recursion at all, just
do sth like this:

template<class Derived>
void incr(const DenseBase<Derived>& mat){

    for(int i=mat.rows(); i>0; --i) {
        ++mat.bottomRightCorner(i,i)(0,0);
    }
}

This version only works for square matrices, I guess you get the idea,
though.


Christoph