Simple program, that only produce zeros, bug? [dup

Question

This question already has an answer here:

Why are the elements of an array formatted as zeros when they are multiplied by 1/2 or 1/3? (1 answer)

Closed 2 years ago.

Why does this fortran program produce only zeros? When I print it out i get -0.00000 everywhere! What have I done wrong? In matlab it runs perfectly. I dont see any reason why its not working to be honest!

It seems like its the fraction that messes it up. if I set x equal to some decimal number it works.

program main


implicit none
  integer iMax, jMax
  double precision, dimension(:,:), allocatable :: T

double precision x, dx,f,L2old,L2norm,y

integer i, j,n,bc

n=10

 allocate(T(1:n+2, 1:n+2))

T=0.0d0

do i=2,n+1
 do j=2,n+1

 x=(j+1)*1/24
 y=(i+1)*1/24


 T(i,j)= -18*(x**2+y**2)**2

 Write(*,*)'T(',i,'',j,'', T(i,j)

end do
end do

Write(*,*)'T(1,1)',T(1,1)


end program main

Answer 1

x=(j+1)*1/24

1/24 is an integer division that rounds down to 0. You should be able to force floating point division by making at least one of the operands floating point, e.g.

x=(j+1)*1.0/24.0

Answer 2

As was indicated by Jim Lewis, the answer to the OP's question was indeed the integer division used.

Nonehteless, I think it is important to point out that one should take care of how the floating point fraction is written down. As the OP's program shows, x was of type DOUBLE PRECISION. Then the correct result should be

x=(j+1)*1.0D0/24.0D0

The difference here is that now you ensure that the division happens with the same precision as x was declared.

To following program demonstrates the problem ::

program test
WRITE(*,'(A43)')   "0.0416666666666666666666666666666666..."
WRITE(*,'(F40.34)') 1/24
WRITE(*,'(F40.34)') 1.0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0
WRITE(*,'(F40.34)') 1.0D0/24.0D0
end program test

which as the output

0.0416666666666666666666666666666666...
0.0000000000000000000000000000000000
0.0416666679084300994873046875000000
0.0416666666666666643537020320309239
0.0416666666666666643537020320309239

You clearly see the differences. The first line is the mathematical correct result. The second line is the integer division leading to zero. The third line, shows the output in case the division is computed as REAL while the fourth and fifth line are in DOUBLE PRECISION. Please take into account that in my case REAL implies a 32bit floating point number and DOUBLE PRECISION a 64 bit version. The precision and representation of both REAL and DOUBLE PRECISION is compiler dependent and not defined in the Standard. It only requires that DOUBLE PRECISION has a higher precision than REAL.

4.4.2.3 Real type

1 The real type has values that approximate the mathematical real numbers. The processor shall provide two or more approximation methods that define sets of values for data of type real. Each such method has a representation method and is characterized by a value for the kind type parameter KIND. The kind type parameter of an approximation method is returned by the intrinsic function KIND (13.7.89).

5 If the type keyword REAL is used without a kind type parameter, the real type with default real kind is specified and the kind value is KIND (0.0). The type specifier DOUBLE PRECISION specifies type real with double precision kind; the kind value is KIND (0.0D0). The decimal precision of the double precision real approximation method shall be greater than that of the default real method.

This actually implies that, if you want to ensure that your computations are done using 32bit, 64bit or 128bit floating point representations, you are advised to use the correct KIND values as defined in the intrinsic module ISO_FORTRAN_ENV.

13.8.2.21 REAL32, REAL64, and REAL128

1 The values of these default integer scalar named constants shall be those of the kind type parameters that specify a REAL type whose storage size expressed in bits is 32, 64, and 128 respectively. If, for any of these constants, the processor supports more than one kind of that size, it is processor dependent which kind value is provided. If the processor supports no kind of a particular size, that constant shall be equal to −2 if the processor supports kinds of a larger size and −1 otherwise.

So this would lead to the following code

PROGRAM main
  USE iso_fortran_env, ONLY : DP => REAL64
  IMPLICIT NONE
  ...
  REAL(DP) :: x
  ...
  x = (j+1)*1.0_DP/24.0_DP
  ...
END PROGRAM main