share/minpack/fortran/fdjac1.f

   1       subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
   2      *                  wa1,wa2)
   3       integer n,ldfjac,iflag,ml,mu
   4       double precision epsfcn
   5       double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
   6 c     **********
   7 c
   8 c     subroutine fdjac1
   9 c
  10 c     this subroutine computes a forward-difference approximation
  11 c     to the n by n jacobian matrix associated with a specified
  12 c     problem of n functions in n variables. if the jacobian has
  13 c     a banded form, then function evaluations are saved by only
  14 c     approximating the nonzero terms.
  15 c
  16 c     the subroutine statement is
  17 c
  18 c       subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
  19 c                         wa1,wa2)
  20 c
  21 c     where
  22 c
  23 c       fcn is the name of the user-supplied subroutine which
  24 c         calculates the functions. fcn must be declared
  25 c         in an external statement in the user calling
  26 c         program, and should be written as follows.
  27 c
  28 c         subroutine fcn(n,x,fvec,iflag)
  29 c         integer n,iflag
  30 c         double precision x(n),fvec(n)
  31 c         ----------
  32 c         calculate the functions at x and
  33 c         return this vector in fvec.
  34 c         ----------
  35 c         return
  36 c         end
  37 c
  38 c         the value of iflag should not be changed by fcn unless
  39 c         the user wants to terminate execution of fdjac1.
  40 c         in this case set iflag to a negative integer.
  41 c
  42 c       n is a positive integer input variable set to the number
  43 c         of functions and variables.
  44 c
  45 c       x is an input array of length n.
  46 c
  47 c       fvec is an input array of length n which must contain the
  48 c         functions evaluated at x.
  49 c
  50 c       fjac is an output n by n array which contains the
  51 c         approximation to the jacobian matrix evaluated at x.
  52 c
  53 c       ldfjac is a positive integer input variable not less than n
  54 c         which specifies the leading dimension of the array fjac.
  55 c
  56 c       iflag is an integer variable which can be used to terminate
  57 c         the execution of fdjac1. see description of fcn.
  58 c
  59 c       ml is a nonnegative integer input variable which specifies
  60 c         the number of subdiagonals within the band of the
  61 c         jacobian matrix. if the jacobian is not banded, set
  62 c         ml to at least n - 1.
  63 c
  64 c       epsfcn is an input variable used in determining a suitable
  65 c         step length for the forward-difference approximation. this
  66 c         approximation assumes that the relative errors in the
  67 c         functions are of the order of epsfcn. if epsfcn is less
  68 c         than the machine precision, it is assumed that the relative
  69 c         errors in the functions are of the order of the machine
  70 c         precision.
  71 c
  72 c       mu is a nonnegative integer input variable which specifies
  73 c         the number of superdiagonals within the band of the
  74 c         jacobian matrix. if the jacobian is not banded, set
  75 c         mu to at least n - 1.
  76 c
  77 c       wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
  78 c         least n, then the jacobian is considered dense, and wa2 is
  79 c         not referenced.
  80 c
  81 c     subprograms called
  82 c
  83 c       minpack-supplied ... dpmpar
  84 c
  85 c       fortran-supplied ... dabs,dmax1,dsqrt
  86 c
  87 c     argonne national laboratory. minpack project. march 1980.
  88 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  89 c
  90 c     **********
  91       integer i,j,k,msum
  92       double precision eps,epsmch,h,temp,zero
  93       double precision dpmpar
  94       data zero /0.0d0/
  95 c
  96 c     epsmch is the machine precision.
  97 c
  98       epsmch = dpmpar(1)
  99 c
 100       eps = dsqrt(dmax1(epsfcn,epsmch))
 101       msum = ml + mu + 1
 102       if (msum .lt. n) go to 40
 103 c
 104 c        computation of dense approximate jacobian.
 105 c
 106          do 20 j = 1, n
 107             temp = x(j)
 108             h = eps*dabs(temp)
 109             if (h .eq. zero) h = eps
 110             x(j) = temp + h
 111             call fcn(n,x,wa1,iflag)
 112             if (iflag .lt. 0) go to 30
 113             x(j) = temp
 114             do 10 i = 1, n
 115                fjac(i,j) = (wa1(i) - fvec(i))/h
 116    10          continue
 117    20       continue
 118    30    continue
 119          go to 110
 120    40 continue
 121 c
 122 c        computation of banded approximate jacobian.
 123 c
 124          do 90 k = 1, msum
 125             do 60 j = k, n, msum
 126                wa2(j) = x(j)
 127                h = eps*dabs(wa2(j))
 128                if (h .eq. zero) h = eps
 129                x(j) = wa2(j) + h
 130    60          continue
 131             call fcn(n,x,wa1,iflag)
 132             if (iflag .lt. 0) go to 100
 133             do 80 j = k, n, msum
 134                x(j) = wa2(j)
 135                h = eps*dabs(wa2(j))
 136                if (h .eq. zero) h = eps
 137                do 70 i = 1, n
 138                   fjac(i,j) = zero
 139                   if (i .ge. j - mu .and. i .le. j + ml)
 140      *               fjac(i,j) = (wa1(i) - fvec(i))/h
 141    70             continue
 142    80          continue
 143    90       continue
 144   100    continue
 145   110 continue
 146       return
 147 c
 148 c     last card of subroutine fdjac1.
 149 c
 150       end
 151