Re: [eigen] TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other)

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


It seems that you also need to specialize:

template<typename T> struct ScalarBinaryOpTraits<MyType< std::complex<T> >,MyType<T> > { typedef MyType< std::complex<T> > ReturnType; };
template<typename T> struct ScalarBinaryOpTraits<MyType<T>,MyType< std::complex<T> > > { typedef MyType< std::complex<T> > ReturnType; };

then proper comparisons for reals, a few mixed real and complex operations (-, *, /), and some -=, /=, *= operators, and everything compiles fine use Eigen's head.

See attached file.

Cheers,
Gael

On Mon, Sep 5, 2016 at 8:23 PM, Peter <list@xxxxxxxxxxxxxxxxx> wrote:
Dear Gael,



Am 05.09.2016 um 18:17 schrieb Gael Guennebaud:
Hi,

the types TpFloat and NumTraits<TpFloat>::Real should be fully compatible, just like std::complex<T> and T. Currently, we cannot even do:

TpFloat c;
NumTraits<TpFloat>::Real r;
c = r;

 After adding a constructor

template< typename Q = T, class = typename boost::enable_if< boost::is_complex<Q>, typename Q::value_type  > >
MyType( const MyType < typename Q::value_type>& a ) {};

your mentioned assignment actually compiles.
Including eigensolverH.compute() produces then > 15000 lines of error, so this will take me some time.

So the first step is to make them compatible and/or define NumTraits<TpFloat>::Real to be TpFloatBase ???

 Well, that would be in contradiction to < https://eigen.tuxfamily.org/dox/structEigen_1_1NumTraits.html > :
" A typedef Real, giving the "real part" type of T. If T is already real, then Real is just a typedef to T. If T is std::complex<U> then Real is a typedef to U. "

This should be enough to compile with the devel branch.

 O.k., I'll try harder.

Best regards,
Peter



#include <stdlib.h>

#include <cmath>

#include <complex>

#include <Eigen/Dense>
#include <Eigen/Sparse>

#include <boost/utility/enable_if.hpp> //needed for enable_if
#include <boost/type_traits/is_complex.hpp>

template<typename T>
class MyType
{
  public:

    typedef T value_type;
	
    MyType( const T a, const T b ) : za(a), zb(b) {};
    MyType( const MyType& a) : za(a.za), zb(a.zb) {};
    MyType( const T& z) : za(z), zb(0) {};
    MyType( const int z) : za(z), zb(0) {}; //needed for eigen
    MyType() {};

    template< typename Q = T, class = typename boost::enable_if< boost::is_complex<Q>, typename Q::value_type  > >
    MyType( const MyType < typename Q::value_type>& a ) {};

    const T& z1() const { return za; };
    T& z1() { return za; };

    const T& z2() const { return zb; };
    T& z2() { return zb; };

    MyType& operator+=(const T a)      { za += a; return *this; }
    MyType& operator+=(const MyType a) { za += a.za; zb += a.zb; return *this; }

    MyType& operator-=(const T a)      { za -= a; return *this; }
    MyType& operator-=(const MyType a) { za -= a.za; zb -= a.zb; return *this; }

    MyType& operator/=(const T a)      { za /= a; return *this; }
    MyType& operator/=(const MyType a) { za /= a.za; zb /= a.zb; return *this; }

    MyType& operator*=(const T a)      { za *= a; return *this; }
    MyType& operator*=(const MyType a) { za *= a.za; zb /= a.zb; return *this; }

    bool operator!=( const MyType& b ) { return (z1() != b.z1() ) || (z2() != b.z2() ); }

protected:
    T za, zb;
};


typedef   double TpFloatBase;
typedef   std::complex<TpFloatBase> TpComplex;

typedef   MyType < TpFloatBase > TpRFloat;
typedef   MyType < TpComplex   > TpFloat;


///=====================================================
/// Eigen typedefs
///=====================================================


typedef Eigen::Matrix< TpFloat,   Eigen::Dynamic, 1> TpVector;
typedef Eigen::Matrix< TpFloat,   Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TpMatrix;

typedef Eigen::Array< TpFloat  , Eigen::Dynamic, 1>  TpArray;

template<class T> MyType<T> operator+(const  T a , const MyType<T> b)           { return MyType<T>( a + b.z1() , b.z2() );}
template<class T> MyType<T> operator+(const  MyType<T> a , const T b)           { return MyType<T>( a.z1() + b , a.z2() );}
template<class T> MyType<T> operator+(const  MyType<T> a , const MyType<T> b)   { return MyType<T>( a.z1() + b.z1() , a.z2() + b.z2() );}

template<class T> MyType<T> operator-(const  T a , const MyType<T> b)           { return MyType<T>( a - b.z1() , - b.z2() );}
template<class T> MyType<T> operator-(const  MyType<T> a ,const  T b)           { return MyType<T>( a.z1() - b , - a.z2() );}
template<class T> MyType<T> operator-(const  MyType<T> a ,const  MyType<T> b)   { return MyType<T>( a.z1() - b.z1() , a.z2() - b.z2() );}

template<class T> MyType<T> operator*( const T a , const MyType<T> b)           { return MyType<T>( a * b.z1() , a*b.z2() ); }
template<class T> MyType<T> operator*( const MyType<T> a , const  T b)          { return MyType<T>( a.z1() * b , a.z2()*b ); }
template<class T> MyType<T> operator*( const MyType<T> a , const MyType<T> b)   { return MyType<T>( a.z1() * b.z1() , a.z2()*b.z2() ) ; }

template<class T> MyType<T> operator/( const MyType<T> a , const MyType<T> b)   { return MyType<T>( a.z1() / b.z1() , a.z2() / b.z2() ) ; }

template<class T>
MyType< std::complex<T> > operator-( const MyType< std::complex<T> >&  a , const MyType<T>& b)
{
  return MyType< std::complex<T> >(a.z1()-b.z1(), a.z2()-b.z2());
}

template<class T>
MyType< std::complex<T> > operator-( const MyType< T >&  a , const MyType< std::complex<T> >& b)
{
  return MyType< std::complex<T> >(a.z1()-b.z1(), a.z2()-b.z2());
}

template<class T> 
MyType< std::complex<T> > operator*( const MyType< std::complex<T> >&  a , const MyType<T>& b)
{ 
	return MyType< std::complex<T> >(a.z1()*b.z1(), a.z2()*b.z2());
}

template<class T>
MyType< std::complex<T> > operator*( const MyType< T >&  a , const MyType< std::complex<T> >& b)
{
  return MyType< std::complex<T> >(a.z1()*b.z1(), a.z2()*b.z2());
}

template<class T> 
MyType< std::complex<T> > operator/( const MyType< std::complex<T> >&  z , const MyType< std::complex<T> > & q)  
{ 
	return MyType< std::complex<T > >( z.z1() / q.z1() , z.z2() / q.z2() ); 
}

template<class T> 
MyType< std::complex<T> > operator/( const MyType< std::complex<T> >&  z , const MyType<T>& b)  
{ 
	return z / MyType< std::complex<T > >( b.z1() , b.z2()  ); ///< Could be optimized
}

template<class T> 
MyType< std::complex<T> > operator/( const MyType<T>& a, const MyType< std::complex<T> >&  z)  
{ 
	return MyType< std::complex<T > >( a.z1() , a.z2()  ) / z ; ///< Could be optimized
}

template<class Q> 
typename boost::enable_if< boost::is_complex<Q>, MyType<Q> >::type conj( const MyType< Q >& a) { return MyType<Q>( conj( a.z1() ), conj( a.z2() ));  };

template<class Q> 
typename boost::enable_if< boost::is_complex<Q>, MyType<typename Q::value_type> >::type real( const MyType< Q >& a) { return MyType<typename Q::value_type>( real( a.z1() ), real( a.z2() ));  };

template<class Q> 
typename boost::enable_if< boost::is_complex<Q>, MyType<typename Q::value_type> >::type imag( const MyType< Q >& a) { return MyType<typename Q::value_type>( imag( a.z1() ), imag( a.z2() ));  };

template<class Q> 
typename boost::disable_if< boost::is_complex<Q>, MyType<Q> >::type abs(const MyType<Q> a)
{
  return (a.z1() == 0 ) ? MyType<Q>( (a.z1() , Q(0) )) :  (a.z1() > 0 ) ?  a : -a;
};

template<class Q> 
typename boost::enable_if< boost::is_complex<Q>, MyType<typename Q::value_type> >::type abs( const MyType< Q >& a) { return MyType<typename Q::value_type>( abs( a.z1() ), abs( a.z2() ));  };


template<class T>
bool operator==( const MyType<T>& a, const MyType<T>& b )
{
  return (a.z1() == b.z1() ) && (a.z2() == b.z2() );
};


template<class T>
bool operator<( const MyType<T>& a, const MyType<T>& b)
{
  return (a.z1() == b.z1()) ? (a.z2()<b.z2()) : (a.z1()<b.z1());
};

template<class T>
bool operator>( const MyType<T>& a, const MyType<T>& b)
{
  return (a.z1() == b.z1()) ? (a.z2()>b.z2()) : (a.z1()>b.z1());
};

template<class T>
bool operator<=( const MyType<T>& a, const MyType<T>& b)
{
  return (a.z1() == b.z1()) ? (a.z2()<=b.z2()) : (a.z1()<=b.z1());
};

template<class T>
bool operator>=( const MyType<T>& a, const MyType<T>& b)
{
  return (a.z1() == b.z1()) ? (a.z2()>=b.z2()) : (a.z1()>=b.z1());
};

template<typename T>
MyType<T> operator-( const MyType<T>& x)
{
  return MyType<T>(-(x.z1()), -(x.z2()));
};

template<typename T>
MyType<T> sqrt( const MyType<T>& x)
{
  using std::sqrt;
  return MyType<T>(sqrt(x.z1()), sqrt(x.z2()));
};


///=====================================================
/// Needed for Eigen IO
///=====================================================

namespace Eigen 
{

  template<typename T> struct ScalarBinaryOpTraits<MyType< std::complex<T> >,MyType<T> > { typedef MyType< std::complex<T> > ReturnType; };
  template<typename T> struct ScalarBinaryOpTraits<MyType<T>,MyType< std::complex<T> > > { typedef MyType< std::complex<T> > ReturnType; };

	template<> struct NumTraits<TpFloat> : NumTraits<TpComplex> // permits to get the epsilon, dummy_precision, lowest, highest functions
	{
			typedef TpRFloat Real;
			typedef TpFloat  NonInteger;
			typedef TpFloat  Nested;

			static inline Real epsilon()         { return TpRFloat( NumTraits<TpFloatBase>::epsilon()          ); }
			static inline Real dummy_precision() { return TpRFloat( NumTraits<TpFloatBase>::dummy_precision()  ); }
	};

}

int main( const int argc, const char *argv[] )
{
	const int M=2;
	
	TpMatrix H(M,M);
	
	H << 1,2,2,4;
	
	TpMatrix rh0  = H.adjoint();      // o.k.
	TpMatrix rho1 =  H * H;           // o.k.
	TpMatrix rho2 =  H + H.adjoint(); // o.k.

	TpMatrix rho3 =  H * H.adjoint(); // o.k., fixed in development branch of  eigen.
	
	///-----------------------------------------------
	
	TpRFloat r(1);
	TpFloat  z(2);
	TpFloat  q = z / r; // o.k.

	TpMatrix A;
  TpFloat ca;
  TpRFloat a;
  ca = a;

  A.triangularView<Eigen::Upper>() /= a;
	
	Eigen::SelfAdjointEigenSolver<TpMatrix> eigensolverH;
	eigensolverH.compute(H);

	return 0;
}

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