share/minpack/fortran/chkder.f

   1       subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
   2       integer m,n,ldfjac,mode
   3       double precision x(n),fvec(m),fjac(ldfjac,n),xp(n),fvecp(m),
   4      *                 err(m)
   5 c     **********
   6 c
   7 c     subroutine chkder
   8 c
   9 c     this subroutine checks the gradients of m nonlinear functions
  10 c     in n variables, evaluated at a point x, for consistency with
  11 c     the functions themselves. the user must call chkder twice,
  12 c     first with mode = 1 and then with mode = 2.
  13 c
  14 c     mode = 1. on input, x must contain the point of evaluation.
  15 c               on output, xp is set to a neighboring point.
  16 c
  17 c     mode = 2. on input, fvec must contain the functions and the
  18 c                         rows of fjac must contain the gradients
  19 c                         of the respective functions each evaluated
  20 c                         at x, and fvecp must contain the functions
  21 c                         evaluated at xp.
  22 c               on output, err contains measures of correctness of
  23 c                          the respective gradients.
  24 c
  25 c     the subroutine does not perform reliably if cancellation or
  26 c     rounding errors cause a severe loss of significance in the
  27 c     evaluation of a function. therefore, none of the components
  28 c     of x should be unusually small (in particular, zero) or any
  29 c     other value which may cause loss of significance.
  30 c
  31 c     the subroutine statement is
  32 c
  33 c       subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err)
  34 c
  35 c     where
  36 c
  37 c       m is a positive integer input variable set to the number
  38 c         of functions.
  39 c
  40 c       n is a positive integer input variable set to the number
  41 c         of variables.
  42 c
  43 c       x is an input array of length n.
  44 c
  45 c       fvec is an array of length m. on input when mode = 2,
  46 c         fvec must contain the functions evaluated at x.
  47 c
  48 c       fjac is an m by n array. on input when mode = 2,
  49 c         the rows of fjac must contain the gradients of
  50 c         the respective functions evaluated at x.
  51 c
  52 c       ldfjac is a positive integer input parameter not less than m
  53 c         which specifies the leading dimension of the array fjac.
  54 c
  55 c       xp is an array of length n. on output when mode = 1,
  56 c         xp is set to a neighboring point of x.
  57 c
  58 c       fvecp is an array of length m. on input when mode = 2,
  59 c         fvecp must contain the functions evaluated at xp.
  60 c
  61 c       mode is an integer input variable set to 1 on the first call
  62 c         and 2 on the second. other values of mode are equivalent
  63 c         to mode = 1.
  64 c
  65 c       err is an array of length m. on output when mode = 2,
  66 c         err contains measures of correctness of the respective
  67 c         gradients. if there is no severe loss of significance,
  68 c         then if err(i) is 1.0 the i-th gradient is correct,
  69 c         while if err(i) is 0.0 the i-th gradient is incorrect.
  70 c         for values of err between 0.0 and 1.0, the categorization
  71 c         is less certain. in general, a value of err(i) greater
  72 c         than 0.5 indicates that the i-th gradient is probably
  73 c         correct, while a value of err(i) less than 0.5 indicates
  74 c         that the i-th gradient is probably incorrect.
  75 c
  76 c     subprograms called
  77 c
  78 c       minpack supplied ... dpmpar
  79 c
  80 c       fortran supplied ... dabs,dlog10,dsqrt
  81 c
  82 c     argonne national laboratory. minpack project. march 1980.
  83 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  84 c
  85 c     **********
  86       integer i,j
  87       double precision eps,epsf,epslog,epsmch,factor,one,temp,zero
  88       double precision dpmpar
  89       data factor,one,zero /1.0d2,1.0d0,0.0d0/
  90 c
  91 c     epsmch is the machine precision.
  92 c
  93       epsmch = dpmpar(1)
  94 c
  95       eps = dsqrt(epsmch)
  96 c
  97       if (mode .eq. 2) go to 20
  98 c
  99 c        mode = 1.
 100 c
 101          do 10 j = 1, n
 102             temp = eps*dabs(x(j))
 103             if (temp .eq. zero) temp = eps
 104             xp(j) = x(j) + temp
 105    10       continue
 106          go to 70
 107    20 continue
 108 c
 109 c        mode = 2.
 110 c
 111          epsf = factor*epsmch
 112          epslog = dlog10(eps)
 113          do 30 i = 1, m
 114             err(i) = zero
 115    30       continue
 116          do 50 j = 1, n
 117             temp = dabs(x(j))
 118             if (temp .eq. zero) temp = one
 119             do 40 i = 1, m
 120                err(i) = err(i) + temp*fjac(i,j)
 121    40          continue
 122    50       continue
 123          do 60 i = 1, m
 124             temp = one
 125             if (fvec(i) .ne. zero .and. fvecp(i) .ne. zero
 126      *          .and. dabs(fvecp(i)-fvec(i)) .ge. epsf*dabs(fvec(i)))
 127      *         temp = eps*dabs((fvecp(i)-fvec(i))/eps-err(i))
 128      *                /(dabs(fvec(i)) + dabs(fvecp(i)))
 129             err(i) = one
 130             if (temp .gt. epsmch .and. temp .lt. eps)
 131      *         err(i) = (dlog10(temp) - epslog)/epslog
 132             if (temp .ge. eps) err(i) = zero
 133    60       continue
 134    70 continue
 135 c
 136       return
 137 c
 138 c     last card of subroutine chkder.
 139 c
 140       end