New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 408 - block() does not allow recursion
Summary: block() does not allow recursion
Status: REVIEWNEEDED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Unknown
Assignee: Gael Guennebaud
URL:
Whiteboard:
Keywords:
Depends on: 58 ExprEval
Blocks: 3.x
  Show dependency treegraph
 
Reported: 2012-01-16 00:19 UTC by Jonathan Schmidt-Dominé
Modified: 2018-10-09 20:09 UTC (History)
4 users (show)



Attachments
Patch to fuse Block<Block<X>> as Block<X> (14.62 KB, patch)
2016-02-23 10:16 UTC, Gael Guennebaud
no flags Details | Diff
A second approach to merge Block<Block<A>> as Block<A> (34.54 KB, patch)
2016-02-29 09:57 UTC, Gael Guennebaud
no flags Details | Diff
New version of the patch merging both previous ones (41.29 KB, patch)
2016-03-07 14:55 UTC, Gael Guennebaud
no flags Details | Diff

Description Jonathan Schmidt-Dominé 2012-01-16 00:19:52 UTC
Hi,

when calling block() on a matrix as returned by block() I get an object of a different type. Thus I cannot call block() recursively (when using templates I get infinite many instantiations). I think it would be better if block-objects would have their own block()-implementation generating an object of the same type.

Regards
Comment 1 Christoph Hertzberg 2012-01-16 00:26:59 UTC
That's basically bug 58.
As a work-around you can manually generate Map objects with the corresponding stride, for now.
Comment 2 Gael Guennebaud 2012-01-18 19:56:40 UTC
"evaluators" should solve this problem too.
Comment 3 Christoph Hertzberg 2014-11-02 17:31:51 UTC
I'm not sure how exactly evaluators will solve this, unless they already do some kind of tree-optimization when determining the type.

E.g., Block<Block<Matrix> > is a Block<Matrix> (with adapted indexes) or things like Block<Product<A, B> > are simplified to Product<Block<A>, Block<B> > before passing them further (the same for all kinds of Cwise operations).
Comment 4 Gael Guennebaud 2014-11-04 09:26:34 UTC
yes I was probably thinking about tree-optimizers, but indeed, for Block<Block<>> it is better to change the Block return type at the cost of many boilerplate lines of code because there are numerous methods returning a Block.
Comment 5 Gael Guennebaud 2016-02-23 10:16:20 UTC
Created attachment 655 [details]
Patch to fuse Block<Block<X>> as Block<X>

Here is a patch to fuse Block<Block<X>> as Block<X>. This implementation simply add specializations of Block<Block<X>> which inherits Block<X>. The actual type is thus still a Block<Block<X>> but the Derived class onto which any method apply is a Block<X>. This is similar to VectorBlock which is only used as a helper to construct sub-vectors.

Please share your thoughts before I push. I might have overseen some aspects.
Comment 6 Jonathan Schmidt-Dominé 2016-02-23 13:05:33 UTC
This works well as long as people use block() and do not use Block<> manually (but auto). Thus most cases are covered. However, for the second case it would be cleaner to use a template typedef in case that C++11 is used. However, there would be some boiler plate code like that:
template<typename A>
struct BlockImpl;

template<typename A>
struct BlockHelper
{
	typedef BlockImpl<A> Block;
};

template<typename A>
using Block = typename BlockHelper<A>::Block;

template<typename A>
struct BlockHelper<BlockImpl<A>>
{
	typedef typename BlockHelper<A>::Block Block;
};

template<typename A>
struct BlockImpl
{
// implementation goes here
};
Comment 7 Christoph Hertzberg 2016-02-25 10:29:57 UTC
The patch looks good in solving the immediate issue, i.e., it allows things like

  template<class MatrixType>
  void recursion(const MatrixType& m)
  {
    if(m.rows()>1)
      recursion(m.topRows(m.rows()/2));
    else { /* do other stuff */ }
  }

So, it allows doing basic divide-and-conquer implementations on matrices.


Anything more sophisticated, like the following will still fail:

      recursion(m.topRows(m.rows()-1).transpose());

On my system, eventually it fails with
    virtual memory exhausted: Cannot allocate memory
To reproduce this in a finite amount of time, I recommend limiting the virtual memory (e.g., using ulimit -v 500000 on linux)

On the other hand tranpose alone has the same issue:
  recursion(m.transpose());

i.e., m.transpose().transpose() is not simplified to m
Comment 8 Gael Guennebaud 2016-02-27 10:09:32 UTC
I'm not sure to see how the template typedef work work unless adding ctor Block(Block<A> &) to class Block<A>... that looks rather dangerous. We already provide all the proper typedefs so that users do not have to directly deal with Block<>, e.g.:

A::ConstBlockXpr


Regarding the current implementation, in the following example:

template<typename MatrixType>
void block_rec(const MatrixType& a)
{
  if(a.cols()>1)
  {
    block_rec(a.leftCols(a.cols()/2));
    block_rec(a.rightCols(a.cols()/2));
  }
}

block_rec(MatrixXd(..));

three versions of block_rec will be instantiated: one for MatrixXd, one for Block<MatrixXd> and one for Block<Block<MatrixXd> >.

The last one could be avoided by the user by calling .derived():

    block_rec(a.leftCols(a.cols()/2).derived());


On Eigen's side we could make sure that the returned types of block/leftCols/etc. are the a Block<MatrixXd> and not a Block<Block<MatrixXd> >. I see two options for that:

1 - make them return typename "current-return-type"::Derived, but then the following won't work:

    Block<A>::BlockXpr blA = a..block(...).block(...);

    because the type returned by the second block won't be exactly Block<A>::BlockXpr.


2 - declare the public return types as:
    typedef typename "current-return-type"::Derived new-return-type;
    but then the following won't work:

    Block<A>::BlockXpr blA(a.block(),...);

    because Block<A>::BlockXpr won't have the proper constructor. This version is similar to the template typedef proposed by Jonathan.
Comment 9 Gael Guennebaud 2016-02-29 09:57:55 UTC
Created attachment 660 [details]
A second approach to merge Block<Block<A>> as Block<A>

Here is another working attempt. The main advantage of this version is that:

A a;
a.block(..).block(..);

truly returns a Block<A>. Block<A>::BlockXpr is also a Block<A>.

On the other hand, compared to the previous patch, if the user explicitly build a Block<Block<A> >, then no merging will occur. This behavior could be achieved by combining this new patch with the first one, or, in c++11, using the approach suggested by Jonathan. However, the user might have a good reason to nest blocks, and it might be confusing that Block<Block<A> >::NestedExpression is not a Block<A>. So after all, it's probably better for API compatibility to keep such explicit constructions as is.

So the main drawback of this new implementation are the novel constructors of Block for which I envision a risk of unpredicted and wrong behavior.

Another risk, is if someone write code like:

template <typename T>
void foo1(const T &x) {
  typename T::ConstBlockXpr blx = x.block(...); // won't compile anymore
}

template <typename D>
void foo2(const MatrixBase<D> &x) {
  foo1(x);
}

foo2(a.block(...));

The user-side fix is to call foo1(x.derived()), which is what the user should do anyway.
Comment 10 Christoph Hertzberg 2016-02-29 12:34:05 UTC
I guess the best solution for the user who wants to avoid multiple instantiations is to declare the method as:

  void recursive(const Ref<const MatrixXd>& a);

Unfortunately, that would require declaring methods (or at least wrapper functions) for every scalar type. The following won't work, since templated functions do not auto-convert their arguments:
  template<typename Scalar>
  void recursive(const Ref<const Matrix<Scalar, Dynamic, Dynamic> >& a)

I don't see a clean solution which avoids any duplicate instantiation (for Matrix and Block<Matrix> of the same scalar type) but still allows calling with arbitrary types (and no unnecessary copying).
Comment 11 Gael Guennebaud 2016-03-02 12:16:06 UTC
Unifying Matrix<> and Block<Matrix<>> in a templated way is another topic.

Regarding this Block<Block<> > topic, I'm still not sure about which implementation to chose if any. Indeed, the second one is more appealing because it truly remove Block<Block<>> types, but it also breaks existing code like:

  A a;
  typedef Block<A> B;
  B b = a.block(...);
  /*... later ...*/
  Block<B> c = b.block(...); // won't compile anymore

Of course, the user should rather write:

  Block<B> c(b,...);

or better:

  auto c = b.block();
  B::BlockXpr c = b.block(...);

but I've seen the first version is several user codes.

On our side, it seems that this issue could be solved by combining both approaches, that is, by adding to the second version specializations of Block<Block<X> > inheriting Block<X>::BlockXpr with additional implicit conversion constructors taking a Block<X>::BlockXpr as argument and passing it to its base class.

Nonetheless, all these "tricks" looks rather fragile...

Perhaps, it is safer to defer its resolution to 3.4?
Comment 12 Christoph Hertzberg 2016-03-02 12:55:53 UTC
I'd be ok with breaking code that explicitly writes
  Block<B> c = b.block(...);
rather than:
  B::BlockXpr c = b.block();

If B already was a Block<Block<...> >, wouldn't it work with your first implementation, either?

I think we should not encourage to explicitly use Block<B> anywhere -- if we wanted to keep this legal, we would also lose the opportunity for optimizations like Matrix::BlockXpr==Map<...>, etc. 
I'm not sure if it really makes a difference if we postpone that decision to 3.4 (sure, we could explicitly discourage using Block<Xpr> for 3.3 and really change it in 3.4).
Comment 13 Gael Guennebaud 2016-03-02 21:12:26 UTC
Actually, there is no B::BlockXpr in 3.2, not even in 3.3-beta1. I've added this missing typedef only recently.
Comment 14 Christoph Hertzberg 2016-03-03 17:27:07 UTC
(In reply to Gael Guennebaud from comment #13)
> Actually, there is no B::BlockXpr in 3.2, not even in 3.3-beta1. I've added
> this missing typedef only recently.

Ok, I actually did not check that (just looked at the current dox-devel page).

So essentially, what you suggest in comment 11 is (omitting some parts for brevity):

  Block<A,r1,c1>::block<r2,c2>(...) returns a Block<A,r2,c2>

  Block<Block<A,r1,c1>,r2,c2> inherits from Block<A,r2,c2> and has an implicit constructor from type Block<A,r2,c2> and explicit constructors taking (Block<A,r1,c1>, startRow, startCol, rows, cols) [and some other variants].

All constructors of Block<Block<A,...>,...> could (later) be marked deprecated.

On first glance this looks sound, though I'm not sure either if there are any loopholes introduces by this ...
Comment 15 Christoph Hertzberg 2016-03-03 17:31:49 UTC
Actually, one thing that might break: If you have
  MatrixXd A;
  auto X = A.block<4,4>(...).block<Dynamic,Dynamic>(...);
Is  traits<decltype(X)>::MaxRowsAtCompileTime==4 ?

Admittedly, a rather pathological example and it's not unlikely that I'm not fully understanding your implementation.

So perhaps too complicated for 3.3 indeed -- if someone wants recursion, Ref<Matrix> can be used as an (even better) workaround. 
But the BlockXpr typedefs would be a good step in making future code more portable/easier to refactor.
Comment 16 Gael Guennebaud 2016-03-07 14:41:21 UTC
Yes for both. Comment 14 well summaries the ideas, and in comment 15, you are right that the max-sizes information will be lost.

The only way to fix this issue would be to extend the Block class to take max-sizes as template parameters.

So as it is still unclear how to properly fix that entry, I again suggest to postpone it for 3.4. Releasing 3.3 asap is key priority for me.
Comment 17 Gael Guennebaud 2016-03-07 14:55:45 UTC
Created attachment 664 [details]
New version of the patch merging both previous ones

For the record, here is a new version of the patch combining both approaches with extended regression unit tests. The ones checking for the preservation of max-sizes are failing.

Note You need to log in before you can comment on or make changes to this bug.