[Users] best way to maintain CCTK_COMPLEX precision with Fortran?

[Bernard Kelly]

Thanks, Erik.

I wasn't taking any magnitudes, but it looks like my problem was just
blindly using the "cmplx" function to form the complex variable
without a KIND argument, and it defaults to KIND=1 (i.e.
single-precision) even if the real & imaginary components enter as
double-precision.

For now, I've switched to the GNU "complex" method, which seems to
have fixed my problem on my Macbook.

Bernard

On 19 July 2015 at 14:52, Erik Schnetter <schnetter at cct.lsu.edu> wrote:
> Bernard
>
> In Fortran, real(z) (where z is complex) does maintain precision. The
> problem is likely caused by your using "-Ofast", "-ffast-math", or an Intel
> compiler that uses the equivalent by default (so you have to explicitly
> switch it off). With such an optimization, e.g. abs(z) is implemented as
> sqrt(real(z)**2 + aimag(z)**2), which overflows the intermediate squares.
> Without -ffast-math, compilers use a more expensive implementation that does
> not overflow.
>
> In your case, I'd assume that the problem is no the call to real, but other
> complex or real functions. Can you show us a larger part of the code, and
> the compiler vendor, version, and options you are using?
>
> The C functions CCTK_Cmplx* are deprecated. These days, C compilers
> implement complex numbers as defined by the standard, and you'd use the
> respective C standard function instead.
>
> -erik
>
>
> On Sun, Jul 19, 2015 at 2:16 PM, Bernard Kelly <physicsbeany at gmail.com>
> wrote:
>>
>> Hi all.
>>
>> The Cactus documentation shows a set of C routines CCTK_CmplxXXXXXX
>> that presumably preserve precision between CCTK_REAL and CCTK_COMPLEX
>> types. Is there a "right" way to handle this in Cactus with Fortran
>> (that isn't "recode in C")?
>>
>> Context: I have some Fortran code in a Cactus thorn that works with
>> real & complex fields. I've noticed that when the numbers are very
>> large in magnitude, taking a simple real(z) drops precision --- a
>> finite but large real part ends up as +/- Inf instead. I suspect it's
>> dropping to single precision at some stage. Switching to "dreal"
>> doesn't help, and avoiding the use of CCTK_COMPLEX altogether ---
>> while possible --- leads to longer, less readable code.
>>
>> Thanks,
>>
>> Bernard
>> _______________________________________________
>> Users mailing list
>> Users at einsteintoolkit.org
>> http://lists.einsteintoolkit.org/mailman/listinfo/users
>
>
>
>
> --
> Erik Schnetter <schnetter at cct.lsu.edu>
> http://www.perimeterinstitute.ca/personal/eschnetter/