subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,stpmin, + stpmax) logical brackt double precision stx, fx, dx, sty, fy, dy, stp, fp, dp, stpmin, + stpmax c ********** c c Subroutine dcstep c c This subroutine computes a safeguarded step for a search c procedure and updates an interval that contains a step that c satisfies a sufficient decrease and a curvature condition. c c The parameter stx contains the step with the least function c value. If brackt is set to .true. then a minimizer has c been bracketed in an interval with endpoints stx and sty. c The parameter stp contains the current step. c The subroutine assumes that if brackt is set to .true. then c c min(stx,sty) < stp < max(stx,sty), c c and that the derivative at stx is negative in the direction c of the step. c c The subroutine statement is c c subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, c stpmin,stpmax) c c where c c stx is a double precision variable. c On entry stx is the best step obtained so far and is an c endpoint of the interval that contains the minimizer. c On exit stx is the updated best step. c c fx is a double precision variable. c On entry fx is the function at stx. c On exit fx is the function at stx. c c dx is a double precision variable. c On entry dx is the derivative of the function at c stx. The derivative must be negative in the direction of c the step, that is, dx and stp - stx must have opposite c signs. c On exit dx is the derivative of the function at stx. c c sty is a double precision variable. c On entry sty is the second endpoint of the interval that c contains the minimizer. c On exit sty is the updated endpoint of the interval that c contains the minimizer. c c fy is a double precision variable. c On entry fy is the function at sty. c On exit fy is the function at sty. c c dy is a double precision variable. c On entry dy is the derivative of the function at sty. c On exit dy is the derivative of the function at the exit sty. c c stp is a double precision variable. c On entry stp is the current step. If brackt is set to .true. c then on input stp must be between stx and sty. c On exit stp is a new trial step. c c fp is a double precision variable. c On entry fp is the function at stp c On exit fp is unchanged. c c dp is a double precision variable. c On entry dp is the the derivative of the function at stp. c On exit dp is unchanged. c c brackt is an logical variable. c On entry brackt specifies if a minimizer has been bracketed. c Initially brackt must be set to .false. c On exit brackt specifies if a minimizer has been bracketed. c When a minimizer is bracketed brackt is set to .true. c c stpmin is a double precision variable. c On entry stpmin is a lower bound for the step. c On exit stpmin is unchanged. c c stpmax is a double precision variable. c On entry stpmax is an upper bound for the step. c On exit stpmax is unchanged. c c MINPACK-1 Project. June 1983 c Argonne National Laboratory. c Jorge J. More' and David J. Thuente. c c MINPACK-2 Project. November 1993. c Argonne National Laboratory and University of Minnesota. c Brett M. Averick and Jorge J. More'. c c ********** double precision zero, p66, two, three parameter (zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0) double precision gamma, p, q, r, s, sgnd, stpc, stpf, stpq, theta sgnd = dp*(dx/abs(dx)) c First case: A higher function value. The minimum is bracketed. c If the cubic step is closer to stx than the quadratic step, the c cubic step is taken, otherwise the average of the cubic and c quadratic steps is taken. if (fp .gt. fx) then theta = three*(fx-fp)/(stp-stx) + dx + dp s = max(abs(theta),abs(dx),abs(dp)) gamma = s*sqrt((theta/s)**2-(dx/s)*(dp/s)) if (stp .lt. stx) gamma = -gamma p = (gamma-dx) + theta q = ((gamma-dx)+gamma) + dp r = p/q stpc = stx + r*(stp-stx) stpq = stx + ((dx/((fx-fp)/(stp-stx)+dx))/two)*(stp-stx) if (abs(stpc-stx) .lt. abs(stpq-stx)) then stpf = stpc else stpf = stpc + (stpq-stpc)/two end if brackt = .true. c Second case: A lower function value and derivatives of opposite c sign. The minimum is bracketed. If the cubic step is farther from c stp than the secant step, the cubic step is taken, otherwise the c secant step is taken. else if (sgnd .lt. zero) then theta = three*(fx-fp)/(stp-stx) + dx + dp s = max(abs(theta),abs(dx),abs(dp)) gamma = s*sqrt((theta/s)**2-(dx/s)*(dp/s)) if (stp .gt. stx) gamma = -gamma p = (gamma-dp) + theta q = ((gamma-dp)+gamma) + dx r = p/q stpc = stp + r*(stx-stp) stpq = stp + (dp/(dp-dx))*(stx-stp) if (abs(stpc-stp) .gt. abs(stpq-stp)) then stpf = stpc else stpf = stpq end if brackt = .true. c Third case: A lower function value, derivatives of the same sign, c and the magnitude of the derivative decreases. else if (abs(dp) .lt. abs(dx)) then c The cubic step is computed only if the cubic tends to infinity c in the direction of the step or if the minimum of the cubic c is beyond stp. Otherwise the cubic step is defined to be the c secant step. theta = three*(fx-fp)/(stp-stx) + dx + dp s = max(abs(theta),abs(dx),abs(dp)) c The case gamma = 0 only arises if the cubic does not tend c to infinity in the direction of the step. gamma = s*sqrt(max(zero,(theta/s)**2-(dx/s)*(dp/s))) if (stp .gt. stx) gamma = -gamma p = (gamma-dp) + theta q = (gamma+(dx-dp)) + gamma r = p/q if (r .lt. zero .and. gamma .ne. zero) then stpc = stp + r*(stx-stp) else if (stp .gt. stx) then stpc = stpmax else stpc = stpmin end if stpq = stp + (dp/(dp-dx))*(stx-stp) if (brackt) then c A minimizer has been bracketed. If the cubic step is c closer to stp than the secant step, the cubic step is c taken, otherwise the secant step is taken. if (abs(stpc-stp) .lt. abs(stpq-stp)) then stpf = stpc else stpf = stpq end if if (stp .gt. stx) then stpf = min(stp+p66*(sty-stp),stpf) else stpf = max(stp+p66*(sty-stp),stpf) end if else c A minimizer has not been bracketed. If the cubic step is c farther from stp than the secant step, the cubic step is c taken, otherwise the secant step is taken. if (abs(stpc-stp) .gt. abs(stpq-stp)) then stpf = stpc else stpf = stpq end if stpf = min(stpmax,stpf) stpf = max(stpmin,stpf) end if c Fourth case: A lower function value, derivatives of the same sign, c and the magnitude of the derivative does not decrease. If the c minimum is not bracketed, the step is either stpmin or stpmax, c otherwise the cubic step is taken. else if (brackt) then theta = three*(fp-fy)/(sty-stp) + dy + dp s = max(abs(theta),abs(dy),abs(dp)) gamma = s*sqrt((theta/s)**2-(dy/s)*(dp/s)) if (stp .gt. sty) gamma = -gamma p = (gamma-dp) + theta q = ((gamma-dp)+gamma) + dy r = p/q stpc = stp + r*(sty-stp) stpf = stpc else if (stp .gt. stx) then stpf = stpmax else stpf = stpmin end if end if c Update the interval which contains a minimizer. if (fp .gt. fx) then sty = stp fy = fp dy = dp else if (sgnd .lt. zero) then sty = stx fy = fx dy = dx end if stx = stp fx = fp dx = dp end if c Compute the new step. stp = stpf end