SUBROUTINE zero_search(bb_0, cc_0, ans, fn_id, n_iuser0, iuser0, n_ruser0, ruser0, ifail) !************************************************************************************* ! ! Linear evaluation of zero based upon the Bus-Dekker algorithm ! ! Bus, J.C.P. & Dekker, T.J. (1975), "Two Efficient Algorithms with Guaranteed ! Convergenece for finding a Zero of a Fuction", ACM Trans. on Math. Software (TOMS), ! 1, pp. 330-345. ! as in NAG library function C05ADF ! ! INPUTS: ! bb_0 = lower bound ! cc_0 = upper bound ! fn_id = internal identification number for function to evaluate zero for ! n_iuser0 = dimension of vector of integer parameter inputs fed into routine ! iuser0 = vector of integer parameter inputs fed into routine ! n_ruser0 = dimension of vector of real parameter inputs fed into routine ! ruser0 = vector of real parameter inputs fed into routine ! OUTPUTS: ! ans = solution ! ifail = flag to indicate errors ! = 0 solution found ! = 1 function doesn't cross zero in interval [bb_0,cc_0] ! = 2 too much accuracy ! = 3 potential pole ! ! subroutine fn_zero passes test point to function for evaluation ! ! code developed for the SIDD microsimulation model ! ! Author: J. van de Ven ! Date: 02 Aug 2016 ! !************************************************************************************* USE global_ParamStore USE global_CommonAll IMPLICIT none integer(4), intent(in) :: fn_id, n_iuser0, n_ruser0 integer(4) :: iuser0(n_iuser0), ifail real(8), intent(in) :: bb_0, cc_0 real(8) :: ruser0(n_ruser0) real(8), intent(out) :: ans ! local variables integer(4) :: stepType real(8) :: tol, tol2, bb, cc, b2, c2, pp, qq, fb, fc, bm1, fbm1, bm2, midPoint, deltaB, temp, fbm2, fbm2b, fbm2bm1 logical :: lgcal, lgcal2 !************************************** ! begin code !************************************** !!!!! The methods used are a combination of two types of rational extrapolation, !!!!! linear interpolation and a bisection algorithm: !!!! stepType=1 bisect !!!! stepType=2 linear1 !!!! stepType=3 linear2 !!!! stepType=4 rational stepType = 1_4 tol = 1.0E-7 ! initialise search values bb = bb_0 cc = cc_0 call fn_zero(bb, fb, fn_id, n_iuser0, iuser0, n_ruser0, ruser0) call fn_zero(cc, fc, fn_id, n_iuser0, iuser0, n_ruser0, ruser0) ! bm1 is previous value of b; this is used in extrapolation ! as the second point (along with b); on entry we want to ! interpolate bm1 = cc; fbm1 = fc; ! Interpolation is just a special case of extrapolation ! where the initial endpoints are used as the extrapolands; ! hence interpolation will be described as LINEAR1; if ( (fb*fc).gt.0.0 ) then tol2 = tol tol = 2*tol*abs(bb) midPoint = bb lgcal2 = .true. ifail = 1 else tol2 = tol tol = 2*tol*abs(bb) midPoint = (cc+bb)/2.0 deltaB = midPoint - bb lgcal2 = .false. ifail = 0 end if do while( abs(midPoint-bb).gt.tol ) if ( stepType.ne.1 ) then ! Extrapolate routines form a sequence of values of the best ! endpoint which aim to converge while keeping the opposite sign ! endpoint fixed. if ( abs(fc).lt.abs(fb) ) then ! We are not close enough to the root and ! have encountered a minimum in the residual ! we mustn't extrapolate as this would take us away from the root ! so we interpolate using the interval endpoints. temp = bb bb = cc cc = temp temp = fb fb = fc fc = temp deltaB = midPoint - bb ! This only matters in RATIONAL_EXTRAPOLATE if ( bm1.ne.bb ) then ! we need to keep a distinct values of b, bm1 and bm2 bm2 = bm1 fbm2 = fbm1 end if ! This matters in LINEAR2 or RATIONAL_EXTRAPOLATE bm1 = cc fbm1 = fc end if if ( deltaB.ne.0.0 ) then tol = tol*deltaB/abs(deltaB) ! so b+tol is inside interval [b,c] else tol = 0.0 end if pp = (bb-bm1)*fb if ( stepType.le.3 ) then ! Not enough points to do a ! rational extrapolation qq = fbm1 - fb ! so b+p/q= b +(b-bm1)*fb/(fb-fbm1) else ! rational extrapolation from ! pts b, bm1 and bm2 fbm2b = (fbm2-fb)/(bm2-bb) fbm2bm1 = (fbm2-fbm1)/(bm2-bm1) pp = fbm2bm1*pp qq = fbm2b*fbm1-fbm2bm1*fb end if ! Decide if we want to move only sign of p/ ! matters so make p>0 to make tests simpler if ( pp.lt.0.0 ) then pp = -pp qq = -qq end if ! default here is bisection value ! change if we prefer extrapolate/interpolated ! value if ( ( pp.eq.0.0 ).or.( pp.le.(qq*tol) ) ) then ! For stability we must make sure that ! the point moves by at least tol each time deltaB = tol elseif ( pp.lt.(qq*deltaB) ) then ! must land in half closer to best value deltaB = pp / qq else ! reject and use bisection , but EXTRAPOLATE NEXT time stepType = 1 end if end if ! end switch ! Shift old b values down stack bm2 = bm1 fbm2 = fbm1 bm1 = bb fbm1 = fb bb = bb + deltaB call fn_zero(bb, fb, fn_id, n_iuser0, iuser0, n_ruser0, ruser0) ! Decide what solution method to use for next step if ( fc.ge.0.0 ) then lgcal = (fb.ge.0.0) else lgcal = (fb.le.0.0) end if if ( lgcal.eq..true.) then ! better than bisection we have reduced the ! interval by more than 0.5 and we replace c ! by old b endpoint cc = bm1 fc = fbm1 ! don't advance state if ( stepType.eq.1 ) then stepType = 2 end if else ! worse than bisect so advance state if ( stepType.eq.4 ) then stepType = 1 else stepType = stepType + 1 end if end if !tol2=tol+tol*abs(b) midPoint = 0.5*(bb+cc) deltaB = midPoint-bb tol = 2*tol2*abs(bb) !max(tol2,tol+tol*abs(b)) end do if ( lgcal2.eq..true. ) then if ( abs(fc).lt.abs(fb) ) then ans = cc else ans = bb end if else ans = 0.5*(bb+cc) end if if (ifail.ne.1_4) then if ( abs(cc-bb).gt.tol ) then ! too much accuracy requested ifail=2 elseif ( abs(fc-fb).gt.(250.0*tol) ) then ! solution at probable pole ifail=3 if ( abs(fc).lt.abs(fb) ) then ! would like to take c, but need to be careful over boundaries if ( abs(cc-c2).lt.tol ) then ! at upper bound, take positive value function evaluation if ( fc.ge.0.0 ) then ! can take fc ans = cc elseif ( fb.ge.0.0 ) then ! can take fb ans = bb else ! we have a problem Houston !$OMP CRITICAL(error_1) print *, 'problem in finding acceptable root' pause 'press enter to exit' stop continue !$OMP END CRITICAL(error_1) end if elseif ( abs(cc-b2).lt.tol ) then ! at lower bound, take negative value function evaluation if ( fc.le.0.0 ) then ! can take fc ans = cc elseif ( fb.le.0.0 ) then ! can take fb ans = bb else ! we have a problem Houston !$OMP CRITICAL(error_2) print *, 'problem in finding acceptable root' pause 'press enter to exit' stop continue !$OMP END CRITICAL(error_2) end if else ans = cc end if else ! would like to take b, but need to be careful over boundaries if ( abs(bb-c2).lt.tol ) then ! at upper bound, take positive value function evaluation if ( fb.ge.0.0 ) then ! can take fb ans = bb elseif ( fc.ge.0.0 ) then ! can take fc ans = cc else ! we have a problem Houston !$OMP CRITICAL(error_4) print *, 'problem in finding acceptable root' pause 'press enter to exit' stop continue !$OMP END CRITICAL(error_4) end if elseif ( abs(bb-b2).lt.tol ) then ! at lower bound, take negative value function evaluation if ( fb.le.0.0 ) then ! can take fb ans = bb elseif ( fc.le.0.0 ) then ! can take fc ans = cc else ! we have a problem Houston !$OMP CRITICAL(error_5) print *, 'problem in finding acceptable root' pause 'press enter to exit' stop continue !$OMP END CRITICAL(error_5) end if else ans = bb end if end if end if end if END SUBROUTINE zero_search