MODULE Line_Root_Finders USE Precision USE Error_Handling USE System_Monitors IMPLICIT NONE PRIVATE PUBLIC::Line_Root_Finder INTEGER,PARAMETER,PUBLIC::& root_finding_error=-1,root_finding_success=0,root_finding_failure=1 INTEGER,PARAMETER,PUBLIC::& free_root_finder=0,inactive_root_finder=-1,active_root_finder=1 TYPE Line_Root_Finder INTEGER::PIN=1 INTEGER::status=free_root_finder,convergence=root_finding_failure INTEGER::n_iterations=0,max_iterations=-1 REAL(KIND=r_wp),DIMENSION(2)::interval=(/0.0_r_wp,1.0_r_wp/) REAL(KIND=r_wp)::root=1.0_r_wp,tolerance=-1.00_r_wp CHARACTER::method='N' ENDTYPE END MODULE Line_Root_Finders MODULE Dual_Line_Minimizers USE Precision USE Error_Handling USE System_Monitors USE Graph_Algorithms USE Vector_Operations USE Network_Matrix_Operations USE Network_Data_Types USE Line_Root_Finders IMPLICIT NONE PUBLIC::InitializeDualLineSearch,DestroyDualLineSearch,PerformDualSearch_Root,& DualLineDerivative,DualCostDerivatives,ExcessPower,FindRoot PUBLIC::Dual_Search_Timers,Dual_Line_Search,Dual_Search_Handle PUBLIC::root_finding_error,root_finding_success,root_finding_failure PUBLIC::free_root_finder,inactive_root_finder,active_root_finder LOGICAL,PUBLIC::debug_dual_search=.FALSE. TYPE Dual_Search_Timers INTEGER::cost_evaluation_timer=-1,bracketing_timer=-1,root_finding_timer=-1 ENDTYPE INTEGER,PARAMETER,PUBLIC::dual_lagrangian=1,primal_lagrangian=2 TYPE Dual_Line_Search INTEGER::PIN=1 INTEGER::lagrangian=dual_lagrangian TYPE(Line_Root_Finder)::root_finder TYPE(Directed_Graph),POINTER::graph=>NULL() TYPE(Network_SC_Cost),POINTER::cost_function=>NULL() TYPE(Network_Arrays),POINTER::arrays=>NULL() TYPE(Dual_Search_Timers)::timers REAL(KIND=r_wp),DIMENSION(:),POINTER::test_steps REAL(KIND=r_wp)::excess_power=0.0_r_wp REAL(KIND=r_wp)::step_size=1.0_r_wp,tolerance=0.1_r_wp LOGICAL::warm_starts=.FALSE. ENDTYPE TYPE Dual_Search_Handle TYPE(Dual_Line_Search),POINTER::search=>NULL() ENDTYPE TYPE(Dual_Search_Handle),DIMENSION(:),ALLOCATABLE,PUBLIC::dual_searches CONTAINS SUBROUTINE InitializeDualLineSearch(line_search) IMPLICIT NONE TYPE(Dual_Line_Search),INTENT(INOUT),TARGET::line_search TYPE(Directed_Graph),POINTER::graph TYPE(Network_Arrays),POINTER::arrays INTEGER::PIN,alloc_status LOGICAL::valid_problem graph=>line_search%graph arrays=>line_search%arrays valid_problem=.TRUE. IF(.NOT.ASSOCIATED(graph))THEN valid_problem=.FALSE. ELSE IF(.NOT.(ASSOCIATED(graph%heads_tails)))valid_problem=.FALSE. END IF IF(.NOT.ASSOCIATED(arrays))THEN valid_problem=.FALSE. ELSE IF(.NOT.(ASSOCIATED(arrays%arcs_conductances).AND.ASSOCIATED(arrays%nodes_potentials).AND.& ASSOCIATED(arrays%arcs_flows).AND.ASSOCIATED(arrays%arcs_voltages)))& valid_problem=.FALSE. IF(line_search%lagrangian==dual_lagrangian)THEN IF(.NOT.(ASSOCIATED(arrays%nodes_excess_potentials).AND.ASSOCIATED(arrays%supplies_demands).AND.& ASSOCIATED(arrays%arcs_excess_voltages)))& valid_problem=.FALSE. ELSE IF(line_search%lagrangian==primal_lagrangian)THEN IF(.NOT.(ASSOCIATED(arrays%arcs_excess_flows)))valid_problem=.FALSE. END IF END IF IF(.NOT.ASSOCIATED(line_search%cost_function))THEN valid_problem=.FALSE. ELSE END IF IF(.NOT.valid_problem)THEN CALL CriticalError(message="Network system not fully specified",& caller="InitializeDualLineSearch") RETURN END IF IF(.NOT.ASSOCIATED(line_search%test_steps))THEN ALLOCATE(line_search%test_steps(0:10)) line_search%test_steps=REAL((/0.0,1.0,2.0,3.5,5.0,7.50,10.0,15.0,20.0,50.0/),r_wp) END IF IF(.NOT.ALLOCATED(dual_searches))THEN ALLOCATE(dual_searches(1:MAX(10,line_search%PIN))) END IF END SUBROUTINE InitializeDualLineSearch SUBROUTINE DestroyDualLineSearch(line_search) IMPLICIT NONE TYPE(Dual_Line_Search),INTENT(INOUT),TARGET::line_search END SUBROUTINE DestroyDualLineSearch SUBROUTINE PerformDualSearch_Root(line_search,ElementalCosts,initialize) IMPLICIT NONE TYPE(Dual_Line_Search),INTENT(INOUT),TARGET::line_search INTERFACE SUBROUTINE ElementalCosts(cost_function,arguments_status,tolerance,& arcs_indices,arcs_flows,arcs_voltages,arcs_resistances,arcs_costs) USE Precision USE Network_Data_Types TYPE(Network_SC_Cost),INTENT(INOUT),OPTIONAL::cost_function CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::arcs_indices REAL(KIND=r_wp),DIMENSION(arcs_indices(1):),INTENT(INOUT),OPTIONAL::& arcs_flows,arcs_voltages,arcs_resistances,arcs_costs END SUBROUTINE ElementalCosts END INTERFACE LOGICAL,INTENT(IN),OPTIONAL::initialize TYPE(Directed_Graph),POINTER::graph TYPE(Network_Arrays),POINTER::arrays TYPE(Line_Root_Finder),POINTER::root_finder REAL(KIND=r_wp)::alpha,line_derivative INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::arc,node INTEGER::test_try LOGICAL::pre_calculate,bracketed IF((line_search%PININT(UBOUND(dual_searches,DIM=1),KIND=i_wp)))THEN CALL CriticalError(message="Not enough dual systems available or PIN out of range",& caller="InitializeDualNewton") RETURN ELSE IF(ASSOCIATED(dual_searches(line_search%PIN)%search))THEN CALL CriticalError(message="Dual system already taken",caller="InitializeDualNewton") RETURN ELSE dual_searches(line_search%PIN)%search=>line_search END IF END IF line_search%root_finder%PIN=line_search%PIN graph=>line_search%graph arrays=>line_search%arrays root_finder=>line_search%root_finder n_special_nodes=graph%n_special_nodes n_nodes=graph%n_nodes n_special_arcs=graph%n_special_arcs n_arcs=graph%n_arcs pre_calculate=.FALSE. IF(PRESENT(initialize))pre_calculate=initialize IF(pre_calculate)THEN IF(line_search%lagrangian==dual_lagrangian)THEN CALL ArcsVoltages(heads_tails=graph%heads_tails,node_offset=n_special_nodes,& nodes_potentials=arrays%nodes_potentials,arcs_voltages=arrays%arcs_voltages) CALL ArcsVoltages(heads_tails=graph%heads_tails,node_offset=n_special_nodes,& nodes_potentials=arrays%nodes_excess_potentials,& arcs_voltages=arrays%arcs_excess_voltages) line_search%excess_power=DOT_PRODUCT(arrays%supplies_demands,arrays%nodes_excess_potentials) ELSE CALL ArcsVoltages(heads_tails=graph%heads_tails,node_offset=n_special_nodes,& nodes_potentials=arrays%nodes_potentials,arcs_voltages=arrays%arcs_excess_voltages) line_search%excess_power=DOT_PRODUCT(arrays%arcs_excess_voltages,arrays%arcs_excess_flows) END IF END IF CALL StartTimer(line_search%timers%bracketing_timer) bracketed=.FALSE. alpha=2.0_r_wp Bracket:DO test_try=LBOUND(line_search%test_steps,DIM=1),UBOUND(line_search%test_steps,DIM=1) root_finder%root=line_search%test_steps(test_try) CALL DualLineDerivative(root_finder=root_finder,ElementalCosts=ElementalCosts,& first_derivative=line_derivative) IF(test_try==0)THEN IF(line_derivative>=0.0_r_wp)THEN bracketed=.FALSE. CALL Warning("Non descent search direction in dual line search",& caller="PerformDualSearch_Root") EXIT Bracket END IF ELSE IF(line_derivative>=0.0_r_wp)THEN bracketed=.TRUE. alpha=root_finder%root EXIT Bracket END IF END IF END DO Bracket IF(bracketed)THEN root_finder%interval(1)=line_search%test_steps(test_try-1) root_finder%interval(2)=line_search%test_steps(test_try) END IF CALL StopTimer(line_search%timers%bracketing_timer) IF(bracketed)THEN root_finder%tolerance=line_search%tolerance*(root_finder%interval(2)-root_finder%interval(1)) WRITE(message_log_unit,"(A, 2F5.2, A, F6.4)")& "-> Starting dual line search in [",root_finder%interval,& "] with tolerance ",root_finder%tolerance CALL StartTimer(line_search%timers%root_finding_timer) CALL FindRoot(ZeroOf=DualLineDerivative,root_finder=root_finder,& AuxilliaryRoutine=ElementalCosts) CALL StopTimer(line_search%timers%root_finding_timer) line_search%step_size=MAX(root_finder%tolerance,root_finder%root) IF(root_finder%convergence==root_finding_success)THEN WRITE(message_log_unit,*)"--> Minimum of dual Lagrangian found for step ",line_search%step_size ELSE IF(root_finder%convergence==root_finding_failure)THEN WRITE(message_log_unit,*)"--> Dual line search did not converge.",& " Last step size is ",line_search%step_size ELSE CALL Warning("Line-search failed to bracket root",caller="PerformDualSearch_Root") line_search%step_size=1.0_r_wp END IF ELSE line_search%step_size=1.0_r_wp CALL Warning("Dual Lagrangian minimum could not be bracketed",& caller="PerformDualSearch_Root") END IF NULLIFY(dual_searches(line_search%PIN)%search) END SUBROUTINE PerformDualSearch_Root SUBROUTINE DualLineDerivative(root_finder,first_derivative,second_derivative,ElementalCosts) IMPLICIT NONE TYPE(Line_Root_Finder),INTENT(INOUT)::root_finder REAL(KIND=r_wp),INTENT(OUT)::first_derivative REAL(KIND=r_wp),INTENT(OUT),OPTIONAL::second_derivative OPTIONAL::ElementalCosts INTERFACE SUBROUTINE ElementalCosts(cost_function,arguments_status,tolerance,& arcs_indices,arcs_flows,arcs_voltages,arcs_resistances,arcs_costs) USE Precision USE Network_Data_Types TYPE(Network_SC_Cost),INTENT(INOUT),OPTIONAL::cost_function CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::arcs_indices REAL(KIND=r_wp),DIMENSION(arcs_indices(1):),INTENT(INOUT),OPTIONAL::& arcs_flows,arcs_voltages,arcs_resistances,arcs_costs END SUBROUTINE ElementalCosts END INTERFACE TYPE(Dual_Line_Search),POINTER::line_search TYPE(Directed_Graph),POINTER::graph TYPE(Network_Arrays),POINTER::arrays INTEGER(KIND=i_wp)::arc,node REAL(KIND=r_wp),PARAMETER::eps=EPSILON(1.0_r_wp) REAL(KIND=r_wp)::dummy,arc_resistance,arc_voltage,arc_flow IF(.NOT.PRESENT(ElementalCosts))RETURN line_search=>dual_searches(root_finder%PIN)%search graph=>line_search%graph arrays=>line_search%arrays PrimalOrDual:IF(line_search%lagrangian==dual_lagrangian)THEN CALL StartTimer(line_search%timers%cost_evaluation_timer) IF(PRESENT(second_derivative))THEN CALL DualCostDerivatives(ElementalCosts=ElementalCosts,& cost_function=line_search%cost_function,lagrangian_function=line_search%lagrangian,& arc_offset=graph%n_special_arcs,& arcs_voltages=arrays%arcs_voltages,arcs_excess_voltages=arrays%arcs_excess_voltages,& step_size=root_finder%root,arcs_flows=arrays%arcs_excess_flows,& arcs_conductances=arrays%arcs_conductances,warm_starts=line_search%warm_starts,& conductances_interval=(/line_search%cost_function%low_conductance,line_search%cost_function%high_conductance/)) ELSE CALL DualCostDerivatives(ElementalCosts=ElementalCosts,& cost_function=line_search%cost_function,lagrangian_function=line_search%lagrangian,& arc_offset=graph%n_special_arcs,& arcs_voltages=arrays%arcs_voltages,arcs_excess_voltages=arrays%arcs_excess_voltages,& step_size=root_finder%root,arcs_flows=arrays%arcs_excess_flows,warm_starts=line_search%warm_starts) END IF CALL StopTimer(line_search%timers%cost_evaluation_timer) first_derivative=-line_search%excess_power+& ExcessPower(arcs_voltages=arrays%arcs_excess_voltages,arcs_flows=arrays%arcs_excess_flows) IF(PRESENT(second_derivative))THEN second_derivative=ExcessPower(arcs_voltages=arrays%arcs_excess_voltages,& arcs_conductances=arrays%arcs_conductances) END IF ELSE IF(line_search%lagrangian==primal_lagrangian)THEN PrimalOrDual CALL StartTimer(line_search%timers%cost_evaluation_timer) IF(PRESENT(second_derivative))THEN CALL DualCostDerivatives(ElementalCosts=ElementalCosts,& cost_function=line_search%cost_function,lagrangian_function=line_search%lagrangian,& arc_offset=graph%n_special_arcs,& arcs_flows=arrays%arcs_flows,arcs_excess_flows=arrays%arcs_excess_flows,& step_size=root_finder%root,arcs_voltages=arrays%arcs_excess_voltages,& arcs_resistances=arrays%arcs_conductances,warm_starts=line_search%warm_starts,& resistances_interval=& (/1.0_r_wp/line_search%cost_function%high_conductance,1.0_r_wp/line_search%cost_function%low_conductance/)) ELSE CALL DualCostDerivatives(ElementalCosts=ElementalCosts,& cost_function=line_search%cost_function,lagrangian_function=line_search%lagrangian,& arc_offset=graph%n_special_arcs,& arcs_flows=arrays%arcs_flows,arcs_excess_flows=arrays%arcs_excess_flows,& step_size=root_finder%root,arcs_voltages=arrays%arcs_excess_voltages,warm_starts=line_search%warm_starts) END IF CALL StopTimer(line_search%timers%cost_evaluation_timer) first_derivative=-line_search%excess_power+& ExcessPower(arcs_voltages=arrays%arcs_excess_voltages,arcs_flows=arrays%arcs_excess_flows) IF(PRESENT(second_derivative))THEN second_derivative=ExcessPower(arcs_flows=arrays%arcs_excess_flows,& arcs_resistances=arrays%arcs_conductances) END IF CALL VectorExponentiate(vector=arrays%arcs_conductances) END IF PrimalOrDual END SUBROUTINE DualLineDerivative FUNCTION ExcessPower(arcs_flows,arcs_voltages,& arcs_conductances,arcs_resistances)RESULT(power) REAL(KIND=r_wp),DIMENSION(:),INTENT(IN),OPTIONAL::& arcs_flows,arcs_voltages,& arcs_conductances,arcs_resistances REAL(KIND=r_wp)::power INTEGER::arc power=0.0_r_wp InputChoices:IF(PRESENT(arcs_flows).AND.PRESENT(arcs_voltages))THEN power=DOT_PRODUCT(arcs_flows,arcs_voltages) ELSE IF(PRESENT(arcs_voltages))THEN InputChoices IF(PRESENT(arcs_conductances))THEN DO arc=1,INT(SIZE(arcs_voltages),KIND=i_wp) power=power+arcs_voltages(arc)**2*arcs_conductances(arc) END DO ELSE IF(PRESENT(arcs_resistances))THEN DO arc=1,INT(SIZE(arcs_voltages),KIND=i_wp) power=power+arcs_voltages(arc)**2/arcs_resistances(arc) END DO END IF ELSE IF(PRESENT(arcs_flows).AND.PRESENT(arcs_conductances))THEN InputChoices IF(PRESENT(arcs_conductances))THEN DO arc=1,INT(SIZE(arcs_voltages),KIND=i_wp) power=power+arcs_voltages(arc)**2*arcs_conductances(arc) END DO ELSE IF(PRESENT(arcs_resistances))THEN DO arc=1,INT(SIZE(arcs_flows),KIND=i_wp) power=power+arcs_flows(arc)**2*arcs_resistances(arc) END DO END IF END IF InputChoices END FUNCTION ExcessPower SUBROUTINE DualCostDerivatives(ElementalCosts,cost_function,lagrangian_function,& arc_offset,arcs_voltages,arcs_flows,& arcs_excess_voltages,arcs_excess_flows,step_size,& arcs_conductances,arcs_resistances,& conductances_interval,resistances_interval,& warm_starts) IMPLICIT NONE INTERFACE SUBROUTINE ElementalCosts(cost_function,arguments_status,tolerance,& arcs_indices,arcs_flows,arcs_voltages,arcs_resistances,arcs_costs) USE Precision USE Network_Data_Types TYPE(Network_SC_Cost),INTENT(INOUT),OPTIONAL::cost_function CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::arcs_indices REAL(KIND=r_wp),DIMENSION(arcs_indices(1):),INTENT(INOUT),OPTIONAL::& arcs_flows,arcs_voltages,arcs_resistances,arcs_costs END SUBROUTINE ElementalCosts END INTERFACE TYPE(Network_SC_Cost),INTENT(INOUT)::cost_function INTEGER,OPTIONAL::lagrangian_function INTEGER(KIND=i_wp),INTENT(IN)::arc_offset REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(INOUT),OPTIONAL::arcs_voltages,arcs_flows REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::& arcs_excess_voltages,arcs_excess_flows REAL(KIND=r_wp),INTENT(IN),OPTIONAL::step_size REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(INOUT),OPTIONAL::& arcs_conductances,arcs_resistances REAL(KIND=r_wp),DIMENSION(2),INTENT(IN),OPTIONAL::& conductances_interval,resistances_interval LOGICAL,INTENT(IN),OPTIONAL::warm_starts REAL(KIND=r_wp)::alpha,arc_voltage,arc_flow,arc_resistance,arc_conductance,dummy(1) REAL(KIND=r_wp),PARAMETER::eps=EPSILON(1.0_r_wp) INTEGER(KIND=i_wp)::arc,n_special_arcs,n_arcs INTEGER::lagrangian CHARACTER(LEN=4)::arguments_status CHARACTER::warm_or_cold LOGICAL::calculate_flows,calculate_voltages,& calculate_resistances,calculate_conductances,& bound_conductances,bound_resistances n_special_arcs=-INT(LBOUND(arcs_voltages,DIM=1),KIND=i_wp) n_arcs=INT(UBOUND(arcs_voltages,DIM=1),KIND=i_wp) lagrangian=dual_lagrangian IF(PRESENT(lagrangian_function))lagrangian=lagrangian_function IF(PRESENT(step_size))THEN alpha=step_size ELSE alpha=1.0_r_wp END IF IF(PRESENT(warm_starts))THEN IF(warm_starts)THEN warm_or_cold='W' ELSE warm_or_cold='C' END IF ELSE warm_or_cold='C' END IF IF(lagrangian==dual_lagrangian)THEN IF(.NOT.PRESENT(arcs_voltages))& CALL CriticalError(message="Arcs voltages needed but not supplied",& caller="DualCostDerivatives") arguments_status="DFDD" calculate_voltages=PRESENT(arcs_excess_voltages) calculate_flows=PRESENT(arcs_flows) IF(calculate_flows)arguments_status(1:1)=warm_or_cold ELSE IF(lagrangian==primal_lagrangian)THEN IF(.NOT.PRESENT(arcs_flows))& CALL CriticalError(message="Arcs flows needed but not supplied",& caller="DualCostDerivatives") arguments_status="FDDD" calculate_flows=PRESENT(arcs_excess_flows) calculate_voltages=PRESENT(arcs_voltages) IF(calculate_voltages)arguments_status(2:2)=warm_or_cold END IF calculate_conductances=PRESENT(arcs_conductances) calculate_resistances=PRESENT(arcs_resistances) IF(calculate_conductances.OR.calculate_resistances)arguments_status(3:3)=warm_or_cold bound_conductances=PRESENT(conductances_interval) bound_resistances=PRESENT(resistances_interval) IF(lagrangian==dual_lagrangian)THEN IF(calculate_voltages)THEN arcs_voltages=arcs_voltages+alpha*arcs_excess_voltages END IF ELSE IF(lagrangian==primal_lagrangian)THEN IF(calculate_flows)THEN arcs_flows=arcs_flows+alpha*arcs_excess_flows END IF END IF CallElementalCosts:IF(calculate_resistances)THEN CALL ElementalCosts(cost_function=cost_function,arcs_flows=arcs_flows,& arcs_voltages=arcs_voltages,arcs_resistances=arcs_resistances,& arcs_costs=dummy,arcs_indices=(/-n_special_arcs,n_arcs/),& arguments_status=arguments_status,tolerance=cost_function%eps) IF(bound_resistances)THEN DO arc=-n_special_arcs,n_arcs arcs_resistances(arc)=& MIN(resistances_interval(2),MAX(arcs_resistances(arc),resistances_interval(1))) END DO END IF IF(calculate_conductances)arcs_conductances=1.0_r_wp/(arcs_resistances+eps) ELSE IF(calculate_conductances)THEN CallElementalCosts CALL ElementalCosts(cost_function=cost_function,arcs_flows=arcs_flows,& arcs_voltages=arcs_voltages,arcs_resistances=arcs_conductances,& arcs_costs=dummy,arcs_indices=(/-n_special_arcs,n_arcs/),& arguments_status=arguments_status,tolerance=cost_function%eps) arcs_conductances=1.0_r_wp/(arcs_conductances+eps) ELSE CallElementalCosts CALL ElementalCosts(cost_function=cost_function,arcs_flows=arcs_flows,& arcs_voltages=arcs_voltages,arcs_resistances=dummy,& arcs_costs=dummy,arcs_indices=(/-n_special_arcs,n_arcs/),& arguments_status=arguments_status,tolerance=cost_function%eps) END IF CallElementalCosts IF(calculate_conductances.AND.bound_conductances)THEN DO arc=-n_special_arcs,n_arcs arcs_conductances(arc)=& MIN(conductances_interval(2),MAX(arcs_conductances(arc),conductances_interval(1))) END DO END IF IF(lagrangian==dual_lagrangian)THEN IF(calculate_voltages)THEN arcs_voltages=arcs_voltages-alpha*arcs_excess_voltages END IF ELSE IF(lagrangian==primal_lagrangian)THEN IF(calculate_flows)THEN arcs_flows=arcs_flows-alpha*arcs_excess_flows END IF END IF END SUBROUTINE DualCostDerivatives SUBROUTINE FindRoot(ZeroOf,root_finder,AuxilliaryRoutine) IMPLICIT NONE INTERFACE SUBROUTINE ZeroOf(root_finder,value,derivative,Auxilliary) USE Precision USE Line_Root_Finders TYPE(Line_Root_Finder),INTENT(INOUT)::root_finder REAL(KIND=r_wp),INTENT(OUT)::value REAL(KIND=r_wp),INTENT(OUT),OPTIONAL::derivative EXTERNAL::Auxilliary OPTIONAL::Auxilliary END SUBROUTINE ZeroOf END INTERFACE TYPE(Line_Root_Finder),INTENT(INOUT)::root_finder EXTERNAL::AuxilliaryRoutine OPTIONAL::AuxilliaryRoutine INTEGER::iteration REAL(KIND=r_wp)::x,x_new,x_old,dx,dx_old,& a,b,midpoint,f,fa,fb,df,f_left,f_right LOGICAL::use_derivative root_finder%status=active_root_finder root_finder%convergence=root_finding_failure root_finder%n_iterations=0 use_derivative=((root_finder%method(1:1)=='N').OR.(root_finder%method(1:1)=='n')) a=root_finder%interval(1) b=root_finder%interval(2) IF(PRESENT(AuxilliaryRoutine))THEN CALL ZeroOf(root_finder,value=fb,Auxilliary=AuxilliaryRoutine) ELSE CALL ZeroOf(root_finder,value=fb) END IF root_finder%root=a IF(use_derivative)THEN IF(PRESENT(AuxilliaryRoutine))THEN CALL ZeroOf(root_finder=root_finder,value=fa,derivative=df,& Auxilliary=AuxilliaryRoutine) ELSE CALL ZeroOf(root_finder=root_finder,value=fa,derivative=df) END IF ELSE IF(PRESENT(AuxilliaryRoutine))THEN CALL ZeroOf(root_finder=root_finder,value=fa,Auxilliary=AuxilliaryRoutine) ELSE CALL ZeroOf(root_finder=root_finder,value=fa) END IF END IF IF(.NOT.use_derivative)df=(fa-fb)/(a-b) x=a f=fa IF(fa*fb>0.0_r_wp)THEN CALL Warning(message="Root bracketing interval incorrect input",caller="FindRoot") root_finder%convergence=root_finding_error root_finder%status=free_root_finder RETURN END IF dx_old=0.5_r_wp*ABS(a-b) iteration=0 CutAndBisect:DO IF(iteration>=root_finder%max_iterations)THEN root_finder%convergence=root_finding_failure EXIT CutAndBisect END IF iteration=iteration+1 IF(ABS(a-b)=MAX(a,b)).OR.& (ABS(2*dx)>ABS(dx_old)))THEN midpoint=0.5_r_wp*(a+b) dx=midpoint-x x=midpoint ELSE x=x_new END IF dx_old=dx IF(ABS(dx)