Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness*From*: Manoj Rajagopalan <rmanoj@xxxxxxxxx>*Date*: Wed, 14 Jul 2010 19:41:58 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

Attached compilable C++ program. thanks, Manoj On Wednesday 14 July 2010 07:14:02 pm Benoit Jacob wrote: > Hi, > > It would help if you could make a compilable test program, so we can > quickly run stuff. > > Benoit > > 2010/7/14 Manoj Rajagopalan <rmanoj@xxxxxxxxx>: > > Hi all, > > > > Could someone please help me figure out the source of the following > > discrepancy? > > > > Given: > > int const N = 5; > > MatrixXi A(N,N); > > A.setRandom(); > > > > TriangularView<MatrixXi,Lower> const A_l(A); > > TriangularView<MatrixXi,Lower|ZeroDiag> const A_zl(A); > > > > the following must be true: > > > > A.selfadjointView<Lower>() == A_l + A_zl.transpose() > > > > Basically, I am adding the lower triangle of A to its own transpose, > > without adding the diagonal a second time. So the following is calculated > > correctly by Eigen: > > > > MatrixXi C1(N,N), C2(N,N); > > C1 = A.selfadjointView<Lower>(); > > C2 = A_l.toDenseMatrix() + A_zl.transpose().toDenseMatrix(); > > assert(C1 == C2); > > > > But, the following code produces a discrepancy: > > > > MatrixXi B(N,N); > > B.setRandom(); > > C1 = A.selfadjointView<Lower>() * B; > > C2 = (A_l * B) + (A_zl.transpose() * B); > > assert(C1 == C2); // <-------- FALSE! > > > > After a long day my brain is unable to figure out if I am doing > > something extremely silly - it seems too trivial to be a bug. Can someone > > offer comments? > > > > thanks, > > Manoj

#include <Eigen/Core> #include <iostream> using namespace Eigen; using namespace std; int main(void) { int const N = 5; MatrixXi A(N,N), B(N,N), C1(N,N), C2(N,N), C_diff(N,N); A.setRandom(); B.setRandom(); C1 = A.selfadjointView<Lower>(); // C1 = A.selfadjointView<Lower>() * B; TriangularView<MatrixXi,Lower> const A_l(A); TriangularView<MatrixXi,Lower|ZeroDiag> const A_zl(A); // C2 = A_l * B; // C2 += A_zl.transpose() * B; C2 = A_l.toDenseMatrix() + A_zl.transpose().toDenseMatrix(); C_diff = C1 - C2; cout << "A =\n" << A << endl; cout << "\nA.triangularView<Lower>() =\n" << A_l.toDenseMatrix() << endl; cout << "\nA.triangularView<Lower|ZeroDiag>() =\n" << A_zl.toDenseMatrix() << endl; cout << "C1 =\n" << C1 << endl; cout << "C2 =\n" << C2 << endl; cout << "C_diff =\n" << C_diff << endl; assert(C1 == C2); return 0; }

**References**:**[eigen] discrepancy in (triangular view + transpose) and self-adjointness***From:*Manoj Rajagopalan

**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness** - Next by Date:
**Re: [eigen] Implementing a custom-matrix product** - Previous by thread:
**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness** - Next by thread:
**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |