conjugate() and adjoint() bug example (Re: [eigen] Lazy evaluation and adjoint())

[Peter <list@xxxxxxxxxxxxxxxxx>]

Dear All,

I've composed a small example that shows highly non-expected behavior for conjugate and adjoint() for user defined Type,
see attached source file, which demonstrates that  conj() is not called in these cases for user defined types,
while it is called for std::complex.

//  typedef std::complex<double> MyUsedType;   ///< Use this for std complex type
typedef MyType MyUsedType;                   ///< Use this for user defined complex type


//typedef Eigen::Matrix< MyUsedType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TpMatrix;
typedef Eigen::Matrix< MyUsedType, 2,2 > TpMatrix;


int main()
{
	MyUsedType z =  TpComplex(0, 1);
	MyUsedType u =  TpComplex(0.5, 0.5);
	
	TpMatrix  A(2,2), U(2,2);
	
	A << z, z,z,z;
	U << u, u, u, u;
	
	std::cout <<"# ------- Mark 1: ------- \n";
	MyType q = conj( z );
	
	std::cout <<"# ------- Mark 2: ------- \n";
	TpMatrix  R1 = U.adjoint() * A;
	
	std::cout <<"# ------- Mark 3: ------- \n";
	TpMatrix  R2 = U.adjoint();

	std::cout <<"# ------- Mark 4: ------- \n";
	A.adjointInPlace();

	std::cout <<"# ------- Mark 5: ------- \n";
	TpMatrix  R3 = A.conjugate();

	std::cout <<"# ------- Mark 6: ------- \n";

	std::cout << U(0,0) << " : " << R2(0,0) << std::endl;
	std::cout << A(0,0) << " : " << R3(0,0) << std::endl;
	
  return 0;
}

reports only one usage of conj() for MyType, for fixed size and dynamical matrices.
While the last 2 lines of output show, that for std::complex conjugation is performed.

Best regards,
Peter

#include <complex>
#include <iostream>

#include <Eigen/Dense>

// g++   -std=c++11 -g  AdjTest.C  -I $HOME/Eigen/eigen -o AdjTest
//  g++ -DEIGEN_DONT_VECTORIZE  -DEIGEN_UNROLLING_LIMIT=0  -std=c++11 -g  AdjTest.C  -I $HOME/Eigen/eigen -o AdjTest && ./AdjTest 


typedef std::complex<double> TpComplex;

class MyType
{
public:
	
	MyType() {};
	MyType( const TpComplex& z) : value(z) {};
	MyType( const double   & z) : value(z) {};
	MyType( const int      & z) : value(z) {};
	
	const TpComplex& z() const { return value; }
	TpComplex& z() { return value; }
	
	MyType& operator+=( const MyType& a ) { value += a.z(); return *this; }
	MyType& operator-=( const MyType& a ) { value -= a.z(); return *this; }
	MyType& operator*=( const MyType& a ) { value *= a.z(); return *this; }
	MyType& operator/=( const MyType& a ) { value /= a.z(); return *this; }
	
  private:
    TpComplex value;
};

MyType conj( const MyType & z )
{	static int counter(0);
	std::cout << "#diff conjugate called " << ++counter << " : " << z.z() << std::endl;
	return conj( z.z() );
}

MyType operator+( const MyType& a, const MyType& b ) { return a.z() + b.z() ; }
MyType operator-( const MyType& a, const MyType& b ) { return a.z() - b.z() ; }
MyType operator*( const MyType& a, const MyType& b ) { return a.z() * b.z() ; }
MyType operator/( const MyType& a, const MyType& b ) { return a.z() / b.z() ; }

std::ostream& operator<<(std::ostream &s,const MyType& a)
{
  s << a.z();
  return s;
};

//typedef std::complex<double> MyUsedType;   ///< Use this for std complex type
typedef MyType MyUsedType;                   ///< Use this for user defined complex type


//typedef Eigen::Matrix< MyUsedType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TpMatrix;
typedef Eigen::Matrix< MyUsedType, 2,2 > TpMatrix;


int main()
{
	MyUsedType z =  TpComplex(0, 1);
	MyUsedType u =  TpComplex(0.5, 0.5);
	
	TpMatrix  A(2,2), U(2,2);
	
	A << z, z,z,z;
	U << u, u, u, u;
	
	std::cout <<"# ------- Mark 1: ------- \n";
	MyType q = conj( z );
	
	std::cout <<"# ------- Mark 2: ------- \n";
	TpMatrix  R1 = U.adjoint() * A;
	
	std::cout <<"# ------- Mark 3: ------- \n";
	TpMatrix  R2 = U.adjoint();

	std::cout <<"# ------- Mark 4: ------- \n";
	A.adjointInPlace();

	std::cout <<"# ------- Mark 5: ------- \n";
	TpMatrix  R3 = A.conjugate();

	std::cout <<"# ------- Mark 6: ------- \n";

	std::cout << U(0,0) << " : " << R2(0,0) << std::endl;
	std::cout << A(0,0) << " : " << R3(0,0) << std::endl;
	
  return 0;
}