Commit d574f95f authored by Sébastien Villemot's avatar Sébastien Villemot

Merge tag 'upstream/19961126+dfsg1'

Upstream version 19961126+dfsg1
parents 60808eae 2ca2ecb2
# Don't require all the GNU mandated files
# Remove comment from no-dependencies for distributions compatable with
# traditional 'make'
AUTOMAKE_OPTIONS = 1.2 foreign # no-dependencies
lib_LTLIBRARIES = libminpack.la
VERSION_CURRENT=@VERSION_CURRENT@
VERSION_REVISION=@VERSION_REVISION@
VERSION_AGE=@VERSION_AGE@
libminpack_la_LDFLAGS = -version-info $(VERSION_CURRENT):$(VERSION_REVISION):$(VERSION_AGE) -D_REENTRANT,Wl,-z,defs
# Library sources
libminpack_la_SOURCES = chkder.f dogleg.f dpmpar.f enorm.f fdjac1.f \
fdjac2.f hybrd.f hybrd1.f hybrj.f hybrj1.f lmder.f lmder1.f lmdif.f \
lmdif1.f lmpar.f lmstr.f lmstr1.f qform.f qrfac.f qrsolv.f r1mpyq.f \
r1updt.f rwupdt.f
# Source files not used for Unix library but which are distributed
#EXTRA_libminpack_la_SOURCES = mac.c nt.c vms.c
# Headers which are installed
include_HEADERS = minpack.h
# Headers which are not installed but which are distributed
#noinst_HEADERS = Colorlist.h animate.h display.h formats.h logo.h mac.h \
# nt.h pict.h version.h vms.h
# Additional files to distribute
#EXTRA_DIST = Imakefile Make.com
# Ensure that configuration header is not distributed
#dist-hook:
# $(RM) $(distdir)/config.h
INCLUDES = #$(X_CFLAGS)
# Since we are building a library, no need for LIBS
#LIBS =
This diff is collapsed.
This diff is collapsed.
subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
integer m,n,ldfjac,mode
double precision x(n),fvec(m),fjac(ldfjac,n),xp(n),fvecp(m),
* err(m)
c **********
c
c subroutine chkder
c
c this subroutine checks the gradients of m nonlinear functions
c in n variables, evaluated at a point x, for consistency with
c the functions themselves. the user must call chkder twice,
c first with mode = 1 and then with mode = 2.
c
c mode = 1. on input, x must contain the point of evaluation.
c on output, xp is set to a neighboring point.
c
c mode = 2. on input, fvec must contain the functions and the
c rows of fjac must contain the gradients
c of the respective functions each evaluated
c at x, and fvecp must contain the functions
c evaluated at xp.
c on output, err contains measures of correctness of
c the respective gradients.
c
c the subroutine does not perform reliably if cancellation or
c rounding errors cause a severe loss of significance in the
c evaluation of a function. therefore, none of the components
c of x should be unusually small (in particular, zero) or any
c other value which may cause loss of significance.
c
c the subroutine statement is
c
c subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
c
c where
c
c m is a positive integer input variable set to the number
c of functions.
c
c n is a positive integer input variable set to the number
c of variables.
c
c x is an input array of length n.
c
c fvec is an array of length m. on input when mode = 2,
c fvec must contain the functions evaluated at x.
c
c fjac is an m by n array. on input when mode = 2,
c the rows of fjac must contain the gradients of
c the respective functions evaluated at x.
c
c ldfjac is a positive integer input parameter not less than m
c which specifies the leading dimension of the array fjac.
c
c xp is an array of length n. on output when mode = 1,
c xp is set to a neighboring point of x.
c
c fvecp is an array of length m. on input when mode = 2,
c fvecp must contain the functions evaluated at xp.
c
c mode is an integer input variable set to 1 on the first call
c and 2 on the second. other values of mode are equivalent
c to mode = 1.
c
c err is an array of length m. on output when mode = 2,
c err contains measures of correctness of the respective
c gradients. if there is no severe loss of significance,
c then if err(i) is 1.0 the i-th gradient is correct,
c while if err(i) is 0.0 the i-th gradient is incorrect.
c for values of err between 0.0 and 1.0, the categorization
c is less certain. in general, a value of err(i) greater
c than 0.5 indicates that the i-th gradient is probably
c correct, while a value of err(i) less than 0.5 indicates
c that the i-th gradient is probably incorrect.
c
c subprograms called
c
c minpack supplied ... dpmpar
c
c fortran supplied ... dabs,dlog10,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j
double precision eps,epsf,epslog,epsmch,factor,one,temp,zero
double precision dpmpar
data factor,one,zero /1.0d2,1.0d0,0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
eps = dsqrt(epsmch)
c
if (mode .eq. 2) go to 20
c
c mode = 1.
c
do 10 j = 1, n
temp = eps*dabs(x(j))
if (temp .eq. zero) temp = eps
xp(j) = x(j) + temp
10 continue
go to 70
20 continue
c
c mode = 2.
c
epsf = factor*epsmch
epslog = dlog10(eps)
do 30 i = 1, m
err(i) = zero
30 continue
do 50 j = 1, n
temp = dabs(x(j))
if (temp .eq. zero) temp = one
do 40 i = 1, m
err(i) = err(i) + temp*fjac(i,j)
40 continue
50 continue
do 60 i = 1, m
temp = one
if (fvec(i) .ne. zero .and. fvecp(i) .ne. zero
* .and. dabs(fvecp(i)-fvec(i)) .ge. epsf*dabs(fvec(i)))
* temp = eps*dabs((fvecp(i)-fvec(i))/eps-err(i))
* /(dabs(fvec(i)) + dabs(fvecp(i)))
err(i) = one
if (temp .gt. epsmch .and. temp .lt. eps)
* err(i) = (dlog10(temp) - epslog)/epslog
if (temp .ge. eps) err(i) = zero
60 continue
70 continue
c
return
c
c last card of subroutine chkder.
c
end
This diff is collapsed.
/* config.h.in. Generated automatically from configure.in by autoheader 2.13. */
/* Define if you have the <dlfcn.h> header file. */
#undef HAVE_DLFCN_H
/* Name of package */
#undef PACKAGE
/* Version number of package */
#undef VERSION
This diff is collapsed.
This diff is collapsed.
dnl Process this file with autoconf to produce a configure script.
AC_INIT(lmdif.f)
dnl Every other copy of the package version number gets its value from here
AM_INIT_AUTOMAKE(minpack, 19961126)
# shared library version control
VERSION_CURRENT=1
VERSION_REVISION=0
VERSION_AGE=0
AC_SUBST(VERSION)
AC_SUBST(VERSION_CURRENT)
AC_SUBST(VERSION_REVISION)
AC_SUBST(VERSION_AGE)
dnl Checks for programs.
AC_PROG_INSTALL
AC_PROG_F77
AC_PROG_CC
AM_PROG_LIBTOOL
dnl maybe bash is not in /bin on this system
AC_PATH_PROG(bash, bash, FAIL)
if test "$bash" = "FAIL"; then
AC_MSG_ERROR(Cannot continue: bash not found)
fi
dnl Checks for libraries.
dnl Checks for header files.
dnl Checks for typedefs, structures, and compiler characteristics.
dnl Checks for library functions.
AC_OUTPUT(Makefile)
subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
integer n,lr
double precision delta
double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
c **********
c
c subroutine dogleg
c
c given an m by n matrix a, an n by n nonsingular diagonal
c matrix d, an m-vector b, and a positive number delta, the
c problem is to determine the convex combination x of the
c gauss-newton and scaled gradient directions that minimizes
c (a*x - b) in the least squares sense, subject to the
c restriction that the euclidean norm of d*x be at most delta.
c
c this subroutine completes the solution of the problem
c if it is provided with the necessary information from the
c qr factorization of a. that is, if a = q*r, where q has
c orthogonal columns and r is an upper triangular matrix,
c then dogleg expects the full upper triangle of r and
c the first n components of (q transpose)*b.
c
c the subroutine statement is
c
c subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c
c where
c
c n is a positive integer input variable set to the order of r.
c
c r is an input array of length lr which must contain the upper
c triangular matrix r stored by rows.
c
c lr is a positive integer input variable not less than
c (n*(n+1))/2.
c
c diag is an input array of length n which must contain the
c diagonal elements of the matrix d.
c
c qtb is an input array of length n which must contain the first
c n elements of the vector (q transpose)*b.
c
c delta is a positive input variable which specifies an upper
c bound on the euclidean norm of d*x.
c
c x is an output array of length n which contains the desired
c convex combination of the gauss-newton direction and the
c scaled gradient direction.
c
c wa1 and wa2 are work arrays of length n.
c
c subprograms called
c
c minpack-supplied ... dpmpar,enorm
c
c fortran-supplied ... dabs,dmax1,dmin1,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jj,jp1,k,l
double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
* temp,zero
double precision dpmpar,enorm
data one,zero /1.0d0,0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
c first, calculate the gauss-newton direction.
c
jj = (n*(n + 1))/2 + 1
do 50 k = 1, n
j = n - k + 1
jp1 = j + 1
jj = jj - k
l = jj + 1
sum = zero
if (n .lt. jp1) go to 20
do 10 i = jp1, n
sum = sum + r(l)*x(i)
l = l + 1
10 continue
20 continue
temp = r(jj)
if (temp .ne. zero) go to 40
l = j
do 30 i = 1, j
temp = dmax1(temp,dabs(r(l)))
l = l + n - i
30 continue
temp = epsmch*temp
if (temp .eq. zero) temp = epsmch
40 continue
x(j) = (qtb(j) - sum)/temp
50 continue
c
c test whether the gauss-newton direction is acceptable.
c
do 60 j = 1, n
wa1(j) = zero
wa2(j) = diag(j)*x(j)
60 continue
qnorm = enorm(n,wa2)
if (qnorm .le. delta) go to 140
c
c the gauss-newton direction is not acceptable.
c next, calculate the scaled gradient direction.
c
l = 1
do 80 j = 1, n
temp = qtb(j)
do 70 i = j, n
wa1(i) = wa1(i) + r(l)*temp
l = l + 1
70 continue
wa1(j) = wa1(j)/diag(j)
80 continue
c
c calculate the norm of the scaled gradient and test for
c the special case in which the scaled gradient is zero.
c
gnorm = enorm(n,wa1)
sgnorm = zero
alpha = delta/qnorm
if (gnorm .eq. zero) go to 120
c
c calculate the point along the scaled gradient
c at which the quadratic is minimized.
c
do 90 j = 1, n
wa1(j) = (wa1(j)/gnorm)/diag(j)
90 continue
l = 1
do 110 j = 1, n
sum = zero
do 100 i = j, n
sum = sum + r(l)*wa1(i)
l = l + 1
100 continue
wa2(j) = sum
110 continue
temp = enorm(n,wa2)
sgnorm = (gnorm/temp)/temp
c
c test whether the scaled gradient direction is acceptable.
c
alpha = zero
if (sgnorm .ge. delta) go to 120
c
c the scaled gradient direction is not acceptable.
c finally, calculate the point along the dogleg
c at which the quadratic is minimized.
c
bnorm = enorm(n,qtb)
temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
temp = temp - (delta/qnorm)*(sgnorm/delta)**2
* + dsqrt((temp-(delta/qnorm))**2
* +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
120 continue
c
c form appropriate convex combination of the gauss-newton
c direction and the scaled gradient direction.
c
temp = (one - alpha)*dmin1(sgnorm,delta)
do 130 j = 1, n
x(j) = temp*wa1(j) + alpha*x(j)
130 continue
140 continue
return
c
c last card of subroutine dogleg.
c
end
double precision function dpmpar(i)
integer i
c **********
c
c Function dpmpar
c
c This function provides double precision machine parameters
c when the appropriate set of data statements is activated (by
c removing the c from column 1) and all other data statements are
c rendered inactive. Most of the parameter values were obtained
c from the corresponding Bell Laboratories Port Library function.
c
c The function statement is
c
c double precision function dpmpar(i)
c
c where
c
c i is an integer input variable set to 1, 2, or 3 which
c selects the desired machine parameter. If the machine has
c t base b digits and its smallest and largest exponents are
c emin and emax, respectively, then these parameters are
c
c dpmpar(1) = b**(1 - t), the machine precision,
c
c dpmpar(2) = b**(emin - 1), the smallest magnitude,
c
c dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
c
c Argonne National Laboratory. MINPACK Project. November 1996.
c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
c
c **********
integer mcheps(4)
integer minmag(4)
integer maxmag(4)
double precision dmach(3)
equivalence (dmach(1),mcheps(1))
equivalence (dmach(2),minmag(1))
equivalence (dmach(3),maxmag(1))
c
c Machine constants for the IBM 360/370 series,
c the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
c the Xerox Sigma 5/7/9 and the Sel systems 85/86.
c
c data mcheps(1),mcheps(2) / z34100000, z00000000 /
c data minmag(1),minmag(2) / z00100000, z00000000 /
c data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
c
c Machine constants for the Honeywell 600/6000 series.
c
c data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
c data minmag(1),minmag(2) / o402400000000, o000000000000 /
c data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
c
c Machine constants for the CDC 6000/7000 series.
c
c data mcheps(1) / 15614000000000000000b /
c data mcheps(2) / 15010000000000000000b /
c
c data minmag(1) / 00604000000000000000b /
c data minmag(2) / 00000000000000000000b /
c
c data maxmag(1) / 37767777777777777777b /
c data maxmag(2) / 37167777777777777777b /
c
c Machine constants for the PDP-10 (KA processor).
c
c data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
c data minmag(1),minmag(2) / "033400000000, "000000000000 /
c data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
c
c Machine constants for the PDP-10 (KI processor).
c
c data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
c data minmag(1),minmag(2) / "000400000000, "000000000000 /
c data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
c
c Machine constants for the PDP-11.
c
c data mcheps(1),mcheps(2) / 9472, 0 /
c data mcheps(3),mcheps(4) / 0, 0 /
c
c data minmag(1),minmag(2) / 128, 0 /
c data minmag(3),minmag(4) / 0, 0 /
c
c data maxmag(1),maxmag(2) / 32767, -1 /
c data maxmag(3),maxmag(4) / -1, -1 /
c
c Machine constants for the Burroughs 6700/7700 systems.
c
c data mcheps(1) / o1451000000000000 /
c data mcheps(2) / o0000000000000000 /
c
c data minmag(1) / o1771000000000000 /
c data minmag(2) / o7770000000000000 /
c
c data maxmag(1) / o0777777777777777 /
c data maxmag(2) / o7777777777777777 /
c
c Machine constants for the Burroughs 5700 system.
c
c data mcheps(1) / o1451000000000000 /
c data mcheps(2) / o0000000000000000 /
c
c data minmag(1) / o1771000000000000 /
c data minmag(2) / o0000000000000000 /
c
c data maxmag(1) / o0777777777777777 /
c data maxmag(2) / o0007777777777777 /
c
c Machine constants for the Burroughs 1700 system.
c
c data mcheps(1) / zcc6800000 /
c data mcheps(2) / z000000000 /
c
c data minmag(1) / zc00800000 /
c data minmag(2) / z000000000 /
c
c data maxmag(1) / zdffffffff /
c data maxmag(2) / zfffffffff /
c
c Machine constants for the Univac 1100 series.
c
c data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
c data minmag(1),minmag(2) / o000040000000, o000000000000 /
c data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
c
c Machine constants for the Data General Eclipse S/200.
c
c Note - it may be appropriate to include the following card -
c static dmach(3)
c
c data minmag/20k,3*0/,maxmag/77777k,3*177777k/
c data mcheps/32020k,3*0/
c
c Machine constants for the Harris 220.
c
c data mcheps(1),mcheps(2) / '20000000, '00000334 /
c data minmag(1),minmag(2) / '20000000, '00000201 /
c data maxmag(1),maxmag(2) / '37777777, '37777577 /
c
c Machine constants for the Cray-1.
c
c data mcheps(1) / 0376424000000000000000b /
c data mcheps(2) / 0000000000000000000000b /
c
c data minmag(1) / 0200034000000000000000b /
c data minmag(2) / 0000000000000000000000b /
c
c data maxmag(1) / 0577777777777777777777b /
c data maxmag(2) / 0000007777777777777776b /
c
c Machine constants for the Prime 400.
c
c data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
c data minmag(1),minmag(2) / :10000000000, :00000100000 /
c data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
c
c Machine constants for the VAX-11.
c
c data mcheps(1),mcheps(2) / 9472, 0 /
c data minmag(1),minmag(2) / 128, 0 /
c data maxmag(1),maxmag(2) / -32769, -1 /
c
c Machine constants for IEEE machines.
c
data dmach(1) /2.22044604926d-16/
data dmach(2) /2.22507385852d-308/
data dmach(3) /1.79769313485d+308/
c
dpmpar = dmach(i)
return
c
c Last card of function dpmpar.
c
end
double precision function enorm(n,x)
integer n
double precision x(n)
c **********
c
c function enorm
c
c given an n-vector x, this function calculates the
c euclidean norm of x.
c
c the euclidean norm is computed by accumulating the sum of
c squares in three different sums. the sums of squares for the
c small and large components are scaled so that no overflows
c occur. non-destructive underflows are permitted. underflows
c and overflows do not occur in the computation of the unscaled
c sum of squares for the intermediate components.
c the definitions of small, intermediate and large components
c depend on two constants, rdwarf and rgiant. the main
c restrictions on these constants are that rdwarf**2 not
c underflow and rgiant**2 not overflow. the constants
c given here are suitable for every known computer.
c
c the function statement is
c
c double precision function enorm(n,x)
c
c where
c
c n is a positive integer input variable.
c
c x is an input array of length n.
c
c subprograms called
c
c fortran-supplied ... dabs,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i
double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,
* x1max,x3max,zero
data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
s1 = zero
s2 = zero
s3 = zero
x1max = zero
x3max = zero
floatn = n
agiant = rgiant/floatn
do 90 i = 1, n
xabs = dabs(x(i))
if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
if (xabs .le. rdwarf) go to 30
c
c sum for large components.
c
if (xabs .le. x1max) go to 10
s1 = one + s1*(x1max/xabs)**2
x1max = xabs
go to 20
10 continue
s1 = s1 + (xabs/x1max)**2
20 continue
go to 60
30 continue
c
c sum for small components.
c
if (xabs .le. x3max) go to 40
s3 = one + s3*(x3max/xabs)**2
x3max = xabs
go to 50
40 continue
if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
50 continue
60 continue
go to 80
70 continue
c
c sum for intermediate components.
c
s2 = s2 + xabs**2
80 continue
90 continue
c