| [eigen] TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
Dear All,
I'm trying to extend my simple test code using eigen for a user defined type
including the diagonalization. Sadly, I got lost in the template machinery.
In short
template<typename T>
class MyType
{
public:
typedef T value_type;
// ...
protected:
T za, zb;
};
typedef double TpFloatBase;
typedef std::complex<TpFloatBase> TpComplex;
typedef MyType < TpFloatBase > TpRFloat;
typedef MyType < TpComplex > TpFloat;
Eigen::SelfAdjointEigenSolver<TpMatrix> eigensolverH;
eigensolverH.compute(H);
triggers
/home/peter/Eigen/eigen/Eigen/src/Core/TriangularMatrix.h:386:26: note: Eigen::TriangularViewImpl<_MatrixType, _Mode, Eigen::Dense>::TriangularViewType& Eigen::TriangularViewImpl<_MatrixType, _Mode, Eigen::Dense>::operator/=(const
typename Eigen::internal::traits<T>::Scalar&) [with _MatrixType = Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1>; unsigned int _Mode = 1u; Eigen::TriangularViewImpl<_MatrixType, _Mode,
Eigen::Dense>::TriangularViewType = Eigen::TriangularView<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1>, 1u>; typename Eigen::internal::traits<T>::Scalar = MyType<std::complex<double> >]
TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
^
/home/peter/Eigen/eigen/Eigen/src/Core/TriangularMatrix.h:386:26: note: no known conversion for argument 1 from ‘Eigen::SelfAdjointEigenSolver<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 1> >::RealScalar {aka
MyType<double>}’ to ‘const Scalar& {aka const MyType<std::complex<double> >&}’
I thought it would be sufficient to define
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
}
And indeed
TpRFloat r(1);
TpFloat z(2);
TpFloat q = z / r; // o.k.
compiles fine.
Could you please give be hint, which operator I'm missing?
I've attached the complete source file that triggers the error, including the error messages I get using
"g++ -std=c++11 -g Test.C -I $HOME/Eigen/eigen -o Test " .
Note, that some the operator definition doesn't make sense mathematically, I just put something to get it compiled.
Thanks in advance,
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>
//#include <Eigen/KroneckerProduct>
// g++ -std=c++11 -g Test.C -I $HOME/Eigen/eigen -o Test 2>&1 | tee Err.txt | more
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() {};
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; }
//template<class Q = T>
//typename boost::enable_if< boost::is_complex<Q>, MyType<Q> >::type& operator/=( const MyType< typename Q::value_type>& a) { return *this; };
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< 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< 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() )); };
///=====================================================
/// Needed for Eigen diagonalizaiton, e.g. we need a sorting
///=====================================================
TpComplex operator +( const TpComplex& r1, const int& r2)
{
return TpComplex(r1.real() + r2, r1.imag() ) ;
}
int operator ==( const TpComplex& r1, const int& r2)
{
return r1 == TpComplex(r2) ;
}
int operator <( const TpComplex& r1, const TpComplex& r2)
{
// return std::abs(r1) < std::abs(r2);
return r1.real() == r2.real() ? r1.imag() < r2.imag() : r1.real() < r2.real();
}
int operator <=( const TpComplex& r1, const TpComplex& r2)
{
// return std::abs(r1) <= std::abs(r2);
return r1.real() == r2.real() ? r1.imag() <= r2.imag() : r1.real() <= r2.real();
}
int operator >( const TpComplex& r1, const TpComplex& r2)
{
return r1.real() == r2.real() ? r1.imag() > r2.imag() : r1.real() > r2.real();
// return std::abs(r1) > std::abs(r2);
}
int operator >=( const TpComplex& r1, const TpComplex& r2)
{
// return std::abs(r1) >= std::abs(r2);
return r1.real() == r2.real() ? r1.imag() >= r2.imag() : r1.real() >= r2.real();
}
template<class T> MyType<T>
operator+(const MyType<T> a , const int b) { return MyType<T>( a.z1() + b , a.z2() );}
template<class T>
int operator>( const MyType<T>& a, const int n )// needed for eigen
{
return (a.z1() == T(0) ) ? ( a.z2()>T(0) ) : ( a.z1() > T(n) );
}
template<class T>
int operator==( const MyType<T>& a, const int n) // needed for eigen
{
return (a.z1() == T(n) ) && (a.z2() == T(0) );
};
template<class T>
int operator==( const MyType<T>& a, const MyType<T>& b ) // needed for eigen
{
return (a.z1() == b.z1() ) && (a.z2() == b.z2() );
};
template<class T>
int 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>
int operator!=( const MyType<T>& a, const int n) // needed for eigen
{
return (a.z1() != T(n) ) || (a.z2() != T(0) );
};
///=====================================================
/// Needed for Eigen IO
///=====================================================
namespace Eigen
{
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() ); }
};
namespace internal
{
template<> struct significant_decimals_impl<TpFloat> : public significant_decimals_impl<TpFloatBase>
{
};
}
}
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.
Eigen::SelfAdjointEigenSolver<TpMatrix> eigensolverH;
eigensolverH.compute(H);
return 0;
}In file included from /home/peter/Eigen/eigen/Eigen/Eigenvalues:39:0,
from /home/peter/Eigen/eigen/Eigen/Dense:7,
from Test.C:7:
/home/peter/Eigen/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h: In instantiation of ‘Eigen::SelfAdjointEigenSolver<MatrixType>& Eigen::SelfAdjointEigenSolver<_MatrixType>::compute(const Eigen::EigenBase<OtherDerived>&, int) [with InputType = Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 1>; _MatrixType = Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 1>]’:
Test.C:239:24: required from here
/home/peter/Eigen/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h:434:40: error: no match for ‘operator/=’ (operand types are ‘Eigen::MatrixBase<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1> >::TriangularViewReturnType<1u>::Type {aka Eigen::TriangularView<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1>, 1u>}’ and ‘Eigen::SelfAdjointEigenSolver<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 1> >::RealScalar {aka MyType<double>}’)
mat.template triangularView<Lower>() /= scale;
^
/home/peter/Eigen/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h:434:40: note: candidate is:
In file included from /home/peter/Eigen/eigen/Eigen/Core:437:0,
from /home/peter/Eigen/eigen/Eigen/Dense:1,
from Test.C:7:
/home/peter/Eigen/eigen/Eigen/src/Core/TriangularMatrix.h:386:26: note: Eigen::TriangularViewImpl<_MatrixType, _Mode, Eigen::Dense>::TriangularViewType& Eigen::TriangularViewImpl<_MatrixType, _Mode, Eigen::Dense>::operator/=(const typename Eigen::internal::traits<T>::Scalar&) [with _MatrixType = Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1>; unsigned int _Mode = 1u; Eigen::TriangularViewImpl<_MatrixType, _Mode, Eigen::Dense>::TriangularViewType = Eigen::TriangularView<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 0, -1, -1>, 1u>; typename Eigen::internal::traits<T>::Scalar = MyType<std::complex<double> >]
TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
^
/home/peter/Eigen/eigen/Eigen/src/Core/TriangularMatrix.h:386:26: note: no known conversion for argument 1 from ‘Eigen::SelfAdjointEigenSolver<Eigen::Matrix<MyType<std::complex<double> >, -1, -1, 1> >::RealScalar {aka MyType<double>}’ to ‘const Scalar& {aka const MyType<std::complex<double> >&}’