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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

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.

 #ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD

     EIGEN_UNUSED_VARIABLE(tmp);

-    c = pmadd(a,b,c);

+    c = cj.pmadd(a,b,c);

 #else

-    tmp = b; tmp = pmul(a,tmp); c = padd(c,tmp);

+    tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);

 #endif

   }

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

   {

     r = pmadd(c,alpha,r);

   }

gael.

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,
Peter