share/minpack/fortran/lmstr1.f

   1       subroutine lmstr1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,ipvt,wa,
   2      *                  lwa)
   3       integer m,n,ldfjac,info,lwa
   4       integer ipvt(n)
   5       double precision tol
   6       double precision x(n),fvec(m),fjac(ldfjac,n),wa(lwa)
   7       external fcn
   8 c     **********
   9 c
  10 c     subroutine lmstr1
  11 c
  12 c     the purpose of lmstr1 is to minimize the sum of the squares of
  13 c     m nonlinear functions in n variables by a modification of
  14 c     the levenberg-marquardt algorithm which uses minimal storage.
  15 c     this is done by using the more general least-squares solver
  16 c     lmstr. the user must provide a subroutine which calculates
  17 c     the functions and the rows of the jacobian.
  18 c
  19 c     the subroutine statement is
  20 c
  21 c       subroutine lmstr1(fcn,m,n,x,fvec,fjac,ldfjac,tol,info,
  22 c                         ipvt,wa,lwa)
  23 c
  24 c     where
  25 c
  26 c       fcn is the name of the user-supplied subroutine which
  27 c         calculates the functions and the rows of the jacobian.
  28 c         fcn must be declared in an external statement in the
  29 c         user calling program, and should be written as follows.
  30 c
  31 c         subroutine fcn(m,n,x,fvec,fjrow,iflag)
  32 c         integer m,n,iflag
  33 c         double precision x(n),fvec(m),fjrow(n)
  34 c         ----------
  35 c         if iflag = 1 calculate the functions at x and
  36 c         return this vector in fvec.
  37 c         if iflag = i calculate the (i-1)-st row of the
  38 c         jacobian at x and return this vector in fjrow.
  39 c         ----------
  40 c         return
  41 c         end
  42 c
  43 c         the value of iflag should not be changed by fcn unless
  44 c         the user wants to terminate execution of lmstr1.
  45 c         in this case set iflag to a negative integer.
  46 c
  47 c       m is a positive integer input variable set to the number
  48 c         of functions.
  49 c
  50 c       n is a positive integer input variable set to the number
  51 c         of variables. n must not exceed m.
  52 c
  53 c       x is an array of length n. on input x must contain
  54 c         an initial estimate of the solution vector. on output x
  55 c         contains the final estimate of the solution vector.
  56 c
  57 c       fvec is an output array of length m which contains
  58 c         the functions evaluated at the output x.
  59 c
  60 c       fjac is an output n by n array. the upper triangle of fjac
  61 c         contains an upper triangular matrix r such that
  62 c
  63 c                t     t           t
  64 c               p *(jac *jac)*p = r *r,
  65 c
  66 c         where p is a permutation matrix and jac is the final
  67 c         calculated jacobian. column j of p is column ipvt(j)
  68 c         (see below) of the identity matrix. the lower triangular
  69 c         part of fjac contains information generated during
  70 c         the computation of r.
  71 c
  72 c       ldfjac is a positive integer input variable not less than n
  73 c         which specifies the leading dimension of the array fjac.
  74 c
  75 c       tol is a nonnegative input variable. termination occurs
  76 c         when the algorithm estimates either that the relative
  77 c         error in the sum of squares is at most tol or that
  78 c         the relative error between x and the solution is at
  79 c         most tol.
  80 c
  81 c       info is an integer output variable. if the user has
  82 c         terminated execution, info is set to the (negative)
  83 c         value of iflag. see description of fcn. otherwise,
  84 c         info is set as follows.
  85 c
  86 c         info = 0  improper input parameters.
  87 c
  88 c         info = 1  algorithm estimates that the relative error
  89 c                   in the sum of squares is at most tol.
  90 c
  91 c         info = 2  algorithm estimates that the relative error
  92 c                   between x and the solution is at most tol.
  93 c
  94 c         info = 3  conditions for info = 1 and info = 2 both hold.
  95 c
  96 c         info = 4  fvec is orthogonal to the columns of the
  97 c                   jacobian to machine precision.
  98 c
  99 c         info = 5  number of calls to fcn with iflag = 1 has
 100 c                   reached 100*(n+1).
 101 c
 102 c         info = 6  tol is too small. no further reduction in
 103 c                   the sum of squares is possible.
 104 c
 105 c         info = 7  tol is too small. no further improvement in
 106 c                   the approximate solution x is possible.
 107 c
 108 c       ipvt is an integer output array of length n. ipvt
 109 c         defines a permutation matrix p such that jac*p = q*r,
 110 c         where jac is the final calculated jacobian, q is
 111 c         orthogonal (not stored), and r is upper triangular.
 112 c         column j of p is column ipvt(j) of the identity matrix.
 113 c
 114 c       wa is a work array of length lwa.
 115 c
 116 c       lwa is a positive integer input variable not less than 5*n+m.
 117 c
 118 c     subprograms called
 119 c
 120 c       user-supplied ...... fcn
 121 c
 122 c       minpack-supplied ... lmstr
 123 c
 124 c     argonne national laboratory. minpack project. march 1980.
 125 c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
 126 c     jorge j. more
 127 c
 128 c     **********
 129       integer maxfev,mode,nfev,njev,nprint
 130       double precision factor,ftol,gtol,xtol,zero
 131       data factor,zero /1.0d2,0.0d0/
 132       info = 0
 133 c
 134 c     check the input parameters for errors.
 135 c
 136       if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. n .or. tol .lt. zero
 137      *    .or. lwa .lt. 5*n + m) go to 10
 138 c
 139 c     call lmstr.
 140 c
 141       maxfev = 100*(n + 1)
 142       ftol = tol
 143       xtol = tol
 144       gtol = zero
 145       mode = 1
 146       nprint = 0
 147       call lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,maxfev,
 148      *           wa(1),mode,factor,nprint,info,nfev,njev,ipvt,wa(n+1),
 149      *           wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
 150       if (info .eq. 8) info = 4
 151    10 continue
 152       return
 153 c
 154 c     last card of subroutine lmstr1.
 155 c
 156       end