search

roystgnr / MetaPhysicL

Metaprogramming and operator-overloaded classes for numerical simulations
Other
22 stars 12 forks source link

Make std::complex<DualNumber> work #66

Open dschwen opened 3 years ago

dschwen commented 3 years ago

We're looking into using the Eigen library to obtain eigen decompositions of matrices of dual numbers. This is of course because LAPACK precludes us from using dual numbers entirely while Eigen is beautifully templated.

I seem to have hit a snag though as there are a bunch of functions used in <complex> that are not overloaded (or specialized) for dual numbers. This results in errors like

/Users/schwd/miniconda3/envs/moose/bin/../include/c++/v1/complex:685:15: error: no matching function for call to 'scalbn'
    _Tp __y = scalbn((__b * __c - __a * __d) / __denom, -__ilogbw);

and

/Users/schwd/miniconda3/envs/moose/bin/../include/c++/v1/complex:676:29: error: no matching function for call to 'fabs'
    _Tp __logbw = logb(fmax(fabs(__c), fabs(__d)));
dschwen commented 3 years ago

Hm, I guess complex numbers are not necessary if we can assume the matices are symmetrical.

roystgnr commented 3 years ago

I do hope you've got another way out of this. DualNumber<complex> should be workable in theory, but std::complex<DualNumber> isn't; C++ says that Specializing the template std::complex for any type other than float, double, and long double is unspecified.

In my experience those specializations are often workable in practice, but not always; I must have hit a compiler that agreed with the standard instead of with me at some point, or I wouldn't have remembered the disagreement. To be fully portable you wouldn't be able to just specialize std::complex; you'd have to create MetaPhysicL::complex, as well as versions of all the functions you'd want to act on that.

roystgnr commented 3 years ago

Looks like even boost might have run into this problem? They at least try to specialize complex<boost::multiprecision::float128>, but when you can help it they recommend you use boost::multiprecision::complex128 instead, and I'd bet that compatibility issues with the former overloads are why.

cticenhour commented 3 years ago

I'm running into this issue more generally while converting some MOOSE objects that utilize std::complex into AD objects. I often utilize complex to make field calculations easier in my electromagnetics code. I think I'll be able to find a workaround, but I just wanted to bump this thread so you know that others are interested in any improvements here (even if I don't have the time to implement it myself right now 😢).

dschwen commented 3 years ago

So the answer would be to roll our own moose::complex number type?!

roystgnr commented 3 years ago

So the answer would be to roll our own moose::complex number type?!

I absolutely hate this, but: yes, probably. It would be a shim to std::complex for float/double/long double, and to a new MetaPhysicL::complex<T> for MetaPhysicL types, and to boost::multiprecision::complex128 for quad precision builds ...

And it would be incredibly annoying but I fear it would be the only way to guarantee things not breaking with some new compiler years down the road.

We could probably make std::complex<MetaPhysicL::foo> a shim to MetaPhysicL::complex<foo>, so that on compiler+library combinations which support it we'd get additional backwards compatibility with generic codes that explicitly reference std::complex rather than using context-dependent namespace lookup, but I'd want to rely on that as little as possible.

dschwen commented 3 years ago

Yeah, I was thinking fallback to std::complex for supported types would be best. I'm guessing there might be compiler optimizations hinging on the restrictive templating of std::complex

On Wed, Mar 10, 2021 at 8:29 AM roystgnr @.***> wrote:

So the answer would be to roll our own moose::complex number type?!

I absolutely hate this, but: yes, probably. It would be a shim to std::complex for float/double/long double, and to a new MetaPhysicL::complex for MetaPhysicL types, and to boost::multiprecision::complex128 for quad precision builds ...

And it would be incredibly annoying but I fear it would be the only way to guarantee things not breaking with some new compiler years down the road.

We could probably make std::complex a shim to MetaPhysicL::complex, so that on compiler+library combinations which support it we'd get additional backwards compatibility with generic codes that explicitly reference std::complex rather than using context-dependent namespace lookup, but I'd want to rely on that as little as possible.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/roystgnr/MetaPhysicL/issues/66#issuecomment-795605913, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABRMPWM25P2INU2IHIRFHTTC564JANCNFSM4T5CM4SA .