#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# invmf
# This archive created: Mon Jan 23 10:51:04 1995
export PATH; PATH=/bin:$PATH
if test ! -d 'invmf'
then
mkdir 'invmf'
fi
cd 'invmf'
if test -f 'INSTALL'
then
echo shar: will not over-write existing file "'INSTALL'"
else
cat << \SHAR_EOF > 'INSTALL'
INVMF - Invert a Monotone Function
Version 1.0 - October 1994
INTRODUCTION
The invmf function finds the value of x such that f(x) = y.
IMPROVEMENTS over previous version, which was called zero.
(1) Arbitrary parameters can now be passed to the function of which
the inverse value is sought.
(2) The packaging has been changed slightly so that a search on an
interval can be made efficiently.
INSTALLATION
1. FORTRAN and C SOURCE
The files cdinvr.c, devlsf.c, dinvr.f and dzror.f contain the FORTRAN
and C source code is neccessary to use invmf.
The source files must be compiled into an object file. I have
included makefiles which may be helpful compiling the source code.
Check to see if the options in the makefile are appropriate to your
system. If you are using a SUN OS type
% make -f makefile.sun
If you are using any other kind of UNIX OS type
% make
Before using the invmf S functions, the object file (called invmf.o)
must be dynamically loaded into S. The method for dynamically loading
object files in S will depend on your system and version of S. On our
system we use the S command `dyn.load':
> dyn.load("invmf.o")
This should be done only once at the beginning of each S session which
will use the invmf function.
2. S CODE
The two files, invmf.s and invmf.f.define.s contain the S code for the
three integration functions. They can be loaded by using the S
`source' command.
For example,
> source("invmf.s")
> source("invmf.f.define.s")
3. HELPFILE
The file invmf.d contains the help file for the invmf function. This
file should be copied into the .Data/.Help directory without the `.d'
extension. This assumes your system reads help files using some form
of nroff.
LEGALITIES
The reference for the Fortran routine, dzror, which finds the zero of
a function is:
Algorithm R of the paper 'Two Efficient Algorithms with Guaranteed
Convergence for Finding a Zero of a Function' by J. C. P. Bus and
T. J. Dekker in ACM Transactions on Mathematical Software, Volume 1,
no. 4 page 330 (Dec. '75)
The original code was written in Algol; we transliterated to Fortran.
Code that appeared in an ACM publication is subject to their
algorithms policy:
Submittal of an algorithm for publication in one of the ACM
Transactions implies that unrestricted use of the algorithm within a
computer is permissible. General permission to copy and distribute
the algorithm without fee is granted provided that the copies are not
made or distributed for direct commercial advantage. The ACM
copyright notice and the title of the publication and its date appear,
and notice is given that copying is by permission of the Association
for Computing Machinery. To copy otherwise, or to republish, requires
a fee and/or specific permission.
Krogh, F. Algorithms Policy. ACM Tran. Math. Softw. 13(1987),
183-186.
We place the code that we have written in the public domain.
DISCLAIMER: Although care was taken in the creation and testing of
this system, any user employs it at his own risk. Neither the authors
nor the University of Texas assumes any legal responsibility for the
correct functioning of the system.
CONTACTING THE AUTHORS
Please send comments, etc. to
kruss@odin.mda.uth.tmc.edu.
OBTAINING INVMF
The file is available by anonymous ftp from odin.mda.uth.tmc.edu
(129.106.3.17) in
/pub/S/invmf.shar.Z
It will also be posted to statlib.
Barry W. Brown
Department of Biomathematics,
Box 237
University of Texas M. D.
Anderson Hospital
1515 Holcombe Blvd
Houston, TX 77030
bwb@odin.mda.uth.tmc.edu
internet address is (129.106.3.17)
SHAR_EOF
fi # end of overwriting check
if test -f 'invmf.d'
then
echo shar: will not over-write existing file "'invmf.d'"
else
cat << \SHAR_EOF > 'invmf.d'
.BG
.FN invmf
.TL
Invert a monotone function
.DN
Finds the value of x such that f(x) = fequals for a provided
function monotone in x, f, and value, fequals.
.CS
invmf(f, init, fequals, ..., small=-10^(300), big=10^300,
name="x", absstp=10^(-5), relstp=10^(-3), stpmul=2,
abstol=10^(-8), reltol=10^(-6))
.RA
.AG f
The S function to be inverted. f is a function of x and possibly
other arguments.
.AG init
Initial guess as to the value of x such that f(x) = fequals.
.AG fequals
Value of function for which x is sought.
.OA
.AG ...
Other arguments to f. These must be specified by name not by
position. For example, use 'a=10' instead of 10 for the argument 'a'
of f. The arguments provided here are not varied when f is evaluated.
.AG small
Low value of interval in which solutions are sought. f is evaluated
at small so the value should be within the range of definition of f.
.AG big
High value of interval in which solutions are sought. f is evaluated
at big so the value should be within the range of definition of f.
.AG name
Name of the variable on which f is to be inverted.
.AG absstp
Minimum initial step size. If absstp >= (big-small) then algorithm R
of the reference is invoked immediately on the interval small to big.
See DETAILS.
.AG relstp
Relative minimum initial step size. See DETAILS.
.AG stpmul
Multiplier of step size when current step does not bound x. See DETAILS.
.AG abstol
Absolute tolerance of answer. See DETAILS.
.AG reltol
Relative tolerance of answer. See DETAILS.
.RT
The value of x approximately satisfying f(x) = fequals if such a value
is found. If a value is not found, NA is returned; this
happens only when f(small) - fequals has the same sign as f(big) -
fequals.
.SH NOTES
invmf can not be called recursively. For example,
.PP
f <- function(x){return(invmf(sqrt,x,x,small=0))}
invmf(f,16,16,small=0,big=100)
.PP
will not give the desired result.
.PP
Unless you are familiar with the arcana of S argument matching, use the
full names of all arguments in '...' and in the optional arguments following '...'.
.DT
Step 0: Determine whether f is increasing or decreasing in x.
.PP
f is evaluated at small and big. If f(small)-fequals has the same
sign as f(big)-fequals return NA. If absstp is greater than or equal
to (big-small) go to Step 3. Determine whether f is increasing or
decreasing in x.
.PP
Step 1: Bound x.
.PP
The initial step size is calculated.
.ce
.br
step <- max( absstp, relstep*init).
Evaluate f at a distance step from init -- the direction is determined
by whether f increases or decreases with x. If the answer is between init
and this evaluation point go to Step 3.
.PP
Whenever x is not bounded by a step, the step size is multiplied by
stpmul and f is evaluated at the resulting new point. If the desired
x is bounded by the new point and the previous one, go to step 3.
.PP
Step 3. Refine the value of x.
.PP
Algorithm R of the REFERENCE is invoked to find x on an interval by a
combination of polynomial approximation and interval bisection. The
interval searched has enpoints small and big if this Step 3 was
arrived at from Step 0, otherwise the bounding values determined by
Step 1 are the interval endpoints.
.PP
To determine whenter a value of x is acceptable, define
.ce
.br
tol <- max( abstol. reltol*x ).
The value is acceptable if the current bounding interval is small
enough to guarantee that the answer lies between x-tol and x+tol.
.SH REFERENCES
Bus, J. C. P. and Dekker, T. J. (1975) 'Two Efficient Algorithms
with Guaranteed Convergence for Finding a Zero of a Function'
ACM Transactions on Mathematical Software: 1, 330-343.
.EX
invmf(sqrt,4,4,small=0,big=100) # solves sqrt(x) = 4
# Find inverse of f(p) = 0.8 with status = 0
f <- function(status,p){if (status == 0) return(p) else return(1-p)}
invmf(f,0.5,0.8,status=0,small=0,big=1,name="p")
.KW math
.KW optimize
.WR
SHAR_EOF
fi # end of overwriting check
if test -f 'invmf.f.define.s'
then
echo shar: will not over-write existing file "'invmf.f.define.s'"
else
cat << \SHAR_EOF > 'invmf.f.define.s'
invmf.f.define <- function(f,name,f.extra)
{
#
# Assign Variables in f.extra to frame 1
n.extra <- length(f.extra)
call.extra <- ""
if (n.extra > 0) {
names.extra <- names(f.extra)
for (i in 1:n.extra) {
temp <- paste("invmf.param.",i,sep="")
call.extra <- paste(call.extra,",",names.extra[i],"=",temp,sep="")
assign(temp,f.extra[[i]],frame=1)
}
}
#
# Assign function f to frame 1
call.func <- "invmf.f"
assign(call.func,f,frame=1)
#
# Find parameters of the function, f
f.param <- names(f)
f.param <- f.param[f.param != ""]
#
# Make sure that the integration variable name is a parameter of f
if (!any(f.param == name))
stop(paste(name,"is not a parameter of f"))
#
# Define a function call to f and assign it to frame 1
f.call <- paste(call.func,"(",name,"=x9999.z1111.q13",call.extra,sep="")
f.call <- paste(f.call,")",sep="")
assign("invmf.f.call",f.call,frame=1)
#
# Define a new function, f2, with a single variable
f2 <- function(x9999.z1111.q13) {
ans <- as.double(eval(parse(text=invmf.f.call)))
if ((!is.numeric(ans)) || any(is.na(ans)) || any(is.inf(ans)))
stop("function, `f' must return a numeric value")
return(ans)
}
return(f2)
}
SHAR_EOF
fi # end of overwriting check
if test -f 'invmf.s'
then
echo shar: will not over-write existing file "'invmf.s'"
else
cat << \SHAR_EOF > 'invmf.s'
invmf <- function(f,init,fequals,...,small=-10^(300),big=10^300,
name="x",absstp=10^(-5),relstp=10^(-3),stpmul=2,
abstol=10^(-8),reltol=10^(-6))
{
#####################################################################
#
# invmf
#
# Function:
# Find X such that F(X)=FEQUALS for S function F.
#
# Arguments:
#
# f - Name of S function which is being inverted in
# terms of the variable defined in name.
#
# init - Initial guess to solution.
#
# fequals - Value to which f is set.
#
# small - Low value in interval of possible solutions.
#
# big - Hi value in interval of possible solutions.
#
# name - a character string containing the name of the variable
# to be solved for. (default is "x")
#
# absstp - Absolute step size. (default is 10^-5)
#
# relstp - Relative step size. (default is 10^-3)
#
# stpmul - Stop mulitplier. (default is 2)
#
# abstol - Absolute tolerance. (default is 10^-8)
#
# reltol - Relative tolerance. (default is 10^-6)
#
# ... : Any additional arguments to the function, f.
#
#####################################################################
#
# check type of arguments input
#
if( missing(f) || !is.function(f) )
stop("Arguement `f' must be a function")
if( missing(init) || !is.numeric(init) )
stop("Arguement `init' must be numeric (initial guess)")
if( missing(fequals) || !is.numeric(fequals) )
stop("Arguement `fequals' must be numeric (F(x)=FEQUALS)")
if( !is.numeric(small) )
stop("Arguement `small' must be numeric (low bnd for solution interval)")
if( !is.numeric(big) )
stop("Arguement `big' must be numeric (hi bnd for solution interval)")
if(!is.numeric(absstp) )
stop("Argument `absstp' must be numeric")
if(!is.numeric(relstp) )
stop("Argument `relstp' must be numeric")
if(!is.numeric(abstol) )
stop("Argument `abstol' must be numeric")
if(!is.numeric(reltol) )
stop("Argument `reltol' must be numeric")
if((!is.numeric(stpmul)) || (stpmul <= 1))
stop("Argument `stpmul' must be greater than 1")
if ( !is.character(name))
stop ("Argument `name' should be a character string")
#
# check absstp and relstp for legality
#
if( ( absstp < 0 ) && ( relstp < 1.0e-12 ) )
stop("\Bound from absstp and relstp too small")
#
# check to find out if small < init < big
#
if( init < small || init > big )
stop("\Initial solution intervals do not include initial guess")
#
# Create, f2, a function of the single variable
#
f.extra <- list(...)
f2 <- invmf.f.define(f,name,f.extra)
#
# check to find out if function is non-increasing or no-decreasing
#
fs <- f2(small)
fb <- f2(big)
fx <- f2(init)
if( is.na(fs) || is.na(fb) || is.na(fx) )
stop("\nBad values were entered for init,small, or big")
if( ((fs < fx) && (fb < fx)) || ((fs > fx) && (fb > fx)) )
stop("\nThis function is not a monotone function between small and big")
#
# invoke c driver
#
rzero.find <- .C("cdinvr",
list(f2),
x=as.double(init),
as.double(fequals),
status=as.integer(0),
qleft=as.integer(0),
qhi=as.integer(0),
as.double(small),
as.double(big),
as.double(absstp),
as.double(relstp),
as.double(stpmul),
as.double(abstol),
as.double(reltol))
#
# check exit conditions to see if all was ok
#
if (rzero.find$status) {
if( rzero.find$qleft == 1 )
warning("\nstepping search terminated at small")
else
warning("\nstepping search terminated at big")
if( rzero.find$qhi == 1 )
warning("\nF(X) is greater than FEQUALS at stopping point")
else
warning("\nF(X) is less than FEQUALS at stopping point")
stop("\nInput conditions were not met with small and big") }
#
# return answer
#
return(rzero.find$x) }
SHAR_EOF
fi # end of overwriting check
if test -f 'makefile'
then
echo shar: will not over-write existing file "'makefile'"
else
cat << \SHAR_EOF > 'makefile'
# Makefile for invmf source
# FC is name of Fortran Compiler
FC = f77
# FCFLAGS are flags to be set for the compile
# Note the -c flag is hard wired in
FCFLAGS =
# CC is name of C Compiler
CC = cc
# CCFLAGS are flags to be set for the compile
# Note the -c flag is hard wired in
CCFLAGS =
# Libraries needed in the load process go here.
# There may be none, or perhaps Fortran or C Libraries
# This is very system specific
LIBS =
#----------------------------------------------------------------------
# The user probably doesn't want to change anything below here
TARGET = invmf.o
OBJ = \
cdinvr.o \
devlsf.o \
dinvr.o \
dzror.o
# Set up suffixes for created implicit rules
.SUFFIXES:
.SUFFIXES: .f .c .o $(SUFFIXES)
# Rule for .f to .o
.f.o:
$(FC) $(FCFLAGS) -c $*.f
# Rule for .c to .o
.c.o:
$(CC) $(CCFLAGS) -c $*.c
# Create relocatable load file
all:$(TARGET)
$(TARGET):$(OBJ)
ld -r -d -o $(TARGET) $(OBJ) $(LIBS)
clean:
rm $(OBJ)
SHAR_EOF
fi # end of overwriting check
if test -f 'dinvr.f'
then
echo shar: will not over-write existing file "'dinvr.f'"
else
cat << \SHAR_EOF > 'dinvr.f'
SUBROUTINE dinvr(func,x,y,status,qleft,qhi,small,big,absstp,
+ relstp,stpmul,abstol,reltol)
C**********************************************************************
C
C SUBROUTINE DINVR(FUNC, X, Y, STATUS, QLEFT, QHI, SMALL, BIG,
C ABSSTP, RELSTP, STPMUL, ABSTOL, RELTOL)
C Double precision
C bounds the zero of the function and invokes zror
C
C
C Function
C
C
C Concise Description - Given a monotone function F finds X
C such that F(X) = Y.
C
C More Precise Description of INVR -
C
C F must be a monotone function, the results are
C otherwise undefined. Let QINCR be .TRUE. if F is non-
C decreasing and .FALSE. if F is non-increasing.
C
C F(SMALL) and F(BIG) must bracket Y, i. e.,
C QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
C QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
C
C The X returned satisfies the following condition. let
C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
C then if QINCR is .TRUE.,
C F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
C and if QINCR is .FALSE.
C F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
C
C Method: Bounds the function and invokes ZROR to perform the zero
C finding.
C
C
C Arguments
C
C
C FUNC --> External function which is being inverted.
C Should be a double precision function of
C one double precision variable.
C DOUBLE PRECISION FUNC
C EXTERNAL FUNC
C
C X <-- The value of X at which F(X) is to be evaluated.
C DOUBLE PRECISION X
C
C FX --> F(X) = Y
C DOUBLE PRECISION Y
C
C STATUS <-- When INVR has finished without error, it will return
C with STATUS 0. In that case X is approximately a root
C of F(X).
C
C If INVR cannot bound the function, it returns status
C 1 and sets QLEFT and QHI.
C INTEGER STATUS
C
C QLEFT <-- Defined only if STATUS is 1. In that
C case it is .TRUE. If the stepping search terminated
C unsucessfully at SMALL. If it is .FALSE. the search
C terminated unsucessfully at BIG.
C QLEFT is LOGICAL
C
C QHI <-- Defined only if STATUS is 1. In that
C case it is .TRUE. if F(X) .GT. Y at the termination
C of the search and .FALSE. if F(X) .LT. Y at the
C termination of the search.
C QHI is LOGICAL
C
C SMALL --> The left endpoint of the interval to be
C searched for a solution.
C SMALL is DOUBLE PRECISION
C
C BIG --> The right endpoint of the interval to be
C searched for a solution.
C BIG is DOUBLE PRECISION
C
C ABSSTP, RELSTP --> The initial step size in the search
C is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
C ABSSTP is DOUBLE PRECISION
C RELSTP is DOUBLE PRECISION
C
C STPMUL --> When a step doesn't bound the zero, the step
C size is multiplied by STPMUL and another step
C taken. A popular value is 2.0
C DOUBLE PRECISION STPMUL
C
C ABSTOL, RELTOL --> Two numbers that determine the accuracy
C of the solution. See function for a precise definition.
C ABSTOL is DOUBLE PRECISION
C RELTOL is DOUBLE PRECISION
C
C**********************************************************************
C .. Scalar Arguments ..
DOUBLE PRECISION y,x,absstp,relstp,stpmul,abstol,reltol,small,big,
+ ans
INTEGER status
LOGICAL qhi,qleft
C ..
C .. Local Scalars ..
DOUBLE PRECISION fx,fbig,fsmall,step,xhi,xlb,xlo,xsave,xub,yy
C + ,ZX,ZY,ZZ
LOGICAL qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup
C ..
C .. External Subroutines ..
EXTERNAL func,dzror
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,max,min
C ..
C .. Statement Functions ..
C LOGICAL QXMON
C ..
C .. Save statement ..
SAVE
C ..
C .. Statement Function definitions ..
C QXMON(ZX,ZY,ZZ) = ZX .LE. ZY .AND. ZY .LE. ZZ
C ..
C .. Executable Statements ..
C
C Initalize Values
C
qleft = .FALSE.
qhi = .FALSE.
status = 0
C
C Check to see that (SMALL <= X) && (X <= BIG)
C
C QCOND = .NOT.QXMON(SMALL,X,BIG)
C IF (QCOND) STOP ' SMALL, X, BIG not monotone in INVR'
xsave = x
C
C See that SMALL and BIG bound the zero and set QINCR
C
x = small
C GET-FUNCTION-VALUE
ASSIGN 10 TO i99999
GO TO 280
10 fsmall = fx
x = big
C GET-FUNCTION-VALUE
ASSIGN 20 TO i99999
GO TO 280
20 fbig = fx
qincr = fbig .GT. fsmall
IF (.NOT. (qincr)) GO TO 50
IF (fsmall.LE.0.0D0) GO TO 30
status = 1
qleft = .TRUE.
qhi = .TRUE.
RETURN
30 IF (fbig.GE.0.0D0) GO TO 40
status = 1
qleft = .FALSE.
qhi = .FALSE.
RETURN
40 GO TO 80
50 IF (fsmall.GE.0.0D0) GO TO 60
status = 1
qleft = .TRUE.
qhi = .FALSE.
RETURN
60 IF (fbig.LE.0.0D0) GO TO 70
status = 1
qleft = .FALSE.
qhi = .TRUE.
RETURN
70 CONTINUE
80 x = xsave
step = max(absstp,relstp*abs(x))
C GET-FUNCTION-VALUE
ASSIGN 90 TO i99999
GO TO 280
90 yy = fx
IF (.NOT. (yy.EQ.0.0D0)) GO TO 100
status = 0
qok = .TRUE.
RETURN
100 qup = (qincr .AND. (yy.LT.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.GT.0.0D0))
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C HANDLE CASE IN WHICH ABSSTP >= (BIG-SMALL)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
IF (.NOT. (absstp.GE.abs(big-small))) GO TO 120
x = xsave
xlb = small
xub = big
C CALL-ZROR
ASSIGN 110 TO i99987
GO TO 290
110 RETURN
120 IF (.NOT. (qup)) GO TO 190
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C HANDLE CASE IN WHICH WE MUST STEP HIGHER
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xlb = xsave
xub = min(xlb+step,big)
GO TO 140
130 IF (qcond) GO TO 170
140 x = xub
C GET-FUNCTION-VALUE
ASSIGN 150 TO i99999
GO TO 280
150 yy = fx
qbdd = (qincr .AND. (yy.GE.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.LE.0.0D0))
qlim = xub .GE. big
qcond = qbdd .OR. qlim
IF (qcond) GO TO 160
step = stpmul*step
xlb = xub
xub = min(xlb+step,big)
160 GO TO 130
170 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 180
status = 1
qleft = .FALSE.
qhi = .NOT. qincr
x = big
RETURN
180 GO TO 260
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C HANDLE CASE IN WHICH WE MUST STEP LOWER
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
190 xub = xsave
xlb = max(xub-step,small)
GO TO 210
200 IF (qcond) GO TO 240
210 x = xlb
C GET-FUNCTION-VALUE
ASSIGN 220 TO i99999
GO TO 280
220 yy = fx
qbdd = (qincr .AND. (yy.LE.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.GE.0.0D0))
qlim = xlb .LE. small
qcond = qbdd .OR. qlim
IF (qcond) GO TO 230
step = stpmul*step
xub = xlb
xlb = max(xub-step,small)
230 GO TO 200
240 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 250
status = 1
qleft = .TRUE.
qhi = qincr
x = small
RETURN
250 CONTINUE
C CALL-ZROR
260 ASSIGN 270 TO i99987
GO TO 290
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
270 RETURN
C ..
C .. Flecs In-Line Procedure
C TO GET-FUNCTION-VALUE
280 CALL devlsf(func,x,ans)
fx = ans - y
GO TO i99999
C TO CALL-ZROR
290 CALL dstzr(xlb,xub,abstol,reltol)
status = 0
GO TO 310
300 IF (.NOT. (status.EQ.1)) GO TO 340
310 CALL dzror(status,x,fx,xlo,xhi,qdum1,qdum2)
IF (.NOT. (status.EQ.1)) GO TO 330
C GET-FUNCTION-VALUE
ASSIGN 320 TO i99999
GO TO 280
320 CONTINUE
330 GO TO 300
340 x = xlo
status = 0
GO TO i99987
END
SHAR_EOF
fi # end of overwriting check
if test -f 'dzror.f'
then
echo shar: will not over-write existing file "'dzror.f'"
else
cat << \SHAR_EOF > 'dzror.f'
SUBROUTINE dzror(status,x,fx,xlo,xhi,qleft,qhi)
C**********************************************************************
C
C SUBROUTINE DZROR(STATUS, X, FX, XLO, XHI, QLEFT, QHI)
C Double precision ZeRo of a function -- Reverse Communication
C
C
C Function
C
C
C Performs the zero finding. STZROR must have been called before
C this routine in order to set its parameters.
C
C
C Arguments
C
C
C STATUS <--> At the beginning of a zero finding problem, STATUS
C should be set to 0 and ZROR invoked. (The value
C of other parameters will be ignored on this call.)
C
C When ZROR needs the function evaluated, it will set
C STATUS to 1 and return. The value of the function
C should be set in FX and ZROR again called without
C changing any of its other parameters.
C
C When ZROR has finished without error, it will return
C with STATUS 0. In that case (XLO,XHI) bound the answe
C
C If ZROR finds an error (which implies that F(XLO)-Y an
C F(XHI)-Y have the same sign, it returns STATUS -1. In
C this case, XLO and XHI are undefined.
C INTEGER STATUS
C
C X <-- The value of X at which F(X) is to be evaluated.
C DOUBLE PRECISION X
C
C FX --> The value of F(X) calculated when ZROR returns with
C STATUS = 1.
C DOUBLE PRECISION FX
C
C XLO <-- When ZROR returns with STATUS = 0, XLO bounds the
C inverval in X containing the solution below.
C DOUBLE PRECISION XLO
C
C XHI <-- When ZROR returns with STATUS = 0, XHI bounds the
C inverval in X containing the solution above.
C DOUBLE PRECISION XHI
C
C QLEFT <-- .TRUE. if the stepping search terminated unsucessfully
C at XLO. If it is .FALSE. the search terminated
C unsucessfully at XHI.
C QLEFT is LOGICAL
C
C QHI <-- .TRUE. if F(X) .GT. Y at the termination of the
C search and .FALSE. if F(X) .LT. Y at the
C termination of the search.
C QHI is LOGICAL
C
C**********************************************************************
C .. Scalar Arguments ..
DOUBLE PRECISION fx,x,xhi,xlo,zabstl,zreltl,zxhi,zxlo
INTEGER status
LOGICAL qhi,qleft
C ..
C .. Save statement ..
SAVE
C ..
C .. Local Scalars ..
DOUBLE PRECISION a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q,
+ reltol,tol,w,xxhi,xxlo,zx
INTEGER ext,i99999
LOGICAL first,qrzero
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,max,sign
C ..
C .. Statement Functions ..
DOUBLE PRECISION ftol
C ..
C .. Statement Function definitions ..
ftol(zx) = 0.5D0*max(abstol,reltol*abs(zx))
C ..
C .. Executable Statements ..
IF (status.GT.0) GO TO 280
xlo = xxlo
xhi = xxhi
b = xlo
x = xlo
C GET-FUNCTION-VALUE
ASSIGN 10 TO i99999
GO TO 270
10 fb = fx
xlo = xhi
a = xlo
x = xlo
C GET-FUNCTION-VALUE
ASSIGN 20 TO i99999
GO TO 270
C
C Check that F(ZXLO) < 0 < F(ZXHI) or
C F(ZXLO) > 0 > F(ZXHI)
C
20 IF (.NOT. (fb.LT.0.0D0)) GO TO 40
IF (.NOT. (fx.LT.0.0D0)) GO TO 30
status = -1
qleft = fx .LT. fb
qhi = .FALSE.
RETURN
30 CONTINUE
40 IF (.NOT. (fb.GT.0.0D0)) GO TO 60
IF (.NOT. (fx.GT.0.0D0)) GO TO 50
status = -1
qleft = fx .GT. fb
qhi = .TRUE.
RETURN
50 CONTINUE
60 fa = fx
C
first = .TRUE.
70 c = a
fc = fa
ext = 0
80 IF (.NOT. (abs(fc).LT.abs(fb))) GO TO 100
IF (.NOT. (c.NE.a)) GO TO 90
d = a
fd = fa
90 a = b
fa = fb
xlo = c
b = xlo
fb = fc
c = a
fc = fa
100 tol = ftol(xlo)
m = (c+b)*.5D0
mb = m - b
IF (.NOT. (abs(mb).GT.tol)) GO TO 240
IF (.NOT. (ext.GT.3)) GO TO 110
w = mb
GO TO 190
110 tol = sign(tol,mb)
p = (b-a)*fb
IF (.NOT. (first)) GO TO 120
q = fa - fb
first = .FALSE.
GO TO 130
120 fdb = (fd-fb)/ (d-b)
fda = (fd-fa)/ (d-a)
p = fda*p
q = fdb*fa - fda*fb
130 IF (.NOT. (p.LT.0.0D0)) GO TO 140
p = -p
q = -q
140 IF (ext.EQ.3) p = p*2.0D0
IF (.NOT. ((p*1.0D0).EQ.0.0D0.OR.p.LE. (q*tol))) GO TO 150
w = tol
GO TO 180
150 IF (.NOT. (p.LT. (mb*q))) GO TO 160
w = p/q
GO TO 170
160 w = mb
170 CONTINUE
180 CONTINUE
190 d = a
fd = fa
a = b
fa = fb
b = b + w
xlo = b
x = xlo
C GET-FUNCTION-VALUE
ASSIGN 200 TO i99999
GO TO 270
200 fb = fx
IF (.NOT. ((fc*fb).GE.0.0D0)) GO TO 210
GO TO 70
210 IF (.NOT. (w.EQ.mb)) GO TO 220
ext = 0
GO TO 230
220 ext = ext + 1
230 GO TO 80
240 xhi = c
qrzero = (fc.GE.0.0D0 .AND. fb.LE.0.0D0) .OR.
+ (fc.LT.0.0D0 .AND. fb.GE.0.0D0)
IF (.NOT. (qrzero)) GO TO 250
status = 0
GO TO 260
250 status = -1
260 RETURN
ENTRY dstzr(zxlo,zxhi,zabstl,zreltl)
C**********************************************************************
C
C SUBROUTINE DSTZR( XLO, XHI, ABSTOL, RELTOL )
C Double precision SeT ZeRo finder - Reverse communication version
C
C
C Function
C
C
C
C Sets quantities needed by ZROR. The function of ZROR
C and the quantities set is given here.
C
C Concise Description - Given a function F
C find XLO such that F(XLO) = 0.
C
C More Precise Description -
C
C Input condition. F is a double precision function of a single
C double precision argument and XLO and XHI are such that
C F(XLO)*F(XHI) .LE. 0.0
C
C If the input condition is met, QRZERO returns .TRUE.
C and output values of XLO and XHI satisfy the following
C F(XLO)*F(XHI) .LE. 0.
C ABS(F(XLO) .LE. ABS(F(XHI)
C ABS(XLO-XHI) .LE. TOL(X)
C where
C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
C
C If this algorithm does not find XLO and XHI satisfying
C these conditions then QRZERO returns .FALSE. This
C implies that the input condition was not met.
C
C
C Arguments
C
C
C XLO --> The left endpoint of the interval to be
C searched for a solution.
C XLO is DOUBLE PRECISION
C
C XHI --> The right endpoint of the interval to be
C for a solution.
C XHI is DOUBLE PRECISION
C
C ABSTOL, RELTOL --> Two numbers that determine the accuracy
C of the solution. See function for a
C precise definition.
C ABSTOL is DOUBLE PRECISION
C RELTOL is DOUBLE PRECISION
C
C
C Method
C
C
C Algorithm R of the paper 'Two Efficient Algorithms with
C Guaranteed Convergence for Finding a Zero of a Function'
C by J. C. P. Bus and T. J. Dekker in ACM Transactions on
C Mathematical Software, Volume 1, no. 4 page 330
C (Dec. '75) is employed to find the zero of F(X)-Y.
C
C**********************************************************************
xxlo = zxlo
xxhi = zxhi
abstol = zabstl
reltol = zreltl
RETURN
C TO GET-FUNCTION-VALUE
270 status = 1
RETURN
280 CONTINUE
GO TO i99999
END
SHAR_EOF
fi # end of overwriting check
if test -f 'cdinvr.c'
then
echo shar: will not over-write existing file "'cdinvr.c'"
else
cat << \SHAR_EOF > 'cdinvr.c'
void cdinvr( f, x, y, status, qleft, qhi, small, big, absstp, relstp,
stpmul,abstol,reltol)
/* Note that variables starting with q are logical in the
sense that 0 is FALSE and 1 is TRUE */
void **f;
float *x, *y, *absstp, *relstp, **reltol, *abstol;
float *small, *big;
int *status, *qleft, *qhi;
{
void dinvr();
dinvr(f,x,y,status,qleft,qhi,small,big,absstp,relstp,stpmul,
abstol,reltol); }
SHAR_EOF
fi # end of overwriting check
if test -f 'devlsf.c'
then
echo shar: will not over-write existing file "'devlsf.c'"
else
cat << \SHAR_EOF > 'devlsf.c'
/***********************************************************************
file devlsf.c
++++++++++++++++++++
Evaluation of an S Function
of One Argument
______________________________________________________________________
This is c code to evaluate an S function f at x and return
the value in ans as type float.
Note that the (void) function name is devlsf_ so that it can be
invoked from Fortran.
f takes one argument which must be type float.
It is suggested that a strategy such as that described on
page 210 of New S be employed to assure that a single number
is returned by the user function.
If this is not done, the user must assure that one
float value is returned by his S function
Arguments
f --> S function of one argument of type double. f is
passed here through a trivial C driver that invokes a
Fortran routine in which all statements FX = F(X) are changed
to CALL DEVLSF( F, X, FX ).
x --> value at which f is to be evaluated.
x is double
ans <-- value of f(x).
ans is double
***********************************************************************/
#include
void devlsf( f, x, ans )
char **f; double *x; double *ans;
{ char *func;
char *args[1]; char *result[1];
char *mode[1]; long length[1];
mode[0] = "double"; length[0] = 1;
args[0] = (char *) x;
func = f[0];
call_S( func, 1L, args, mode, length, 0, 1L, result);
*ans = *( (double *) result[0] );
}
SHAR_EOF
fi # end of overwriting check
if test -f 'devlsfs.c'
then
echo shar: will not over-write existing file "'devlsfs.c'"
else
cat << \SHAR_EOF > 'devlsfs.c'
/***********************************************************************
file devlsf.c
++++++++++++++++++++
Evaluation of an S Function
of One Argument
______________________________________________________________________
This is c code to evaluate an S function f at x and return
the value in ans as type float.
Note that the (void) function name is devlsf_ so that it can be
invoked from Fortran.
f takes one argument which must be type float.
It is suggested that a strategy such as that described on
page 210 of New S be employed to assure that a single number
is returned by the user function.
If this is not done, the user must assure that one
float value is returned by his S function
Arguments
f --> S function of one argument of type double. f is
passed here through a trivial C driver that invokes a
Fortran routine in which all statements FX = F(X) are changed
to CALL DEVLSF( F, X, FX ).
x --> value at which f is to be evaluated.
x is double
ans <-- value of f(x).
ans is double
***********************************************************************/
#include
void devlsf_( f, x, ans )
char **f; double *x; double *ans;
{ char *func;
char *args[1]; char *result[1];
char *mode[1]; long length[1];
mode[0] = "double"; length[0] = 1;
args[0] = (char *) x;
func = f[0];
call_S( func, 1L, args, mode, length, 0, 1L, result);
*ans = *( (double *) result[0] );
}
SHAR_EOF
fi # end of overwriting check
if test -f 'cdinvrs.c'
then
echo shar: will not over-write existing file "'cdinvrs.c'"
else
cat << \SHAR_EOF > 'cdinvrs.c'
void cdinvr( f, x, y, status, qleft, qhi, small, big, absstp, relstp,
stpmul,abstol,reltol)
/* Note that variables starting with q are logical in the
sense that 0 is FALSE and 1 is TRUE */
void **f;
float *x, *y, *absstp, *relstp, **reltol, *abstol;
float *small, *big;
int *status, *qleft, *qhi;
{
void dinvr();
dinvr_(f,x,y,status,qleft,qhi,small,big,absstp,relstp,stpmul,
abstol,reltol); }
SHAR_EOF
fi # end of overwriting check
if test -f 'makefile.sun'
then
echo shar: will not over-write existing file "'makefile.sun'"
else
cat << \SHAR_EOF > 'makefile.sun'
# Makefile for invmf source
# FC is name of Fortran Compiler
FC = f77
# FCFLAGS are flags to be set for the compile
# Note the -c flag is hard wired in
FCFLAGS =
# CC is name of C Compiler
CC = cc
# CCFLAGS are flags to be set for the compile
# Note the -c flag is hard wired in
CCFLAGS =
# Libraries needed in the load process go here.
# There may be none, or perhaps Fortran or C Libraries
# This is very system specific
LIBS =
#----------------------------------------------------------------------
# The user probably doesn't want to change anything below here
TARGET = invmf.o
OBJ = \
cdinvrs.o \
devlsfs.o \
dinvr.o \
dzror.o
# Set up suffixes for created implicit rules
.SUFFIXES:
.SUFFIXES: .f .c .o $(SUFFIXES)
# Rule for .f to .o
.f.o:
$(FC) $(FCFLAGS) -c $*.f
# Rule for .c to .o
.c.o:
$(CC) $(CCFLAGS) -c $*.c
# Create relocatable load file
all:$(TARGET)
$(TARGET):$(OBJ)
ld -r -d -o $(TARGET) $(OBJ) $(LIBS)
clean:
rm $(OBJ)
SHAR_EOF
fi # end of overwriting check
cd ..
# End of shell archive
exit 0