Add current minpack from netlib.org

LICENSE

@@ -0,0 +1,51 @@
+Minpack Copyright Notice (1999) University of Chicago.  All rights reserved
+
+Redistribution and use in source and binary forms, with or
+without modification, are permitted provided that the
+following conditions are met:
+
+1. Redistributions of source code must retain the above
+copyright notice, this list of conditions and the following
+disclaimer.
+
+2. Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following
+disclaimer in the documentation and/or other materials
+provided with the distribution.
+
+3. The end-user documentation included with the
+redistribution, if any, must include the following
+acknowledgment:
+
+   "This product includes software developed by the
+   University of Chicago, as Operator of Argonne National
+   Laboratory.
+
+Alternately, this acknowledgment may appear in the software
+itself, if and wherever such third-party acknowledgments
+normally appear.
+
+4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+BE CORRECTED.
+
+5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+POSSIBILITY OF SUCH LOSS OR DAMAGES.

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

chkdrv.f

@@ -0,0 +1,87 @@
+c     **********
+c
+c     this program tests the ability of chkder to detect
+c     inconsistencies between functions and their first derivatives.
+c     fourteen test function vectors and jacobians are used. eleven of
+c     the tests are false(f), i.e. there are inconsistencies between
+c     the function vectors and the corresponding jacobians. three of
+c     the tests are true(t), i.e. there are no inconsistencies. the
+c     driver reads in data, calls chkder and prints out information
+c     required by and received from chkder.
+c
+c     subprograms called
+c
+c       minpack supplied ... chkder,errjac,initpt,vecfcn
+c
+c     argonne national laboratory. minpack project. march 1980.
+c     burton s. garbow, kenneth e. hillstrom, jorge j. more
+c
+c     **********
+      integer i,ldfjac,lnp,mode,n,nprob,nread,nwrite
+      integer na(14),np(14)
+      logical a(14)
+      double precision cp,one
+      double precision diff(10),err(10),errmax(14),errmin(14),
+     *                 fjac(10,10),fvec1(10),fvec2(10),x1(10),x2(10)
+c
+c     logical input unit is assumed to be number 5.
+c     logical output unit is assumed to be number 6.
+c
+      data nread,nwrite /5,6/
+c
+      data a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9),a(10),a(11),
+     *     a(12),a(13),a(14)
+     *     /.false.,.false.,.false.,.true.,.false.,.false.,.false.,
+     *      .true.,.false.,.false.,.false.,.false.,.true.,.false./
+      data cp,one /1.23d-1,1.0d0/
+      ldfjac = 10
+   10 continue
+         read (nread,60) nprob,n
+         if (nprob .le. 0) go to 40
+         call initpt(n,x1,nprob,one)
+         do 20 i = 1, n
+            x1(i) = x1(i) + cp
+            cp = -cp
+   20       continue
+         write (nwrite,70) nprob,n,a(nprob)
+         mode = 1
+         call chkder(n,n,x1,fvec1,fjac,ldfjac,x2,fvec2,mode,err)
+         mode = 2
+         call vecfcn(n,x1,fvec1,nprob)
+         call errjac(n,x1,fjac,ldfjac,nprob)
+         call vecfcn(n,x2,fvec2,nprob)
+         call chkder(n,n,x1,fvec1,fjac,ldfjac,x2,fvec2,mode,err)
+         errmin(nprob) = err(1)
+         errmax(nprob) = err(1)
+         do 30 i = 1, n
+            diff(i) = fvec2(i) - fvec1(i)
+            if (errmin(nprob) .gt. err(i)) errmin(nprob) = err(i)
+            if (errmax(nprob) .lt. err(i)) errmax(nprob) = err(i)
+   30       continue
+         np(nprob) = nprob
+         lnp = nprob
+         na(nprob) = n
+         write (nwrite,80) (fvec1(i), i = 1, n)
+         write (nwrite,90) (diff(i), i = 1, n)
+         write (nwrite,100) (err(i), i = 1, n)
+         go to 10
+   40 continue
+      write (nwrite,110) lnp
+      write (nwrite,120)
+      do 50 i = 1, lnp
+         write (nwrite,130) np(i),na(i),a(i),errmin(i),errmax(i)
+   50    continue
+      stop
+   60 format (2i5)
+   70 format ( /// 5x, 8h problem, i5, 5x, 15h with dimension, i5, 2x,
+     *         5h is  , l1)
+   80 format ( // 5x, 25h first function vector    // (5x, 5d15.7))
+   90 format ( // 5x, 27h function difference vector // (5x, 5d15.7))
+  100 format ( // 5x, 13h error vector // (5x, 5d15.7))
+  110 format (12h1summary of , i3, 16h tests of chkder /)
+  120 format (46h nprob   n    status     errmin         errmax /)
+  130 format (i4, i6, 6x, l1, 3x, 2d15.7)
+c
+c     last card of derivative check test driver.
+c
+      end