Ceiling function Overflow issue

1

Dear all,

I am translating some Matlab version of matrix exponential solver to Fortran, one of it is the expmv from Higham,

Those solvers contain some large integer operations.
For a minimal example, for some matrix, in the solver it will calculate things like

ceiling(3.06921040378475e+21)

But the result is just wrong since its overflow. I don’t think brute force adding a kind parameter such as int32,64,128 in the ceiling function can really solve the problem, because overflow can still happen.

However, Matlab’s ceil function does not have problem. Its output is just 3.06921040378475e+21.
I know its ceil function seems like a double precision version, it returns the nearest integer in the format of double precision, so no overflow.

I wonder, is there a way to let the ceiling function in Fortran behaves like Matlab?
I mean the ceiling function do not really return a strict integer, it just need a correct value in double precision.

Many thanks in advance!

PS.
I can think of a stupid way, roughly, just define a double precision version of ceiling function to perhaps mimic Matlab’s ceil function, in it, it will check the input variable x. If x is very big, then ceiling function output x. Otherwise ceiling function output ceiling(x), and then convert it to double precision. But should there be better solutions?

  module ceil_mod
  implicit none
  contains
  elemental real function ceil(x)
  implicit none
  !integer, parameter :: i4=selected_int_kind(9)
  integer, parameter :: i8=selected_int_kind(15)
  integer, parameter :: r8=selected_real_kind(15,9)
  real(r8), intent(in) :: x
  integer(i8), parameter :: ceil_limit = 2.0**62
  if (abs(x)<ceil_limit) then
    ceil = ceiling(x,kind=i8)
  else
    ceil = x
  endif
  return
  end function ceil
  end module
2

Perhaps the ANINT() function is a better match? Of course, it returns the NEAREST integer value (as a floaitng-point number), but that could be amended with anint(x+0.5).

1 Like
3

The nearest integer to a float may not be representable for a real variable of a specified kind if the magnitude of the number is too large. For example, with gfortran

use iso_fortran_env, only: real128
print*,spacing(3.06921040378475e+21)
print*,spacing(3.06921040378475d+21)
print*,spacing(3.06921040378475e+21_real128)
end

gives

   2.81474977E+14
   524288.00000000000     
   4.54747350886464118957519531250000000E-0013

so you need quad precision to represent the nearest integer to your float.

2 Likes
4

Thank you very @Arjen @Beliavsky. The anint(x+0.5) is a great solution!

There is just one small annoying issue, and this can cause a big difference in the calculations in those solvers.
Thing is, for an exact integer N
ceiling(N) gives me exactly N. For example N=0, so ceiling(0) gives me 0 which is the expected value.
However, anint(N+0.5) gives N+1. For example N=0, so anint(0.5) gives me 1 which is not the expected value.

Should we add something slightly smaller than 0.5? If so, what should be that value that is slightly smaller than 0.5? :laughing:

5

Hm, yes, that edge case is annoying. How about 0.5 - epsilon(x) ?

1 Like
6

That works, fantastic! Thank you so much @Arjen !
I previously tried things like 0.5-tiny(1.0_r8), but the tiny does not work. The epsilon works!
So perhaps the below ceil function can mimic Matlab’s ceil function,

  module ceil_mod
  implicit none
  contains
  elemental real function ceil(x) ! To mimic Matlab's ceil function
  implicit none
  integer, parameter :: r8=selected_real_kind(15,9)
  real(r8), parameter :: almost_one_half = 0.5_r8-epsilon(1.0_r8)
  real(r8), intent(in) :: x
  ceil=anint(x+almost_one_half)
  return
  end function ceil  
  end module
1 Like
7

Thanks. I had to go back and verify a few procedures were never within spacing(x)/2.0 range of huge(0) after thinking about that; but they were OK.

I do not think I have seen a description of CEILING(3f) or FLOOR(3f) that points out that just not input values clearly outside of the range -huge(0)-1 to huge(0) are undefined but that near those boundaries equivalent real values are not representable in 32-bit and 64-bit values. It’s clear AFTER you think about it!

write(*,*)ceiling(huge(0)+0.0)
end

will set off a red flag from now on; kicking myself for not looking out for something more subtle but like that (that example is likely to get caught at compile time, at least with the right flags on).

1 Like