[3168] in Athena Bugs
RT Fortran compiler
daemon@ATHENA.MIT.EDU (James D. Oliver III)
Mon Sep 11 19:58:48 1989
To: bugs@ATHENA.MIT.EDU
Date: Mon, 11 Sep 89 19:58:30 EDT
From: James D. Oliver III <oliver@ATHENA.MIT.EDU>
The following subroutine does not compile on the RT, giving the following
message. The same occurs when compiled without optimization. The
program is successfully compiled on VAXstations.
-
f77 -O -c rtsafe.f
rtsafe.f:
rtsafe:
Termination code 11
*** Exit 3
-
*********************************************************************
C FUNCTION TO FIND ROOT USING COMBINATION OF NEWTON-RAPHSON AND
* BISECTION
double precision function rtsafe(funcd,x1,x2,xacc)
implicit double precision(a-h,o-z), integer(i-n)
external funcd
parameter(maxit=100)
call funcd(x1,fl,df)
call funcd(x2,fh,df)
if (fl*fh .ge. 0.d0) then
pause 'root in RTSAFE must be bracketed'
stop
endif
if (fl .lt. 0.d0) then
xl=x1
xh=x2
else
xh=x1
xl=x2
swap=fl
fl=fh
fh=swap
endif
call funcd(rtsf,f,df)
do 10 j=1,maxit
if (((rtsf-xh)*df-f)*((rtsf-xl)*df-f) .ge. 0.d0
$ .or. dabs(2.d0*f) .gt. dabs(dxold*df)) then
dxold=dx
dx=0.5d0*(xh-xl)
rtsf=xl+dx
if (xl .eq. rtsf) then
rtsafe=rtsf
return
endif
else
dxold=dx
dx=f/df
temp=rtsf
rtsf=rtsf-dx
if (temp .eq. rtsf) then
rtsafe=rtsf
return
endif
endif
if (dabs(dx) .lt. xacc) then
rtsafe=rtsf
return
endif
call funcd(rtsf,f,df)
if (f .lt. 0.d0) then
xl=rtsf
fl=f
else
xh=rtsf
fl=f
endif
10 continue
pause 'RTSAFE exceeding maximum iterations'
rtsafe=rtsf
return
end
