share/minpack/fortran/hybrd.f

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