Repatch (from linux) of OSX IRAF

23 files changed, 4949 insertions, 0 deletions

diff --git a/math/minpack/chkder.f b/math/minpack/chkder.f
new file mode 100644
index 00000000..29578fc4
--- /dev/null
+++ b/math/minpack/chkder.f
@@ -0,0 +1,140 @@
+      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
diff --git a/math/minpack/dogleg.f b/math/minpack/dogleg.f
new file mode 100644
index 00000000..b812f196
--- /dev/null
+++ b/math/minpack/dogleg.f
@@ -0,0 +1,177 @@
+      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
diff --git a/math/minpack/dpmpar.f b/math/minpack/dpmpar.f
new file mode 100644
index 00000000..776bf3d1
--- /dev/null
+++ b/math/minpack/dpmpar.f
@@ -0,0 +1,171 @@
+      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. June 1983.
+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
+      data mcheps(1),mcheps(2) /   9472,  0 /
+      data minmag(1),minmag(2) /    128,  0 /
+      data maxmag(1),maxmag(2) / -32769, -1 /
+c
+      dpmpar = dmach(i)
+      return
+c
+c     Last card of function dpmpar.
+c
+      end
diff --git a/math/minpack/enorm.f b/math/minpack/enorm.f
new file mode 100644
index 00000000..2cb5b607
--- /dev/null
+++ b/math/minpack/enorm.f
@@ -0,0 +1,108 @@
+      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
+c     calculation of norm.
+c
+      if (s1 .eq. zero) go to 100
+         enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
+         go to 130
+  100 continue
+         if (s2 .eq. zero) go to 110
+            if (s2 .ge. x3max)
+     *         enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
+            if (s2 .lt. x3max)
+     *         enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
+            go to 120
+  110    continue
+            enorm = x3max*dsqrt(s3)
+  120    continue
+  130 continue
+      return
+c
+c     last card of function enorm.
+c
+      end
diff --git a/math/minpack/fdjac1.f b/math/minpack/fdjac1.f
new file mode 100644
index 00000000..aa058010
--- /dev/null
+++ b/math/minpack/fdjac1.f
@@ -0,0 +1,150 @@
+      subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
+     *                  wa1,wa2)
+      integer n,ldfjac,iflag,ml,mu
+      double precision epsfcn
+      double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
+c     **********
+c
+c     subroutine fdjac1
+c
+c     this subroutine computes a forward-difference approximation
+c     to the n by n jacobian matrix associated with a specified
+c     problem of n functions in n variables. if the jacobian has
+c     a banded form, then function evaluations are saved by only
+c     approximating the nonzero terms.
+c
+c     the subroutine statement is
+c
+c       subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
+c                         wa1,wa2)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,iflag)
+c         integer n,iflag
+c         double precision x(n),fvec(n)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of fdjac1.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an input array of length n.
+c
+c       fvec is an input array of length n which must contain the
+c         functions evaluated at x.
+c
+c       fjac is an output n by n array which contains the
+c         approximation to the jacobian matrix evaluated at x.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       iflag is an integer variable which can be used to terminate
+c         the execution of fdjac1. see description of fcn.
+c
+c       ml is a nonnegative integer input variable which specifies
+c         the number of subdiagonals within the band of the
+c         jacobian matrix. if the jacobian is not banded, set
+c         ml to at least n - 1.
+c
+c       epsfcn is an input variable used in determining a suitable
+c         step length for the forward-difference approximation. this
+c         approximation assumes that the relative errors in the
+c         functions are of the order of epsfcn. if epsfcn is less
+c         than the machine precision, it is assumed that the relative
+c         errors in the functions are of the order of the machine
+c         precision.
+c
+c       mu is a nonnegative integer input variable which specifies
+c         the number of superdiagonals within the band of the
+c         jacobian matrix. if the jacobian is not banded, set
+c         mu to at least n - 1.
+c
+c       wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
+c         least n, then the jacobian is considered dense, and wa2 is
+c         not referenced.
+c
+c     subprograms called
+c
+c       minpack-supplied ... dpmpar
+c
+c       fortran-supplied ... dabs,dmax1,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,k,msum
+      double precision eps,epsmch,h,temp,zero
+      double precision dpmpar
+      data zero /0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      eps = dsqrt(dmax1(epsfcn,epsmch))
+      msum = ml + mu + 1
+      if (msum .lt. n) go to 40
+c
+c        computation of dense approximate jacobian.
+c
+         do 20 j = 1, n
+            temp = x(j)
+            h = eps*dabs(temp)
+            if (h .eq. zero) h = eps
+            x(j) = temp + h
+            call fcn(n,x,wa1,iflag)
+            if (iflag .lt. 0) go to 30
+            x(j) = temp
+            do 10 i = 1, n
+               fjac(i,j) = (wa1(i) - fvec(i))/h
+   10          continue
+   20       continue
+   30    continue
+         go to 110
+   40 continue
+c
+c        computation of banded approximate jacobian.
+c
+         do 90 k = 1, msum
+            do 60 j = k, n, msum
+               wa2(j) = x(j)
+               h = eps*dabs(wa2(j))
+               if (h .eq. zero) h = eps
+               x(j) = wa2(j) + h
+   60          continue
+            call fcn(n,x,wa1,iflag)
+            if (iflag .lt. 0) go to 100
+            do 80 j = k, n, msum
+               x(j) = wa2(j)
+               h = eps*dabs(wa2(j))
+               if (h .eq. zero) h = eps
+               do 70 i = 1, n
+                  fjac(i,j) = zero
+                  if (i .ge. j - mu .and. i .le. j + ml)
+     *               fjac(i,j) = (wa1(i) - fvec(i))/h
+   70             continue
+   80          continue
+   90       continue
+  100    continue
+  110 continue
+      return
+c
+c     last card of subroutine fdjac1.
+c
+      end
diff --git a/math/minpack/fdjac2.f b/math/minpack/fdjac2.f
new file mode 100644
index 00000000..218ab94c
--- /dev/null
+++ b/math/minpack/fdjac2.f
@@ -0,0 +1,107 @@
+      subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
+      integer m,n,ldfjac,iflag
+      double precision epsfcn
+      double precision x(n),fvec(m),fjac(ldfjac,n),wa(m)
+c     **********
+c
+c     subroutine fdjac2
+c
+c     this subroutine computes a forward-difference approximation
+c     to the m by n jacobian matrix associated with a specified
+c     problem of m functions in n variables.
+c
+c     the subroutine statement is
+c
+c       subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,iflag)
+c         integer m,n,iflag
+c         double precision x(n),fvec(m)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of fdjac2.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an input array of length n.
+c
+c       fvec is an input array of length m which must contain the
+c         functions evaluated at x.
+c
+c       fjac is an output m by n array which contains the
+c         approximation to the jacobian matrix evaluated at x.
+c
+c       ldfjac is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array fjac.
+c
+c       iflag is an integer variable which can be used to terminate
+c         the execution of fdjac2. see description of fcn.
+c
+c       epsfcn is an input variable used in determining a suitable
+c         step length for the forward-difference approximation. this
+c         approximation assumes that the relative errors in the
+c         functions are of the order of epsfcn. if epsfcn is less
+c         than the machine precision, it is assumed that the relative
+c         errors in the functions are of the order of the machine
+c         precision.
+c
+c       wa is a work array of length m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dpmpar
+c
+c       fortran-supplied ... dabs,dmax1,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,epsmch,h,temp,zero
+      double precision dpmpar
+      data zero /0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      eps = dsqrt(dmax1(epsfcn,epsmch))
+      do 20 j = 1, n
+         temp = x(j)
+         h = eps*dabs(temp)
+         if (h .eq. zero) h = eps
+         x(j) = temp + h
+         call fcn(m,n,x,wa,iflag)
+         if (iflag .lt. 0) go to 30
+         x(j) = temp
+         do 10 i = 1, m
+            fjac(i,j) = (wa(i) - fvec(i))/h
+   10       continue
+   20    continue
+   30 continue
+      return
+c
+c     last card of subroutine fdjac2.
+c
+      end
diff --git a/math/minpack/hybrd.f b/math/minpack/hybrd.f
new file mode 100644
index 00000000..fc0b4c26
--- /dev/null
+++ b/math/minpack/hybrd.f
@@ -0,0 +1,459 @@
+      subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,
+     *                 mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr,
+     *                 qtf,wa1,wa2,wa3,wa4)
+      integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
+      double precision xtol,epsfcn,factor
+      double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),
+     *                 qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
+      external fcn
+c     **********
+c
+c     subroutine hybrd
+c
+c     the purpose of hybrd is to find a zero of a system of
+c     n nonlinear functions in n variables by a modification
+c     of the powell hybrid method. the user must provide a
+c     subroutine which calculates the functions. the jacobian is
+c     then calculated by a forward-difference approximation.
+c
+c     the subroutine statement is
+c
+c       subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,
+c                        diag,mode,factor,nprint,info,nfev,fjac,
+c                        ldfjac,r,lr,qtf,wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,iflag)
+c         integer n,iflag
+c         double precision x(n),fvec(n)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ---------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of hybrd.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length n which contains
+c         the functions evaluated at the output x.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn is at least maxfev
+c         by the end of an iteration.
+c
+c       ml is a nonnegative integer input variable which specifies
+c         the number of subdiagonals within the band of the
+c         jacobian matrix. if the jacobian is not banded, set
+c         ml to at least n - 1.
+c
+c       mu is a nonnegative integer input variable which specifies
+c         the number of superdiagonals within the band of the
+c         jacobian matrix. if the jacobian is not banded, set
+c         mu to at least n - 1.
+c
+c       epsfcn is an input variable used in determining a suitable
+c         step length for the forward-difference approximation. this
+c         approximation assumes that the relative errors in the
+c         functions are of the order of epsfcn. if epsfcn is less
+c         than the machine precision, it is assumed that the relative
+c         errors in the functions are of the order of the machine
+c         precision.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.). 100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x and fvec available
+c         for printing. if nprint is not positive, no special calls
+c         of fcn with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0   improper input parameters.
+c
+c         info = 1   relative error between two consecutive iterates
+c                    is at most xtol.
+c
+c         info = 2   number of calls to fcn has reached or exceeded
+c                    maxfev.
+c
+c         info = 3   xtol is too small. no further improvement in
+c                    the approximate solution x is possible.
+c
+c         info = 4   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    five jacobian evaluations.
+c
+c         info = 5   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    ten iterations.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn.
+c
+c       fjac is an output n by n array which contains the
+c         orthogonal matrix q produced by the qr factorization
+c         of the final approximate jacobian.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       r is an output array of length lr which contains the
+c         upper triangular matrix produced by the qr factorization
+c         of the final approximate jacobian, stored rowwise.
+c
+c       lr is a positive integer input variable not less than
+c         (n*(n+1))/2.
+c
+c       qtf is an output array of length n which contains
+c         the vector (q transpose)*fvec.
+c
+c       wa1, wa2, wa3, and wa4 are work arrays of length n.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dogleg,dpmpar,enorm,fdjac1,
+c                            qform,qrfac,r1mpyq,r1updt
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,min0,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,jm1,l,msum,ncfail,ncsuc,nslow1,nslow2
+      integer iwa(1)
+      logical jeval,sing
+      double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
+     *                 prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
+     *                 zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p001,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0
+     *    .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero
+     *    .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) go to 300
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 300
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(n,x,fvec,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 300
+      fnorm = enorm(n,fvec)
+c
+c     determine the number of calls to fcn needed to compute
+c     the jacobian matrix.
+c
+      msum = min0(ml+mu+1,n)
+c
+c     initialize iteration counter and monitors.
+c
+      iter = 1
+      ncsuc = 0
+      ncfail = 0
+      nslow1 = 0
+      nslow2 = 0
+c
+c     beginning of the outer loop.
+c
+   30 continue
+         jeval = .true.
+c
+c        calculate the jacobian matrix.
+c
+         iflag = 2
+         call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,
+     *               wa2)
+         nfev = nfev + msum
+         if (iflag .lt. 0) go to 300
+c
+c        compute the qr factorization of the jacobian.
+c
+         call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 70
+         if (mode .eq. 2) go to 50
+         do 40 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+   40       continue
+   50    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 60 j = 1, n
+            wa3(j) = diag(j)*x(j)
+   60       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+   70    continue
+c
+c        form (q transpose)*fvec and store in qtf.
+c
+         do 80 i = 1, n
+            qtf(i) = fvec(i)
+   80       continue
+         do 120 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 110
+            sum = zero
+            do 90 i = j, n
+               sum = sum + fjac(i,j)*qtf(i)
+   90          continue
+            temp = -sum/fjac(j,j)
+            do 100 i = j, n
+               qtf(i) = qtf(i) + fjac(i,j)*temp
+  100          continue
+  110       continue
+  120       continue
+c
+c        copy the triangular factor of the qr factorization into r.
+c
+         sing = .false.
+         do 150 j = 1, n
+            l = j
+            jm1 = j - 1
+            if (jm1 .lt. 1) go to 140
+            do 130 i = 1, jm1
+               r(l) = fjac(i,j)
+               l = l + n - i
+  130          continue
+  140       continue
+            r(l) = wa1(j)
+            if (wa1(j) .eq. zero) sing = .true.
+  150       continue
+c
+c        accumulate the orthogonal factor in fjac.
+c
+         call qform(n,n,fjac,ldfjac,wa1)
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 170
+         do 160 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  160       continue
+  170    continue
+c
+c        beginning of the inner loop.
+c
+  180    continue
+c
+c           if requested, call fcn to enable printing of iterates.
+c
+            if (nprint .le. 0) go to 190
+            iflag = 0
+            if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag)
+            if (iflag .lt. 0) go to 300
+  190       continue
+c
+c           determine the direction p.
+c
+            call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 200 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  200          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(n,wa2,wa4,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 300
+            fnorm1 = enorm(n,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction.
+c
+            l = 1
+            do 220 i = 1, n
+               sum = zero
+               do 210 j = i, n
+                  sum = sum + r(l)*wa1(j)
+                  l = l + 1
+  210             continue
+               wa3(i) = qtf(i) + sum
+  220          continue
+            temp = enorm(n,wa3)
+            prered = zero
+            if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .gt. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .ge. p1) go to 230
+               ncsuc = 0
+               ncfail = ncfail + 1
+               delta = p5*delta
+               go to 240
+  230       continue
+               ncfail = 0
+               ncsuc = ncsuc + 1
+               if (ratio .ge. p5 .or. ncsuc .gt. 1)
+     *            delta = dmax1(delta,pnorm/p5)
+               if (dabs(ratio-one) .le. p1) delta = pnorm/p5
+  240       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 260
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 250 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+               fvec(j) = wa4(j)
+  250          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  260       continue
+c
+c           determine the progress of the iteration.
+c
+            nslow1 = nslow1 + 1
+            if (actred .ge. p001) nslow1 = 0
+            if (jeval) nslow2 = nslow2 + 1
+            if (actred .ge. p1) nslow2 = 0
+c
+c           test for convergence.
+c
+            if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
+            if (info .ne. 0) go to 300
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 2
+            if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
+            if (nslow2 .eq. 5) info = 4
+            if (nslow1 .eq. 10) info = 5
+            if (info .ne. 0) go to 300
+c
+c           criterion for recalculating jacobian approximation
+c           by forward differences.
+c
+            if (ncfail .eq. 2) go to 290
+c
+c           calculate the rank one modification to the jacobian
+c           and update qtf if necessary.
+c
+            do 280 j = 1, n
+               sum = zero
+               do 270 i = 1, n
+                  sum = sum + fjac(i,j)*wa4(i)
+  270             continue
+               wa2(j) = (sum - wa3(j))/pnorm
+               wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
+               if (ratio .ge. p0001) qtf(j) = sum
+  280          continue
+c
+c           compute the qr factorization of the updated jacobian.
+c
+            call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
+            call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
+            call r1mpyq(1,n,qtf,1,wa2,wa3)
+c
+c           end of the inner loop.
+c
+            jeval = .false.
+            go to 180
+  290    continue
+c
+c        end of the outer loop.
+c
+         go to 30
+  300 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(n,x,fvec,iflag)
+      return
+c
+c     last card of subroutine hybrd.
+c
+      end
diff --git a/math/minpack/hybrd1.f b/math/minpack/hybrd1.f
new file mode 100644
index 00000000..c0a85927
--- /dev/null
+++ b/math/minpack/hybrd1.f
@@ -0,0 +1,123 @@
+      subroutine hybrd1(fcn,n,x,fvec,tol,info,wa,lwa)
+      integer n,info,lwa
+      double precision tol
+      double precision x(n),fvec(n),wa(lwa)
+      external fcn
+c     **********
+c
+c     subroutine hybrd1
+c
+c     the purpose of hybrd1 is to find a zero of a system of
+c     n nonlinear functions in n variables by a modification
+c     of the powell hybrid method. this is done by using the
+c     more general nonlinear equation solver hybrd. the user
+c     must provide a subroutine which calculates the functions.
+c     the jacobian is then calculated by a forward-difference
+c     approximation.
+c
+c     the subroutine statement is
+c
+c       subroutine hybrd1(fcn,n,x,fvec,tol,info,wa,lwa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,iflag)
+c         integer n,iflag
+c         double precision x(n),fvec(n)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ---------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of hybrd1.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length n which contains
+c         the functions evaluated at the output x.
+c
+c       tol is a nonnegative input variable. termination occurs
+c         when the algorithm estimates that the relative error
+c         between x and the solution is at most tol.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0   improper input parameters.
+c
+c         info = 1   algorithm estimates that the relative error
+c                    between x and the solution is at most tol.
+c
+c         info = 2   number of calls to fcn has reached or exceeded
+c                    200*(n+1).
+c
+c         info = 3   tol is too small. no further improvement in
+c                    the approximate solution x is possible.
+c
+c         info = 4   iteration is not making good progress.
+c
+c       wa is a work array of length lwa.
+c
+c       lwa is a positive integer input variable not less than
+c         (n*(3*n+13))/2.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... hybrd
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer index,j,lr,maxfev,ml,mode,mu,nfev,nprint
+      double precision epsfcn,factor,one,xtol,zero
+      data factor,one,zero /1.0d2,1.0d0,0.0d0/
+      info = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. tol .lt. zero .or. lwa .lt. (n*(3*n + 13))/2)
+     *   go to 20
+c
+c     call hybrd.
+c
+      maxfev = 200*(n + 1)
+      xtol = tol
+      ml = n - 1
+      mu = n - 1
+      epsfcn = zero
+      mode = 2
+      do 10 j = 1, n
+         wa(j) = one
+   10    continue
+      nprint = 0
+      lr = (n*(n + 1))/2
+      index = 6*n + lr
+      call hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,wa(1),mode,
+     *           factor,nprint,info,nfev,wa(index+1),n,wa(6*n+1),lr,
+     *           wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
+      if (info .eq. 5) info = 4
+   20 continue
+      return
+c
+c     last card of subroutine hybrd1.
+c
+      end
diff --git a/math/minpack/hybrj.f b/math/minpack/hybrj.f
new file mode 100644
index 00000000..3070dad3
--- /dev/null
+++ b/math/minpack/hybrj.f
@@ -0,0 +1,440 @@
+      subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode,
+     *                 factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,
+     *                 wa3,wa4)
+      integer n,ldfjac,maxfev,mode,nprint,info,nfev,njev,lr
+      double precision xtol,factor
+      double precision x(n),fvec(n),fjac(ldfjac,n),diag(n),r(lr),
+     *                 qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
+c     **********
+c
+c     subroutine hybrj
+c
+c     the purpose of hybrj is to find a zero of a system of
+c     n nonlinear functions in n variables by a modification
+c     of the powell hybrid method. the user must provide a
+c     subroutine which calculates the functions and the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,
+c                        mode,factor,nprint,info,nfev,njev,r,lr,qtf,
+c                        wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the jacobian. fcn must
+c         be declared in an external statement in the user
+c         calling program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
+c         integer n,ldfjac,iflag
+c         double precision x(n),fvec(n),fjac(ldfjac,n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec. do not alter fjac.
+c         if iflag = 2 calculate the jacobian at x and
+c         return this matrix in fjac. do not alter fvec.
+c         ---------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of hybrj.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length n which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output n by n array which contains the
+c         orthogonal matrix q produced by the qr factorization
+c         of the final approximate jacobian.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn with iflag = 1
+c         has reached maxfev.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.). 100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x and fvec available
+c         for printing. fvec and fjac should not be altered.
+c         if nprint is not positive, no special calls of fcn
+c         with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0   improper input parameters.
+c
+c         info = 1   relative error between two consecutive iterates
+c                    is at most xtol.
+c
+c         info = 2   number of calls to fcn with iflag = 1 has
+c                    reached maxfev.
+c
+c         info = 3   xtol is too small. no further improvement in
+c                    the approximate solution x is possible.
+c
+c         info = 4   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    five jacobian evaluations.
+c
+c         info = 5   iteration is not making good progress, as
+c                    measured by the improvement from the last
+c                    ten iterations.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn with iflag = 1.
+c
+c       njev is an integer output variable set to the number of
+c         calls to fcn with iflag = 2.
+c
+c       r is an output array of length lr which contains the
+c         upper triangular matrix produced by the qr factorization
+c         of the final approximate jacobian, stored rowwise.
+c
+c       lr is a positive integer input variable not less than
+c         (n*(n+1))/2.
+c
+c       qtf is an output array of length n which contains
+c         the vector (q transpose)*fvec.
+c
+c       wa1, wa2, wa3, and wa4 are work arrays of length n.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dogleg,dpmpar,enorm,
+c                            qform,qrfac,r1mpyq,r1updt
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,jm1,l,ncfail,ncsuc,nslow1,nslow2
+      integer iwa(1)
+      logical jeval,sing
+      double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
+     *                 prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
+     *                 zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p001,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+      njev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. ldfjac .lt. n .or. xtol .lt. zero
+     *    .or. maxfev .le. 0 .or. factor .le. zero
+     *    .or. lr .lt. (n*(n + 1))/2) go to 300
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 300
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(n,x,fvec,fjac,ldfjac,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 300
+      fnorm = enorm(n,fvec)
+c
+c     initialize iteration counter and monitors.
+c
+      iter = 1
+      ncsuc = 0
+      ncfail = 0
+      nslow1 = 0
+      nslow2 = 0
+c
+c     beginning of the outer loop.
+c
+   30 continue
+         jeval = .true.
+c
+c        calculate the jacobian matrix.
+c
+         iflag = 2
+         call fcn(n,x,fvec,fjac,ldfjac,iflag)
+         njev = njev + 1
+         if (iflag .lt. 0) go to 300
+c
+c        compute the qr factorization of the jacobian.
+c
+         call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 70
+         if (mode .eq. 2) go to 50
+         do 40 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+   40       continue
+   50    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 60 j = 1, n
+            wa3(j) = diag(j)*x(j)
+   60       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+   70    continue
+c
+c        form (q transpose)*fvec and store in qtf.
+c
+         do 80 i = 1, n
+            qtf(i) = fvec(i)
+   80       continue
+         do 120 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 110
+            sum = zero
+            do 90 i = j, n
+               sum = sum + fjac(i,j)*qtf(i)
+   90          continue
+            temp = -sum/fjac(j,j)
+            do 100 i = j, n
+               qtf(i) = qtf(i) + fjac(i,j)*temp
+  100          continue
+  110       continue
+  120       continue
+c
+c        copy the triangular factor of the qr factorization into r.
+c
+         sing = .false.
+         do 150 j = 1, n
+            l = j
+            jm1 = j - 1
+            if (jm1 .lt. 1) go to 140
+            do 130 i = 1, jm1
+               r(l) = fjac(i,j)
+               l = l + n - i
+  130          continue
+  140       continue
+            r(l) = wa1(j)
+            if (wa1(j) .eq. zero) sing = .true.
+  150       continue
+c
+c        accumulate the orthogonal factor in fjac.
+c
+         call qform(n,n,fjac,ldfjac,wa1)
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 170
+         do 160 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  160       continue
+  170    continue
+c
+c        beginning of the inner loop.
+c
+  180    continue
+c
+c           if requested, call fcn to enable printing of iterates.
+c
+            if (nprint .le. 0) go to 190
+            iflag = 0
+            if (mod(iter-1,nprint) .eq. 0)
+     *         call fcn(n,x,fvec,fjac,ldfjac,iflag)
+            if (iflag .lt. 0) go to 300
+  190       continue
+c
+c           determine the direction p.
+c
+            call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 200 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  200          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(n,wa2,wa4,fjac,ldfjac,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 300
+            fnorm1 = enorm(n,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction.
+c
+            l = 1
+            do 220 i = 1, n
+               sum = zero
+               do 210 j = i, n
+                  sum = sum + r(l)*wa1(j)
+                  l = l + 1
+  210             continue
+               wa3(i) = qtf(i) + sum
+  220          continue
+            temp = enorm(n,wa3)
+            prered = zero
+            if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .gt. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .ge. p1) go to 230
+               ncsuc = 0
+               ncfail = ncfail + 1
+               delta = p5*delta
+               go to 240
+  230       continue
+               ncfail = 0
+               ncsuc = ncsuc + 1
+               if (ratio .ge. p5 .or. ncsuc .gt. 1)
+     *            delta = dmax1(delta,pnorm/p5)
+               if (dabs(ratio-one) .le. p1) delta = pnorm/p5
+  240       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 260
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 250 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+               fvec(j) = wa4(j)
+  250          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  260       continue
+c
+c           determine the progress of the iteration.
+c
+            nslow1 = nslow1 + 1
+            if (actred .ge. p001) nslow1 = 0
+            if (jeval) nslow2 = nslow2 + 1
+            if (actred .ge. p1) nslow2 = 0
+c
+c           test for convergence.
+c
+            if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
+            if (info .ne. 0) go to 300
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 2
+            if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
+            if (nslow2 .eq. 5) info = 4
+            if (nslow1 .eq. 10) info = 5
+            if (info .ne. 0) go to 300
+c
+c           criterion for recalculating jacobian.
+c
+            if (ncfail .eq. 2) go to 290
+c
+c           calculate the rank one modification to the jacobian
+c           and update qtf if necessary.
+c
+            do 280 j = 1, n
+               sum = zero
+               do 270 i = 1, n
+                  sum = sum + fjac(i,j)*wa4(i)
+  270             continue
+               wa2(j) = (sum - wa3(j))/pnorm
+               wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
+               if (ratio .ge. p0001) qtf(j) = sum
+  280          continue
+c
+c           compute the qr factorization of the updated jacobian.
+c
+            call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
+            call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
+            call r1mpyq(1,n,qtf,1,wa2,wa3)
+c
+c           end of the inner loop.
+c
+            jeval = .false.
+            go to 180
+  290    continue
+c
+c        end of the outer loop.
+c
+         go to 30
+  300 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(n,x,fvec,fjac,ldfjac,iflag)
+      return
+c
+c     last card of subroutine hybrj.
+c
+      end
diff --git a/math/minpack/hybrj1.f b/math/minpack/hybrj1.f
new file mode 100644
index 00000000..9f51c496
--- /dev/null
+++ b/math/minpack/hybrj1.f
@@ -0,0 +1,127 @@
+      subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
+      integer n,ldfjac,info,lwa
+      double precision tol
+      double precision x(n),fvec(n),fjac(ldfjac,n),wa(lwa)
+      external fcn
+c     **********
+c
+c     subroutine hybrj1
+c
+c     the purpose of hybrj1 is to find a zero of a system of
+c     n nonlinear functions in n variables by a modification
+c     of the powell hybrid method. this is done by using the
+c     more general nonlinear equation solver hybrj. the user
+c     must provide a subroutine which calculates the functions
+c     and the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the jacobian. fcn must
+c         be declared in an external statement in the user
+c         calling program, and should be written as follows.
+c
+c         subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
+c         integer n,ldfjac,iflag
+c         double precision x(n),fvec(n),fjac(ldfjac,n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec. do not alter fjac.
+c         if iflag = 2 calculate the jacobian at x and
+c         return this matrix in fjac. do not alter fvec.
+c         ---------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of hybrj1.
+c         in this case set iflag to a negative integer.
+c
+c       n is a positive integer input variable set to the number
+c         of functions and variables.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length n which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output n by n array which contains the
+c         orthogonal matrix q produced by the qr factorization
+c         of the final approximate jacobian.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       tol is a nonnegative input variable. termination occurs
+c         when the algorithm estimates that the relative error
+c         between x and the solution is at most tol.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0   improper input parameters.
+c
+c         info = 1   algorithm estimates that the relative error
+c                    between x and the solution is at most tol.
+c
+c         info = 2   number of calls to fcn with iflag = 1 has
+c                    reached 100*(n+1).
+c
+c         info = 3   tol is too small. no further improvement in
+c                    the approximate solution x is possible.
+c
+c         info = 4   iteration is not making good progress.
+c
+c       wa is a work array of length lwa.
+c
+c       lwa is a positive integer input variable not less than
+c         (n*(n+13))/2.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... hybrj
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer j,lr,maxfev,mode,nfev,njev,nprint
+      double precision factor,one,xtol,zero
+      data factor,one,zero /1.0d2,1.0d0,0.0d0/
+      info = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. ldfjac .lt. n .or. tol .lt. zero
+     *    .or. lwa .lt. (n*(n + 13))/2) go to 20
+c
+c     call hybrj.
+c
+      maxfev = 100*(n + 1)
+      xtol = tol
+      mode = 2
+      do 10 j = 1, n
+         wa(j) = one
+   10    continue
+      nprint = 0
+      lr = (n*(n + 1))/2
+      call hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,wa(1),mode,
+     *           factor,nprint,info,nfev,njev,wa(6*n+1),lr,wa(n+1),
+     *           wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
+      if (info .eq. 5) info = 4
+   20 continue
+      return
+c
+c     last card of subroutine hybrj1.
+c
+      end
diff --git a/math/minpack/lmder.f b/math/minpack/lmder.f
new file mode 100644
index 00000000..8797d8be
--- /dev/null
+++ b/math/minpack/lmder.f
@@ -0,0 +1,452 @@
+      subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+     *                 maxfev,diag,mode,factor,nprint,info,nfev,njev,
+     *                 ipvt,qtf,wa1,wa2,wa3,wa4)
+      integer m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev
+      integer ipvt(n)
+      double precision ftol,xtol,gtol,factor
+      double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n),
+     *                 wa1(n),wa2(n),wa3(n),wa4(m)
+c     **********
+c
+c     subroutine lmder
+c
+c     the purpose of lmder is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of
+c     the levenberg-marquardt algorithm. the user must provide a
+c     subroutine which calculates the functions and the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+c                        maxfev,diag,mode,factor,nprint,info,nfev,
+c                        njev,ipvt,qtf,wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the jacobian. fcn must
+c         be declared in an external statement in the user
+c         calling program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+c         integer m,n,ldfjac,iflag
+c         double precision x(n),fvec(m),fjac(ldfjac,n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec. do not alter fjac.
+c         if iflag = 2 calculate the jacobian at x and
+c         return this matrix in fjac. do not alter fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmder.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output m by n array. the upper n by n submatrix
+c         of fjac contains an upper triangular matrix r with
+c         diagonal elements of nonincreasing magnitude such that
+c
+c                t     t           t
+c               p *(jac *jac)*p = r *r,
+c
+c         where p is a permutation matrix and jac is the final
+c         calculated jacobian. column j of p is column ipvt(j)
+c         (see below) of the identity matrix. the lower trapezoidal
+c         part of fjac contains information generated during
+c         the computation of r.
+c
+c       ldfjac is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array fjac.
+c
+c       ftol is a nonnegative input variable. termination
+c         occurs when both the actual and predicted relative
+c         reductions in the sum of squares are at most ftol.
+c         therefore, ftol measures the relative error desired
+c         in the sum of squares.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol. therefore, xtol measures the
+c         relative error desired in the approximate solution.
+c
+c       gtol is a nonnegative input variable. termination
+c         occurs when the cosine of the angle between fvec and
+c         any column of the jacobian is at most gtol in absolute
+c         value. therefore, gtol measures the orthogonality
+c         desired between the function vector and the columns
+c         of the jacobian.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn with iflag = 1
+c         has reached maxfev.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.).100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x, fvec, and fjac
+c         available for printing. fvec and fjac should not be
+c         altered. if nprint is not positive, no special calls
+c         of fcn with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  both actual and predicted relative reductions
+c                   in the sum of squares are at most ftol.
+c
+c         info = 2  relative error between two consecutive iterates
+c                   is at most xtol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  the cosine of the angle between fvec and any
+c                   column of the jacobian is at most gtol in
+c                   absolute value.
+c
+c         info = 5  number of calls to fcn with iflag = 1 has
+c                   reached maxfev.
+c
+c         info = 6  ftol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  xtol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c         info = 8  gtol is too small. fvec is orthogonal to the
+c                   columns of the jacobian to machine precision.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn with iflag = 1.
+c
+c       njev is an integer output variable set to the number of
+c         calls to fcn with iflag = 2.
+c
+c       ipvt is an integer output array of length n. ipvt
+c         defines a permutation matrix p such that jac*p = q*r,
+c         where jac is the final calculated jacobian, q is
+c         orthogonal (not stored), and r is upper triangular
+c         with diagonal elements of nonincreasing magnitude.
+c         column j of p is column ipvt(j) of the identity matrix.
+c
+c       qtf is an output array of length n which contains
+c         the first n elements of the vector (q transpose)*fvec.
+c
+c       wa1, wa2, and wa3 are work arrays of length n.
+c
+c       wa4 is a work array of length m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dpmpar,enorm,lmpar,qrfac
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,l
+      double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
+     *                 one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
+     *                 sum,temp,temp1,temp2,xnorm,zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p25,p75,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+      njev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
+     *    .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
+     *    .or. maxfev .le. 0 .or. factor .le. zero) go to 300
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 300
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 300
+      fnorm = enorm(m,fvec)
+c
+c     initialize levenberg-marquardt parameter and iteration counter.
+c
+      par = zero
+      iter = 1
+c
+c     beginning of the outer loop.
+c
+   30 continue
+c
+c        calculate the jacobian matrix.
+c
+         iflag = 2
+         call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+         njev = njev + 1
+         if (iflag .lt. 0) go to 300
+c
+c        if requested, call fcn to enable printing of iterates.
+c
+         if (nprint .le. 0) go to 40
+         iflag = 0
+         if (mod(iter-1,nprint) .eq. 0)
+     *      call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+         if (iflag .lt. 0) go to 300
+   40    continue
+c
+c        compute the qr factorization of the jacobian.
+c
+         call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 80
+         if (mode .eq. 2) go to 60
+         do 50 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+   50       continue
+   60    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 70 j = 1, n
+            wa3(j) = diag(j)*x(j)
+   70       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+   80    continue
+c
+c        form (q transpose)*fvec and store the first n components in
+c        qtf.
+c
+         do 90 i = 1, m
+            wa4(i) = fvec(i)
+   90       continue
+         do 130 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 120
+            sum = zero
+            do 100 i = j, m
+               sum = sum + fjac(i,j)*wa4(i)
+  100          continue
+            temp = -sum/fjac(j,j)
+            do 110 i = j, m
+               wa4(i) = wa4(i) + fjac(i,j)*temp
+  110          continue
+  120       continue
+            fjac(j,j) = wa1(j)
+            qtf(j) = wa4(j)
+  130       continue
+c
+c        compute the norm of the scaled gradient.
+c
+         gnorm = zero
+         if (fnorm .eq. zero) go to 170
+         do 160 j = 1, n
+            l = ipvt(j)
+            if (wa2(l) .eq. zero) go to 150
+            sum = zero
+            do 140 i = 1, j
+               sum = sum + fjac(i,j)*(qtf(i)/fnorm)
+  140          continue
+            gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
+  150       continue
+  160       continue
+  170    continue
+c
+c        test for convergence of the gradient norm.
+c
+         if (gnorm .le. gtol) info = 4
+         if (info .ne. 0) go to 300
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 190
+         do 180 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  180       continue
+  190    continue
+c
+c        beginning of the inner loop.
+c
+  200    continue
+c
+c           determine the levenberg-marquardt parameter.
+c
+            call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
+     *                 wa3,wa4)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 210 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  210          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(m,n,wa2,wa4,fjac,ldfjac,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 300
+            fnorm1 = enorm(m,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction and
+c           the scaled directional derivative.
+c
+            do 230 j = 1, n
+               wa3(j) = zero
+               l = ipvt(j)
+               temp = wa1(l)
+               do 220 i = 1, j
+                  wa3(i) = wa3(i) + fjac(i,j)*temp
+  220             continue
+  230          continue
+            temp1 = enorm(n,wa3)/fnorm
+            temp2 = (dsqrt(par)*pnorm)/fnorm
+            prered = temp1**2 + temp2**2/p5
+            dirder = -(temp1**2 + temp2**2)
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .ne. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .gt. p25) go to 240
+               if (actred .ge. zero) temp = p5
+               if (actred .lt. zero)
+     *            temp = p5*dirder/(dirder + p5*actred)
+               if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
+               delta = temp*dmin1(delta,pnorm/p1)
+               par = par/temp
+               go to 260
+  240       continue
+               if (par .ne. zero .and. ratio .lt. p75) go to 250
+               delta = pnorm/p5
+               par = p5*par
+  250          continue
+  260       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 290
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 270 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+  270          continue
+            do 280 i = 1, m
+               fvec(i) = wa4(i)
+  280          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  290       continue
+c
+c           tests for convergence.
+c
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one) info = 1
+            if (delta .le. xtol*xnorm) info = 2
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
+            if (info .ne. 0) go to 300
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 5
+            if (dabs(actred) .le. epsmch .and. prered .le. epsmch
+     *          .and. p5*ratio .le. one) info = 6
+            if (delta .le. epsmch*xnorm) info = 7
+            if (gnorm .le. epsmch) info = 8
+            if (info .ne. 0) go to 300
+c
+c           end of the inner loop. repeat if iteration unsuccessful.
+c
+            if (ratio .lt. p0001) go to 200
+c
+c        end of the outer loop.
+c
+         go to 30
+  300 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+      return
+c
+c     last card of subroutine lmder.
+c
+      end
diff --git a/math/minpack/lmder1.f b/math/minpack/lmder1.f
new file mode 100644
index 00000000..d691940f
--- /dev/null
+++ b/math/minpack/lmder1.f
@@ -0,0 +1,156 @@
+      subroutine lmder1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,ipvt,wa,
+     *                  lwa)
+      integer m,n,ldfjac,info,lwa
+      integer ipvt(n)
+      double precision tol
+      double precision x(n),fvec(m),fjac(ldfjac,n),wa(lwa)
+      external fcn
+c     **********
+c
+c     subroutine lmder1
+c
+c     the purpose of lmder1 is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of the
+c     levenberg-marquardt algorithm. this is done by using the more
+c     general least-squares solver lmder. the user must provide a
+c     subroutine which calculates the functions and the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine lmder1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,
+c                         ipvt,wa,lwa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the jacobian. fcn must
+c         be declared in an external statement in the user
+c         calling program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag)
+c         integer m,n,ldfjac,iflag
+c         double precision x(n),fvec(m),fjac(ldfjac,n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec. do not alter fjac.
+c         if iflag = 2 calculate the jacobian at x and
+c         return this matrix in fjac. do not alter fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmder1.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output m by n array. the upper n by n submatrix
+c         of fjac contains an upper triangular matrix r with
+c         diagonal elements of nonincreasing magnitude such that
+c
+c                t     t           t
+c               p *(jac *jac)*p = r *r,
+c
+c         where p is a permutation matrix and jac is the final
+c         calculated jacobian. column j of p is column ipvt(j)
+c         (see below) of the identity matrix. the lower trapezoidal
+c         part of fjac contains information generated during
+c         the computation of r.
+c
+c       ldfjac is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array fjac.
+c
+c       tol is a nonnegative input variable. termination occurs
+c         when the algorithm estimates either that the relative
+c         error in the sum of squares is at most tol or that
+c         the relative error between x and the solution is at
+c         most tol.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  algorithm estimates that the relative error
+c                   in the sum of squares is at most tol.
+c
+c         info = 2  algorithm estimates that the relative error
+c                   between x and the solution is at most tol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  fvec is orthogonal to the columns of the
+c                   jacobian to machine precision.
+c
+c         info = 5  number of calls to fcn with iflag = 1 has
+c                   reached 100*(n+1).
+c
+c         info = 6  tol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  tol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c       ipvt is an integer output array of length n. ipvt
+c         defines a permutation matrix p such that jac*p = q*r,
+c         where jac is the final calculated jacobian, q is
+c         orthogonal (not stored), and r is upper triangular
+c         with diagonal elements of nonincreasing magnitude.
+c         column j of p is column ipvt(j) of the identity matrix.
+c
+c       wa is a work array of length lwa.
+c
+c       lwa is a positive integer input variable not less than 5*n+m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... lmder
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer maxfev,mode,nfev,njev,nprint
+      double precision factor,ftol,gtol,xtol,zero
+      data factor,zero /1.0d2,0.0d0/
+      info = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m .or. tol .lt. zero
+     *    .or. lwa .lt. 5*n + m) go to 10
+c
+c     call lmder.
+c
+      maxfev = 100*(n + 1)
+      ftol = tol
+      xtol = tol
+      gtol = zero
+      mode = 1
+      nprint = 0
+      call lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,
+     *           wa(1),mode,factor,nprint,info,nfev,njev,ipvt,wa(n+1),
+     *           wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
+      if (info .eq. 8) info = 4
+   10 continue
+      return
+c
+c     last card of subroutine lmder1.
+c
+      end
diff --git a/math/minpack/lmdif.f b/math/minpack/lmdif.f
new file mode 100644
index 00000000..dd3d4ee2
--- /dev/null
+++ b/math/minpack/lmdif.f
@@ -0,0 +1,454 @@
+      subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
+     *                 diag,mode,factor,nprint,info,nfev,fjac,ldfjac,
+     *                 ipvt,qtf,wa1,wa2,wa3,wa4)
+      integer m,n,maxfev,mode,nprint,info,nfev,ldfjac
+      integer ipvt(n)
+      double precision ftol,xtol,gtol,epsfcn,factor
+      double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n),
+     *                 wa1(n),wa2(n),wa3(n),wa4(m)
+      external fcn
+c     **********
+c
+c     subroutine lmdif
+c
+c     the purpose of lmdif is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of
+c     the levenberg-marquardt algorithm. the user must provide a
+c     subroutine which calculates the functions. the jacobian is
+c     then calculated by a forward-difference approximation.
+c
+c     the subroutine statement is
+c
+c       subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
+c                        diag,mode,factor,nprint,info,nfev,fjac,
+c                        ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,iflag)
+c         integer m,n,iflag
+c         double precision x(n),fvec(m)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmdif.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       ftol is a nonnegative input variable. termination
+c         occurs when both the actual and predicted relative
+c         reductions in the sum of squares are at most ftol.
+c         therefore, ftol measures the relative error desired
+c         in the sum of squares.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol. therefore, xtol measures the
+c         relative error desired in the approximate solution.
+c
+c       gtol is a nonnegative input variable. termination
+c         occurs when the cosine of the angle between fvec and
+c         any column of the jacobian is at most gtol in absolute
+c         value. therefore, gtol measures the orthogonality
+c         desired between the function vector and the columns
+c         of the jacobian.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn is at least
+c         maxfev by the end of an iteration.
+c
+c       epsfcn is an input variable used in determining a suitable
+c         step length for the forward-difference approximation. this
+c         approximation assumes that the relative errors in the
+c         functions are of the order of epsfcn. if epsfcn is less
+c         than the machine precision, it is assumed that the relative
+c         errors in the functions are of the order of the machine
+c         precision.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.). 100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x and fvec available
+c         for printing. if nprint is not positive, no special calls
+c         of fcn with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  both actual and predicted relative reductions
+c                   in the sum of squares are at most ftol.
+c
+c         info = 2  relative error between two consecutive iterates
+c                   is at most xtol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  the cosine of the angle between fvec and any
+c                   column of the jacobian is at most gtol in
+c                   absolute value.
+c
+c         info = 5  number of calls to fcn has reached or
+c                   exceeded maxfev.
+c
+c         info = 6  ftol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  xtol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c         info = 8  gtol is too small. fvec is orthogonal to the
+c                   columns of the jacobian to machine precision.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn.
+c
+c       fjac is an output m by n array. the upper n by n submatrix
+c         of fjac contains an upper triangular matrix r with
+c         diagonal elements of nonincreasing magnitude such that
+c
+c                t     t           t
+c               p *(jac *jac)*p = r *r,
+c
+c         where p is a permutation matrix and jac is the final
+c         calculated jacobian. column j of p is column ipvt(j)
+c         (see below) of the identity matrix. the lower trapezoidal
+c         part of fjac contains information generated during
+c         the computation of r.
+c
+c       ldfjac is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array fjac.
+c
+c       ipvt is an integer output array of length n. ipvt
+c         defines a permutation matrix p such that jac*p = q*r,
+c         where jac is the final calculated jacobian, q is
+c         orthogonal (not stored), and r is upper triangular
+c         with diagonal elements of nonincreasing magnitude.
+c         column j of p is column ipvt(j) of the identity matrix.
+c
+c       qtf is an output array of length n which contains
+c         the first n elements of the vector (q transpose)*fvec.
+c
+c       wa1, wa2, and wa3 are work arrays of length n.
+c
+c       wa4 is a work array of length m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,l
+      double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
+     *                 one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
+     *                 sum,temp,temp1,temp2,xnorm,zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p25,p75,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
+     *    .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
+     *    .or. maxfev .le. 0 .or. factor .le. zero) go to 300
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 300
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(m,n,x,fvec,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 300
+      fnorm = enorm(m,fvec)
+c
+c     initialize levenberg-marquardt parameter and iteration counter.
+c
+      par = zero
+      iter = 1
+c
+c     beginning of the outer loop.
+c
+   30 continue
+c
+c        calculate the jacobian matrix.
+c
+         iflag = 2
+         call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
+         nfev = nfev + n
+         if (iflag .lt. 0) go to 300
+c
+c        if requested, call fcn to enable printing of iterates.
+c
+         if (nprint .le. 0) go to 40
+         iflag = 0
+         if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag)
+         if (iflag .lt. 0) go to 300
+   40    continue
+c
+c        compute the qr factorization of the jacobian.
+c
+         call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 80
+         if (mode .eq. 2) go to 60
+         do 50 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+   50       continue
+   60    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 70 j = 1, n
+            wa3(j) = diag(j)*x(j)
+   70       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+   80    continue
+c
+c        form (q transpose)*fvec and store the first n components in
+c        qtf.
+c
+         do 90 i = 1, m
+            wa4(i) = fvec(i)
+   90       continue
+         do 130 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 120
+            sum = zero
+            do 100 i = j, m
+               sum = sum + fjac(i,j)*wa4(i)
+  100          continue
+            temp = -sum/fjac(j,j)
+            do 110 i = j, m
+               wa4(i) = wa4(i) + fjac(i,j)*temp
+  110          continue
+  120       continue
+            fjac(j,j) = wa1(j)
+            qtf(j) = wa4(j)
+  130       continue
+c
+c        compute the norm of the scaled gradient.
+c
+         gnorm = zero
+         if (fnorm .eq. zero) go to 170
+         do 160 j = 1, n
+            l = ipvt(j)
+            if (wa2(l) .eq. zero) go to 150
+            sum = zero
+            do 140 i = 1, j
+               sum = sum + fjac(i,j)*(qtf(i)/fnorm)
+  140          continue
+            gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
+  150       continue
+  160       continue
+  170    continue
+c
+c        test for convergence of the gradient norm.
+c
+         if (gnorm .le. gtol) info = 4
+         if (info .ne. 0) go to 300
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 190
+         do 180 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  180       continue
+  190    continue
+c
+c        beginning of the inner loop.
+c
+  200    continue
+c
+c           determine the levenberg-marquardt parameter.
+c
+            call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
+     *                 wa3,wa4)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 210 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  210          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(m,n,wa2,wa4,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 300
+            fnorm1 = enorm(m,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction and
+c           the scaled directional derivative.
+c
+            do 230 j = 1, n
+               wa3(j) = zero
+               l = ipvt(j)
+               temp = wa1(l)
+               do 220 i = 1, j
+                  wa3(i) = wa3(i) + fjac(i,j)*temp
+  220             continue
+  230          continue
+            temp1 = enorm(n,wa3)/fnorm
+            temp2 = (dsqrt(par)*pnorm)/fnorm
+            prered = temp1**2 + temp2**2/p5
+            dirder = -(temp1**2 + temp2**2)
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .ne. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .gt. p25) go to 240
+               if (actred .ge. zero) temp = p5
+               if (actred .lt. zero)
+     *            temp = p5*dirder/(dirder + p5*actred)
+               if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
+               delta = temp*dmin1(delta,pnorm/p1)
+               par = par/temp
+               go to 260
+  240       continue
+               if (par .ne. zero .and. ratio .lt. p75) go to 250
+               delta = pnorm/p5
+               par = p5*par
+  250          continue
+  260       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 290
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 270 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+  270          continue
+            do 280 i = 1, m
+               fvec(i) = wa4(i)
+  280          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  290       continue
+c
+c           tests for convergence.
+c
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one) info = 1
+            if (delta .le. xtol*xnorm) info = 2
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
+            if (info .ne. 0) go to 300
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 5
+            if (dabs(actred) .le. epsmch .and. prered .le. epsmch
+     *          .and. p5*ratio .le. one) info = 6
+            if (delta .le. epsmch*xnorm) info = 7
+            if (gnorm .le. epsmch) info = 8
+            if (info .ne. 0) go to 300
+c
+c           end of the inner loop. repeat if iteration unsuccessful.
+c
+            if (ratio .lt. p0001) go to 200
+c
+c        end of the outer loop.
+c
+         go to 30
+  300 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
+      return
+c
+c     last card of subroutine lmdif.
+c
+      end
diff --git a/math/minpack/lmdif1.f b/math/minpack/lmdif1.f
new file mode 100644
index 00000000..70f8aae0
--- /dev/null
+++ b/math/minpack/lmdif1.f
@@ -0,0 +1,135 @@
+      subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
+      integer m,n,info,lwa
+      integer iwa(n)
+      double precision tol
+      double precision x(n),fvec(m),wa(lwa)
+      external fcn
+c     **********
+c
+c     subroutine lmdif1
+c
+c     the purpose of lmdif1 is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of the
+c     levenberg-marquardt algorithm. this is done by using the more
+c     general least-squares solver lmdif. the user must provide a
+c     subroutine which calculates the functions. the jacobian is
+c     then calculated by a forward-difference approximation.
+c
+c     the subroutine statement is
+c
+c       subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions. fcn must be declared
+c         in an external statement in the user calling
+c         program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,iflag)
+c         integer m,n,iflag
+c         double precision x(n),fvec(m)
+c         ----------
+c         calculate the functions at x and
+c         return this vector in fvec.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmdif1.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       tol is a nonnegative input variable. termination occurs
+c         when the algorithm estimates either that the relative
+c         error in the sum of squares is at most tol or that
+c         the relative error between x and the solution is at
+c         most tol.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  algorithm estimates that the relative error
+c                   in the sum of squares is at most tol.
+c
+c         info = 2  algorithm estimates that the relative error
+c                   between x and the solution is at most tol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  fvec is orthogonal to the columns of the
+c                   jacobian to machine precision.
+c
+c         info = 5  number of calls to fcn has reached or
+c                   exceeded 200*(n+1).
+c
+c         info = 6  tol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  tol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c       iwa is an integer work array of length n.
+c
+c       wa is a work array of length lwa.
+c
+c       lwa is a positive integer input variable not less than
+c         m*n+5*n+m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... lmdif
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer maxfev,mode,mp5n,nfev,nprint
+      double precision epsfcn,factor,ftol,gtol,xtol,zero
+      data factor,zero /1.0d2,0.0d0/
+      info = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. tol .lt. zero
+     *    .or. lwa .lt. m*n + 5*n + m) go to 10
+c
+c     call lmdif.
+c
+      maxfev = 200*(n + 1)
+      ftol = tol
+      xtol = tol
+      gtol = zero
+      epsfcn = zero
+      mode = 1
+      nprint = 0
+      mp5n = m + 5*n
+      call lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,wa(1),
+     *           mode,factor,nprint,info,nfev,wa(mp5n+1),m,iwa,
+     *           wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
+      if (info .eq. 8) info = 4
+   10 continue
+      return
+c
+c     last card of subroutine lmdif1.
+c
+      end
diff --git a/math/minpack/lmpar.f b/math/minpack/lmpar.f
new file mode 100644
index 00000000..26c422a7
--- /dev/null
+++ b/math/minpack/lmpar.f
@@ -0,0 +1,264 @@
+      subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1,
+     *                 wa2)
+      integer n,ldr
+      integer ipvt(n)
+      double precision delta,par
+      double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n),
+     *                 wa2(n)
+c     **********
+c
+c     subroutine lmpar
+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,
+c     the problem is to determine a value for the parameter
+c     par such that if x solves the system
+c
+c           a*x = b ,     sqrt(par)*d*x = 0 ,
+c
+c     in the least squares sense, and dxnorm is the euclidean
+c     norm of d*x, then either par is zero and
+c
+c           (dxnorm-delta) .le. 0.1*delta ,
+c
+c     or par is positive and
+c
+c           abs(dxnorm-delta) .le. 0.1*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, with column pivoting, of a. that is, if
+c     a*p = q*r, where p is a permutation matrix, q has orthogonal
+c     columns, and r is an upper triangular matrix with diagonal
+c     elements of nonincreasing magnitude, then lmpar expects
+c     the full upper triangle of r, the permutation matrix p,
+c     and the first n components of (q transpose)*b. on output
+c     lmpar also provides an upper triangular matrix s such that
+c
+c            t   t                   t
+c           p *(a *a + par*d*d)*p = s *s .
+c
+c     s is employed within lmpar and may be of separate interest.
+c
+c     only a few iterations are generally needed for convergence
+c     of the algorithm. if, however, the limit of 10 iterations
+c     is reached, then the output par will contain the best
+c     value obtained so far.
+c
+c     the subroutine statement is
+c
+c       subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,
+c                        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 n by n array. on input the full upper triangle
+c         must contain the full upper triangle of the matrix r.
+c         on output the full upper triangle is unaltered, and the
+c         strict lower triangle contains the strict upper triangle
+c         (transposed) of the upper triangular matrix s.
+c
+c       ldr is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array r.
+c
+c       ipvt is an integer input array of length n which defines the
+c         permutation matrix p such that a*p = q*r. column j of p
+c         is column ipvt(j) of the identity matrix.
+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       par is a nonnegative variable. on input par contains an
+c         initial estimate of the levenberg-marquardt parameter.
+c         on output par contains the final estimate.
+c
+c       x is an output array of length n which contains the least
+c         squares solution of the system a*x = b, sqrt(par)*d*x = 0,
+c         for the output par.
+c
+c       sdiag is an output array of length n which contains the
+c         diagonal elements of the upper triangular matrix s.
+c
+c       wa1 and wa2 are work arrays of length n.
+c
+c     subprograms called
+c
+c       minpack-supplied ... dpmpar,enorm,qrsolv
+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,iter,j,jm1,jp1,k,l,nsing
+      double precision dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001,
+     *                 sum,temp,zero
+      double precision dpmpar,enorm
+      data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/
+c
+c     dwarf is the smallest positive magnitude.
+c
+      dwarf = dpmpar(2)
+c
+c     compute and store in x the gauss-newton direction. if the
+c     jacobian is rank-deficient, obtain a least squares solution.
+c
+      nsing = n
+      do 10 j = 1, n
+         wa1(j) = qtb(j)
+         if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
+         if (nsing .lt. n) wa1(j) = zero
+   10    continue
+      if (nsing .lt. 1) go to 50
+      do 40 k = 1, nsing
+         j = nsing - k + 1
+         wa1(j) = wa1(j)/r(j,j)
+         temp = wa1(j)
+         jm1 = j - 1
+         if (jm1 .lt. 1) go to 30
+         do 20 i = 1, jm1
+            wa1(i) = wa1(i) - r(i,j)*temp
+   20       continue
+   30    continue
+   40    continue
+   50 continue
+      do 60 j = 1, n
+         l = ipvt(j)
+         x(l) = wa1(j)
+   60    continue
+c
+c     initialize the iteration counter.
+c     evaluate the function at the origin, and test
+c     for acceptance of the gauss-newton direction.
+c
+      iter = 0
+      do 70 j = 1, n
+         wa2(j) = diag(j)*x(j)
+   70    continue
+      dxnorm = enorm(n,wa2)
+      fp = dxnorm - delta
+      if (fp .le. p1*delta) go to 220
+c
+c     if the jacobian is not rank deficient, the newton
+c     step provides a lower bound, parl, for the zero of
+c     the function. otherwise set this bound to zero.
+c
+      parl = zero
+      if (nsing .lt. n) go to 120
+      do 80 j = 1, n
+         l = ipvt(j)
+         wa1(j) = diag(l)*(wa2(l)/dxnorm)
+   80    continue
+      do 110 j = 1, n
+         sum = zero
+         jm1 = j - 1
+         if (jm1 .lt. 1) go to 100
+         do 90 i = 1, jm1
+            sum = sum + r(i,j)*wa1(i)
+   90       continue
+  100    continue
+         wa1(j) = (wa1(j) - sum)/r(j,j)
+  110    continue
+      temp = enorm(n,wa1)
+      parl = ((fp/delta)/temp)/temp
+  120 continue
+c
+c     calculate an upper bound, paru, for the zero of the function.
+c
+      do 140 j = 1, n
+         sum = zero
+         do 130 i = 1, j
+            sum = sum + r(i,j)*qtb(i)
+  130       continue
+         l = ipvt(j)
+         wa1(j) = sum/diag(l)
+  140    continue
+      gnorm = enorm(n,wa1)
+      paru = gnorm/delta
+      if (paru .eq. zero) paru = dwarf/dmin1(delta,p1)
+c
+c     if the input par lies outside of the interval (parl,paru),
+c     set par to the closer endpoint.
+c
+      par = dmax1(par,parl)
+      par = dmin1(par,paru)
+      if (par .eq. zero) par = gnorm/dxnorm
+c
+c     beginning of an iteration.
+c
+  150 continue
+         iter = iter + 1
+c
+c        evaluate the function at the current value of par.
+c
+         if (par .eq. zero) par = dmax1(dwarf,p001*paru)
+         temp = dsqrt(par)
+         do 160 j = 1, n
+            wa1(j) = temp*diag(j)
+  160       continue
+         call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2)
+         do 170 j = 1, n
+            wa2(j) = diag(j)*x(j)
+  170       continue
+         dxnorm = enorm(n,wa2)
+         temp = fp
+         fp = dxnorm - delta
+c
+c        if the function is small enough, accept the current value
+c        of par. also test for the exceptional cases where parl
+c        is zero or the number of iterations has reached 10.
+c
+         if (dabs(fp) .le. p1*delta
+     *       .or. parl .eq. zero .and. fp .le. temp
+     *            .and. temp .lt. zero .or. iter .eq. 10) go to 220
+c
+c        compute the newton correction.
+c
+         do 180 j = 1, n
+            l = ipvt(j)
+            wa1(j) = diag(l)*(wa2(l)/dxnorm)
+  180       continue
+         do 210 j = 1, n
+            wa1(j) = wa1(j)/sdiag(j)
+            temp = wa1(j)
+            jp1 = j + 1
+            if (n .lt. jp1) go to 200
+            do 190 i = jp1, n
+               wa1(i) = wa1(i) - r(i,j)*temp
+  190          continue
+  200       continue
+  210       continue
+         temp = enorm(n,wa1)
+         parc = ((fp/delta)/temp)/temp
+c
+c        depending on the sign of the function, update parl or paru.
+c
+         if (fp .gt. zero) parl = dmax1(parl,par)
+         if (fp .lt. zero) paru = dmin1(paru,par)
+c
+c        compute an improved estimate for par.
+c
+         par = dmax1(parl,par+parc)
+c
+c        end of an iteration.
+c
+         go to 150
+  220 continue
+c
+c     termination.
+c
+      if (iter .eq. 0) par = zero
+      return
+c
+c     last card of subroutine lmpar.
+c
+      end
diff --git a/math/minpack/lmstr.f b/math/minpack/lmstr.f
new file mode 100644
index 00000000..d9a7893f
--- /dev/null
+++ b/math/minpack/lmstr.f
@@ -0,0 +1,466 @@
+      subroutine lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+     *                 maxfev,diag,mode,factor,nprint,info,nfev,njev,
+     *                 ipvt,qtf,wa1,wa2,wa3,wa4)
+      integer m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev
+      integer ipvt(n)
+      logical sing
+      double precision ftol,xtol,gtol,factor
+      double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n),
+     *                 wa1(n),wa2(n),wa3(n),wa4(m)
+c     **********
+c
+c     subroutine lmstr
+c
+c     the purpose of lmstr is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of
+c     the levenberg-marquardt algorithm which uses minimal storage.
+c     the user must provide a subroutine which calculates the
+c     functions and the rows of the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
+c                        maxfev,diag,mode,factor,nprint,info,nfev,
+c                        njev,ipvt,qtf,wa1,wa2,wa3,wa4)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the rows of the jacobian.
+c         fcn must be declared in an external statement in the
+c         user calling program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,fjrow,iflag)
+c         integer m,n,iflag
+c         double precision x(n),fvec(m),fjrow(n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec.
+c         if iflag = i calculate the (i-1)-st row of the
+c         jacobian at x and return this vector in fjrow.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmstr.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output n by n array. the upper triangle of fjac
+c         contains an upper triangular matrix r such that
+c
+c                t     t           t
+c               p *(jac *jac)*p = r *r,
+c
+c         where p is a permutation matrix and jac is the final
+c         calculated jacobian. column j of p is column ipvt(j)
+c         (see below) of the identity matrix. the lower triangular
+c         part of fjac contains information generated during
+c         the computation of r.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       ftol is a nonnegative input variable. termination
+c         occurs when both the actual and predicted relative
+c         reductions in the sum of squares are at most ftol.
+c         therefore, ftol measures the relative error desired
+c         in the sum of squares.
+c
+c       xtol is a nonnegative input variable. termination
+c         occurs when the relative error between two consecutive
+c         iterates is at most xtol. therefore, xtol measures the
+c         relative error desired in the approximate solution.
+c
+c       gtol is a nonnegative input variable. termination
+c         occurs when the cosine of the angle between fvec and
+c         any column of the jacobian is at most gtol in absolute
+c         value. therefore, gtol measures the orthogonality
+c         desired between the function vector and the columns
+c         of the jacobian.
+c
+c       maxfev is a positive integer input variable. termination
+c         occurs when the number of calls to fcn with iflag = 1
+c         has reached maxfev.
+c
+c       diag is an array of length n. if mode = 1 (see
+c         below), diag is internally set. if mode = 2, diag
+c         must contain positive entries that serve as
+c         multiplicative scale factors for the variables.
+c
+c       mode is an integer input variable. if mode = 1, the
+c         variables will be scaled internally. if mode = 2,
+c         the scaling is specified by the input diag. other
+c         values of mode are equivalent to mode = 1.
+c
+c       factor is a positive input variable used in determining the
+c         initial step bound. this bound is set to the product of
+c         factor and the euclidean norm of diag*x if nonzero, or else
+c         to factor itself. in most cases factor should lie in the
+c         interval (.1,100.). 100. is a generally recommended value.
+c
+c       nprint is an integer input variable that enables controlled
+c         printing of iterates if it is positive. in this case,
+c         fcn is called with iflag = 0 at the beginning of the first
+c         iteration and every nprint iterations thereafter and
+c         immediately prior to return, with x and fvec available
+c         for printing. if nprint is not positive, no special calls
+c         of fcn with iflag = 0 are made.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  both actual and predicted relative reductions
+c                   in the sum of squares are at most ftol.
+c
+c         info = 2  relative error between two consecutive iterates
+c                   is at most xtol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  the cosine of the angle between fvec and any
+c                   column of the jacobian is at most gtol in
+c                   absolute value.
+c
+c         info = 5  number of calls to fcn with iflag = 1 has
+c                   reached maxfev.
+c
+c         info = 6  ftol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  xtol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c         info = 8  gtol is too small. fvec is orthogonal to the
+c                   columns of the jacobian to machine precision.
+c
+c       nfev is an integer output variable set to the number of
+c         calls to fcn with iflag = 1.
+c
+c       njev is an integer output variable set to the number of
+c         calls to fcn with iflag = 2.
+c
+c       ipvt is an integer output array of length n. ipvt
+c         defines a permutation matrix p such that jac*p = q*r,
+c         where jac is the final calculated jacobian, q is
+c         orthogonal (not stored), and r is upper triangular.
+c         column j of p is column ipvt(j) of the identity matrix.
+c
+c       qtf is an output array of length n which contains
+c         the first n elements of the vector (q transpose)*fvec.
+c
+c       wa1, wa2, and wa3 are work arrays of length n.
+c
+c       wa4 is a work array of length m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... dpmpar,enorm,lmpar,qrfac,rwupdt
+c
+c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
+c     jorge j. more
+c
+c     **********
+      integer i,iflag,iter,j,l
+      double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
+     *                 one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
+     *                 sum,temp,temp1,temp2,xnorm,zero
+      double precision dpmpar,enorm
+      data one,p1,p5,p25,p75,p0001,zero
+     *     /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+      info = 0
+      iflag = 0
+      nfev = 0
+      njev = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. n
+     *    .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
+     *    .or. maxfev .le. 0 .or. factor .le. zero) go to 340
+      if (mode .ne. 2) go to 20
+      do 10 j = 1, n
+         if (diag(j) .le. zero) go to 340
+   10    continue
+   20 continue
+c
+c     evaluate the function at the starting point
+c     and calculate its norm.
+c
+      iflag = 1
+      call fcn(m,n,x,fvec,wa3,iflag)
+      nfev = 1
+      if (iflag .lt. 0) go to 340
+      fnorm = enorm(m,fvec)
+c
+c     initialize levenberg-marquardt parameter and iteration counter.
+c
+      par = zero
+      iter = 1
+c
+c     beginning of the outer loop.
+c
+   30 continue
+c
+c        if requested, call fcn to enable printing of iterates.
+c
+         if (nprint .le. 0) go to 40
+         iflag = 0
+         if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,wa3,iflag)
+         if (iflag .lt. 0) go to 340
+   40    continue
+c
+c        compute the qr factorization of the jacobian matrix
+c        calculated one row at a time, while simultaneously
+c        forming (q transpose)*fvec and storing the first
+c        n components in qtf.
+c
+         do 60 j = 1, n
+            qtf(j) = zero
+            do 50 i = 1, n
+               fjac(i,j) = zero
+   50          continue
+   60       continue
+         iflag = 2
+         do 70 i = 1, m
+            call fcn(m,n,x,fvec,wa3,iflag)
+            if (iflag .lt. 0) go to 340
+            temp = fvec(i)
+            call rwupdt(n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2)
+            iflag = iflag + 1
+   70       continue
+         njev = njev + 1
+c
+c        if the jacobian is rank deficient, call qrfac to
+c        reorder its columns and update the components of qtf.
+c
+         sing = .false.
+         do 80 j = 1, n
+            if (fjac(j,j) .eq. zero) sing = .true.
+            ipvt(j) = j
+            wa2(j) = enorm(j,fjac(1,j))
+   80       continue
+         if (.not.sing) go to 130
+         call qrfac(n,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
+         do 120 j = 1, n
+            if (fjac(j,j) .eq. zero) go to 110
+            sum = zero
+            do 90 i = j, n
+               sum = sum + fjac(i,j)*qtf(i)
+   90          continue
+            temp = -sum/fjac(j,j)
+            do 100 i = j, n
+               qtf(i) = qtf(i) + fjac(i,j)*temp
+  100          continue
+  110       continue
+            fjac(j,j) = wa1(j)
+  120       continue
+  130    continue
+c
+c        on the first iteration and if mode is 1, scale according
+c        to the norms of the columns of the initial jacobian.
+c
+         if (iter .ne. 1) go to 170
+         if (mode .eq. 2) go to 150
+         do 140 j = 1, n
+            diag(j) = wa2(j)
+            if (wa2(j) .eq. zero) diag(j) = one
+  140       continue
+  150    continue
+c
+c        on the first iteration, calculate the norm of the scaled x
+c        and initialize the step bound delta.
+c
+         do 160 j = 1, n
+            wa3(j) = diag(j)*x(j)
+  160       continue
+         xnorm = enorm(n,wa3)
+         delta = factor*xnorm
+         if (delta .eq. zero) delta = factor
+  170    continue
+c
+c        compute the norm of the scaled gradient.
+c
+         gnorm = zero
+         if (fnorm .eq. zero) go to 210
+         do 200 j = 1, n
+            l = ipvt(j)
+            if (wa2(l) .eq. zero) go to 190
+            sum = zero
+            do 180 i = 1, j
+               sum = sum + fjac(i,j)*(qtf(i)/fnorm)
+  180          continue
+            gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
+  190       continue
+  200       continue
+  210    continue
+c
+c        test for convergence of the gradient norm.
+c
+         if (gnorm .le. gtol) info = 4
+         if (info .ne. 0) go to 340
+c
+c        rescale if necessary.
+c
+         if (mode .eq. 2) go to 230
+         do 220 j = 1, n
+            diag(j) = dmax1(diag(j),wa2(j))
+  220       continue
+  230    continue
+c
+c        beginning of the inner loop.
+c
+  240    continue
+c
+c           determine the levenberg-marquardt parameter.
+c
+            call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
+     *                 wa3,wa4)
+c
+c           store the direction p and x + p. calculate the norm of p.
+c
+            do 250 j = 1, n
+               wa1(j) = -wa1(j)
+               wa2(j) = x(j) + wa1(j)
+               wa3(j) = diag(j)*wa1(j)
+  250          continue
+            pnorm = enorm(n,wa3)
+c
+c           on the first iteration, adjust the initial step bound.
+c
+            if (iter .eq. 1) delta = dmin1(delta,pnorm)
+c
+c           evaluate the function at x + p and calculate its norm.
+c
+            iflag = 1
+            call fcn(m,n,wa2,wa4,wa3,iflag)
+            nfev = nfev + 1
+            if (iflag .lt. 0) go to 340
+            fnorm1 = enorm(m,wa4)
+c
+c           compute the scaled actual reduction.
+c
+            actred = -one
+            if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
+c
+c           compute the scaled predicted reduction and
+c           the scaled directional derivative.
+c
+            do 270 j = 1, n
+               wa3(j) = zero
+               l = ipvt(j)
+               temp = wa1(l)
+               do 260 i = 1, j
+                  wa3(i) = wa3(i) + fjac(i,j)*temp
+  260             continue
+  270          continue
+            temp1 = enorm(n,wa3)/fnorm
+            temp2 = (dsqrt(par)*pnorm)/fnorm
+            prered = temp1**2 + temp2**2/p5
+            dirder = -(temp1**2 + temp2**2)
+c
+c           compute the ratio of the actual to the predicted
+c           reduction.
+c
+            ratio = zero
+            if (prered .ne. zero) ratio = actred/prered
+c
+c           update the step bound.
+c
+            if (ratio .gt. p25) go to 280
+               if (actred .ge. zero) temp = p5
+               if (actred .lt. zero)
+     *            temp = p5*dirder/(dirder + p5*actred)
+               if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
+               delta = temp*dmin1(delta,pnorm/p1)
+               par = par/temp
+               go to 300
+  280       continue
+               if (par .ne. zero .and. ratio .lt. p75) go to 290
+               delta = pnorm/p5
+               par = p5*par
+  290          continue
+  300       continue
+c
+c           test for successful iteration.
+c
+            if (ratio .lt. p0001) go to 330
+c
+c           successful iteration. update x, fvec, and their norms.
+c
+            do 310 j = 1, n
+               x(j) = wa2(j)
+               wa2(j) = diag(j)*x(j)
+  310          continue
+            do 320 i = 1, m
+               fvec(i) = wa4(i)
+  320          continue
+            xnorm = enorm(n,wa2)
+            fnorm = fnorm1
+            iter = iter + 1
+  330       continue
+c
+c           tests for convergence.
+c
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one) info = 1
+            if (delta .le. xtol*xnorm) info = 2
+            if (dabs(actred) .le. ftol .and. prered .le. ftol
+     *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
+            if (info .ne. 0) go to 340
+c
+c           tests for termination and stringent tolerances.
+c
+            if (nfev .ge. maxfev) info = 5
+            if (dabs(actred) .le. epsmch .and. prered .le. epsmch
+     *          .and. p5*ratio .le. one) info = 6
+            if (delta .le. epsmch*xnorm) info = 7
+            if (gnorm .le. epsmch) info = 8
+            if (info .ne. 0) go to 340
+c
+c           end of the inner loop. repeat if iteration unsuccessful.
+c
+            if (ratio .lt. p0001) go to 240
+c
+c        end of the outer loop.
+c
+         go to 30
+  340 continue
+c
+c     termination, either normal or user imposed.
+c
+      if (iflag .lt. 0) info = iflag
+      iflag = 0
+      if (nprint .gt. 0) call fcn(m,n,x,fvec,wa3,iflag)
+      return
+c
+c     last card of subroutine lmstr.
+c
+      end
diff --git a/math/minpack/lmstr1.f b/math/minpack/lmstr1.f
new file mode 100644
index 00000000..2fa8ee1c
--- /dev/null
+++ b/math/minpack/lmstr1.f
@@ -0,0 +1,156 @@
+      subroutine lmstr1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,ipvt,wa,
+     *                  lwa)
+      integer m,n,ldfjac,info,lwa
+      integer ipvt(n)
+      double precision tol
+      double precision x(n),fvec(m),fjac(ldfjac,n),wa(lwa)
+      external fcn
+c     **********
+c
+c     subroutine lmstr1
+c
+c     the purpose of lmstr1 is to minimize the sum of the squares of
+c     m nonlinear functions in n variables by a modification of
+c     the levenberg-marquardt algorithm which uses minimal storage.
+c     this is done by using the more general least-squares solver
+c     lmstr. the user must provide a subroutine which calculates
+c     the functions and the rows of the jacobian.
+c
+c     the subroutine statement is
+c
+c       subroutine lmstr1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,
+c                         ipvt,wa,lwa)
+c
+c     where
+c
+c       fcn is the name of the user-supplied subroutine which
+c         calculates the functions and the rows of the jacobian.
+c         fcn must be declared in an external statement in the
+c         user calling program, and should be written as follows.
+c
+c         subroutine fcn(m,n,x,fvec,fjrow,iflag)
+c         integer m,n,iflag
+c         double precision x(n),fvec(m),fjrow(n)
+c         ----------
+c         if iflag = 1 calculate the functions at x and
+c         return this vector in fvec.
+c         if iflag = i calculate the (i-1)-st row of the
+c         jacobian at x and return this vector in fjrow.
+c         ----------
+c         return
+c         end
+c
+c         the value of iflag should not be changed by fcn unless
+c         the user wants to terminate execution of lmstr1.
+c         in this case set iflag to a negative integer.
+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. n must not exceed m.
+c
+c       x is an array of length n. on input x must contain
+c         an initial estimate of the solution vector. on output x
+c         contains the final estimate of the solution vector.
+c
+c       fvec is an output array of length m which contains
+c         the functions evaluated at the output x.
+c
+c       fjac is an output n by n array. the upper triangle of fjac
+c         contains an upper triangular matrix r such that
+c
+c                t     t           t
+c               p *(jac *jac)*p = r *r,
+c
+c         where p is a permutation matrix and jac is the final
+c         calculated jacobian. column j of p is column ipvt(j)
+c         (see below) of the identity matrix. the lower triangular
+c         part of fjac contains information generated during
+c         the computation of r.
+c
+c       ldfjac is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array fjac.
+c
+c       tol is a nonnegative input variable. termination occurs
+c         when the algorithm estimates either that the relative
+c         error in the sum of squares is at most tol or that
+c         the relative error between x and the solution is at
+c         most tol.
+c
+c       info is an integer output variable. if the user has
+c         terminated execution, info is set to the (negative)
+c         value of iflag. see description of fcn. otherwise,
+c         info is set as follows.
+c
+c         info = 0  improper input parameters.
+c
+c         info = 1  algorithm estimates that the relative error
+c                   in the sum of squares is at most tol.
+c
+c         info = 2  algorithm estimates that the relative error
+c                   between x and the solution is at most tol.
+c
+c         info = 3  conditions for info = 1 and info = 2 both hold.
+c
+c         info = 4  fvec is orthogonal to the columns of the
+c                   jacobian to machine precision.
+c
+c         info = 5  number of calls to fcn with iflag = 1 has
+c                   reached 100*(n+1).
+c
+c         info = 6  tol is too small. no further reduction in
+c                   the sum of squares is possible.
+c
+c         info = 7  tol is too small. no further improvement in
+c                   the approximate solution x is possible.
+c
+c       ipvt is an integer output array of length n. ipvt
+c         defines a permutation matrix p such that jac*p = q*r,
+c         where jac is the final calculated jacobian, q is
+c         orthogonal (not stored), and r is upper triangular.
+c         column j of p is column ipvt(j) of the identity matrix.
+c
+c       wa is a work array of length lwa.
+c
+c       lwa is a positive integer input variable not less than 5*n+m.
+c
+c     subprograms called
+c
+c       user-supplied ...... fcn
+c
+c       minpack-supplied ... lmstr
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
+c     jorge j. more
+c
+c     **********
+      integer maxfev,mode,nfev,njev,nprint
+      double precision factor,ftol,gtol,xtol,zero
+      data factor,zero /1.0d2,0.0d0/
+      info = 0
+c
+c     check the input parameters for errors.
+c
+      if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. n .or. tol .lt. zero
+     *    .or. lwa .lt. 5*n + m) go to 10
+c
+c     call lmstr.
+c
+      maxfev = 100*(n + 1)
+      ftol = tol
+      xtol = tol
+      gtol = zero
+      mode = 1
+      nprint = 0
+      call lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,
+     *           wa(1),mode,factor,nprint,info,nfev,njev,ipvt,wa(n+1),
+     *           wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
+      if (info .eq. 8) info = 4
+   10 continue
+      return
+c
+c     last card of subroutine lmstr1.
+c
+      end
diff --git a/math/minpack/qform.f b/math/minpack/qform.f
new file mode 100644
index 00000000..087b2478
--- /dev/null
+++ b/math/minpack/qform.f
@@ -0,0 +1,95 @@
+      subroutine qform(m,n,q,ldq,wa)
+      integer m,n,ldq
+      double precision q(ldq,m),wa(m)
+c     **********
+c
+c     subroutine qform
+c
+c     this subroutine proceeds from the computed qr factorization of
+c     an m by n matrix a to accumulate the m by m orthogonal matrix
+c     q from its factored form.
+c
+c     the subroutine statement is
+c
+c       subroutine qform(m,n,q,ldq,wa)
+c
+c     where
+c
+c       m is a positive integer input variable set to the number
+c         of rows of a and the order of q.
+c
+c       n is a positive integer input variable set to the number
+c         of columns of a.
+c
+c       q is an m by m array. on input the full lower trapezoid in
+c         the first min(m,n) columns of q contains the factored form.
+c         on output q has been accumulated into a square matrix.
+c
+c       ldq is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array q.
+c
+c       wa is a work array of length m.
+c
+c     subprograms called
+c
+c       fortran-supplied ... min0
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,j,jm1,k,l,minmn,np1
+      double precision one,sum,temp,zero
+      data one,zero /1.0d0,0.0d0/
+c
+c     zero out upper triangle of q in the first min(m,n) columns.
+c
+      minmn = min0(m,n)
+      if (minmn .lt. 2) go to 30
+      do 20 j = 2, minmn
+         jm1 = j - 1
+         do 10 i = 1, jm1
+            q(i,j) = zero
+   10       continue
+   20    continue
+   30 continue
+c
+c     initialize remaining columns to those of the identity matrix.
+c
+      np1 = n + 1
+      if (m .lt. np1) go to 60
+      do 50 j = np1, m
+         do 40 i = 1, m
+            q(i,j) = zero
+   40       continue
+         q(j,j) = one
+   50    continue
+   60 continue
+c
+c     accumulate q from its factored form.
+c
+      do 120 l = 1, minmn
+         k = minmn - l + 1
+         do 70 i = k, m
+            wa(i) = q(i,k)
+            q(i,k) = zero
+   70       continue
+         q(k,k) = one
+         if (wa(k) .eq. zero) go to 110
+         do 100 j = k, m
+            sum = zero
+            do 80 i = k, m
+               sum = sum + q(i,j)*wa(i)
+   80          continue
+            temp = sum/wa(k)
+            do 90 i = k, m
+               q(i,j) = q(i,j) - temp*wa(i)
+   90          continue
+  100       continue
+  110    continue
+  120    continue
+      return
+c
+c     last card of subroutine qform.
+c
+      end
diff --git a/math/minpack/qrfac.f b/math/minpack/qrfac.f
new file mode 100644
index 00000000..cb686086
--- /dev/null
+++ b/math/minpack/qrfac.f
@@ -0,0 +1,164 @@
+      subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
+      integer m,n,lda,lipvt
+      integer ipvt(lipvt)
+      logical pivot
+      double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
+c     **********
+c
+c     subroutine qrfac
+c
+c     this subroutine uses householder transformations with column
+c     pivoting (optional) to compute a qr factorization of the
+c     m by n matrix a. that is, qrfac determines an orthogonal
+c     matrix q, a permutation matrix p, and an upper trapezoidal
+c     matrix r with diagonal elements of nonincreasing magnitude,
+c     such that a*p = q*r. the householder transformation for
+c     column k, k = 1,2,...,min(m,n), is of the form
+c
+c                           t
+c           i - (1/u(k))*u*u
+c
+c     where u has zeros in the first k-1 positions. the form of
+c     this transformation and the method of pivoting first
+c     appeared in the corresponding linpack subroutine.
+c
+c     the subroutine statement is
+c
+c       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
+c
+c     where
+c
+c       m is a positive integer input variable set to the number
+c         of rows of a.
+c
+c       n is a positive integer input variable set to the number
+c         of columns of a.
+c
+c       a is an m by n array. on input a contains the matrix for
+c         which the qr factorization is to be computed. on output
+c         the strict upper trapezoidal part of a contains the strict
+c         upper trapezoidal part of r, and the lower trapezoidal
+c         part of a contains a factored form of q (the non-trivial
+c         elements of the u vectors described above).
+c
+c       lda is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array a.
+c
+c       pivot is a logical input variable. if pivot is set true,
+c         then column pivoting is enforced. if pivot is set false,
+c         then no column pivoting is done.
+c
+c       ipvt is an integer output array of length lipvt. ipvt
+c         defines the permutation matrix p such that a*p = q*r.
+c         column j of p is column ipvt(j) of the identity matrix.
+c         if pivot is false, ipvt is not referenced.
+c
+c       lipvt is a positive integer input variable. if pivot is false,
+c         then lipvt may be as small as 1. if pivot is true, then
+c         lipvt must be at least n.
+c
+c       rdiag is an output array of length n which contains the
+c         diagonal elements of r.
+c
+c       acnorm is an output array of length n which contains the
+c         norms of the corresponding columns of the input matrix a.
+c         if this information is not needed, then acnorm can coincide
+c         with rdiag.
+c
+c       wa is a work array of length n. if pivot is false, then wa
+c         can coincide with rdiag.
+c
+c     subprograms called
+c
+c       minpack-supplied ... dpmpar,enorm
+c
+c       fortran-supplied ... dmax1,dsqrt,min0
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,j,jp1,k,kmax,minmn
+      double precision ajnorm,epsmch,one,p05,sum,temp,zero
+      double precision dpmpar,enorm
+      data one,p05,zero /1.0d0,5.0d-2,0.0d0/
+c
+c     epsmch is the machine precision.
+c
+      epsmch = dpmpar(1)
+c
+c     compute the initial column norms and initialize several arrays.
+c
+      do 10 j = 1, n
+         acnorm(j) = enorm(m,a(1,j))
+         rdiag(j) = acnorm(j)
+         wa(j) = rdiag(j)
+         if (pivot) ipvt(j) = j
+   10    continue
+c
+c     reduce a to r with householder transformations.
+c
+      minmn = min0(m,n)
+      do 110 j = 1, minmn
+         if (.not.pivot) go to 40
+c
+c        bring the column of largest norm into the pivot position.
+c
+         kmax = j
+         do 20 k = j, n
+            if (rdiag(k) .gt. rdiag(kmax)) kmax = k
+   20       continue
+         if (kmax .eq. j) go to 40
+         do 30 i = 1, m
+            temp = a(i,j)
+            a(i,j) = a(i,kmax)
+            a(i,kmax) = temp
+   30       continue
+         rdiag(kmax) = rdiag(j)
+         wa(kmax) = wa(j)
+         k = ipvt(j)
+         ipvt(j) = ipvt(kmax)
+         ipvt(kmax) = k
+   40    continue
+c
+c        compute the householder transformation to reduce the
+c        j-th column of a to a multiple of the j-th unit vector.
+c
+         ajnorm = enorm(m-j+1,a(j,j))
+         if (ajnorm .eq. zero) go to 100
+         if (a(j,j) .lt. zero) ajnorm = -ajnorm
+         do 50 i = j, m
+            a(i,j) = a(i,j)/ajnorm
+   50       continue
+         a(j,j) = a(j,j) + one
+c
+c        apply the transformation to the remaining columns
+c        and update the norms.
+c
+         jp1 = j + 1
+         if (n .lt. jp1) go to 100
+         do 90 k = jp1, n
+            sum = zero
+            do 60 i = j, m
+               sum = sum + a(i,j)*a(i,k)
+   60          continue
+            temp = sum/a(j,j)
+            do 70 i = j, m
+               a(i,k) = a(i,k) - temp*a(i,j)
+   70          continue
+            if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
+            temp = a(j,k)/rdiag(k)
+            rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
+            if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
+            rdiag(k) = enorm(m-j,a(jp1,k))
+            wa(k) = rdiag(k)
+   80       continue
+   90       continue
+  100    continue
+         rdiag(j) = -ajnorm
+  110    continue
+      return
+c
+c     last card of subroutine qrfac.
+c
+      end
diff --git a/math/minpack/qrsolv.f b/math/minpack/qrsolv.f
new file mode 100644
index 00000000..f48954b3
--- /dev/null
+++ b/math/minpack/qrsolv.f
@@ -0,0 +1,193 @@
+      subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
+      integer n,ldr
+      integer ipvt(n)
+      double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
+c     **********
+c
+c     subroutine qrsolv
+c
+c     given an m by n matrix a, an n by n diagonal matrix d,
+c     and an m-vector b, the problem is to determine an x which
+c     solves the system
+c
+c           a*x = b ,     d*x = 0 ,
+c
+c     in the least squares sense.
+c
+c     this subroutine completes the solution of the problem
+c     if it is provided with the necessary information from the
+c     qr factorization, with column pivoting, of a. that is, if
+c     a*p = q*r, where p is a permutation matrix, q has orthogonal
+c     columns, and r is an upper triangular matrix with diagonal
+c     elements of nonincreasing magnitude, then qrsolv expects
+c     the full upper triangle of r, the permutation matrix p,
+c     and the first n components of (q transpose)*b. the system
+c     a*x = b, d*x = 0, is then equivalent to
+c
+c                  t       t
+c           r*z = q *b ,  p *d*p*z = 0 ,
+c
+c     where x = p*z. if this system does not have full rank,
+c     then a least squares solution is obtained. on output qrsolv
+c     also provides an upper triangular matrix s such that
+c
+c            t   t               t
+c           p *(a *a + d*d)*p = s *s .
+c
+c     s is computed within qrsolv and may be of separate interest.
+c
+c     the subroutine statement is
+c
+c       subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
+c
+c     where
+c
+c       n is a positive integer input variable set to the order of r.
+c
+c       r is an n by n array. on input the full upper triangle
+c         must contain the full upper triangle of the matrix r.
+c         on output the full upper triangle is unaltered, and the
+c         strict lower triangle contains the strict upper triangle
+c         (transposed) of the upper triangular matrix s.
+c
+c       ldr is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array r.
+c
+c       ipvt is an integer input array of length n which defines the
+c         permutation matrix p such that a*p = q*r. column j of p
+c         is column ipvt(j) of the identity matrix.
+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       x is an output array of length n which contains the least
+c         squares solution of the system a*x = b, d*x = 0.
+c
+c       sdiag is an output array of length n which contains the
+c         diagonal elements of the upper triangular matrix s.
+c
+c       wa is a work 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,j,jp1,k,kp1,l,nsing
+      double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
+      data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
+c
+c     copy r and (q transpose)*b to preserve input and initialize s.
+c     in particular, save the diagonal elements of r in x.
+c
+      do 20 j = 1, n
+         do 10 i = j, n
+            r(i,j) = r(j,i)
+   10       continue
+         x(j) = r(j,j)
+         wa(j) = qtb(j)
+   20    continue
+c
+c     eliminate the diagonal matrix d using a givens rotation.
+c
+      do 100 j = 1, n
+c
+c        prepare the row of d to be eliminated, locating the
+c        diagonal element using p from the qr factorization.
+c
+         l = ipvt(j)
+         if (diag(l) .eq. zero) go to 90
+         do 30 k = j, n
+            sdiag(k) = zero
+   30       continue
+         sdiag(j) = diag(l)
+c
+c        the transformations to eliminate the row of d
+c        modify only a single element of (q transpose)*b
+c        beyond the first n, which is initially zero.
+c
+         qtbpj = zero
+         do 80 k = j, n
+c
+c           determine a givens rotation which eliminates the
+c           appropriate element in the current row of d.
+c
+            if (sdiag(k) .eq. zero) go to 70
+            if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40
+               cotan = r(k,k)/sdiag(k)
+               sin = p5/dsqrt(p25+p25*cotan**2)
+               cos = sin*cotan
+               go to 50
+   40       continue
+               tan = sdiag(k)/r(k,k)
+               cos = p5/dsqrt(p25+p25*tan**2)
+               sin = cos*tan
+   50       continue
+c
+c           compute the modified diagonal element of r and
+c           the modified element of ((q transpose)*b,0).
+c
+            r(k,k) = cos*r(k,k) + sin*sdiag(k)
+            temp = cos*wa(k) + sin*qtbpj
+            qtbpj = -sin*wa(k) + cos*qtbpj
+            wa(k) = temp
+c
+c           accumulate the tranformation in the row of s.
+c
+            kp1 = k + 1
+            if (n .lt. kp1) go to 70
+            do 60 i = kp1, n
+               temp = cos*r(i,k) + sin*sdiag(i)
+               sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
+               r(i,k) = temp
+   60          continue
+   70       continue
+   80       continue
+   90    continue
+c
+c        store the diagonal element of s and restore
+c        the corresponding diagonal element of r.
+c
+         sdiag(j) = r(j,j)
+         r(j,j) = x(j)
+  100    continue
+c
+c     solve the triangular system for z. if the system is
+c     singular, then obtain a least squares solution.
+c
+      nsing = n
+      do 110 j = 1, n
+         if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
+         if (nsing .lt. n) wa(j) = zero
+  110    continue
+      if (nsing .lt. 1) go to 150
+      do 140 k = 1, nsing
+         j = nsing - k + 1
+         sum = zero
+         jp1 = j + 1
+         if (nsing .lt. jp1) go to 130
+         do 120 i = jp1, nsing
+            sum = sum + r(i,j)*wa(i)
+  120       continue
+  130    continue
+         wa(j) = (wa(j) - sum)/sdiag(j)
+  140    continue
+  150 continue
+c
+c     permute the components of z back to components of x.
+c
+      do 160 j = 1, n
+         l = ipvt(j)
+         x(l) = wa(j)
+  160    continue
+      return
+c
+c     last card of subroutine qrsolv.
+c
+      end
diff --git a/math/minpack/r1mpyq.f b/math/minpack/r1mpyq.f
new file mode 100644
index 00000000..ec99b96c
--- /dev/null
+++ b/math/minpack/r1mpyq.f
@@ -0,0 +1,92 @@
+      subroutine r1mpyq(m,n,a,lda,v,w)
+      integer m,n,lda
+      double precision a(lda,n),v(n),w(n)
+c     **********
+c
+c     subroutine r1mpyq
+c
+c     given an m by n matrix a, this subroutine computes a*q where
+c     q is the product of 2*(n - 1) transformations
+c
+c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
+c
+c     and gv(i), gw(i) are givens rotations in the (i,n) plane which
+c     eliminate elements in the i-th and n-th planes, respectively.
+c     q itself is not given, rather the information to recover the
+c     gv, gw rotations is supplied.
+c
+c     the subroutine statement is
+c
+c       subroutine r1mpyq(m,n,a,lda,v,w)
+c
+c     where
+c
+c       m is a positive integer input variable set to the number
+c         of rows of a.
+c
+c       n is a positive integer input variable set to the number
+c         of columns of a.
+c
+c       a is an m by n array. on input a must contain the matrix
+c         to be postmultiplied by the orthogonal matrix q
+c         described above. on output a*q has replaced a.
+c
+c       lda is a positive integer input variable not less than m
+c         which specifies the leading dimension of the array a.
+c
+c       v is an input array of length n. v(i) must contain the
+c         information necessary to recover the givens rotation gv(i)
+c         described above.
+c
+c       w is an input array of length n. w(i) must contain the
+c         information necessary to recover the givens rotation gw(i)
+c         described above.
+c
+c     subroutines 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,j,nmj,nm1
+      double precision cos,one,sin,temp
+      data one /1.0d0/
+c
+c     apply the first set of givens rotations to a.
+c
+      nm1 = n - 1
+      if (nm1 .lt. 1) go to 50
+      do 20 nmj = 1, nm1
+         j = n - nmj
+         if (dabs(v(j)) .gt. one) cos = one/v(j)
+         if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
+         if (dabs(v(j)) .le. one) sin = v(j)
+         if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
+         do 10 i = 1, m
+            temp = cos*a(i,j) - sin*a(i,n)
+            a(i,n) = sin*a(i,j) + cos*a(i,n)
+            a(i,j) = temp
+   10       continue
+   20    continue
+c
+c     apply the second set of givens rotations to a.
+c
+      do 40 j = 1, nm1
+         if (dabs(w(j)) .gt. one) cos = one/w(j)
+         if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
+         if (dabs(w(j)) .le. one) sin = w(j)
+         if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
+         do 30 i = 1, m
+            temp = cos*a(i,j) + sin*a(i,n)
+            a(i,n) = -sin*a(i,j) + cos*a(i,n)
+            a(i,j) = temp
+   30       continue
+   40    continue
+   50 continue
+      return
+c
+c     last card of subroutine r1mpyq.
+c
+      end
diff --git a/math/minpack/r1updt.f b/math/minpack/r1updt.f
new file mode 100644
index 00000000..e034973d
--- /dev/null
+++ b/math/minpack/r1updt.f
@@ -0,0 +1,207 @@
+      subroutine r1updt(m,n,s,ls,u,v,w,sing)
+      integer m,n,ls
+      logical sing
+      double precision s(ls),u(m),v(n),w(m)
+c     **********
+c
+c     subroutine r1updt
+c
+c     given an m by n lower trapezoidal matrix s, an m-vector u,
+c     and an n-vector v, the problem is to determine an
+c     orthogonal matrix q such that
+c
+c                   t
+c           (s + u*v )*q
+c
+c     is again lower trapezoidal.
+c
+c     this subroutine determines q as the product of 2*(n - 1)
+c     transformations
+c
+c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
+c
+c     where gv(i), gw(i) are givens rotations in the (i,n) plane
+c     which eliminate elements in the i-th and n-th planes,
+c     respectively. q itself is not accumulated, rather the
+c     information to recover the gv, gw rotations is returned.
+c
+c     the subroutine statement is
+c
+c       subroutine r1updt(m,n,s,ls,u,v,w,sing)
+c
+c     where
+c
+c       m is a positive integer input variable set to the number
+c         of rows of s.
+c
+c       n is a positive integer input variable set to the number
+c         of columns of s. n must not exceed m.
+c
+c       s is an array of length ls. on input s must contain the lower
+c         trapezoidal matrix s stored by columns. on output s contains
+c         the lower trapezoidal matrix produced as described above.
+c
+c       ls is a positive integer input variable not less than
+c         (n*(2*m-n+1))/2.
+c
+c       u is an input array of length m which must contain the
+c         vector u.
+c
+c       v is an array of length n. on input v must contain the vector
+c         v. on output v(i) contains the information necessary to
+c         recover the givens rotation gv(i) described above.
+c
+c       w is an output array of length m. w(i) contains information
+c         necessary to recover the givens rotation gw(i) described
+c         above.
+c
+c       sing is a logical output variable. sing is set true if any
+c         of the diagonal elements of the output s are zero. otherwise
+c         sing is set false.
+c
+c     subprograms called
+c
+c       minpack-supplied ... dpmpar
+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     john l. nazareth
+c
+c     **********
+      integer i,j,jj,l,nmj,nm1
+      double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,
+     *                 zero
+      double precision dpmpar
+      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
+c
+c     giant is the largest magnitude.
+c
+      giant = dpmpar(3)
+c
+c     initialize the diagonal element pointer.
+c
+      jj = (n*(2*m - n + 1))/2 - (m - n)
+c
+c     move the nontrivial part of the last column of s into w.
+c
+      l = jj
+      do 10 i = n, m
+         w(i) = s(l)
+         l = l + 1
+   10    continue
+c
+c     rotate the vector v into a multiple of the n-th unit vector
+c     in such a way that a spike is introduced into w.
+c
+      nm1 = n - 1
+      if (nm1 .lt. 1) go to 70
+      do 60 nmj = 1, nm1
+         j = n - nmj
+         jj = jj - (m - j + 1)
+         w(j) = zero
+         if (v(j) .eq. zero) go to 50
+c
+c        determine a givens rotation which eliminates the
+c        j-th element of v.
+c
+         if (dabs(v(n)) .ge. dabs(v(j))) go to 20
+            cotan = v(n)/v(j)
+            sin = p5/dsqrt(p25+p25*cotan**2)
+            cos = sin*cotan
+            tau = one
+            if (dabs(cos)*giant .gt. one) tau = one/cos
+            go to 30
+   20    continue
+            tan = v(j)/v(n)
+            cos = p5/dsqrt(p25+p25*tan**2)
+            sin = cos*tan
+            tau = sin
+   30    continue
+c
+c        apply the transformation to v and store the information
+c        necessary to recover the givens rotation.
+c
+         v(n) = sin*v(j) + cos*v(n)
+         v(j) = tau
+c
+c        apply the transformation to s and extend the spike in w.
+c
+         l = jj
+         do 40 i = j, m
+            temp = cos*s(l) - sin*w(i)
+            w(i) = sin*s(l) + cos*w(i)
+            s(l) = temp
+            l = l + 1
+   40       continue
+   50    continue
+   60    continue
+   70 continue
+c
+c     add the spike from the rank 1 update to w.
+c
+      do 80 i = 1, m
+         w(i) = w(i) + v(n)*u(i)
+   80    continue
+c
+c     eliminate the spike.
+c
+      sing = .false.
+      if (nm1 .lt. 1) go to 140
+      do 130 j = 1, nm1
+         if (w(j) .eq. zero) go to 120
+c
+c        determine a givens rotation which eliminates the
+c        j-th element of the spike.
+c
+         if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
+            cotan = s(jj)/w(j)
+            sin = p5/dsqrt(p25+p25*cotan**2)
+            cos = sin*cotan
+            tau = one
+            if (dabs(cos)*giant .gt. one) tau = one/cos
+            go to 100
+   90    continue
+            tan = w(j)/s(jj)
+            cos = p5/dsqrt(p25+p25*tan**2)
+            sin = cos*tan
+            tau = sin
+  100    continue
+c
+c        apply the transformation to s and reduce the spike in w.
+c
+         l = jj
+         do 110 i = j, m
+            temp = cos*s(l) + sin*w(i)
+            w(i) = -sin*s(l) + cos*w(i)
+            s(l) = temp
+            l = l + 1
+  110       continue
+c
+c        store the information necessary to recover the
+c        givens rotation.
+c
+         w(j) = tau
+  120    continue
+c
+c        test for zero diagonal elements in the output s.
+c
+         if (s(jj) .eq. zero) sing = .true.
+         jj = jj + (m - j + 1)
+  130    continue
+  140 continue
+c
+c     move w back into the last column of the output s.
+c
+      l = jj
+      do 150 i = n, m
+         s(l) = w(i)
+         l = l + 1
+  150    continue
+      if (s(jj) .eq. zero) sing = .true.
+      return
+c
+c     last card of subroutine r1updt.
+c
+      end
diff --git a/math/minpack/rwupdt.f b/math/minpack/rwupdt.f
new file mode 100644
index 00000000..05282b55
--- /dev/null
+++ b/math/minpack/rwupdt.f
@@ -0,0 +1,113 @@
+      subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+      integer n,ldr
+      double precision alpha
+      double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
+c     **********
+c
+c     subroutine rwupdt
+c
+c     given an n by n upper triangular matrix r, this subroutine
+c     computes the qr decomposition of the matrix formed when a row
+c     is added to r. if the row is specified by the vector w, then
+c     rwupdt determines an orthogonal matrix q such that when the
+c     n+1 by n matrix composed of r augmented by w is premultiplied
+c     by (q transpose), the resulting matrix is upper trapezoidal.
+c     the matrix (q transpose) is the product of n transformations
+c
+c           g(n)*g(n-1)* ... *g(1)
+c
+c     where g(i) is a givens rotation in the (i,n+1) plane which
+c     eliminates elements in the (n+1)-st plane. rwupdt also
+c     computes the product (q transpose)*c where c is the
+c     (n+1)-vector (b,alpha). q itself is not accumulated, rather
+c     the information to recover the g rotations is supplied.
+c
+c     the subroutine statement is
+c
+c       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
+c
+c     where
+c
+c       n is a positive integer input variable set to the order of r.
+c
+c       r is an n by n array. on input the upper triangular part of
+c         r must contain the matrix to be updated. on output r
+c         contains the updated triangular matrix.
+c
+c       ldr is a positive integer input variable not less than n
+c         which specifies the leading dimension of the array r.
+c
+c       w is an input array of length n which must contain the row
+c         vector to be added to r.
+c
+c       b is an array of length n. on input b must contain the
+c         first n elements of the vector c. on output b contains
+c         the first n elements of the vector (q transpose)*c.
+c
+c       alpha is a variable. on input alpha must contain the
+c         (n+1)-st element of the vector c. on output alpha contains
+c         the (n+1)-st element of the vector (q transpose)*c.
+c
+c       cos is an output array of length n which contains the
+c         cosines of the transforming givens rotations.
+c
+c       sin is an output array of length n which contains the
+c         sines of the transforming givens rotations.
+c
+c     subprograms called
+c
+c       fortran-supplied ... dabs,dsqrt
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
+c     jorge j. more
+c
+c     **********
+      integer i,j,jm1
+      double precision cotan,one,p5,p25,rowj,tan,temp,zero
+      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
+c
+      do 60 j = 1, n
+         rowj = w(j)
+         jm1 = j - 1
+c
+c        apply the previous transformations to
+c        r(i,j), i=1,2,...,j-1, and to w(j).
+c
+         if (jm1 .lt. 1) go to 20
+         do 10 i = 1, jm1
+            temp = cos(i)*r(i,j) + sin(i)*rowj
+            rowj = -sin(i)*r(i,j) + cos(i)*rowj
+            r(i,j) = temp
+   10       continue
+   20    continue
+c
+c        determine a givens rotation which eliminates w(j).
+c
+         cos(j) = one
+         sin(j) = zero
+         if (rowj .eq. zero) go to 50
+         if (dabs(r(j,j)) .ge. dabs(rowj)) go to 30
+            cotan = r(j,j)/rowj
+            sin(j) = p5/dsqrt(p25+p25*cotan**2)
+            cos(j) = sin(j)*cotan
+            go to 40
+   30    continue
+            tan = rowj/r(j,j)
+            cos(j) = p5/dsqrt(p25+p25*tan**2)
+            sin(j) = cos(j)*tan
+   40    continue
+c
+c        apply the current transformation to r(j,j), b(j), and alpha.
+c
+         r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
+         temp = cos(j)*b(j) + sin(j)*alpha
+         alpha = -sin(j)*b(j) + cos(j)*alpha
+         b(j) = temp
+   50    continue
+   60    continue
+      return
+c
+c     last card of subroutine rwupdt.
+c
+      end