share/minpack/fortran/lmstr.f

   1       subroutine lmstr(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       logical sing
   7       double precision ftol,xtol,gtol,factor
   8       double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n),
   9      *                 wa1(n),wa2(n),wa3(n),wa4(m)
  10 c     **********
  11 c
  12 c     subroutine lmstr
  13 c
  14 c     the purpose of lmstr 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 which uses minimal storage.
  17 c     the user must provide a subroutine which calculates the
  18 c     functions and the rows of the jacobian.
  19 c
  20 c     the subroutine statement is
  21 c
  22 c       subroutine lmstr(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
  23 c                        maxfev,diag,mode,factor,nprint,info,nfev,
  24 c                        njev,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 and the rows of the jacobian.
  30 c         fcn must be declared in an external statement in the
  31 c         user calling program, and should be written as follows.
  32 c
  33 c         subroutine fcn(m,n,x,fvec,fjrow,iflag)
  34 c         integer m,n,iflag
  35 c         double precision x(n),fvec(m),fjrow(n)
  36 c         ----------
  37 c         if iflag = 1 calculate the functions at x and
  38 c         return this vector in fvec.
  39 c         if iflag = i calculate the (i-1)-st row of the
  40 c         jacobian at x and return this vector in fjrow.
  41 c         ----------
  42 c         return
  43 c         end
  44 c
  45 c         the value of iflag should not be changed by fcn unless
  46 c         the user wants to terminate execution of lmstr.
  47 c         in this case set iflag to a negative integer.
  48 c
  49 c       m is a positive integer input variable set to the number
  50 c         of functions.
  51 c
  52 c       n is a positive integer input variable set to the number
  53 c         of variables. n must not exceed m.
  54 c
  55 c       x is an array of length n. on input x must contain
  56 c         an initial estimate of the solution vector. on output x
  57 c         contains the final estimate of the solution vector.
  58 c
  59 c       fvec is an output array of length m which contains
  60 c         the functions evaluated at the output x.
  61 c
  62 c       fjac is an output n by n array. the upper triangle of fjac
  63 c         contains an upper triangular matrix r such that
  64 c
  65 c                t     t           t
  66 c               p *(jac *jac)*p = r *r,
  67 c
  68 c         where p is a permutation matrix and jac is the final
  69 c         calculated jacobian. column j of p is column ipvt(j)
  70 c         (see below) of the identity matrix. the lower triangular
  71 c         part of fjac contains information generated during
  72 c         the computation of r.
  73 c
  74 c       ldfjac is a positive integer input variable not less than n
  75 c         which specifies the leading dimension of the array fjac.
  76 c
  77 c       ftol is a nonnegative input variable. termination
  78 c         occurs when both the actual and predicted relative
  79 c         reductions in the sum of squares are at most ftol.
  80 c         therefore, ftol measures the relative error desired
  81 c         in the sum of squares.
  82 c
  83 c       xtol is a nonnegative input variable. termination
  84 c         occurs when the relative error between two consecutive
  85 c         iterates is at most xtol. therefore, xtol measures the
  86 c         relative error desired in the approximate solution.
  87 c
  88 c       gtol is a nonnegative input variable. termination
  89 c         occurs when the cosine of the angle between fvec and
  90 c         any column of the jacobian is at most gtol in absolute
  91 c         value. therefore, gtol measures the orthogonality
  92 c         desired between the function vector and the columns
  93 c         of the jacobian.
  94 c
  95 c       maxfev is a positive integer input variable. termination
  96 c         occurs when the number of calls to fcn with iflag = 1
  97 c         has reached maxfev.
  98 c
  99 c       diag is an array of length n. if mode = 1 (see
 100 c         below), diag is internally set. if mode = 2, diag
 101 c         must contain positive entries that serve as
 102 c         multiplicative scale factors for the variables.
 103 c
 104 c       mode is an integer input variable. if mode = 1, the
 105 c         variables will be scaled internally. if mode = 2,
 106 c         the scaling is specified by the input diag. other
 107 c         values of mode are equivalent to mode = 1.
 108 c
 109 c       factor is a positive input variable used in determining the
 110 c         initial step bound. this bound is set to the product of
 111 c         factor and the euclidean norm of diag*x if nonzero, or else
 112 c         to factor itself. in most cases factor should lie in the
 113 c         interval (.1,100.). 100. is a generally recommended value.
 114 c
 115 c       nprint is an integer input variable that enables controlled
 116 c         printing of iterates if it is positive. in this case,
 117 c         fcn is called with iflag = 0 at the beginning of the first
 118 c         iteration and every nprint iterations thereafter and
 119 c         immediately prior to return, with x and fvec available
 120 c         for printing. 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         column j of p is column ipvt(j) of the identity matrix.
 165 c
 166 c       qtf is an output array of length n which contains
 167 c         the first n elements of the vector (q transpose)*fvec.
 168 c
 169 c       wa1, wa2, and wa3 are work arrays of length n.
 170 c
 171 c       wa4 is a work array of length m.
 172 c
 173 c     subprograms called
 174 c
 175 c       user-supplied ...... fcn
 176 c
 177 c       minpack-supplied ... dpmpar,enorm,lmpar,qrfac,rwupdt
 178 c
 179 c       fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
 180 c
 181 c     argonne national laboratory. minpack project. march 1980.
 182 c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
 183 c     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. n
 206      *    .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
 207      *    .or. maxfev .le. 0 .or. factor .le. zero) go to 340
 208       if (mode .ne. 2) go to 20
 209       do 10 j = 1, n
 210          if (diag(j) .le. zero) go to 340
 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,wa3,iflag)
 219       nfev = 1
 220       if (iflag .lt. 0) go to 340
 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        if requested, call fcn to enable printing of iterates.
 233 c
 234          if (nprint .le. 0) go to 40
 235          iflag = 0
 236          if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,wa3,iflag)
 237          if (iflag .lt. 0) go to 340
 238    40    continue
 239 c
 240 c        compute the qr factorization of the jacobian matrix
 241 c        calculated one row at a time, while simultaneously
 242 c        forming (q transpose)*fvec and storing the first
 243 c        n components in qtf.
 244 c
 245          do 60 j = 1, n
 246             qtf(j) = zero
 247             do 50 i = 1, n
 248                fjac(i,j) = zero
 249    50          continue
 250    60       continue
 251          iflag = 2
 252          do 70 i = 1, m
 253             call fcn(m,n,x,fvec,wa3,iflag)
 254             if (iflag .lt. 0) go to 340
 255             temp = fvec(i)
 256             call rwupdt(n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2)
 257             iflag = iflag + 1
 258    70       continue
 259          njev = njev + 1
 260 c
 261 c        if the jacobian is rank deficient, call qrfac to
 262 c        reorder its columns and update the components of qtf.
 263 c
 264          sing = .false.
 265          do 80 j = 1, n
 266             if (fjac(j,j) .eq. zero) sing = .true.
 267             ipvt(j) = j
 268             wa2(j) = enorm(j,fjac(1,j))
 269    80       continue
 270          if (.not.sing) go to 130
 271          call qrfac(n,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
 272          do 120 j = 1, n
 273             if (fjac(j,j) .eq. zero) go to 110
 274             sum = zero
 275             do 90 i = j, n
 276                sum = sum + fjac(i,j)*qtf(i)
 277    90          continue
 278             temp = -sum/fjac(j,j)
 279             do 100 i = j, n
 280                qtf(i) = qtf(i) + fjac(i,j)*temp
 281   100          continue
 282   110       continue
 283             fjac(j,j) = wa1(j)
 284   120       continue
 285   130    continue
 286 c
 287 c        on the first iteration and if mode is 1, scale according
 288 c        to the norms of the columns of the initial jacobian.
 289 c
 290          if (iter .ne. 1) go to 170
 291          if (mode .eq. 2) go to 150
 292          do 140 j = 1, n
 293             diag(j) = wa2(j)
 294             if (wa2(j) .eq. zero) diag(j) = one
 295   140       continue
 296   150    continue
 297 c
 298 c        on the first iteration, calculate the norm of the scaled x
 299 c        and initialize the step bound delta.
 300 c
 301          do 160 j = 1, n
 302             wa3(j) = diag(j)*x(j)
 303   160       continue
 304          xnorm = enorm(n,wa3)
 305          delta = factor*xnorm
 306          if (delta .eq. zero) delta = factor
 307   170    continue
 308 c
 309 c        compute the norm of the scaled gradient.
 310 c
 311          gnorm = zero
 312          if (fnorm .eq. zero) go to 210
 313          do 200 j = 1, n
 314             l = ipvt(j)
 315             if (wa2(l) .eq. zero) go to 190
 316             sum = zero
 317             do 180 i = 1, j
 318                sum = sum + fjac(i,j)*(qtf(i)/fnorm)
 319   180          continue
 320             gnorm = dmax1(gnorm,dabs(sum/wa2(l)))
 321   190       continue
 322   200       continue
 323   210    continue
 324 c
 325 c        test for convergence of the gradient norm.
 326 c
 327          if (gnorm .le. gtol) info = 4
 328          if (info .ne. 0) go to 340
 329 c
 330 c        rescale if necessary.
 331 c
 332          if (mode .eq. 2) go to 230
 333          do 220 j = 1, n
 334             diag(j) = dmax1(diag(j),wa2(j))
 335   220       continue
 336   230    continue
 337 c
 338 c        beginning of the inner loop.
 339 c
 340   240    continue
 341 c
 342 c           determine the levenberg-marquardt parameter.
 343 c
 344             call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
 345      *                 wa3,wa4)
 346 c
 347 c           store the direction p and x + p. calculate the norm of p.
 348 c
 349             do 250 j = 1, n
 350                wa1(j) = -wa1(j)
 351                wa2(j) = x(j) + wa1(j)
 352                wa3(j) = diag(j)*wa1(j)
 353   250          continue
 354             pnorm = enorm(n,wa3)
 355 c
 356 c           on the first iteration, adjust the initial step bound.
 357 c
 358             if (iter .eq. 1) delta = dmin1(delta,pnorm)
 359 c
 360 c           evaluate the function at x + p and calculate its norm.
 361 c
 362             iflag = 1
 363             call fcn(m,n,wa2,wa4,wa3,iflag)
 364             nfev = nfev + 1
 365             if (iflag .lt. 0) go to 340
 366             fnorm1 = enorm(m,wa4)
 367 c
 368 c           compute the scaled actual reduction.
 369 c
 370             actred = -one
 371             if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
 372 c
 373 c           compute the scaled predicted reduction and
 374 c           the scaled directional derivative.
 375 c
 376             do 270 j = 1, n
 377                wa3(j) = zero
 378                l = ipvt(j)
 379                temp = wa1(l)
 380                do 260 i = 1, j
 381                   wa3(i) = wa3(i) + fjac(i,j)*temp
 382   260             continue
 383   270          continue
 384             temp1 = enorm(n,wa3)/fnorm
 385             temp2 = (dsqrt(par)*pnorm)/fnorm
 386             prered = temp1**2 + temp2**2/p5
 387             dirder = -(temp1**2 + temp2**2)
 388 c
 389 c           compute the ratio of the actual to the predicted
 390 c           reduction.
 391 c
 392             ratio = zero
 393             if (prered .ne. zero) ratio = actred/prered
 394 c
 395 c           update the step bound.
 396 c
 397             if (ratio .gt. p25) go to 280
 398                if (actred .ge. zero) temp = p5
 399                if (actred .lt. zero)
 400      *            temp = p5*dirder/(dirder + p5*actred)
 401                if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
 402                delta = temp*dmin1(delta,pnorm/p1)
 403                par = par/temp
 404                go to 300
 405   280       continue
 406                if (par .ne. zero .and. ratio .lt. p75) go to 290
 407                delta = pnorm/p5
 408                par = p5*par
 409   290          continue
 410   300       continue
 411 c
 412 c           test for successful iteration.
 413 c
 414             if (ratio .lt. p0001) go to 330
 415 c
 416 c           successful iteration. update x, fvec, and their norms.
 417 c
 418             do 310 j = 1, n
 419                x(j) = wa2(j)
 420                wa2(j) = diag(j)*x(j)
 421   310          continue
 422             do 320 i = 1, m
 423                fvec(i) = wa4(i)
 424   320          continue
 425             xnorm = enorm(n,wa2)
 426             fnorm = fnorm1
 427             iter = iter + 1
 428   330       continue
 429 c
 430 c           tests for convergence.
 431 c
 432             if (dabs(actred) .le. ftol .and. prered .le. ftol
 433      *          .and. p5*ratio .le. one) info = 1
 434             if (delta .le. xtol*xnorm) info = 2
 435             if (dabs(actred) .le. ftol .and. prered .le. ftol
 436      *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
 437             if (info .ne. 0) go to 340
 438 c
 439 c           tests for termination and stringent tolerances.
 440 c
 441             if (nfev .ge. maxfev) info = 5
 442             if (dabs(actred) .le. epsmch .and. prered .le. epsmch
 443      *          .and. p5*ratio .le. one) info = 6
 444             if (delta .le. epsmch*xnorm) info = 7
 445             if (gnorm .le. epsmch) info = 8
 446             if (info .ne. 0) go to 340
 447 c
 448 c           end of the inner loop. repeat if iteration unsuccessful.
 449 c
 450             if (ratio .lt. p0001) go to 240
 451 c
 452 c        end of the outer loop.
 453 c
 454          go to 30
 455   340 continue
 456 c
 457 c     termination, either normal or user imposed.
 458 c
 459       if (iflag .lt. 0) info = iflag
 460       iflag = 0
 461       if (nprint .gt. 0) call fcn(m,n,x,fvec,wa3,iflag)
 462       return
 463 c
 464 c     last card of subroutine lmstr.
 465 c
 466       end