IVF 10 Bug: Wrong branch returned by SQRT for COMPLEX*16 argument

I've been corresponding with a support person for Intel Visual Fortran
for about 3 ws regarding a bug in IVF 10 that did not occur for me
in IVF 9.  I have a complicated Fortran 95 code that contains a SQRT
function applied to a COMPLEX*16 argument that happens to have a
negative real part and zero imaginary part.  When I compile this code
using IVF 10 and with optimization enabled (any level), the result of
the SQRT has zero real part, as expected, but negative imaginary
part.  This, of course, is incorrect.  The SQRT (Fortran) function is
supposed to return values from the principal branch of the complex-
valued (mathematical) square root function, which in this case means
that the imaginary part should be positive.

The support person discovered that if he uses the option /fp:precise
when invoking the compiler, then the correct result is obtained, and
he would like to close the issue.  I maintain that using /fp:precise
is an acceptable temporary workaround, but that the incorrect result
obtained without using it indicates a serious compiler bug that should
be fixed.  According to my understanding of the /fp:precise flag, its
function is to control rounding and order of evaluation, and that it
should have only minor effects on the accuracy of an obtained result.
This bug causes one to obtain zero significant digits of accuracy in
many cases!  Also, if I need to use /fp:precise on this code, I will
not be able to obtain maximum speed, which is important in this
application (it is a very CPU intensive computational electromagnetics
code, requiring hours to run).

I would like to know if my position on this issue is correct, and if
so, how to convince the support person that it is.

Just in case Steve Lionel or someone else from Intel reads this, the
issue number is 443980.

Thanks,
Peter

Peter wrote:
(snip)

> When I compile this code
> using IVF 10 and with optimization enabled (any level), the result of
> the SQRT has zero real part, as expected, but negative imaginary
> part.  This, of course, is incorrect.  The SQRT (Fortran) function is
> supposed to return values from the principal branch of the complex-
> valued (mathematical) square root function, which in this case means
> that the imaginary part should be positive.

From Fortran 2003:

"A result of type complex is the principal value with the real
part greater than or equal to zero. When the real part of the
result is zero, the imaginary part has the same sign as the
imaginary part of X."

Is the imaginary part of X really zero, or just very small (and
negative)?

Note that due to rounding, it could be small and negative even
when you don't expect that it could happen.

-- glen

On Aug 21, 12:31 pm, Peter <petersamsim...@hotmail.com> wrote:
> The SQRT (Fortran) function is
> supposed to return values from the principal branch of the complex-
> valued (mathematical) square root function, which in this case means
> that the imaginary part should be positive.

This statement was true for Fortran 95 but not for Fortran 2003, and
Intel Fortran is following F2003 here.  The new text is:

"When the real part of the result is zero, the imaginary part has the
same sign as the imaginary part of X".

Steve

On Aug 21, 10:54 am, glen herrmannsfeldt <g...@ugcs.caltech.edu>
wrote:
> Peter wrote:
>
> (snip)
> 
>
>  From Fortran 2003:
>
>     "A result of type complex is the principal value with the real
>      part greater than or equal to zero. When the real part of the
>      result is zero, the imaginary part has the same sign as the
>      imaginary part of X."
>
> Is the imaginary part of X really zero, or just very small (and
> negative)?
>
> Note that due to rounding, it could be small and negative even
> when you don't expect that it could happen.
>
> -- glen

It is identically zero.  The quantity being computed is the complex
wavenumber z-component in a lossless medium.  In the absence of
losses, the real part is identically zero. Here is a portion of the
code and the resulting printout when compiled with optimization:

! Obtain wavenumber (Region i) squared (rad^2/m^2):
ksq = k0sq * layers(i)%mur * layers(i)%epsr
gamma = SQRT(betasq - ksq) ! Atten. constant (neper/m).
IF ((REAL(gamma) < 0.0_WP) .OR. &
& REAL(gamma) == 0.0_WP .AND. AIMAG(gamma) < 0.0_WP) THEN
WRITE(0,*) 'Routine: gsm_multi_slab_interface: Bad sqrt: arg = ', &
&  betasq-ksq, ', sqrt(arg) = ', gamma
gamma = -gamma
END IF


Setting Phi to     0.000 deg
Setting Theta to     0.000 deg
Analyzing frequency 2.2000 GHz
Routine: setup_modes: Bad sqrt: arg =
(-2126.00211632318,0.000000000000000E+000) , sqrt(arg) =
(0.000000000000000E+000,-46.1085904829370)


--Peter

On Aug 21, 10:25 am, Peter <petersamsim...@hotmail.com> wrote:
> On Aug 21, 10:54 am, glen herrmannsfeldt <g...@ugcs.caltech.edu>
> wrote:
>
>
> 
> 
> 
> 
> 
> 
> 
> 
>
> It is identically zero.  The quantity being computed is the complex
> wavenumber z-component in a lossless medium.  In the absence of
> losses, the real part is identically zero. Here is a portion of the
> code and the resulting printout when compiled with optimization:
>
> ! Obtain wavenumber (Region i) squared (rad^2/m^2):
> ksq = k0sq * layers(i)%mur * layers(i)%epsr
> gamma = SQRT(betasq - ksq) ! Atten. constant (neper/m).
> IF ((REAL(gamma) < 0.0_WP) .OR. &
>   & REAL(gamma) == 0.0_WP .AND. AIMAG(gamma) < 0.0_WP) THEN
>   WRITE(0,*) 'Routine: gsm_multi_slab_interface: Bad sqrt: arg = ', &
>            &  betasq-ksq, ', sqrt(arg) = ', gamma
>   gamma = -gamma
> END IF
>
> Setting Phi to     0.000 deg
> Setting Theta to     0.000 deg
> Analyzing frequency 2.2000 GHz
>  Routine: setup_modes: Bad sqrt: arg =
>  (-2126.00211632318,0.000000000000000E+000) , sqrt(arg) =
>  (0.000000000000000E+000,-46.1085904829370)
>
> --Peter

I said above that the real part of the argument of SQRT is zero.  Of
course I meant the imaginary part.

glen herrmannsfeldt wrote:
...
>    "A result of type complex is the principal value with the real
>     part greater than or equal to zero. When the real part of the
>     result is zero, the imaginary part has the same sign as the
>     imaginary part of X."
>
> Is the imaginary part of X really zero, or just very small (and
> negative)?

Since I wrote that particular wording, I can tell you
with some confidence that the imaginary part of the
result will be negative if the imaginary part of X is
negative - including minus zero.  So, for example
SQRT(-(1.0,0.0)) is CMPLX(0.0, -1.0), while
SQRT((-1.0, 0.0)) is CMPLX(0.0, 1.0).

Pedantic application of rules is always part of
computing.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies."   --  C. A. R. Hoare

On Aug 21, 10:07 am, Steve Lionel <steve.lio...@intel.com> wrote:
> On Aug 21, 12:31 pm, Peter <petersamsim...@hotmail.com> wrote:
> 
>
> This statement was true for Fortran 95 but not for Fortran 2003, and
> Intel Fortran is following F2003 here.  The new text is:
>
> "When the real part of the result is zero, the imaginary part has the
> same sign as the imaginary part of X".
>
> Steve

Thanks for your prompt reply.  In my case the imaginary part of X is
zero.  What then is the proper interpretation of the phrase "the sign
of the imaginary part of X"?  To my way of thinking, if you have to
assign a sign to 0, it must be positive, not negative!

--Peter

On Tue, 21 Aug 2007 10:54:57 -0700, glen herrmannsfeldt wrote
(in article <ZaadnSwS5-GwhVbbnZ2dnUVZ_sGvnZ2d@comcast.com> ):

> Peter wrote:
> (snip)
> 

> Is the imaginary part of X really zero, or just very small (and
> negative)?
>
> Note that due to rounding, it could be small and negative even
> when you don't expect that it could happen.

That ocurred to me also, particularly given the bit about /fp:precise being 
a
workaround. If you have code that depends on last-bit accuracy of
floating-point computations, in particular on getting a value that is
precisely zero from something like subtraction, then I'd call that a bug in
the code.

If the value is really, truly zero, then it does sound from what Glen quoted
like a compiler violation of the standard. I always have to go check the
standard for those branch cut things because the answers aren't always what
one might expect (and that gets debated here on occasion). But it looks like
Glen did the checking for you.

An alternative workaround, by the way, might be to write your own specific
for sqrt, wrapping around the compiler's version and testing for the
problematic case. The wrapper will slow down the complex sqrt evaluation, bu
t
shouldn't have all the other side effects of /fp:precise.

--
Richard Maine                    | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle           |  -- Mark Twain


 ________________________________________
________

Hogwasher, Premier News and Mail for OS X
http://www.asar.com/cgi-bin/product.../hogwasher.html
 ________________________________________
________

> Thanks for your prompt reply.  In my case the imaginary part of X is
> zero.  What then is the proper interpretation of the phrase "the sign
> of the imaginary part of X"?  To my way of thinking, if you have to
> assign a sign to 0, it must be positive, not negative!

Zero can be positive or negative.

--
FX

Peter wrote:


> Thanks for your prompt reply.  In my case the imaginary part of X is
> zero.  What then is the proper interpretation of the phrase "the sign
> of the imaginary part of X"?  To my way of thinking, if you have to
> assign a sign to 0, it must be positive, not negative!

How sure are you that it is zero?

Is it possible that the value you think is zero isn't the
value that is getting used for the SQRT?

-- glen