Severe Bug in g77 and probably gcc
Frank Kuehnel
kuehnel@pa.msu.edu
Tue Dec 7 19:59:00 GMT 1999
To whom it may concern:
Dear developer of free software, my appreciation is boundless for all
you
people who spend countless hours on the most honorable task to liberate
us
users from the scourge of commercial products. My contribution is, if
any,
only a very small one, and the more, it is only a complaint, a bug so to
speak:
The "bug" occurs in the Fortran program, which I included below. I've to
apologize that I was unable to shorten it, to the very core, but I do
hope
you can make some sense of my annotations! I use REAL numbers in the
program,
instead of DOUBLE PRECISION, and I find that seemingly innocent loops in
the main program produce different results, though they are not supposed
to!
This is not an ACCURACY problem with REAL number!
Here is the system and compiler I used, but it occured also with all
other
versions of gnu compilers, even "f2c" and then compiling reproduced the
bug!
Please read the text in the main program!
My contact address:
Frank Kuehnel
Michigan State University
Department of Physics and Astronomy
kuehnel@pa.msu.edu
g77 version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release) (from
FSF-g77 version 0.5.24-19981002)
Driving: g77 -v -c -xf77-version /dev/null -xnone
Reading specs from /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/specs
gcc version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/cpp -lang-c -v -undef
-D__GNUC__=2
-D__GNUC_MINOR__=91 -D__ELF__ -D__unix__ -D__i386__ -D__i386__
-D__linux__ -D__unix -D__i386 -D__linux
-Asystem(posix) -D_LANGUAGE_FORTRAN -traditional -Asystem(unix)
-Acpu(i386) -Amachine(i386) -Di386
-D__i386 -D__i386__ -D__tune_i386__ /dev/null /dev/null
GNU CPP version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release) (i386
Linux/ELF)
#include "..." search starts here:
#include <...> search starts here:
/usr/local/include
/usr/i386-redhat-linux/include
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include
/usr/include
End of search list.
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/f771 -fnull-version
-quiet -dumpbase g77-version.f -version
-fversion -o /tmp/ccC7Ii9g.s /dev/null
GNU F77 version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
(i386-redhat-linux) compiled by GNU C
version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release).
GNU Fortran Front End version 0.5.24-19981002
as -V -Qy -o /tmp/ccOYUklp.o /tmp/ccC7Ii9g.s
GNU assembler version 2.9-codefusion-990706 (i686-pc-linux-gnulibc2.1)
using BFD version 2.9-codefusion-990706
ld -m elf_i386 -dynamic-linker /lib/ld-linux.so.2 -o /tmp/cc26qZdx
/tmp/ccOYUklp.o /usr/lib/crt1.o
/usr/lib/crti.o
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/crtbegin.o
-L/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66
-L/usr/i386-redhat-linux/lib -lg2c -lm -lgcc -lc -lgcc
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/crtend.o /usr/lib/crtn.o
/tmp/cc26qZdx
__G77_LIBF77_VERSION__: 0.5.24
@(#)LIBF77 VERSION 19970919
__G77_LIBI77_VERSION__: 0.5.24-19981021
@(#) LIBI77 VERSION pjw,dmg-mods 19980617
__G77_LIBU77_VERSION__: 0.5.24-19990305
@(#) LIBU77 VERSION 19980709
Here is the program: Compiled with g77 -Wall -Wimplicit -pedantic
PROGRAM severe_flaw
c ... Task:
c ... given a test function "testf(x)" with a parameter omega
c ... find its minimum position "x0min"
c ... How it works:
c ... 1.) specify parameter omega
c ... 2.) call test(x0min) which returns the minimum position x0min
c ... Observation:
c ... The results depend on prehistory!!!!
c ... Loops with a different start produce different output??! ->
Severe Bug
c ... Cure:
c ... Exchanging DOUBLE PRECISION for REAL throughout the program!
c ...
REAL omega, x0min
COMMON /argument/ omega
INTEGER steps
c ... here we calculate the minimum for the function test with
c ... parameter omega = 0.6 ... 1.0
omega = 0.6
DO steps=1, 5
CALL test(x0min)
PRINT '( "testresult:",F8.2,2(1X,E14.6) )', omega,x0min
omega = omega + 0.1
END DO
c ...
c ... once again calculate the minimum for the function test
c ... with parameter omega 0.7 ... 1.0
c ... The values for 0.7 ... 1.0 should coinced with the previous
c ... calculation! IT DOES NOT! (Something is random here?)
c ...
omega = 0.7
DO steps=1, 4
CALL test(x0min)
PRINT '( "testresult:",F8.2,2(1X,E14.6) )', omega,x0min
omega = omega + 0.1
END DO
omega = 1.0
CALL test(x0min)
PRINT '( "testresult:",F8.2,2(1X,E14.6) )', omega,x0min
STOP
END
c ...
SUBROUTINE test(xmin)
REAL fret,xmin
c ...
REAL TOL
PARAMETER (TOL=1.e-4)
c ...
REAL ax,bx,fa,fb,fx,xx,brent,testf
EXTERNAL testf, brent
c ...
ax = 0.
xx = 1.
CALL mnbrak(ax,xx,bx,fa,fx,fb,testf)
fret=brent(ax,xx,bx,testf,TOL,xmin)
RETURN
END
c ...
FUNCTION testf(x)
REAL testf, x
c ...
REAL PI,LOG8
PARAMETER (PI = 3.141592653589, LOG8 = 2.079441542)
REAL omega
COMMON /argument/ omega
c ...
REAL dirlnG
c ...
dirlnG= 2.*x**2 - LOG((x*omega)**2)+LOG8
testf= PI*(1+EXP(-x**2))*omega**2 + dirlnG
RETURN
END
c ...
c ... the routines mnbrak and brent are from Numerical Recipes
chapter (10.1),
c ... and are found to work reliable!
c ...
SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func)
REAL ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY
EXTERNAL func
PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.e-20)
c ... Given a function func, and given distinct initial points ax
and bx,
c ... this routine searches in the downhill direction (defned by the
function
c ... as evaluated at the initial points) and returns new points ax,
bx, cx
c ... that bracket a minimum of the function. Also returned are the
function
c ... values at the three points, fa, fb, and fc. Parameters: GOLD
is the default
c ... ratio by which successive intervals are magnified; GLIMIT is
the maximum
c ... magnifcation allowed for a parabolic-fit step.
REAL dum,fu,q,r,u,ulim
fa=func(ax)
fb=func(bx)
IF (fb.GT.fa) THEN
c ... Switch roles of a and b so that we can go downhill in the
direction from a to b.
dum=ax
ax=bx
bx=dum
dum=fb
fb=fa
fa=dum
ENDIF
cx=bx+GOLD*(bx-ax)
c ... First guess for c.
fc=func(cx)
1 IF (fb.GE.fc) THEN
c ... "do while": keep returning here until we bracket.
r=(bx-ax)*(fb-fc)
c ... Compute u by parabolic extrapolation from a; b; c. TINY is
used to prevent
c ... any possible division by zero.
q=(bx-cx)*(fb-fa)
u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*SIGN(MAX(ABS(q-r),TINY),q-r))
ulim=bx+GLIMIT*(cx-bx)
c ... We won't go farther than this. Test various possibilities:
IF ((bx-u)*(u-cx).GT.0.) THEN
c ... Parabolic u is between b and c: try it.
fu=func(u)
IF (fu.LT.fc) THEN
c ... Got a minimum between b and c.
ax=bx
fa=fb
bx=u
fb=fu
RETURN
ELSE IF (fu.GT.fb) THEN
c ... Got a minimum between between a and u.
cx=u
fc=fu
RETURN
ENDIF
u=cx+GOLD*(cx-bx)
c ... Parabolic t was no use. Use default magni cation.
fu=func(u)
ELSE IF ((cx-u)*(u-ulim).GT.0.) THEN
c ... Parabolic t is between c and its allowed limit.
fu=func(u)
IF (fu.LT.fc) THEN
bx=cx
cx=u
u=cx+GOLD*(cx-bx)
fb=fc
fc=fu
fu=func(u)
ENDIF
ELSE IF ((u-ulim)*(ulim-cx).GE.0.) THEN
c ... Limit parabolic u to maximum allowed value.
u=ulim
fu=func(u)
ELSE
c ... Reject parabolic u, use default magnifcation.
u=cx+GOLD*(cx-bx)
fu=func(u)
ENDIF
ax=bx
c ... Eliminate oldest point and continue.
bx=cx
cx=u
fa=fb
fb=fc
fc=fu
GOTO 1
ENDIF
RETURN
END
c ...
FUNCTION brent(ax,bx,cx,f,tol,xmin)
INTEGER ITMAX
REAL brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS
EXTERNAL f
PARAMETER (ITMAX=50,CGOLD=.3819660,ZEPS=1.0e-10)
c ... Given a function f, and given a bracketing triplet of
abscissas ax, bx, cx
c ... (such that bx is between ax and cx, and f(bx) is less than
both f(ax) and f(cx)),
c ... this routine isolates the minimum to a fractional precision of
about tol using Brent's method.
c ... The abscissa of the minimum is returned as xmin, and the
minimum function value is
c ... returned as brent, the returned function value. Parameters:
Maximum allowed number of
c ... iterations;
c ... golden ratio; and a small number that protects against trying
to achieve fractional accuracy
c ... for a minimum that happens to be exactly zero.
INTEGER iter
REAL a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
a=MIN(ax,cx)
c a and b must be in ascending order, though the input abscissas
need not be.
b=MAX(ax,cx)
v=bx
c ... Initializations...
w=v
x=v
e=0.
c ,,, This will be the distance moved on the step before last.
fx=f(x)
fv=fx
fw=fx
DO iter=1, ITMAX
c ... Main program loop.
xm=0.5*(a+b)
tol1=tol*ABS(x)+ZEPS
tol2=2.*tol1
IF (ABS(x-xm).LE.(tol2-.5*(b-a))) GOTO 3
c ... Test for done here.
IF (ABS(e).GT.tol1) THEN
c ... Construct a trial parabolic-fit.
r=(x-w)*(fx-fv)
q=(x-v)*(fx-fw)
p=(x-v)*q-(x-w)*r
q=2.*(q-r)
IF (q.GT.0.) p=-p
q=ABS(q)
etemp=e
e=d
IF (ABS(p).GE.ABS(.5*q*etemp).OR.p.LE.q*(a-x).OR.
& p.GE.q*(b-x)) GOTO 1
c ... The above conditions determine the acceptability of the
parabolic fit. Here it is o.k.:
d=p/q
c ... Take the parabolic step.
u=x+d
IF (u-a.LT.tol2 .OR. b-u.LT.tol2) d=SIGN(tol1,xm-x)
GOTO 2
c ... Skip over the golden section step.
ENDIF
1 IF (x.GE.xm) THEN
c ... We arrive here for a golden section step, which we take into
the larger of the two segments.
e=a-x
ELSE
e=b-x
ENDIF
d=CGOLD * e
c ... Take the golden section step.
2 IF (ABS(d).GE.tol1) THEN
c ... Arrive here with d computed either from parabolic fit, or else
from golden section.
u=x+d
ELSE
u=x+SIGN(tol1,d)
ENDIF
fu=f(u)
c ... This is the one function evaluation per iteration,
IF (fu.LE.fx) THEN
c ... and now we have to decide what to do with our function
evaluation. Housekeeping follows:
IF (u.GE.x) THEN
a=x
ELSE
b=x
ENDIF
v=w
fv=fw
w=x
fw=fx
x=u
fx=fu
ELSE
IF (u.LT.x) THEN
a=u
ELSE
b=u
ENDIF
IF (fu.LE.fw .OR. w.EQ.x) THEN
v=w
fv=fw
w=u
fw=fu
ELSE IF (fu.LE.fv .OR. v.EQ.x .OR. v.EQ.w) THEN
v=u
fv=fu
ENDIF
ENDIF
c ... Done with housekeeping. Back for another iteration.
END DO
PAUSE 'brent exceed maximum iterations'
3 xmin=x
c ... Arrive here ready to exit with best values.
brent=fx
RETURN
END
More information about the Gcc-bugs
mailing list