This forum has been archived. All content is frozen. Please use KDE Discuss instead.

Infinite loop in general_matrix_matrix_product

Tags: None
(comma "," separated)
To_N
Registered Member
Posts
6
Karma
0
Hello,
I am using Eigen 3.1.0-alpha1

Using the code below I am able to cause an infinite loop in general_matrix_matrix_product.

JacN is 2x3 fixed size
JacD is 2x2 dynamic size
JacND(2, JacN.cols()+JacD.cols()-2) is 2x3 dynamic size
Code: Select all
JacND.block(0, 0, 2, JacN.cols()) = JacD.block(0, 0, 2, 2) * JacN;


Note:
- The infinite loop is triggered only in rare situation (e.g. in 1 out of 10 runs)
- The code itself is executed in a loop parallelized using OpenMP
- EIGEN_HAS_OPENMP is not defined
- Compiler used is Visual Studio 2010

The infinite loop is caused in GeneralMatrixMatrix.h row 188
for(Index i2=0; i2<rows; i2+=mc)
because mc = 0 (!), rows = 2

In line 81 Index mc is calculated as follows
Index mc = (std::min)(rows,blocking.mc());

level3_blocking<LhsScalar,RhsScalar>& blocking has members
{m_blockA=0x0000000000000000 m_blockB=0x0000000000000000 m_blockW=0x0000000000000000 m_mc=0 m_nc=3 m_kc=2}

The question is if I am doing something wrong here or if this is a bug? I also tried:
Code: Select all
JacND.block(0, 0, 2, JacN.cols()) = JacD.block(0, 0, 2, 2).eval() * JacN;

But it had no effect.

Regards

Call stack is if that is of use
geometry.dll!Eigen::internal::general_matrix_matrix_product<__int64,double,0,0,double,0,0,0>::run(__int64 rows=2, __int64 cols=3, __int64 depth=2, const double * _lhs=0x00000000031cb5e0, __int64 lhsStride=2, const double * _rhs=0x0000000003f2eba0, __int64 rhsStride=2, double * res=0x00000000031cb840, __int64 resStride=2, double alpha=1.0000000000000000, Eigen::internal::level3_blocking<double,double> & blocking={...}, Eigen::internal::GemmParallelInfo<__int64> * info=0x0000000000000000) Row 190 + 0x2b Bytes C++
geometry.dll!Eigen::internal::gemm_functor<double,__int64,Eigen::internal::general_matrix_matrix_product<__int64,double,0,0,double,0,0,0>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,Eigen::Matrix<double,-1,3,0,-1,3>,Eigen::internal::gemm_blocking_space<0,double,double,-1,3,2,0> >::operator()(__int64 row=0, __int64 rows=2, __int64 col=0, __int64 cols=3, Eigen::internal::GemmParallelInfo<__int64> * info=0x0000000000000000) Row 240 C++
geometry.dll!Eigen::internal::parallelize_gemm<1,Eigen::internal::gemm_functor<double,__int64,Eigen::internal::general_matrix_matrix_product<__int64,double,0,0,double,0,0,0>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,Eigen::Matrix<double,-1,3,0,-1,3>,Eigen::internal::gemm_blocking_space<0,double,double,-1,3,2,0> >,__int64>(const Eigen::internal::gemm_functor<double,__int64,Eigen::internal::general_matrix_matrix_product<__int64,double,0,0,double,0,0,0>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,Eigen::Matrix<double,-1,3,0,-1,3>,Eigen::internal::gemm_blocking_space<0,double,double,-1,3,2,0> > & func={...}, __int64 rows=2, __int64 cols=3, bool transpose=false) Row 106 + 0x30 Bytes C++
geometry.dll!Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>::scaleAndAddTo<Eigen::Matrix<double,-1,3,0,-1,3> >(Eigen::Matrix<double,-1,3,0,-1,3> & dst={...}, double alpha=1.0000000000000000) Row 435 + 0x80 Bytes C++
geometry.dll!Eigen::ProductBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3> >::scaleAndAddTo<Eigen::Matrix<double,-1,3,0,-1,3> >(Eigen::Matrix<double,-1,3,0,-1,3> & dst={...}, double alpha=1.0000000000000000) Row 124 + 0x46 Bytes C++
> geometry.dll!Eigen::ProductBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3> >::evalTo<Eigen::Matrix<double,-1,3,0,-1,3> >(Eigen::Matrix<double,-1,3,0,-1,3> & dst={...}) Row 115 + 0x44 Bytes C++
geometry.dll!Eigen::MatrixBase<Eigen::Matrix<double,-1,3,0,-1,3> >::lazyAssign<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3> >(const Eigen::ProductBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3> > & other={...}) Row 284 C++
geometry.dll!Eigen::PlainObjectBase<Eigen::Matrix<double,-1,3,0,-1,3> >::lazyAssign<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> >(const Eigen::DenseBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> > & other={...}) Row 386 C++
geometry.dll!Eigen::internal::assign_selector<Eigen::Matrix<double,-1,3,0,-1,3>,Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,0,0>::run(Eigen::Matrix<double,-1,3,0,-1,3> & dst={...}, const Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> & other={...}) Row 535 + 0x3a Bytes C++
geometry.dll!Eigen::PlainObjectBase<Eigen::Matrix<double,-1,3,0,-1,3> >::_set_noalias<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> >(const Eigen::DenseBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> > & other={...}) Row 592 C++
geometry.dll!Eigen::Matrix<double,-1,3,0,-1,3>::Matrix<double,-1,3,0,-1,3><Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> >(const Eigen::MatrixBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> > & other={...}) Row 294 + 0xf Bytes C++
geometry.dll!Eigen::DenseBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> >::eval() Row 389 + 0x17 Bytes C++
geometry.dll!Eigen::internal::assign_selector<Eigen::Block<Eigen::Matrix<double,-1,-1,0,-1,-1>,-1,-1,0,1>,Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5>,1,0>::run(Eigen::Block<Eigen::Matrix<double,-1,-1,0,-1,-1>,-1,-1,0,1> & dst={...}, const Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> & other={...}) Row 539 + 0x3b Bytes C++
geometry.dll!Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double,-1,-1,0,-1,-1>,-1,-1,0,1> >::operator=<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> >(const Eigen::DenseBase<Eigen::GeneralProduct<Eigen::Matrix<double,-1,-1,0,-1,-1>,Eigen::Matrix<double,2,3,0,2,3>,5> > & other={...}) Row 576 C++
geometry.dll!slSDK::geometry::Intrinsics::getOptimizerJacobian(Eigen::Matrix<double,-1,-1,0,-1,-1> & jacobian={...}, bool includeFocalLength=false, bool fixAspect=false, bool includeC=false, bool includeSkew=false, bool includeDist=false, bool includeTangentDist=false, bool includeDistRad3=false, const Eigen::Matrix<double,3,1,0,3,1> & xC={...}, bool derivePoint=true) Row 565 + 0xd0 Bytes C++
To_N
Registered Member
Posts
6
Karma
0
Note, I am able to reproduce this using the following code. Then in about 1 of 10 execution the code blocks.

Code: Select all
#pragma omp parallel for
for(int i = 0; i < 10000; i++)
{
   Eigen::Matrix<double, 2, 3> JacN;
   JacN.setOnes();
   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> JacD(2, 2);
   JacD.setOnes();
   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> JacND(2, JacN.cols()+JacD.cols()-2);
   JacND.block(0, 0, 2, JacN.cols()) = JacD.block(0, 0, 2, 2) * JacN;
}


The code works well if not run in parallel (i.e. without OpenMP).
So I guess something is not completely thread safe in this configuration of the matrix product (maybe because block() is used).
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
I remember I fixed a similar issue not long ago. Perhaps you re-try with alpha2? Also you might want to disable Eigen's own parallelization using either -DEIGEN_DONT_PARALLELIZE, or at runtime: Eigen::setNbThreads(1); If alpha-2 does not work, can you try to perform a dummy product before the openmp part:

MatrixXf a(1,1); a << 1;
a = a*a;
To_N
Registered Member
Posts
6
Karma
0
Thx for the answer,
I see in the changelog of alpha-2 that bug #406 was fixed, which is explained here:
http://stackoverflow.com/questions/8828 ... ck/8840839

I looks exactly like what I discovered. I will try alpha-2 and if that does not work post again.
Thx for the help!


Bookmarks



Who is online

Registered users: Bing [Bot], claydoh, Google [Bot], rblackwell, Yahoo [Bot]