Fast round up/down double in fortran?

Question

Is there fast way to do round up/down in Fortran?

Because of linear order of bit-representation of positive double numbers it's possible to implement rounding as below.

pinf and ninf are global constants which are +/- infinity respectively

function roundup(x)
    double precision ,intent(in) :: x
    double precision :: roundup

    if (isnan(x))then 
        roundup = pinf
        return
    end if 
    if (x==pinf)then
        roundup = pinf
        return 
    end if
    if (x==ninf)then
        roundup = ninf
        return 
    end if
    if (x>0)then
        roundup = transfer((transfer(x,1_8)+1_8),1d0)
    else if (x<0) then
        roundup = transfer((transfer(x,1_8)-1_8),1d0)
    else
        if (transfer(x,1_8)==Z'0000000000000000')then
            roundup = transfer((transfer(x,1_8)+1_8),1d0)
        else
            roundup = transfer((transfer(-x,1_8)+1_8),1d0)
        end if 
    end if 
end function roundup

I feel it's not the best way to do that because it's slow, but it uses almost only bit-operations.

Another way is using multiplication and some epsilon eps = epsilon (1d0)

function roundup2(x)
    double precision ,intent(in) :: x
    double precision :: roundup2
    if (isnan(x)) then 
        roundup2 = pinf
        return 
    else if (x>=eps) then 
        roundup2 = x*(1d0+eps)
    else if (x<=-eps) then 
        roundup2 = x*(1d0-eps)
    else 
        roundup2 = eps
    end if 
end function roundup2

For some x both functions returns the same result (1d0, 158d0), for some don't (0.1d0, 15d0).

The first function is more accurate, but it's about 3.6 times slower than second (11.1 vs 3.0 seconds on 10^9 rounds test)

    print * ,x,y,abs(x-y)
    do i = 1, 1000000000
        x = roundup(x)
        !y = roundup2(y)
    end do 
    print * ,x,y,abs(x-y)

With no checks for NaN/Infinities first function test takes 8.5 seconds (-20%).

I use round function really hard and it takes a lot of time in profile of program. Is there cross-platform way to round faster with no loose of precision?

Update

The question suspects calls of roundup and rounddown at the time with no ability to reorder them. I didn't mention rounddown to keep topic short.

Hint: First function uses two transfer function and one adding. And it's slower than one multiplication and one adding in the second case. Why transfer cost so much when it doesn't do any with the number's bits? Is it possible to replace transfer by faster function(s) or avoid addition calls at all?

Answer 1

If I'm understanding correctly what you want to do, doesn't the "nearest" intrinsic do what you want, if you feed it +/- infinity as the arguments?

http://gcc.gnu.org/onlinedocs/gfortran/NEAREST.html#NEAREST

This might work, if the compiler implements this with decent performance. If you want NaN to round to Inf, you'll have to add that in a wrapper.

As for why roundup2 is faster, I can't tell for certain what's going on on your machine, but I can say two things:

1. The addition in roundup2 is probably optimized out (if eps is a parameter?) , so there's really just a multiplication.
2. If the transfer really does anything at all, that could easily slow the function down noticeably, since the function itself is so short. That might even be true if the transfer is just making superfluous copies of x.

Answer 2

I would recommend that you look at the Fortran standard IEEE floating point intrinsic modules (IEEE_ARITHMETIC, IEEE_FEATURES, IEEE_EXCEPTIONS). These provide IEEE_SET_ROUNDING_MODE where you can set the rounding mode for subsequent operations. Ideally you'd use IEEE_GET_ROUNDING_MODE to get the current mode and save it, set the new one, do your operations, then restore the mode.

Some caveats - changing the processor rounding mode is itself a slow operation, but if you do it once and then do lots of rounds, that will be a win. Not all current Fortran compilers support the IEEE intrinsic modules, but most reasonable ones should. You might need to tell the compiler you are playing with the IEEE environment - for Intel Fortran, use "-fp-model strict".