Re: [eigen] Custom complex type

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Mon, Jun 20, 2011 at 11:38 AM, Pavel Holoborodko
<pavel@xxxxxxxxxxxxxxx> wrote:
> 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>

hm, I was afraid of that. This part of Eigen's internals requires some
love. In theory, you should not have to do such specializations.

> 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?

I really have to clean up all this mess with complexes!  Could you
check whether m_matrixV * m_matrixU.adjoint() is equivalent to
m_matrixV * m_matrixU.transpose()?

gael

> 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.
>>> >
>>> >
>>> >
>>> >
>>>
>>>
>>
>
>