diff options
Diffstat (limited to 'math/minpack/fdjac1.f')
| -rw-r--r-- | math/minpack/fdjac1.f | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/math/minpack/fdjac1.f b/math/minpack/fdjac1.f new file mode 100644 index 00000000..aa058010 --- /dev/null +++ b/math/minpack/fdjac1.f @@ -0,0 +1,150 @@ + subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn, + * wa1,wa2) + integer n,ldfjac,iflag,ml,mu + double precision epsfcn + double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n) +c ********** +c +c subroutine fdjac1 +c +c this subroutine computes a forward-difference approximation +c to the n by n jacobian matrix associated with a specified +c problem of n functions in n variables. if the jacobian has +c a banded form, then function evaluations are saved by only +c approximating the nonzero terms. +c +c the subroutine statement is +c +c subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn, +c wa1,wa2) +c +c where +c +c fcn is the name of the user-supplied subroutine which +c calculates the functions. fcn must be declared +c in an external statement in the user calling +c program, and should be written as follows. +c +c subroutine fcn(n,x,fvec,iflag) +c integer n,iflag +c double precision x(n),fvec(n) +c ---------- +c calculate the functions at x and +c return this vector in fvec. +c ---------- +c return +c end +c +c the value of iflag should not be changed by fcn unless +c the user wants to terminate execution of fdjac1. +c in this case set iflag to a negative integer. +c +c n is a positive integer input variable set to the number +c of functions and variables. +c +c x is an input array of length n. +c +c fvec is an input array of length n which must contain the +c functions evaluated at x. +c +c fjac is an output n by n array which contains the +c approximation to the jacobian matrix evaluated at x. +c +c ldfjac is a positive integer input variable not less than n +c which specifies the leading dimension of the array fjac. +c +c iflag is an integer variable which can be used to terminate +c the execution of fdjac1. see description of fcn. +c +c ml is a nonnegative integer input variable which specifies +c the number of subdiagonals within the band of the +c jacobian matrix. if the jacobian is not banded, set +c ml to at least n - 1. +c +c epsfcn is an input variable used in determining a suitable +c step length for the forward-difference approximation. this +c approximation assumes that the relative errors in the +c functions are of the order of epsfcn. if epsfcn is less +c than the machine precision, it is assumed that the relative +c errors in the functions are of the order of the machine +c precision. +c +c mu is a nonnegative integer input variable which specifies +c the number of superdiagonals within the band of the +c jacobian matrix. if the jacobian is not banded, set +c mu to at least n - 1. +c +c wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at +c least n, then the jacobian is considered dense, and wa2 is +c not referenced. +c +c subprograms called +c +c minpack-supplied ... dpmpar +c +c fortran-supplied ... dabs,dmax1,dsqrt +c +c argonne national laboratory. minpack project. march 1980. +c burton s. garbow, kenneth e. hillstrom, jorge j. more +c +c ********** + integer i,j,k,msum + double precision eps,epsmch,h,temp,zero + double precision dpmpar + data zero /0.0d0/ +c +c epsmch is the machine precision. +c + epsmch = dpmpar(1) +c + eps = dsqrt(dmax1(epsfcn,epsmch)) + msum = ml + mu + 1 + if (msum .lt. n) go to 40 +c +c computation of dense approximate jacobian. +c + do 20 j = 1, n + temp = x(j) + h = eps*dabs(temp) + if (h .eq. zero) h = eps + x(j) = temp + h + call fcn(n,x,wa1,iflag) + if (iflag .lt. 0) go to 30 + x(j) = temp + do 10 i = 1, n + fjac(i,j) = (wa1(i) - fvec(i))/h + 10 continue + 20 continue + 30 continue + go to 110 + 40 continue +c +c computation of banded approximate jacobian. +c + do 90 k = 1, msum + do 60 j = k, n, msum + wa2(j) = x(j) + h = eps*dabs(wa2(j)) + if (h .eq. zero) h = eps + x(j) = wa2(j) + h + 60 continue + call fcn(n,x,wa1,iflag) + if (iflag .lt. 0) go to 100 + do 80 j = k, n, msum + x(j) = wa2(j) + h = eps*dabs(wa2(j)) + if (h .eq. zero) h = eps + do 70 i = 1, n + fjac(i,j) = zero + if (i .ge. j - mu .and. i .le. j + ml) + * fjac(i,j) = (wa1(i) - fvec(i))/h + 70 continue + 80 continue + 90 continue + 100 continue + 110 continue + return +c +c last card of subroutine fdjac1. +c + end |