share/minpack/fortran/hybrj.f

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