9

In Fortran 90 (using gfortran on Mac OS X) if I assign a value to a double-precision variable without explicitly tacking on a kind, the precision doesn't "take." What I mean is, if I run the following program:

program sample_dp

implicit none

integer, parameter :: sp = kind(1.0)
integer, parameter :: dp = kind(1.0d0)

real(sp) :: a = 0.
real(dp) :: b = 0., c = 0., d = 0.0_dp, e = 0_dp

! assign values
a = 0.12345678901234567890
b = 0.12345678901234567890
c = DBLE(0.12345678901234567890)
d = 0.12345678901234567890_dp

write(*,101) a, b, c, d
101 format(1x, 'Single precision: ',  T27, F17.15, / &
           1x, 'Double precisison: ', T27, F17.15, / &
           1x, 'Double precision (DBLE): ', T27, F17.15, / &
           1x, 'Double precision (_dp): ',  T27, F17.15)

end program

I get the result:

Single precision:        0.123456791043282
Double precision:        0.123456791043282
Double precision (DBLE): 0.123456791043282
Double precision (_dp):  0.123456789012346

The single precision result starts rounding off at the 8th decimal as expected, but only the double precision variable I assigned explicitly with _dp keeps all 16 digits of precision. This seems odd, as I would expect (I'm relatively new to Fortran) that a double precision variable would automatically be double-precision. Is there a better way to assign double precision variables, or do I have to explicitly type them as above?

tomshafer
  • 795
  • 3
  • 8
  • 12
  • 1
    Note that the intitialization of the variable `e` is done with an integer constant; `0_dp` is of type integer with kind number dp. It's an unfortunate reality that on many compilers the kind numbers for reals and integers overlap, so you won't get a compile time error for a mistake for this (I assume it was unintended). – eriktous May 26 '11 at 23:52

2 Answers2

11

A real which isn't marked as double precision will be assumed to be single precision. Just because sometime later you assign it to a double precision variable, or convert it to double precision, that doesn't mean that the value will 'magically' be double precision. It doesn't look ahead to see how the value will be used.

MRAB
  • 20,356
  • 6
  • 40
  • 33
  • 1
    So you mean the fact that I'm assigning the value to a double precision "container" doesn't matter unless I explicitly say it's a double precision value? – tomshafer May 26 '11 at 23:11
  • 2
    Correct. As I said, it doesn't look ahead to see how the value will be used later on. It sees a real literal, it's not marked as double precision, therefore it assumes it's single precision. – MRAB May 26 '11 at 23:45
  • 7
    @tomshafer: it's important to get a clear mental picture of the way an assignment statement is evaluated. First, the RHS is evaluated fully, according to the rules of precedence; at each step, if it's a mixed-mode expression, the required implicit type conversion is performed. Only after the RHS is evaluated fully, the result is assigned to the LHS, again with implicit type conversion, if required. – eriktous May 27 '11 at 00:01
  • Thanks guys, that makes much more sense. – tomshafer May 27 '11 at 00:03
10

There are several questions linking here so it is good to state some details more explicitly with examples, especially for beginners.

As stated by MRAB in his correct answer, an expression is always evaluated without any context, so

 0.12345678901234567890

is a default (single) precision floating literal, no matter where does it appear. The same holds to floating point numbers in the exponential form

 0.12345678901234567890E0

it is also a default precision number.

If one want to use a double precision constant, one can use D instead of E in the above form. Even if such a double precision constant is assigned to a default precision variable, it is first treated as a double precision number and then it is converted to default precision.

The way you are using in your question (employing the kind notation and several kind constants) is more general and more modern, but the principle is the same.

 0.12345678901234567890_sp

is a number of kind sp and 

 0.12345678901234567890_dp

is a number of kind dp and it does not matter where do they appear.

As your example shows, it is not only about assignment. In the line

 c = DBLE(0.12345678901234567890)

first the number 0.12345678901234567890 is default precision. Then it is converted to double precision by DBLE, but that is done after some of the digits are already lost. Then this new double precision number is assigned to c.

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
  • 1
    This question is a duplicate target for many similar questions and this additional explanation is useful. – Ross Mar 08 '17 at 18:43