| Re: [eigen] Recursion and block matrices |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen 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
Before:
10 20 30 40
50 60 70 80
90 100 110 120
After:
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==
==32578== All heap blocks were freed -- no leaks are possible
==32578==
==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>,
0,
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,
Norm
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,
Norm
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>?
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?
Christoph
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
...
}
// EIGEN_MAPBASE_WRITE_PLUGIN
typedef Derived View;
View subm( Index startRow2,
Index startCol2,
Index blockRows2,
Index blockCols2 )
{
return
Derived( &this->operator()(startRow2, startCol2),
blockRows2, blockCols2,
Stride< Dynamic, Dynamic >( this->outerStride(),
this->innerStride() ) );
}// subm
// EIGEN_PLAINOBJECTBASE_PLUGIN
typedef typename StridedMapType< Stride<Dynamic,Dynamic> >::type View;
View subm( Index startRow2,
Index startCol2,
Index blockRows2,
Index blockCols2 )
{
return
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 )
{
++matV(0,0);
if( ( 1 == matV.rows() ) ||
( 1 == matV.cols() ) )
{
return;
}
auto sub0 = matV.subm( 1,
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