Re: [eigen] Custom complex type

[Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx>]

I see, thank you for your prompt response.

Yes, std::complex<mpreal> specialization will be much easier than modification of Eigen for custom complex type.

I wanted mpcomplex to be based on low level C library, which uses POD types for real and imag parts (C pointers actually).  

It would be the most efficient to use them directly without wrapping to C++ numeric class. 

Since it is impossible to specialize std::complex with POD - I had to write my own mpcomplex. 

Now, as Eigen heavily based on std::complex I have no other way than to make std::complex<mpreal> class specialization.

I think it is acceptable and maybe even better solution in a long run (e.g. compatibility with other libraries).  

Thank you for your help.

Pavel Holoborodko.

On Tue, Jun 21, 2011 at 4:29 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

ok, got it. The pb is in gebp_traits
(src/Core/products/GeneralBlockPanelKernel.h) where I have
specializations on std::complex to identify complexes.

Just a remark: cannot you simply specialize std::complex<mpreal> and
implement your optimization there ?

gael

On Mon, Jun 20, 2011 at 9:21 PM, Gael Guennebaud

<gael.guennebaud@xxxxxxxxx> wrote:
> 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.
>>>> >
>>>> >
>>>> >
>>>> >
>>>>
>>>>
>>>
>>
>>
>