Description: Link against minpack from Debian package instead of embedded copy
 - remove all src/*.f files
 - add -lminpack to PKG_LIBS in Makevars
Author: Sébastien Villemot <sebastien@debian.org>
Forwarded: not-needed
Last-Update: 2018-01-30
---
This patch header follows DEP-3: http://dep.debian.net/deps/dep3/
--- a/src/chkder.f
+++ /dev/null
@@ -1,140 +0,0 @@
-      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
--- a/src/dogleg.f
+++ /dev/null
@@ -1,177 +0,0 @@
-      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
--- a/src/dpmpar.f
+++ /dev/null
@@ -1,177 +0,0 @@
-      double precision function dpmpar(i)
-      integer i
-c     **********
-c
-c     Function dpmpar
-c
-c     This function provides double precision machine parameters
-c     when the appropriate set of data statements is activated (by
-c     removing the c from column 1) and all other data statements are
-c     rendered inactive. Most of the parameter values were obtained
-c     from the corresponding Bell Laboratories Port Library function.
-c
-c     The function statement is
-c
-c       double precision function dpmpar(i)
-c
-c     where
-c
-c       i is an integer input variable set to 1, 2, or 3 which
-c         selects the desired machine parameter. If the machine has
-c         t base b digits and its smallest and largest exponents are
-c         emin and emax, respectively, then these parameters are
-c
-c         dpmpar(1) = b**(1 - t), the machine precision,
-c
-c         dpmpar(2) = b**(emin - 1), the smallest magnitude,
-c
-c         dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
-c
-c     Argonne National Laboratory. MINPACK Project. November 1996.
-c     Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
-c
-c     **********
-      integer mcheps(4)
-      integer minmag(4)
-      integer maxmag(4)
-      double precision dmach(3)
-      equivalence (dmach(1),mcheps(1))
-      equivalence (dmach(2),minmag(1))
-      equivalence (dmach(3),maxmag(1))
-c
-c     Machine constants for the IBM 360/370 series,
-c     the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
-c     the Xerox Sigma 5/7/9 and the Sel systems 85/86.
-c
-c     data mcheps(1),mcheps(2) / z34100000, z00000000 /
-c     data minmag(1),minmag(2) / z00100000, z00000000 /
-c     data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
-c
-c     Machine constants for the Honeywell 600/6000 series.
-c
-c     data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
-c     data minmag(1),minmag(2) / o402400000000, o000000000000 /
-c     data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
-c
-c     Machine constants for the CDC 6000/7000 series.
-c
-c     data mcheps(1) / 15614000000000000000b /
-c     data mcheps(2) / 15010000000000000000b /
-c
-c     data minmag(1) / 00604000000000000000b /
-c     data minmag(2) / 00000000000000000000b /
-c
-c     data maxmag(1) / 37767777777777777777b /
-c     data maxmag(2) / 37167777777777777777b /
-c
-c     Machine constants for the PDP-10 (KA processor).
-c
-c     data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
-c     data minmag(1),minmag(2) / "033400000000, "000000000000 /
-c     data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
-c
-c     Machine constants for the PDP-10 (KI processor).
-c
-c     data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
-c     data minmag(1),minmag(2) / "000400000000, "000000000000 /
-c     data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
-c
-c     Machine constants for the PDP-11. 
-c
-c     data mcheps(1),mcheps(2) /   9472,      0 /
-c     data mcheps(3),mcheps(4) /      0,      0 /
-c
-c     data minmag(1),minmag(2) /    128,      0 /
-c     data minmag(3),minmag(4) /      0,      0 /
-c
-c     data maxmag(1),maxmag(2) /  32767,     -1 /
-c     data maxmag(3),maxmag(4) /     -1,     -1 /
-c
-c     Machine constants for the Burroughs 6700/7700 systems.
-c
-c     data mcheps(1) / o1451000000000000 /
-c     data mcheps(2) / o0000000000000000 /
-c
-c     data minmag(1) / o1771000000000000 /
-c     data minmag(2) / o7770000000000000 /
-c
-c     data maxmag(1) / o0777777777777777 /
-c     data maxmag(2) / o7777777777777777 /
-c
-c     Machine constants for the Burroughs 5700 system.
-c
-c     data mcheps(1) / o1451000000000000 /
-c     data mcheps(2) / o0000000000000000 /
-c
-c     data minmag(1) / o1771000000000000 /
-c     data minmag(2) / o0000000000000000 /
-c
-c     data maxmag(1) / o0777777777777777 /
-c     data maxmag(2) / o0007777777777777 /
-c
-c     Machine constants for the Burroughs 1700 system.
-c
-c     data mcheps(1) / zcc6800000 /
-c     data mcheps(2) / z000000000 /
-c
-c     data minmag(1) / zc00800000 /
-c     data minmag(2) / z000000000 /
-c
-c     data maxmag(1) / zdffffffff /
-c     data maxmag(2) / zfffffffff /
-c
-c     Machine constants for the Univac 1100 series.
-c
-c     data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
-c     data minmag(1),minmag(2) / o000040000000, o000000000000 /
-c     data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
-c
-c     Machine constants for the Data General Eclipse S/200.
-c
-c     Note - it may be appropriate to include the following card -
-c     static dmach(3)
-c
-c     data minmag/20k,3*0/,maxmag/77777k,3*177777k/
-c     data mcheps/32020k,3*0/
-c
-c     Machine constants for the Harris 220.
-c
-c     data mcheps(1),mcheps(2) / '20000000, '00000334 /
-c     data minmag(1),minmag(2) / '20000000, '00000201 /
-c     data maxmag(1),maxmag(2) / '37777777, '37777577 /
-c
-c     Machine constants for the Cray-1.
-c
-c     data mcheps(1) / 0376424000000000000000b /
-c     data mcheps(2) / 0000000000000000000000b /
-c
-c     data minmag(1) / 0200034000000000000000b /
-c     data minmag(2) / 0000000000000000000000b /
-c
-c     data maxmag(1) / 0577777777777777777777b /
-c     data maxmag(2) / 0000007777777777777776b /
-c
-c     Machine constants for the Prime 400.
-c
-c     data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
-c     data minmag(1),minmag(2) / :10000000000, :00000100000 /
-c     data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
-c
-c     Machine constants for the VAX-11.
-c
-c     data mcheps(1),mcheps(2) /   9472,  0 /
-c     data minmag(1),minmag(2) /    128,  0 /
-c     data maxmag(1),maxmag(2) / -32769, -1 /
-c
-c     Machine constants for IEEE machines.
-c
-      data dmach(1) /2.22044604926d-16/
-      data dmach(2) /2.22507385852d-308/
-      data dmach(3) /1.79769313485d+308/
-c
-      dpmpar = dmach(i)
-      return
-c
-c     Last card of function dpmpar.
-c
-      end
--- a/src/enorm.f
+++ /dev/null
@@ -1,108 +0,0 @@
-      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
--- a/src/fdjac2.f
+++ /dev/null
@@ -1,107 +0,0 @@
-      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
--- a/src/lmder.f
+++ /dev/null
@@ -1,452 +0,0 @@
-      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
--- a/src/lmdif.f
+++ /dev/null
@@ -1,454 +0,0 @@
-      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
--- a/src/lmpar.f
+++ /dev/null
@@ -1,264 +0,0 @@
-      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
--- a/src/qform.f
+++ /dev/null
@@ -1,95 +0,0 @@
-      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
--- a/src/qrfac.f
+++ /dev/null
@@ -1,164 +0,0 @@
-      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
--- a/src/qrsolv.f
+++ /dev/null
@@ -1,193 +0,0 @@
-      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
--- a/src/r1mpyq.f
+++ /dev/null
@@ -1,92 +0,0 @@
-      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
--- a/src/r1updt.f
+++ /dev/null
@@ -1,207 +0,0 @@
-      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
--- a/src/rwupdt.f
+++ /dev/null
@@ -1,113 +0,0 @@
-      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
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,2 @@
-PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)  
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lminpack