Re: [eigen] Lazy evaluation and adjoint(): Bug?

[ Thread Index | Date Index | More Archives ]


is the result correct? If yes, then this means Eigen is following a custom conjugation path. If not, then could you try the following patch (this code path seems to assume that IsComplex <=> std::complex):

diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
--- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h
+++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
@@ -429,25 +429,26 @@ public:
   EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const
     dest = ploadu<LhsPacketType>(a);

   template<typename LhsPacketType, typename RhsPacketType, typename AccPacketType>
   EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, AccPacketType& tmp) const
+    conj_helper<LhsPacketType,RhsPacketType,ConjLhs,ConjRhs> cj;
     // It would be a lot cleaner to call pmadd all the time. Unfortunately if we
     // let gcc allocate the register in which to store the result of the pmul
     // (in the case where there is no FMA) gcc fails to figure out how to avoid
     // spilling register.
-    c = pmadd(a,b,c);
+    c = cj.pmadd(a,b,c);
-    tmp = b; tmp = pmul(a,tmp); c = padd(c,tmp);
+    tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);

   EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
     r = pmadd(c,alpha,r);


On Wed, Sep 14, 2016 at 2:46 PM, Peter <list@xxxxxxxxxxxxxxxxx> wrote:
Dear All,

while  desperately investigating my current "Learning Eigen" project, I came across the following problem.

I'm using Matrices for self defined complex types, where conj() is defined as

template<class Q>
typename boost::enable_if< boost::is_complex<Q>, MyType<Q> >::type conj( const MyType< Q >& a)
    static int counter(0);
    std::cout << "#MyType conjugate called " << ++counter << " : " << a << std::endl;
    return MyType<Q>( conj( a.z1() ), conj( a.z2() ));

Now, somewhere deep in my code I have the following rotation:

std::cout << "#------Lazy----------\n";
TpMatrix TT0  =  U.adjoint() * T * U;

std::cout << "#------helper-----------\n";
TpMatrix Uadj = U.adjoint()  ;

std::cout << "#-----------------\n";
TpMatrix TT   = Uadj * T * U;
std::cout << "#-------Diff------\n\n" << TT -TT0 << std::endl;;

where the involved matrices have dimension 10x10.

Now, in the evaluation of TT0 MyType::conj gets called precisely one time,
while in the evaluation of Uadj it gets called 100 == 10*10 times, which I would naively expect.

Am I'm missing something?

I also just tested that
TpMatrix TT1  =  (U.adjoint() * T).eval()  * U;

calls MyType::conj only once.

Is this intended behaviour? Am I'm missing something?
Since TT is different from TT0 I assume that Eigen is too lazy here.
Sadly, I still don't know enough of Eigens internals to make sense of the code paths shown by ddd/gdb.

Any help is highly appreciated,
best regards,

Mail converted by MHonArc 2.6.19+