double precision numbers

[hannahse]

I have a program with several thousand lines of code. By default, all variables are double precision except for the usual I through N integers. I have noticed that equations like y = ( x - 1 ) / 2 amy give different results than y = ( x - 1.D0 ) / 2.D0. I assume the latter version is more accurate for double precision math.

The code is way too long to change all of the numbers by hand. Is there a compiler flag that would cause the numbers to be treated as double precision?

The program gives iterative solutions so I need as much consistency and precision as I can get.

[rasa]

I will leave details on what is correct and what's not to the experts.

I personally think that you should start using Kind definitions so that your code becomes portable. You can define kind as:

[fortran]module kinds
integer, parameter, public :: sp_k = kind(1.0) ! Default Type of Real
integer, parameter, public :: dp_k = selected_real_kind(2*precision(1.0_sp_k)) 
! The following statement will also work for Double Precision
! INTEGER, PARAMETER :: dp_k = SELECTED_REAL_KIND(0.0D0)  
! If one has trouble setting the sp_k on a 64-bit machine use the following
! kind definition that will be valid per IEEE arithmetic
! integer, parameter, public :: dp_k = selected_real_kind(p =  15, r = 300)

!  Once you have the Kinds defined, you can define the number constants as
real(kind = dp_k), parameter, public   :: one_quarter = 0.25_dp_k
real(kind = dp_k), parameter, public   :: one_half = 0.5_dp_k
real(kind = dp_k), parameter, public   :: one = 1.0_dp_k
real(kind = dp_k), parameter, public   :: two = 2.0_dp_k
real(kind = dp_k), parameter, public   :: three = 3.0_dp_k

end module kinds[/fortran]

In the code you can say

[fortran]use kinds
! the following is even better
! use kinds, only: dp_k, one, two, three[/fortran]

Then it is a matter of global search and replace. replace 1 by one and 2 by two.

This approach will make your code easier to port in the future.

integer, parameter, public :: sp_k = kind(1.0) ! Default Type of Real
integer, parameter, public :: dp_k = selected_real_kind(2*precision(1.0_sp_k))
! The following statement will also work for Double Precision
! INTEGER, PARAMETER :: dp_k = SELECTED_REAL_KIND(0.0D0)
! If one has trouble setting the sp_k on a 64-bit machine use the following
! kind definition that will be valid per IEEE arithmetic
! integer, parameter, public :: dp_k = selected_real_kind(p = 15, r = 300)

[hannahse]

I can't easily use global replace because
1. sometimes 1 is integer 1
2. sometimes I have 1 for double precision 1
3. sometimes I have 1. for double precision 1
4. sometimes I have 1.0 for double precision 1
5. I might hae a number 1.2325

where 1 is just an example to represent any digit

WhatI need is for any real number (not variable) to be treated as double precision (or kind(2))

[mecej4]

Ragu's advice is worth using in the future. However, since hannahse stated that all variables except those starting with I-N are double precision, it does not address the question raised.

The example given, y = (x - 1)/2, is not a valid example because 1 and 2 are exactly representable in double precision. My recommendation is to go through the code and change the constants as necessary.

That said, you may consider whether the /fpconstant and/or /realsize:8 options will serve your needs.

Note, however, that declaring

double precision, PI = 3.14159d0

does not give you a double precision value of PI, but PI = 2d0*acos(0d0), etc. will.

[DavidWhite]

My understanding is that in a mixed mode operation (X-1), the lower precision number is promoted to the higher level (1 -> 1D0). This means that the results should be the same, although extra time is required to do the conversion.

[TimP]

If x is default typed, there are cases where promotion to single precision won't give as much accuracy as double. An optimizing compiler would change the constants to the required type at compile time, so you wouldn't see a difference in performance. The promotion of /2 to *.5 with ifort occurs only when /Qprec-div is not set (a somewhat undesirable default), so writing in that conversion could be useful.

[nvaneck]

3., 4.,& 5. can be handled by /fpconstant to extend single precision constantsto docuble precision.

[hannahse]

Thanks to all for taking the time to respond. Based on documentation, I believe the fpconstant will be my best option. I will try it and and then compare results from a version where one of my subroutines has the
D0 attached to the numbers.

Although I have been changing numbers by hand as I make updates to the program, It is not practical for me to change all of the numbers by hand. And the more I change by hand, the more likely I am to introduce errors via mis-keying. It kind of falls into the category of things I would like to do but do not have time for ( like changing to implicit none, replacing commons with modules, getting rid of goto(s), etc).

THe specific example I quoted may or may not have different results. I know that some constants are represented exactly. The example was just meant to illustrate the concept. I have changed some of the constants by hand and I know that at least some of those produce different results. Since I'm never sure which do or do not make a difference, I have been trying to change them all one routine at a time.

[mecej4]

Your plan is well-thought out and sensible. I have found file-difference tools such as Ultracompare and Nelson Beebe's utility ndiff
very useful in similar work myself.

[jimdempseyatthecove]

This reply is targeted at Intel (Steve).

Is there, or in the lack thereof,what if there were, a -warn_mixed_double_and_single_precision_expression option (abbreviated of course), such that when enabled, would announce where it the code the occurrences were. Then using the IDE, you could quickly locate the problematic statements.

As a further (and subsequent) add-on, an IDEmacro could be provided that works in conjunction with the warning, such that for literals, it would auto-correct the expression (leaving old statement in annotated comment). e.g. click on warning in compiler message window (which positions to errant expression), then examine line to see if change required and if so, Alt-(macro hot key) to affect change.

This would be handy not only for old code conversonsbut for newly written code. This would make for quick work in making edits and would be (almost) goof proof.

Jim Dempsey

[Steven_L_Intel1]

There is no option to warn of mixed-mode arithmetic. Such a warning, optional, is on a list of requests. I know some other compilers do this.

The Fortran standard has explicit rules allowing and defining mixed-mode arithmetic. In most cases, it does what you want. In the case above, I see no reason it should cause different results. I'd like to see a small but complete program that shows otherwise.

[hannahse]

Here's a small program demonstrating that there is a difference with and without the D0. ( I don't seem to be smart enough to paste without messing up spacing even though I am pasting from a text editor (textpad)

Program test_mixed
implicit none
double precision :: x, y, Z, z1

x = 3.1415967 * 2.5
y = 3.1415967D0 * 2.5D0
z1 = (y - x) * 1.0D6
z = (y - x)
write(*,1) "X= ", X, "Y= ", Y, "Difference= ", z, "large diff= ", z1
1 format(3(2X,A10,2XF12.10/)2X,A10,2XF12.6)

end

I compile this with 8.1 compiler using "ifort /debug:full test2.f90" and get the following results when i run the executable:

X= 7.8539919853
Y= 7.8539917500
Difference -.0000002353
large diff -0.235321

[Steven_L_Intel1]

This is very different than your original question.

The expression 3.1415967 * 2.5 is single-precision. It is evaluated as single-precision and then converted to double by adding 32 binary zero fraction bits. The expression 3.1415967D0 * 2.5D0 is double-precision and has more useful fraction bits. The value 3.1415967 is not exactly representable in either single or double-precision, so the more bits of double gets you closer to the decimal value.

Since you are using 8.1, there is the additional complication of the compiler using the extended-precision X87 instructions by default, which can lead to inconsistent appearance of "extra" precision, though not in this case.

[mecej4]

One of the first things one learns about doing arithmetic on a computer is exactly what you just described.

Here are a few references:

1. float.htm

2. What Every Computer Scientist Should Know About Floating-Point Arithmetic

and the Chauvenet Prize paper of Wilkinson, "The perfidious polynomial", as quoted in

3. Higham's book.