Numerical Error and Ceiling Function

[frank_m1]

Say I have a variable x,

REAL(8) :: x
INTEGER :: y

x = 0.0_8
.
.
much numerical processing
.
.

y =CEILING(x)

y will frequently round up to NINT(x)+1 even when allvisible digits behind decimalare zero.

Is there a way to zero the decinal part of x when the decimal part is very close to zero?

[mecej4 · Accepted solution]

You wrote: y will frequently round up to NINT(x)+1 even when allvisible digits behind decimalare zero.

That is what CEILING is supposed to do, unless x is an exact integer.

What do you mean by "visible digits"? The value of a variable is not affected by whether or how many digits are displayed. If x = 2.00000000001, and I display x using, say, F10.8, I see 2.00000000, but CEILING(x) should still be 3, not 2.

Next, your wish to zero the "decimal part of x when that is very close to zero": what do you wish to do with 1.9999999 and 2.0000001? Will NINT(x-SPACING(x)) fill your needs?

[frank_m1]

Yeah, I was a bit vague. Here is the function I was working on.I thought about doing this as a linear programming problem but this method was recomended by number theorist. Additionally, I tried to speed it up with openMP and found that with even very large lengths of data, the current method is quickest.

FUNCTION FGET_BUFFER_SIZE(datalen)

IMPLICIT NONE

INTEGER :: FGET_BUFFER_SIZE
INTEGER, INTENT(IN) :: datalen
REAL(8) :: ln2, ln3, ln5, ln7, ln11
REAL(8) :: lndata

ln2 = DLOG(2.0_8)
ln3 = DLOG(3.0_8)
ln5 = DLOG(5.0_8)
ln7 = DLOG(7.0_8)
ln11 = DLOG(11.0_8)
lndata = DLOG(DBLE(datalen))

IF (MOD(datalen,2) >= 1) THEN
FGET_BUFFER_SIZE = ODD_BUFFER()
ELSE
FGET_BUFFER_SIZE = EVEN_BUFFER()
END IF

CONTAINS

FUNCTION ODD_BUFFER()
INTEGER :: ODD_BUFFER
INTEGER :: buffsize, buffest
REAL(8) :: j ,k, l, m, N
REAL(8) :: jtmp, ktmp, ltmp, mtmp

buffsize = HUGE(buffsize)
N = DLOG(DBLE(buffsize))
j = 0.0_8
k = 0.0_8
l = 0.0_8
m = 0.0_8
DO WHILE(m*ln11 mtmp=m*ln11
DO WHILE(l*ln7+mtmp < N)
ltmp = l*ln7+mtmp
DO WHILE(k*ln5+ltmp < N)
ktmp = k*ln5+ltmp
j = DNINT((lndata-ktmp)/ln3+0.45_8)
DO WHILE(j*ln3+ktmp < N)
jtmp = j*ln3+ktmp
buffest=NINT(DEXP(jtmp))
IF (buffest .GE. datalen) THEN
N = jtmp
buffsize = buffest
IF(buffsize == datalen) GOTO 200
ENDIF
j =j+1.0_8
ENDDO
j =0.0_8
k =k+1.0_8
ENDDO
k =0.0_8
l =l+1.0_8
ENDDO
l = 0.0_8
m = m+1.0_8
ENDDO
200 CONTINUE
ODD_BUFFER = buffsize
END FUNCTION ODD_BUFFER

FUNCTION EVEN_BUFFER()
INTEGER :: EVEN_BUFFER
INTEGER :: buffsize, buffest
REAL(8) :: i ,j , k, l, m, N
REAL(8) :: itmp, jtmp, ktmp, ltmp, mtmp

buffsize = HUGE(buffsize)
N = DLOG(DBLE(buffsize))
j = 0.0_8
k = 0.0_8
l = 0.0_8
m = 0.0_8
DO WHILE(ln2+m*ln11< N)
mtmp=ln2+m*ln11
DO WHILE(l*ln7+mtmp < N)
ltmp = l*ln7+mtmp
DO WHILE(k*ln5+ltmp < N)
ktmp = k*ln5+ltmp
DO WHILE(j*ln3+ktmp < N)
jtmp = j*ln3+ktmp-ln2
i = DNINT((lndata-jtmp)/ln2+0.45_8)
If(i < 1.0_8) i = 1.0_8
DO WHILE(i*ln2+jtmp < N)
itmp=i*ln2+jtmp
buffest=NINT(DEXP(itmp))
IF (buffest .GE. datalen) THEN
N = itmp
buffsize = buffest
IF(buffsize == datalen) GOTO 300
ENDIF
i =i+1.0_8
ENDDO
i =1.0_8
j =j+1.0_8
ENDDO
j =0.0_8
k =k+1.0_8
ENDDO
k =0.0_8
l =l+1.0_8
ENDDO
l = 0.0_8
m = m+1.0_8
ENDDO
300 CONTINUE
EVEN_BUFFER = buffsize
END FUNCTION EVEN_BUFFER

END FUNCTION FGET_BUFFER_SIZE

[jimdempseyatthecove]

The two contains functions return an unmodified "buffsize".
IOW the code above the 200 CONTINUE and 300 CONTINUE appear to be a busy NO-OP.

Also, check with your theorist. It looks like you are attempting a convergence function where you start at coarse grained convergence then step into finer and finer grains. If this is your intention then your loops are all wrong. You have nesting of loops coarse grainto fine grain. Whereas you would want a chain of loops coarse grain to fine grain

do coarsest
..
end do
do next to coarsest
...
end do
...
do finest
...
end do

IOW similar functionality to a binary search (or Newton-like method of convergence)

The nesting of the loops will result in the correct (same) result, however you will be perfroming unnecessary work.

Jim Dempsey

[mecej4]

Jim Dempsey wrote "The two contains functions return an unmodified "buffsize".
IOW the code above the 200 CONTINUE and 300 CONTINUE appear to be a busy NO-OP."

Not quite so, Jim. You perhaps overlooked his redefinition of "buffsize" in the innermost loops of the two functions.

Your observations as to the loops, etc., hit the nail on the head. The problem that the code is trying to solve appears to be something along the lines of:

'Find the smallest integer values of i,j,k,l such that f(i,j,k,l) >= x for a given value of x.'

Frank: Please post a verbal statement of the problem that you wish to solve along the lines of the preceding statement. Virtuosity in writing code is no match for mathematical skill and insight as far as obtaining an economical solution is concerned.