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

[Peter <list@xxxxxxxxxxxxxxxxx>]

Dear All,

here's a "corrected" example.
I forgot to add the NumTraits specialization.

I've now added

template<> struct NumTraits<MyType> : NumTraits<TpComplex> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
		typedef TpComplex::value_type Real;

};

and the test

std::cout << "# Complex: " <<  Eigen::NumTraits<MyUsedType>::IsComplex << std::endl;

prints a "1".

The behaviour is not changed.

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;

template< class TpR, class TpM >
void Transpose( TpR& R, const TpM & M )
{
	for( int x=0; x<M.rows(); ++x) 
		for( int y=0; y<M.cols(); ++y ) 
			R(x,y) = conj( M(y,x) );

}

namespace Eigen
{
template<> struct NumTraits<MyType> : NumTraits<TpComplex> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
		typedef TpComplex::value_type Real;

};
}

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 << "# Complex: " <<  Eigen::NumTraits<MyUsedType>::IsComplex << std::endl;
	
	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";
	TpMatrix  R4(2,2);
	Transpose( R4, U );

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

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