Re: [eigen] Custom complex type

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Basically you need a functor and to add it to the API. See for
instance changeset 3829 for an example of the addition of asin and
acos. (It contains some other minor modifications and the
factorization of the declaration of standard function via a macro).
The vectorization parts, with the p* functions are optional.

Gael

On Tue, Jun 21, 2011 at 9:44 AM, Pavel Holoborodko
<pavel@xxxxxxxxxxxxxxx> wrote:
> I can and already did - now it works fine with Eigen.
> In previous post I've just explained my reasons why I wanted
> non-std-non-mpreal-based mpcomplex.
> I already switched to std::complex<mpreal> - since it is the only way to
> integrate custom complex type in Eigen.
> BTW, could you give me hints on how to add coeff-wise math functions to
> Eigen?
> Like hyperbolic trigonometric functions or others? So I can use them as
> usual, e.g. m.sqrt().
> Would it be enough to just add them to MathFunctions.h or do I have to
> change something else?
> On Tue, Jun 21, 2011 at 4:24 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
> wrote:
>>
>> Hi,
>>
>> I don't really understand why you cannot specialize std::complex for
>> mpreal, e.g:
>>
>> namespace std {
>> template<> class complex<mpreal> {
>> // copy paste here the content of your current mpcomplex
>> };
>> }
>>
>> gael
>>
>> On Tue, Jun 21, 2011 at 4:06 AM, Pavel Holoborodko
>> <pavel@xxxxxxxxxxxxxxx> wrote:
>> > 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.
>> >> >>>> >
>> >> >>>> >
>> >> >>>> >
>> >> >>>> >
>> >> >>>>
>> >> >>>>
>> >> >>>
>> >> >>
>> >> >>
>> >> >
>> >>
>> >>
>> >
>> >
>>
>>
>
>