share/minpack/fortran/lmder.f

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