Home > Archive > Fortran > August 2007 > IVF 10 Bug: Wrong branch returned by SQRT for COMPLEX*16 argument
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
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 w s 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
| |
| glen herrmannsfeldt 2007-08-21, 7:08 pm |
| 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
| |
| Steve Lionel 2007-08-21, 7:08 pm |
| 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.
| |
| James Giles 2007-08-21, 7:08 pm |
| 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
| |
| Richard Maine 2007-08-21, 7:08 pm |
| On Tue, 21 Aug 2007 10:54:57 -0700, glen herrmannsfeldt wrote
(in article <ZaadnSwS5-GwhVbbnZ2dnUVZ_sGvnZ2d@comcast.com> ):
> Peter wrote:
> (snip)
>
[color=darkred]
> 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, but
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
| |
| glen herrmannsfeldt 2007-08-21, 7:08 pm |
| 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
| |
| glen herrmannsfeldt 2007-08-21, 7:08 pm |
| James Giles wrote:
(snip, I wrote)
[color=darkred]
> 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).
Even though this was a big discussion not so long ago,
I didn't think of it at the time. A negative zero would
give a negative imaginary part for the result and, likely,
print out as a positive zero. This sounds like another
reason for printing negative zeros with a minus sign.
(The first post didn't say how he knew it was zero.)
> Pedantic application of rules is always part of
> computing.
Yep!
-- glen
| |
|
| On Aug 21, 10:46 am, Richard Maine <nos...@see.signature> wrote:
> On Tue, 21 Aug 2007 10:54:57 -0700, glen herrmannsfeldt wrote
> (in article <ZaadnSwS5-GwhVbbnZ2dnUVZ_sGvn...@comcast.com> ):
>
>
>
>
> 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, but
> 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 again to all for the useful replies. I had forgotten that
there is such a thing as a signed zero in Fortran. In my sample code
above, ksq has positive real part and (positive) zero imaginary part.
betasq is a REAL*8 quantity that is set to 0.0D0. Consider (betasq -
ksq), the argument to the SQRT call. Before the subtraction is
performed, betasq is promoted to COMPLEX*16 from REAL*8, presumably
receiving positive zero for its imaginary part. The imaginary part of
the difference is then of the form zero - zero, where both zeros have
positive signs. What is the sign attached to the difference by the
Fortran compiler? If it receives a negative sign, then I guess IVF is
returning the correct answer according to the f2003 standard. If this
is the case, then is it true that using /fp:precise is returning the
wrong answer?? Also, when I print out the value of the argument in my
sample code above, why doesn't the sign appear in the printout?
--Peter
| |
|
| On Aug 21, 12:04 pm, glen herrmannsfeldt <g...@ugcs.caltech.edu>
wrote:
> Peter wrote:
>
> 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
Yes, I am sure it is zero, because it is the result of multiplying,
adding, and subtracting COMPLEX*16 variables that have been previously
set equal to pure real values.
--Peter
| |
| James Van Buskirk 2007-08-21, 7:08 pm |
| "Peter" <petersamsimon2@hotmail.com> wrote in message
news:1187717724.205978.255920@o80g2000hse.googlegroups.com...
> 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!
I think you can check for negative zero via:
logical is_negative_zero
is_negative_zero = aimag(x) == 0 .AND. sign(real(1,kind(x)),aimag(x)) < 0
This could inform you as to whether your problem originates here in
any case.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
| |
| glen herrmannsfeldt 2007-08-22, 4:19 am |
| James Giles wrote:
(snip)
> Not clear why. All languages other than Fortran
> unambiguously require negative zero to print with
> a minus sign. (So does Fortran, but some people
> insist that would lead to a contradiction with some
> rarely used other feature.) Certainly the C language
> definition does. Since most I/O libraries are written
> in C these days, one would suspect that most I/O
> libraries *would* correctly print a minus sign when
> negative zero is output.
As far as I know, C89 doesn't require it and the compilers
I tried didn't do it.
I don't know C99 well enough to say.
-- glen
|
|
|
|
|