share/minpack/fortran/lmdif.f

   1       subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
   2      *                 diag,mode,factor,nprint,info,nfev,fjac,ldfjac,
   3      *                 ipvt,qtf,wa1,wa2,wa3,wa4)
   4       integer m,n,maxfev,mode,nprint,info,nfev,ldfjac
   5       integer ipvt(n)
   6       double precision ftol,xtol,gtol,epsfcn,factor
   7       double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n),
   8      *                 wa1(n),wa2(n),wa3(n),wa4(m)
   9       external fcn
  10 c     **********
  11 c
  12 c     subroutine lmdif
  13 c
  14 c     the purpose of lmdif is to minimize the sum of the squares of
  15 c     m nonlinear functions in n variables by a modification of
  16 c     the levenberg-marquardt algorithm. the user must provide a
  17 c     subroutine which calculates the functions. the jacobian is
  18 c     then calculated by a forward-difference approximation.
  19 c
  20 c     the subroutine statement is
  21 c
  22 c       subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
  23 c                        diag,mode,factor,nprint,info,nfev,fjac,
  24 c                        ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
  25 c
  26 c     where
  27 c
  28 c       fcn is the name of the user-supplied subroutine which
  29 c         calculates the functions. fcn must be declared
  30 c         in an external statement in the user calling
  31 c         program, and should be written as follows.
  32 c
  33 c         subroutine fcn(m,n,x,fvec,iflag)
  34 c         integer m,n,iflag
  35 c         double precision x(n),fvec(m)
  36 c         ----------
  37 c         calculate the functions at x and
  38 c         return this vector in fvec.
  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 lmdif.
  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       ftol is a nonnegative input variable. termination
  61 c         occurs when both the actual and predicted relative
  62 c         reductions in the sum of squares are at most ftol.
  63 c         therefore, ftol measures the relative error desired
  64 c         in the sum of squares.
  65 c
  66 c       xtol is a nonnegative input variable. termination
  67 c         occurs when the relative error between two consecutive
  68 c         iterates is at most xtol. therefore, xtol measures the
  69 c         relative error desired in the approximate solution.
  70 c
  71 c       gtol is a nonnegative input variable. termination
  72 c         occurs when the cosine of the angle between fvec and
  73 c         any column of the jacobian is at most gtol in absolute
  74 c         value. therefore, gtol measures the orthogonality
  75 c         desired between the function vector and the columns
  76 c         of the jacobian.
  77 c
  78 c       maxfev is a positive integer input variable. termination
  79 c         occurs when the number of calls to fcn is at least
  80 c         maxfev by the end of an iteration.
  81 c
  82 c       epsfcn is an input variable used in determining a suitable
  83 c         step length for the forward-difference approximation. this
  84 c         approximation assumes that the relative errors in the
  85 c         functions are of the order of epsfcn. if epsfcn is less
  86 c         than the machine precision, it is assumed that the relative
  87 c         errors in the functions are of the order of the machine
  88 c         precision.
  89 c
  90 c       diag is an array of length n. if mode = 1 (see
  91 c         below), diag is internally set. if mode = 2, diag
  92 c         must contain positive entries that serve as
  93 c         multiplicative scale factors for the variables.
  94 c
  95 c       mode is an integer input variable. if mode = 1, the
  96 c         variables will be scaled internally. if mode = 2,
  97 c         the scaling is specified by the input diag. other
  98 c         values of mode are equivalent to mode = 1.
  99 c
 100 c       factor is a positive input variable used in determining the
 101 c         initial step bound. this bound is set to the product of
 102 c         factor and the euclidean norm of diag*x if nonzero, or else
 103 c         to factor itself. in most cases factor should lie in the
 104 c         interval (.1,100.). 100. is a generally recommended value.
 105 c
 106 c       nprint is an integer input variable that enables controlled
 107 c         printing of iterates if it is positive. in this case,
 108 c         fcn is called with iflag = 0 at the beginning of the first
 109 c         iteration and every nprint iterations thereafter and
 110 c         immediately prior to return, with x and fvec available
 111 c         for printing. if nprint is not positive, no special calls
 112 c         of fcn with iflag = 0 are made.
 113 c
 114 c       info is an integer output variable. if the user has
 115 c         terminated execution, info is set to the (negative)
 116 c         value of iflag. see description of fcn. otherwise,
 117 c         info is set as follows.
 118 c
 119 c         info = 0  improper input parameters.
 120 c
 121 c         info = 1  both actual and predicted relative reductions
 122 c                   in the sum of squares are at most ftol.
 123 c
 124 c         info = 2  relative error between two consecutive iterates
 125 c                   is at most xtol.
 126 c
 127 c         info = 3  conditions for info = 1 and info = 2 both hold.
 128 c
 129 c         info = 4  the cosine of the angle between fvec and any
 130 c                   column of the jacobian is at most gtol in
 131 c                   absolute value.
 132 c
 133 c         info = 5  number of calls to fcn has reached or
 134 c                   exceeded maxfev.
 135 c
 136 c         info = 6  ftol is too small. no further reduction in
 137 c                   the sum of squares is possible.
 138 c
 139 c         info = 7  xtol is too small. no further improvement in
 140 c                   the approximate solution x is possible.
 141 c
 142 c         info = 8  gtol is too small. fvec is orthogonal to the
 143 c                   columns of the jacobian to machine precision.
 144 c
 145 c       nfev is an integer output variable set to the number of
 146 c         calls to fcn.
 147 c
 148 c       fjac is an output m by n array. the upper n by n submatrix
 149 c         of fjac contains an upper triangular matrix r with
 150 c         diagonal elements of nonincreasing magnitude such that
 151 c
 152 c                t     t           t
 153 c               p *(jac *jac)*p = r *r,
 154 c
 155 c         where p is a permutation matrix and jac is the final
 156 c         calculated jacobian. column j of p is column ipvt(j)
 157 c         (see below) of the identity matrix. the lower trapezoidal
 158 c         part of fjac contains information generated during
 159 c         the computation of r.
 160 c
 161 c       ldfjac is a positive integer input variable not less than m
 162 c         which specifies the leading dimension of the array fjac.
 163 c
 164 c       ipvt is an integer output array of length n. ipvt
 165 c         defines a permutation matrix p such that jac*p = q*r,
 166 c         where jac is the final calculated jacobian, q is
 167 c         orthogonal (not stored), and r is upper triangular
 168 c         with diagonal elements of nonincreasing magnitude.
 169 c         column j of p is column ipvt(j) of the identity matrix.
 170 c
 171 c       qtf is an output array of length n which contains
 172 c         the first n elements of the vector (q transpose)*fvec.
 173 c
 174 c       wa1, wa2, and wa3 are work arrays of length n.
 175 c
 176 c       wa4 is a work array of length m.
 177 c
 178 c     subprograms called
 179 c
 180 c       user-supplied ...... fcn
 181 c
 182 c       minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
 183 c
 184 c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
 185 c
 186 c     argonne national laboratory. minpack project. march 1980.
 187 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
 188 c
 189 c     **********
 190       integer i,iflag,iter,j,l
 191       double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,
 192      *                 one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio,
 193      *                 sum,temp,temp1,temp2,xnorm,zero
 194       double precision dpmpar,enorm
 195       data one,p1,p5,p25,p75,p0001,zero
 196      *     /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
 197 c
 198 c     epsmch is the machine precision.
 199 c
 200       epsmch = dpmpar(1)
 201 c
 202       info = 0
 203       iflag = 0
 204       nfev = 0
 205 c
 206 c     check the input parameters for errors.
 207 c
 208       if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
 209      *    .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
 210      *    .or. maxfev .le. 0 .or. factor .le. zero) go to 300
 211       if (mode .ne. 2) go to 20
 212       do 10 j = 1, n
 213          if (diag(j) .le. zero) go to 300
 214    10    continue
 215    20 continue
 216 c
 217 c     evaluate the function at the starting point
 218 c     and calculate its norm.
 219 c
 220       iflag = 1
 221       call fcn(m,n,x,fvec,iflag)
 222       nfev = 1
 223       if (iflag .lt. 0) go to 300
 224       fnorm = enorm(m,fvec)
 225 c
 226 c     initialize levenberg-marquardt parameter and iteration counter.
 227 c
 228       par = zero
 229       iter = 1
 230 c
 231 c     beginning of the outer loop.
 232 c
 233    30 continue
 234 c
 235 c        calculate the jacobian matrix.
 236 c
 237          iflag = 2
 238          call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
 239          nfev = nfev + n
 240          if (iflag .lt. 0) go to 300
 241 c
 242 c        if requested, call fcn to enable printing of iterates.
 243 c
 244          if (nprint .le. 0) go to 40
 245          iflag = 0
 246          if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag)
 247          if (iflag .lt. 0) go to 300
 248    40    continue
 249 c
 250 c        compute the qr factorization of the jacobian.
 251 c
 252          call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
 253 c
 254 c        on the first iteration and if mode is 1, scale according
 255 c        to the norms of the columns of the initial jacobian.
 256 c
 257          if (iter .ne. 1) go to 80
 258          if (mode .eq. 2) go to 60
 259          do 50 j = 1, n
 260             diag(j) = wa2(j)
 261             if (wa2(j) .eq. zero) diag(j) = one
 262    50       continue
 263    60    continue
 264 c
 265 c        on the first iteration, calculate the norm of the scaled x
 266 c        and initialize the step bound delta.
 267 c
 268          do 70 j = 1, n
 269             wa3(j) = diag(j)*x(j)
 270    70       continue
 271          xnorm = enorm(n,wa3)
 272          delta = factor*xnorm
 273          if (delta .eq. zero) delta = factor
 274    80    continue
 275 c
 276 c        form (q transpose)*fvec and store the first n components in
 277 c        qtf.
 278 c
 279          do 90 i = 1, m
 280             wa4(i) = fvec(i)
 281    90       continue
 282          do 130 j = 1, n
 283             if (fjac(j,j) .eq. zero) go to 120
 284             sum = zero
 285             do 100 i = j, m
 286                sum = sum + fjac(i,j)*wa4(i)
 287   100          continue
 288             temp = -sum/fjac(j,j)
 289             do 110 i = j, m
 290                wa4(i) = wa4(i) + fjac(i,j)*temp
 291   110          continue
 292   120       continue
 293             fjac(j,j) = wa1(j)
 294             qtf(j) = wa4(j)
 295   130       continue
 296 c
 297 c        compute the norm of the scaled gradient.
 298 c
 299          gnorm = zero
 300          if (fnorm .eq. zero) go to 170
 301          do 160 j = 1, n
 302             l = ipvt(j)
 303             if (wa2(l) .eq. zero) go to 150
 304             sum = zero
 305             do 140 i = 1, j
 306                sum = sum + fjac(i,j)*(qtf(i)/fnorm)
 307   140          continue
 308             gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
 309   150       continue
 310   160       continue
 311   170    continue
 312 c
 313 c        test for convergence of the gradient norm.
 314 c
 315          if (gnorm .le. gtol) info = 4
 316          if (info .ne. 0) go to 300
 317 c
 318 c        rescale if necessary.
 319 c
 320          if (mode .eq. 2) go to 190
 321          do 180 j = 1, n
 322             diag(j) = dmax1(diag(j),wa2(j))
 323   180       continue
 324   190    continue
 325 c
 326 c        beginning of the inner loop.
 327 c
 328   200    continue
 329 c
 330 c           determine the levenberg-marquardt parameter.
 331 c
 332             call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
 333      *                 wa3,wa4)
 334 c
 335 c           store the direction p and x + p. calculate the norm of p.
 336 c
 337             do 210 j = 1, n
 338                wa1(j) = -wa1(j)
 339                wa2(j) = x(j) + wa1(j)
 340                wa3(j) = diag(j)*wa1(j)
 341   210          continue
 342             pnorm = enorm(n,wa3)
 343 c
 344 c           on the first iteration, adjust the initial step bound.
 345 c
 346             if (iter .eq. 1) delta = dmin1(delta,pnorm)
 347 c
 348 c           evaluate the function at x + p and calculate its norm.
 349 c
 350             iflag = 1
 351             call fcn(m,n,wa2,wa4,iflag)
 352             nfev = nfev + 1
 353             if (iflag .lt. 0) go to 300
 354             fnorm1 = enorm(m,wa4)
 355 c
 356 c           compute the scaled actual reduction.
 357 c
 358             actred = -one
 359             if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
 360 c
 361 c           compute the scaled predicted reduction and
 362 c           the scaled directional derivative.
 363 c
 364             do 230 j = 1, n
 365                wa3(j) = zero
 366                l = ipvt(j)
 367                temp = wa1(l)
 368                do 220 i = 1, j
 369                   wa3(i) = wa3(i) + fjac(i,j)*temp
 370   220             continue
 371   230          continue
 372             temp1 = enorm(n,wa3)/fnorm
 373             temp2 = (dsqrt(par)*pnorm)/fnorm
 374             prered = temp1**2 + temp2**2/p5
 375             dirder = -(temp1**2 + temp2**2)
 376 c
 377 c           compute the ratio of the actual to the predicted
 378 c           reduction.
 379 c
 380             ratio = zero
 381             if (prered .ne. zero) ratio = actred/prered
 382 c
 383 c           update the step bound.
 384 c
 385             if (ratio .gt. p25) go to 240
 386                if (actred .ge. zero) temp = p5
 387                if (actred .lt. zero)
 388      *            temp = p5*dirder/(dirder + p5*actred)
 389                if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
 390                delta = temp*dmin1(delta,pnorm/p1)
 391                par = par/temp
 392                go to 260
 393   240       continue
 394                if (par .ne. zero .and. ratio .lt. p75) go to 250
 395                delta = pnorm/p5
 396                par = p5*par
 397   250          continue
 398   260       continue
 399 c
 400 c           test for successful iteration.
 401 c
 402             if (ratio .lt. p0001) go to 290
 403 c
 404 c           successful iteration. update x, fvec, and their norms.
 405 c
 406             do 270 j = 1, n
 407                x(j) = wa2(j)
 408                wa2(j) = diag(j)*x(j)
 409   270          continue
 410             do 280 i = 1, m
 411                fvec(i) = wa4(i)
 412   280          continue
 413             xnorm = enorm(n,wa2)
 414             fnorm = fnorm1
 415             iter = iter + 1
 416   290       continue
 417 c
 418 c           tests for convergence.
 419 c
 420             if (dabs(actred) .le. ftol .and. prered .le. ftol
 421      *          .and. p5*ratio .le. one) info = 1
 422             if (delta .le. xtol*xnorm) info = 2
 423             if (dabs(actred) .le. ftol .and. prered .le. ftol
 424      *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
 425             if (info .ne. 0) go to 300
 426 c
 427 c           tests for termination and stringent tolerances.
 428 c
 429             if (nfev .ge. maxfev) info = 5
 430             if (dabs(actred) .le. epsmch .and. prered .le. epsmch
 431      *          .and. p5*ratio .le. one) info = 6
 432             if (delta .le. epsmch*xnorm) info = 7
 433             if (gnorm .le. epsmch) info = 8
 434             if (info .ne. 0) go to 300
 435 c
 436 c           end of the inner loop. repeat if iteration unsuccessful.
 437 c
 438             if (ratio .lt. p0001) go to 200
 439 c
 440 c        end of the outer loop.
 441 c
 442          go to 30
 443   300 continue
 444 c
 445 c     termination, either normal or user imposed.
 446 c
 447       if (iflag .lt. 0) info = iflag
 448       iflag = 0
 449       if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
 450       return
 451 c
 452 c     last card of subroutine lmdif.
 453 c
 454       end