CODEWITHSUNDEEP
fortranprecisiondivisionfortran90
Gogrik
Gogrik

Reputation: 48

Division of double-precision by an integer in Fortran

Right now I am dealing with Fortran codes of the man who worked on the project before me. The code is relatively huge so I would not provide it here. In a lot of places in this codes there are division of double precision value over integer. I.e.:

double precision :: a,c
integer :: b
a = 1.0d0
b = 3
c = a/b

This sometimes leads to a mistakes much worse than a regular floating point number mistake. I will copypaste one of such mistakes for you: 3.00000032899483. For his purposes this was OK, but for me this mistake is awful. It seems like fortran is dividing real number over integer and then enlarge number to double precision, adding some random stuff. It can be treated of course if I will write:

c = a/(real(b,8))

So, I can do it, but it will take very long time to search through all his codes. That's why, I tried to overload (/) operator in the next way:

MODULE divisionoverload
!-----------------------------------------------------
   INTERFACE OPERATOR (/)
      MODULE PROCEDURE real_divoverinteger
   END INTERFACE OPERATOR ( / )
CONTAINS
!-----------------------------------------------------

DOUBLE PRECISION FUNCTION real_divoverinteger(a,b)
   IMPLICIT NONE
   DOUBLE PRECISION, INTENT (IN) :: a 
   INTEGER, INTENT(IN) :: b 

   real_divoverinteger = a/real(b,8)
END FUNCTION real_divoverinteger

END MODULE divisionoverload

But the mistake gfortran4.8 gives me clearly indicates that this conflicts with intrinsic division operator:

 MODULE PROCEDURE real_divoverinteger
                                          1
Error: Operator interface at (1) conflicts with intrinsic interface

So what can I do to solve this problem without manually handling all the divisions by integer?

Code which generates such numbers with printing format: 

 hx = (gridx(Nx1+1)-gridx(Nx1))/modx)
 do i = 2,lenx
    sub%subgridx(i)=sub%subgridx(i-1) + hx
 end do
 !*WHERE*
 !hx,gridx[1d array], sub%subgridx[1d array] - double precision
 !modx, lenx - integer

code for printing

subroutine PrintGrid(grid,N)
    integer, intent (in) :: N       
    double precision, dimension(N), intent (in) :: grid
    integer :: i
    print"(15(f20.14))", (grid(i), i=1,N)
        print*, ""
end subroutine PrintGrid

Solution: I apologize for this question. I should read code more careful. I found a mistake. The reproduction for the problem you get get in the next way:

program main
   double precision :: a
   a = 1.0/3
   print*, a
end program main

So basically, he was performing real/integer division and assigned this value to real*8. I changed everything with codes like this: 

 program main
   double precision :: a
   a = 1.0d0/3
   print*, a
end program main

And it worked! Thank you very much for help, without your answers I would be dealing with this problem much more.

Upvotes: 2

Views: 11386

Answers (2)

agentp
agentp

Reputation: 6989

as noted your title claim is incorrect, integer division is automatically upcast and is fine.

Your first example, however does have a problem:

program main
   double precision :: a
   a = 1.0/3
   print*, a
end program main

here we have a single precision 1.. The 3 gets promoted to single precision, so you have a single precision approximation to 1/3 assigned to your double precision a. The fix is to make the 1.0 1.d0.

Monkeying around overloading operators feels like a bad idea to me. If there is really to much code to manually fix see if your compiler has an option to change the default literal type to double.

Upvotes: 4

Gilles
Gilles

Reputation: 9489

I believe the initial assertion of your question, which is (if I understand you well) that dividing a real by an integer leads to lesser precision than dividing by a real, is wrong. Fortran promotes the weakest type of the two operands of an operation (if needed) to the strongest type, prior to compute the operation. Therefore, doing the promotion yourself is unnecessary.

See the following code for example:

program promotion
    implicit none
    integer, parameter :: dp = kind( 1.d0 )
    integer :: i
    real( kind = dp ) :: a

    i = 3
    a = 1
    print *, "division by a real:", a/real(i,dp), "division by an integer:", a/i
    print *, "are they equal?", a/real(i,dp) == a/i
end program promotion

The actual printing of the results may differ from one compiler to the next, but the final test will always evaluate to .true.

~/tmp$ gfortran promotion.f90
~/tmp$ ./a.out 
 division by a real:  0.33333333333333331      division by an integer:  0.33333333333333331     
 are they equal? T

Upvotes: 8

Related Questions