Re: [eigen] Custom complex type

[Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx>]

Hi.

I've finished custom complex scalar type development last week. And now I'm working on its integration to Eigen 3..0.1.

Although it seems to be smooth I encountered some difficulties which I cannot solve by myself. 

I would appreciate any advice on how to proceed. 

Just for recollection, I develop custom multiprecision scalar class mpreal (based on MPFR library) which is already supported by Eigen. 

Now I'm trying to create analogous type for multiprecision complex numbers based on low level MPC library. 

The simplest way would be to use  std::complex<mpreal>. But this approach disables some possibilities for optimization and I've created custom class

mpcomplex based on mpreal with API compatible with std::complex<> for tests. 

In order to integrate it with Eigen I wrote NumTraits<mpc::mpcomplex> specialization with basic stuff:

template<> struct NumTraits<mpc::mpcomplex>

: GenericNumTraits<mpc::mpcomplex>

{

typedef mpfr::mpreal Real;

typedef mpfr::mpreal NonInteger;

enum {

IsComplex = 1,

IsInteger = 0,

RequireInitialization = NumTraits<Real>::RequireInitialization,

ReadCost = 2 * NumTraits<Real>::ReadCost,

AddCost = 2 * NumTraits<Real>::AddCost,

MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost

};

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

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

};

Plus many specializations, like:

template<>  struct conj_impl<mpc::mpcomplex> {...}

template<>  struct real_impl<mpc::mpcomplex> {...}

.....

This kind of NumTraits was enough for mpreal to work with Eigen. But not for mpcomplex.

In order to compile my Eigen-based code I had to add definitions for various kinds of conj_helper<>:

template<> struct conj_helper<mpc::mpcomplex, mpc::mpcomplex,...>

And 

template<> struct scalar_product_traits<mpc::mpcomplex, mpc::mpcomplex>

Code compiled well after that, but after my tests I noticed that in _expression_ (pseudo inverse calculation in SVD):

m_matrixV * invertedSingVals.asDiagonal() * m_matrixU.adjoint();

m_matrixU.adjoint() works just as plain transpose() despite that m_matrixU is complex.  invertedSingVals is of mpreal type.

If I use m_matrixU.adjoint().eval() instead everything is all right (as well as adjointInPlace()) .

I guess some mpcomplex-specialization is missing for proper matrix _expression_ template evaluation.

Could you advise what it could be? 

NumTraits<mpcomplex> is in attach. mpcomplex class as simple as: 

class mpcomplex{

private:

mpreal re;

mpreal im;

// Standard API

.....

Thanks in advance.

Pavel Holoborodko.

P.S.

std::complex<mpreal> works nicely without any shaman dances.

On Wed, Jun 8, 2011 at 12:13 PM, Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx> wrote:

Hi Gael,

Thanks for answers everyone.

Gael, I'll try this way (custom mpcomplex + std comp. API) with hope that it will work with Eigen "in practice" too.

I'll let you know about the results.

On Wed, Jun 8, 2011 at 3:43 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Hi Pavel,

in theory you should be able to write your own mpcomplex class, and
define NumTraits<mpcomplex>::IsComplex to true as long as you
implement a STD compatible API.

gael

On Tue, Jun 7, 2011 at 8:04 AM, Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx> wrote:
> Hello,
> Eigen's magicians, please advise on how to implement custom complex type for
> smooth integration with Eigen.
> I plan to create multi-precision C++ complex scalar class based on
> MPC: http://www.multiprecision.org/index.php?prog=mpc.
> MPC is cousin of MPFR from the same authors - http://mpfr.org.
> There are two possible ways to proceed:
> (1).  Add template specialization for std::complex<mpfr::mpreal>.
> (2).  Create own scalar class something like mpcomplex based on low-level
> MPC.
> Second option seems attractive from optimization point of view - minimum
> overhead, no conversion mpreal <-> MPC low level stuff is required.
> There is only one problem - will it be possible to integrate such complex
> type with Eigen?
> Documentation on NumTraits says:
> "An enum value IsComplex. It is equal to 1 if T is a std::complex type, and
> to 0 otherwise."
> Does it mean that only std::complex<T> custom complex types are supported in
> Eigen?
> If yes, do I need to provide something besides
> NumTraits<std::complex<mpreal> > specialization?
> Thanks in advance,
> Pavel.
>
>
>
>

#ifndef EIGEN_MPCOMPLEXSUPPORT_MODULE_H
#define EIGEN_MPCOMPLEXSUPPORT_MODULE_H

#include <Eigen/Core>
#include <mpcomplex.h>
#include "MPRealSupport.h"

namespace Eigen {
	template<> struct NumTraits<mpc::mpcomplex>
	: GenericNumTraits<mpc::mpcomplex>
	{
		typedef mpfr::mpreal Real;
		typedef mpfr::mpreal NonInteger;

		enum {
			IsComplex = 1,
			IsInteger = 0,
			RequireInitialization = NumTraits<Real>::RequireInitialization,
			ReadCost = 2 * NumTraits<Real>::ReadCost,
			AddCost = 2 * NumTraits<Real>::AddCost,
			MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
		};
		
		inline static Real epsilon() { return NumTraits<Real>::epsilon(); }
		inline static Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
	};

	namespace internal {

		template<>	struct conj_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x) 
			{ return mpc::conj(x); } 
		};

		template<> struct real_impl<mpc::mpcomplex>
		{ 
			inline static const mpfr::mpreal run(const mpc::mpcomplex& x) 
			{ return mpc::real(x); } 
		};

		template<> struct imag_impl<mpc::mpcomplex> 
		{ 
			inline static const mpfr::mpreal run(const mpc::mpcomplex& x)   
			{ return mpc::imag(x); } 
		};

		template<> struct abs_impl<mpc::mpcomplex> 
		{ 
			inline static const mpfr::mpreal run(const mpc::mpcomplex& x)  
			{ return mpc::abs(x); } 
		};

		template<> struct abs2_impl<mpc::mpcomplex> 
		{ 
			inline static const mpfr::mpreal run(const mpc::mpcomplex& x) 
			{ return mpc::norm(x); } 
		};

		template<> struct sqrt_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x) 
			{ return mpc::sqrt(x); } 
		};

		template<> struct exp_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::exp(x); } 
		};

		template<> struct log_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::log(x); } 
		};

		template<> struct sin_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::sin(x); } 
		};

		template<> struct cos_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::cos(x); } 
		};

		template<> struct pow_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x, const mpc::mpcomplex& y)  
			{ return mpc::pow(x, y); } 
		};

		template<> struct tan_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::tan(x); } 
		};

		template<> struct acos_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::acos(x); } 
		};

		template<> struct asin_impl<mpc::mpcomplex> 
		{ 
			inline static const mpc::mpcomplex run(const mpc::mpcomplex& x)  
			{ return mpc::asin(x); } 
		};

		template<> mpc::mpcomplex random<mpc::mpcomplex>()
		{
			return mpc::random();
		}

		template<> mpc::mpcomplex random<mpc::mpcomplex>(const mpc::mpcomplex& x, const mpc::mpcomplex& y)
		{
			// Call RNG from NumTraits<mpfr::mpreal>
			return mpc::mpcomplex(random(real(x), real(y)),
								  random(imag(x), imag(y)));
		}

		template<> struct conj_helper<mpc::mpcomplex, mpc::mpcomplex, false, false>
		{
			typedef mpc::mpcomplex Scalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
			{ return x*y; }
		};

		template<> struct conj_helper<mpc::mpcomplex, mpc::mpcomplex, true, false>
		{
			typedef mpc::mpcomplex Scalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
			{ return mpc::conj(x)*y; }
		};

		template<> struct conj_helper<mpc::mpcomplex, mpc::mpcomplex, false, true>
		{
			typedef mpc::mpcomplex Scalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
			{ return x*mpc::conj(y); }
		};

		template<> struct conj_helper<mpc::mpcomplex, mpc::mpcomplex, true, true>
		{
			typedef mpc::mpcomplex Scalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
			{ return mpc::conj(x)*mpc::conj(y); }
		};

		template<> struct conj_helper<mpc::mpcomplex, mpfr::mpreal, true, false>
		{
			typedef mpc::mpcomplex Scalar;
			typedef mpfr::mpreal RealScalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
			{ return mpc::conj(x)*y; }
		};

		template<> struct conj_helper<mpc::mpcomplex, mpfr::mpreal, false, false>
		{
			typedef mpc::mpcomplex Scalar;
			typedef mpfr::mpreal RealScalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
			{ return x*y; }
		};

		template<> struct conj_helper<mpfr::mpreal, mpc::mpcomplex, false, true>
		{
			typedef mpc::mpcomplex Scalar;
			typedef mpfr::mpreal RealScalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
			{ return x*mpc::conj(y); }
		};

		template<> struct conj_helper<mpfr::mpreal, mpc::mpcomplex, false, false>
		{
			typedef mpc::mpcomplex Scalar;
			typedef mpfr::mpreal RealScalar;

			EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
			{ return c + pmul(x,y); }

			EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
			{ return x*y; }
		};

		template<typename T> struct scalar_product_traits<T, mpc::mpcomplex>
		{
			typedef mpc::mpcomplex ReturnType;
		};

		template<typename T> struct scalar_product_traits<mpc::mpcomplex, T>
		{
			typedef mpc::mpcomplex ReturnType;
		};

		template<> struct scalar_product_traits<mpc::mpcomplex, mpc::mpcomplex>
		{
			typedef mpc::mpcomplex ReturnType;
		};

	} // namespace internal
	
} // namespace Eigen
#endif // EIGEN_MPCOMPLEXSUPPORT_MODULE_H