MODULE Dual_Newton_SSCNO
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Graph_Algorithms
      USE Vector_Operations
      USE Network_Matrix_Operations
      USE Network_Data_Types
      USE Conjugate_Gradient
      USE Dual_Network_Solvers
      USE Dual_Line_Minimizers
      IMPLICIT NONE
      PUBLIC::InitializeDualSSCNO,DestroyDualSSCNO,SolveDualSSCNO_TCGN
      
      PUBLIC::SSCNO_Timers,SSCNO_Error_Tolerances,&
      Dual_Network_Problem,Dual_Problem_Handle
      
      INTEGER,PARAMETER,PUBLIC::&
      SSCNO_unconverged=-1,SSCNO_converged=0,SSCNO_failed=1
      
      TYPE SSCNO_Timers
      INTEGER::initialization_timer=-1,array_update_timer=-1,cost_evaluation_timer=-1,&
      solver_timer=-1,line_search_timer=-1
      END TYPE
      
      TYPE SSCNO_Error_Tolerances
      REAL(KIND=r_wp)::relative_error=1.0_r_wp
      REAL(KIND=r_wp)::voltages_infeasibility_norm=-1.0_r_wp,flows_infeasibility_norm=-1.0_r_wp
      
      REAL(KIND=r_wp)::excess_potentials_norm=-1.0_r_wp,excess_flows_norm=-1.0_r_wp
      
      ENDTYPE
      
      
      INTEGER,PARAMETER,PUBLIC::SSCNO_TDN=1,SSCNO_TSQP=2
      TYPE Dual_Network_Problem
      INTEGER::PIN=1
      INTEGER::method=SSCNO_TDN
      INTEGER(KIND=i_wp)::n_fixed_potentials=0,n_free_special_nodes=0
      TYPE(Network_Problem)::problem
      TYPE(Dual_Network_System)::system
      TYPE(Dual_Line_Search)::line_search
      TYPE(SSCNO_Timers)::timers
      TYPE(SSCNO_Error_Tolerances)::error_tolerances
      REAL(KIND=r_wp)::forcing_prefactor=0.1_r_wp,forcing_exponent=1.0_r_wp
      REAL(KIND=r_wp)::ls_tolerance_prefactor=0.1_r_wp,ls_tolerance_exponent=1.0_r_wp
      INTEGER::convergence=SSCNO_unconverged
      LOGICAL,DIMENSION(4)::monitor_indicators=.TRUE.
      
      INTEGER::n_iterations=0,max_iterations=100
      LOGICAL::log_SSCNO=.TRUE.,warm_starts=.FALSE.,&
      perform_line_search=.TRUE.,force_feasible=.FALSE.
      LOGICAL,DIMENSION(6)::allocated_arrays=.FALSE.
      ENDTYPE
      
      TYPE Dual_Problem_Handle
      TYPE(Dual_Network_Problem),POINTER::dual_problem
      ENDTYPE
      TYPE(Dual_Problem_Handle),DIMENSION(:),ALLOCATABLE,PUBLIC::&
      dual_problems
      
      CONTAINS
      
      
      SUBROUTINE InitializeDualSSCNO(dual_problem)
      IMPLICIT NONE
      TYPE(Dual_Network_Problem),INTENT(INOUT),TARGET::dual_problem
      
      TYPE(Directed_Graph),POINTER::graph,subgraph,special_subgraph
      TYPE(Network_SC_Cost),POINTER::cost_function
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(Dual_Line_Search),POINTER::line_search
      
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs
      INTEGER(KIND=i_wp)::arc,node,special_arc,special_node
      INTEGER::alloc_status
      LOGICAL::valid_problem,fixed_head,fixed_tail
      
      dual_system=>dual_problem%system
      line_search=>dual_problem%line_search
      graph=>dual_problem%problem%graph
      arrays=>dual_problem%problem%arrays
      cost_function=>dual_problem%problem%cost_function
      
      
      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%nodes_potentials).AND.ASSOCIATED(arrays%arcs_flows).AND.&
      ASSOCIATED(arrays%supplies_demands)))THEN
      valid_problem=.FALSE.
      END IF
      END IF
      IF(.NOT.ASSOCIATED(cost_function))THEN
      valid_problem=.FALSE.
      ELSE
      
      END IF
      IF(.NOT.valid_problem)THEN
      CALL CriticalError(message="SSCNO dual problem not fully specified",&
      caller="InitializeDualSSCNO")
      RETURN
      END IF
      
      CALL StartTimer(dual_problem%timers%initialization_timer)
      
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      
      IF(.NOT.ASSOCIATED(arrays%arcs_voltages))THEN
      dual_problem%allocated_arrays(1)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%arcs_voltages))THEN
      ALLOCATE(arrays%arcs_voltages(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      IF(.NOT.ASSOCIATED(arrays%arcs_conductances))THEN
      dual_problem%allocated_arrays(2)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%arcs_conductances))THEN
      ALLOCATE(arrays%arcs_conductances(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      IF(.NOT.ASSOCIATED(arrays%arcs_excess_voltages))THEN
      dual_problem%allocated_arrays(3)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%arcs_excess_voltages))THEN
      ALLOCATE(arrays%arcs_excess_voltages(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      IF(.NOT.ASSOCIATED(arrays%arcs_excess_flows))THEN
      dual_problem%allocated_arrays(4)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%arcs_excess_flows))THEN
      ALLOCATE(arrays%arcs_excess_flows(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      IF(.NOT.ASSOCIATED(arrays%nodes_excess_potentials))THEN
      dual_problem%allocated_arrays(5)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%nodes_excess_potentials))THEN
      ALLOCATE(arrays%nodes_excess_potentials(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      IF(.NOT.ASSOCIATED(arrays%nodes_excess_flows))THEN
      dual_problem%allocated_arrays(6)=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(arrays%nodes_excess_flows))THEN
      ALLOCATE(arrays%nodes_excess_flows(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      END IF
      END IF
      
      
      subgraph=>dual_problem%problem%subgraph
      subgraph%n_nodes=n_nodes
      subgraph%n_arcs=n_arcs
      subgraph%n_special_nodes=0
      subgraph%n_special_arcs=-1
      subgraph%heads_tails=>graph%heads_tails(:,1:)
      
      special_subgraph=>dual_problem%problem%special_subgraph
      special_subgraph%n_nodes=n_nodes
      special_subgraph%n_arcs=n_special_arcs
      special_subgraph%n_special_nodes=n_special_nodes
      special_subgraph%n_special_arcs=-1
      
      
      FixedPotentials:IF(dual_problem%n_fixed_potentials>0)THEN
      
      
      ALLOCATE(special_subgraph%heads_tails(1:2,-special_subgraph%n_special_arcs:special_subgraph%n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=INT(SIZE(special_subgraph%heads_tails),KIND=i_wp),mold=1_i_wp,&
      caller="InitializeDualSSCNO",alloc_status=alloc_status)
      special_subgraph%heads_tails(1:2,1:graph%n_special_arcs)=graph%heads_tails(1:2,(-1):(-n_special_arcs):(-1))
      
      IF(dual_system%method==dual_network_pcg)THEN
      dual_system%ground_dummy=.FALSE.
      ELSE
      dual_system%ground_dummy=.TRUE.
      END IF
      
      RerouteArcs:DO special_arc=1,n_special_arcs
      IF(graph%heads_tails(1,-special_arc)<-dual_problem%n_free_special_nodes)&
      graph%heads_tails(1,-special_arc)=0
      IF(graph%heads_tails(2,-special_arc)<-dual_problem%n_free_special_nodes)&
      graph%heads_tails(2,-special_arc)=0
      END DO RerouteArcs
      
      dual_problem%n_free_special_nodes=graph%n_special_nodes-dual_problem%n_fixed_potentials
      ALLOCATE(dual_system%graph)
      dual_system%graph%heads_tails=>graph%heads_tails
      dual_system%graph%n_special_nodes=dual_problem%n_free_special_nodes
      dual_system%graph%n_special_arcs=n_special_arcs
      dual_system%graph%n_nodes=n_nodes
      dual_system%graph%n_arcs=n_arcs
      
      ALLOCATE(dual_system%arrays)
      dual_system%arrays%arcs_conductances=>arrays%arcs_conductances
      dual_system%arrays%supplies_demands=>SubarrayPointer(offset=-dual_problem%n_free_special_nodes,&
      array=arrays%supplies_demands(-dual_problem%n_free_special_nodes:))
      dual_system%arrays%nodes_excess_potentials=>SubarrayPointer(offset=-dual_problem%n_free_special_nodes,&
      array=arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:))
      dual_system%arrays%nodes_excess_flows=>SubarrayPointer(offset=-dual_problem%n_free_special_nodes,&
      array=arrays%nodes_excess_flows(-dual_problem%n_free_special_nodes:))
      IF(ASSOCIATED(arrays%nodes_coordinates))dual_system%arrays%nodes_coordinates=>&
      SubarrayPointer(row_offset=1,column_offset=-dual_problem%n_free_special_nodes,&
      array=arrays%nodes_coordinates(:,-dual_problem%n_free_special_nodes:))
      
      ELSE FixedPotentials
      
      special_subgraph%heads_tails=>graph%heads_tails(1:2,(-1):(-n_special_arcs):(-1))
      dual_system%ground_dummy=.FALSE.
      
      dual_problem%n_free_special_nodes=n_special_nodes
      graph%heads_tails(0,1)=n_nodes
      dual_system%graph=>graph
      dual_system%arrays=>arrays
      
      END IF FixedPotentials
      
      CALL InitializeDualSolver(dual_system=dual_system)
      
      IF(dual_system%method==dual_network_pcg)THEN
      CALL InitializeDualPreconditioner(dual_system=dual_system)
      END IF
      
      line_search%graph=>graph
      line_search%arrays=>arrays
      line_search%cost_function=>cost_function
      IF(dual_problem%method==SSCNO_TDN)THEN
      line_search%lagrangian=dual_lagrangian
      ELSE
      line_search%lagrangian=primal_lagrangian
      END IF
      CALL InitializeDualLineSearch(line_search=line_search)
      
      CALL StopTimer(dual_problem%timers%initialization_timer)
      
      END SUBROUTINE InitializeDualSSCNO
      
      
      
      SUBROUTINE DestroyDualSSCNO(dual_problem)
      IMPLICIT NONE
      TYPE(Dual_Network_Problem),INTENT(INOUT),TARGET::dual_problem
      
      TYPE(Directed_Graph),POINTER::graph,special_subgraph
      TYPE(Network_Arrays),POINTER::arrays
      
      INTEGER(KIND=i_wp)::arc,node,special_arc,special_node
      INTEGER::alloc_status
      LOGICAL::valid_problem,fixed_head,fixed_tail
      
      graph=>dual_problem%problem%graph
      special_subgraph=>dual_problem%problem%special_subgraph
      arrays=>dual_problem%problem%arrays
      
      IF(dual_problem%n_fixed_potentials>0)THEN
      
      
      graph%heads_tails(1:2,(-1):(-graph%n_special_arcs):(-1))=special_subgraph%heads_tails(1:2,1:graph%n_special_arcs)
      
      CALL RecordAllocation(n_elements=-INT(SIZE(special_subgraph%heads_tails),KIND=i_wp),&
      mold=1_i_wp,caller="DestroyDualSSCNO")
      DEALLOCATE(special_subgraph%heads_tails,STAT=alloc_status)
      DEALLOCATE(dual_problem%system%graph)
      DEALLOCATE(dual_problem%system%arrays)
      
      ELSE
      
      graph%heads_tails(0,1)=0
      
      END IF
      
      IF(dual_problem%allocated_arrays(1))THEN
      
      IF(ASSOCIATED(arrays%arcs_voltages))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%arcs_voltages),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%arcs_voltages,STAT=alloc_status)
      END IF
      END IF
      IF(dual_problem%allocated_arrays(2))THEN
      
      IF(ASSOCIATED(arrays%arcs_conductances))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%arcs_conductances),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%arcs_conductances,STAT=alloc_status)
      END IF
      END IF
      IF(dual_problem%allocated_arrays(3))THEN
      
      IF(ASSOCIATED(arrays%arcs_excess_voltages))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%arcs_excess_voltages),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%arcs_excess_voltages,STAT=alloc_status)
      END IF
      END IF
      IF(dual_problem%allocated_arrays(4))THEN
      
      IF(ASSOCIATED(arrays%arcs_excess_flows))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%arcs_excess_flows),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%arcs_excess_flows,STAT=alloc_status)
      END IF
      END IF
      IF(dual_problem%allocated_arrays(5))THEN
      
      IF(ASSOCIATED(arrays%nodes_excess_potentials))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%nodes_excess_potentials),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%nodes_excess_potentials,STAT=alloc_status)
      END IF
      END IF
      IF(dual_problem%allocated_arrays(6))THEN
      
      IF(ASSOCIATED(arrays%nodes_excess_flows))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(arrays%nodes_excess_flows),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(arrays%nodes_excess_flows,STAT=alloc_status)
      END IF
      END IF
      
      IF(dual_problem%system%method==dual_network_pcg)THEN
      CALL DestroyDualPreconditioner(dual_system=dual_problem%system)
      END IF
      
      CALL DestroyDualSolver(dual_system=dual_problem%system)
      CALL DestroyDualLineSearch(line_search=dual_problem%line_search)
      
      END SUBROUTINE DestroyDualSSCNO
      
      
      
      SUBROUTINE SolveDualSSCNO_TCGN(dual_problem,ElementalCosts)
      IMPLICIT NONE
      TYPE(Dual_Network_Problem),INTENT(INOUT),TARGET::dual_problem
      
      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(Directed_Graph),POINTER::graph,subgraph,special_subgraph
      TYPE(Network_SC_Cost),POINTER::cost_function
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(CG_Solver),POINTER::solver
      TYPE(Dual_Line_Search),POINTER::line_search
      
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs
      INTEGER(KIND=i_wp)::arc,node,special_arc,special_node
      INTEGER::alloc_status,iteration
      
      REAL(KIND=r_wp),DIMENSION(4)::error_tolerances
      REAL(KIND=r_wp)::relative_tolerance,initial_flows_inf_norm,initial_voltages_inf_norm,eta
      LOGICAL,DIMENSION(4)::convergence_indicators
      LOGICAL::relative_convergence
      
      dual_system=>dual_problem%system
      line_search=>dual_problem%line_search
      graph=>dual_problem%problem%graph
      subgraph=>dual_problem%problem%subgraph
      special_subgraph=>dual_problem%problem%special_subgraph
      cost_function=>dual_problem%problem%cost_function
      arrays=>dual_problem%problem%arrays
      IF(dual_system%method==dual_network_pcg)solver=>dual_problem%system%system%solver
      
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      
      dual_system%PIN=dual_problem%PIN
      line_search%PIN=dual_problem%PIN
      line_search%warm_starts=dual_problem%warm_starts
      
      
      IF(cost_function%low_conductance<0.0_r_wp)&
      cost_function%low_conductance=MAX(1.0E-3_r_wp/cost_function%typical_resistance,&
10    .0_r_wp*EPSILON(1.0_r_wp))
      IF(cost_function%high_conductance<0.0_r_wp)&
      cost_function%high_conductance=MIN(1.0E+3_r_wp/cost_function%typical_resistance,&
0     .1_r_wp/EPSILON(1.0_r_wp))
      dual_system%grounding_conductance=1.0_r_wp/cost_function%typical_resistance
      
      
      IF(dual_problem%error_tolerances%relative_error>0.0_r_wp)THEN
      relative_tolerance=dual_problem%error_tolerances%relative_error
      ELSE
      relative_tolerance=1.0_r_wp
      END IF
      error_tolerances=-1.0_r_wp
      
      IF(dual_problem%error_tolerances%voltages_infeasibility_norm>0.0_r_wp)THEN
      dual_problem%monitor_indicators(1)=.TRUE.
      error_tolerances(1)=dual_problem%error_tolerances%voltages_infeasibility_norm
      END IF
      IF(dual_problem%error_tolerances%flows_infeasibility_norm>0.0_r_wp)THEN
      dual_problem%monitor_indicators(2)=.TRUE.
      error_tolerances(2)=dual_problem%error_tolerances%flows_infeasibility_norm
      END IF
      IF(dual_problem%error_tolerances%excess_potentials_norm>0.0_r_wp)THEN
      dual_problem%monitor_indicators(3)=.TRUE.
      error_tolerances(3)=dual_problem%error_tolerances%excess_potentials_norm
      END IF
      IF(dual_problem%error_tolerances%excess_flows_norm>0.0_r_wp)THEN
      dual_problem%monitor_indicators(4)=.TRUE.
      error_tolerances(4)=dual_problem%error_tolerances%excess_flows_norm
      END IF
      convergence_indicators=(error_tolerances<0.0_r_wp)
      relative_convergence=.FALSE.
      
      initial_flows_inf_norm=1.0_r_wp
      initial_voltages_inf_norm=1.0_r_wp
      
      PrintHeader:IF(dual_problem%log_SSCNO)THEN
      WRITE(message_print_unit,*)
      WRITE(message_print_unit,"(A)")&
      "_______________________________________Start of TSQP statistics"
      IF(dual_system%method==dual_network_pcg)THEN
      WRITE(message_print_unit,"(A10,A15,4A15,A10,A5,A10,A5,A10)")&
      "Iter.#","Rel.Err.","Inf.Volt.","Inf.Flow","Ex.Pot.","Ex.Flow","LS step","#LS","CG #iter","RC","Conv."
      ELSE
      WRITE(message_print_unit,"(A10,A15,4A15,A10,A5,A10,A5,A10)")&
      "Iter.#","Rel.Err.","Inf.Volt.","Inf.Flow","Ex.Pot.","Ex.Flow","LS step","#LS","Fill-in","RC","Conv."
      END IF
      WRITE(message_print_unit,"(A)")&
      "_______________________________________________________________"
      WRITE(message_print_unit,*)
      END IF PrintHeader
      
      iteration=0
      TCGN:DO
      
      IF(iteration>=dual_problem%max_iterations)THEN
      dual_problem%convergence=SSCNO_unconverged
      EXIT TCGN
      END IF
      iteration=iteration+1
      
      MonitorStatistics:IF(dual_problem%log_SSCNO)THEN
      PrintStatistics:IF(iteration>1)THEN
      
      IF(dual_system%method==dual_network_pcg)THEN
      WRITE(message_print_unit,"(I10,F15.10,4E15.3,F10.3,I5,I10,L5,4L2)")&
      iteration-1,dual_problem%error_tolerances%relative_error,&
      dual_problem%error_tolerances%voltages_infeasibility_norm,dual_problem%error_tolerances%flows_infeasibility_norm,&
      dual_problem%error_tolerances%excess_potentials_norm,dual_problem%error_tolerances%excess_flows_norm,&
      line_search%step_size,line_search%root_finder%n_iterations,&
      solver%n_iterations,relative_convergence,convergence_indicators
      ELSE
      WRITE(message_print_unit,"(I10,F15.10,4E15.3,F10.3,I5,A10,L5,4L2)")&
      iteration-1,dual_problem%error_tolerances%relative_error,&
      dual_problem%error_tolerances%voltages_infeasibility_norm,dual_problem%error_tolerances%flows_infeasibility_norm,&
      dual_problem%error_tolerances%excess_potentials_norm,dual_problem%error_tolerances%excess_flows_norm,&
      line_search%step_size,line_search%root_finder%n_iterations,&
      "",relative_convergence,convergence_indicators
      END IF
      END IF PrintStatistics
      
      line_search%step_size=-1.0_r_wp
      line_search%root_finder%n_iterations=-1
      IF(dual_system%method==dual_network_pcg)solver%n_iterations=-1
      END IF MonitorStatistics
      
      
      IF(ALL(convergence_indicators).AND.relative_convergence)THEN
      
      dual_problem%convergence=SSCNO_converged
      EXIT TCGN
      END IF
      
      
      DualGradient:IF(dual_problem%method==SSCNO_TDN)THEN
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL ArcsVoltages(heads_tails=subgraph%heads_tails,node_offset=subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-subgraph%n_special_nodes:subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_voltages(-subgraph%n_special_arcs:subgraph%n_arcs))
      CALL ArcsVoltages(heads_tails=special_subgraph%heads_tails,node_offset=special_subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-special_subgraph%n_special_nodes:special_subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_voltages(special_subgraph%n_special_arcs:(-special_subgraph%n_arcs):(-1)))
      arrays%arcs_voltages(0)=0.0_r_wp
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      IF(dual_problem%monitor_indicators(4))THEN
      
      CALL VectorCopy(source=arrays%arcs_flows,target=arrays%arcs_excess_flows)
      END IF
      
      CALL StartTimer(dual_problem%timers%cost_evaluation_timer)
      CALL DualCostDerivatives(ElementalCosts=ElementalCosts,&
      cost_function=cost_function,lagrangian_function=dual_lagrangian,&
      arc_offset=n_special_arcs,arcs_voltages=arrays%arcs_voltages,arcs_flows=arrays%arcs_flows,&
      arcs_conductances=arrays%arcs_conductances,warm_starts=dual_problem%warm_starts,&
      conductances_interval=(/cost_function%low_conductance,cost_function%high_conductance/))
      arrays%arcs_flows(0)=0.0_r_wp
      arrays%arcs_conductances(0)=1.0_r_wp/cost_function%typical_resistance
      CALL StopTimer(dual_problem%timers%cost_evaluation_timer)
      
      IF(dual_problem%monitor_indicators(4))THEN
      CALL VectorSubtraction(from=arrays%arcs_flows,what=arrays%arcs_excess_flows,&
      reverse=.TRUE.)
      
      dual_problem%error_tolerances%excess_flows_norm=&
      SQRT(DOT_PRODUCT(arrays%arcs_excess_flows,arrays%arcs_excess_flows))
      IF(error_tolerances(4)>0.0_r_wp)&
      convergence_indicators(4)=(dual_problem%error_tolerances%excess_flows_norm<error_tolerances(4))
      END IF
      
      IF(dual_problem%monitor_indicators(1))THEN
      dual_problem%error_tolerances%voltages_infeasibility_norm=0.0_r_wp
      IF(error_tolerances(1)>0.0_r_wp)&
      convergence_indicators(1)=(dual_problem%error_tolerances%flows_infeasibility_norm<error_tolerances(1))
      END IF
      
      
      END IF DualGradient
      
      IF((dual_problem%method==SSCNO_TDN).OR.dual_problem%monitor_indicators(2))THEN
      
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL NodesExcessFlows(heads_tails=graph%heads_tails,node_offset=n_special_nodes,&
      arcs_flows=arrays%arcs_flows,excess_flows=arrays%nodes_excess_flows)
      CALL VectorSubtraction(from=arrays%supplies_demands(-dual_problem%n_free_special_nodes:),&
      what=arrays%nodes_excess_flows(-dual_problem%n_free_special_nodes:),&
      reverse=.TRUE.)
      
      IF(dual_problem%method/=SSCNO_TDN)CALL StopTimer(dual_problem%timers%array_update_timer)
      dual_problem%error_tolerances%flows_infeasibility_norm=SQRT(DOT_PRODUCT(arrays%nodes_excess_flows(-dual_problem%n_free_speci&
     &al_nodes:),&
      arrays%nodes_excess_flows(-dual_problem%n_free_special_nodes:)))
      IF(dual_problem%method==SSCNO_TDN)CALL StopTimer(dual_problem%timers%array_update_timer)
      
      IF(error_tolerances(2)>0.0_r_wp)THEN
      convergence_indicators(2)=(dual_problem%error_tolerances%flows_infeasibility_norm<error_tolerances(2))
      END IF
      IF(iteration==1)initial_flows_inf_norm=dual_problem%error_tolerances%flows_infeasibility_norm
      
      END IF
      
      LagrangianGradient:IF(dual_problem%method==SSCNO_TDN)THEN
      
      dual_problem%error_tolerances%relative_error=dual_problem%error_tolerances%flows_infeasibility_norm/initial_flows_inf_norm
      relative_convergence=(dual_problem%error_tolerances%flows_infeasibility_norm<=relative_tolerance*initial_flows_inf_norm)
      IF(ALL(convergence_indicators).AND.relative_convergence)THEN
      
      CYCLE TCGN
      END IF
      
      eta=MAX(10.0_r_wp*EPSILON(1.0_r_wp),&
      MIN(dual_problem%forcing_prefactor/REAL(iteration,r_wp)**dual_problem%forcing_exponent,&
      dual_problem%error_tolerances%flows_infeasibility_norm))
      WRITE(message_log_unit,"(A,I5,A,F10.5,A,E10.3)")"DTN iteration ",iteration,&
      " has relative error ",dual_problem%error_tolerances%relative_error," and forcing sequence",eta
      
      ELSE LagrangianGradient
      
      CALL StartTimer(dual_problem%timers%cost_evaluation_timer)
      arrays%arcs_flows(0)=0.0_r_wp
      CALL DualCostDerivatives(ElementalCosts=ElementalCosts,&
      cost_function=cost_function,lagrangian_function=primal_lagrangian,&
      arc_offset=n_special_arcs,arcs_voltages=arrays%arcs_voltages,arcs_flows=arrays%arcs_flows,&
      arcs_conductances=arrays%arcs_conductances,warm_starts=dual_problem%warm_starts,&
      conductances_interval=(/cost_function%low_conductance,cost_function%high_conductance/))
      arrays%arcs_voltages(0)=0.0_r_wp
      arrays%arcs_conductances(0)=1.0_r_wp/cost_function%typical_resistance
      CALL StopTimer(dual_problem%timers%cost_evaluation_timer)
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL ArcsVoltages(heads_tails=subgraph%heads_tails,node_offset=subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-subgraph%n_special_nodes:subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(-subgraph%n_special_arcs:subgraph%n_arcs))
      CALL ArcsVoltages(heads_tails=special_subgraph%heads_tails,node_offset=special_subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-special_subgraph%n_special_nodes:special_subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(special_subgraph%n_special_arcs:(-special_subgraph%n_arcs):(-1)))
      arrays%arcs_excess_voltages(0)=0.0_r_wp
      CALL VectorSubtraction(from=arrays%arcs_voltages,what=arrays%arcs_excess_voltages,reverse=.TRUE.)
      
      CALL VectorMultiplication(first=arrays%arcs_conductances,second=arrays%arcs_excess_voltages,&
      product=arrays%arcs_excess_flows)
      CALL VectorSubtraction(from=arrays%arcs_excess_flows,what=arrays%arcs_flows,reverse=.FALSE.)
      
      CALL NodesExcessFlows(heads_tails=graph%heads_tails,node_offset=n_special_nodes,&
      arcs_flows=arrays%arcs_excess_flows,excess_flows=arrays%nodes_excess_flows)
      
      CALL VectorAddition(first=arrays%nodes_excess_flows(-dual_problem%n_free_special_nodes:),&
      second=arrays%supplies_demands(-dual_problem%n_free_special_nodes:))
      
      
      CALL StopTimer(dual_problem%timers%array_update_timer)
      IF(iteration==1)dual_problem%error_tolerances%excess_flows_norm=&
      SQRT(DOT_PRODUCT(arrays%arcs_flows,arrays%arcs_flows))
      eta=MAX(10.0_r_wp*EPSILON(1.0_r_wp),&
      MIN(dual_problem%forcing_prefactor/REAL(iteration,r_wp)**dual_problem%forcing_exponent,&
      dual_problem%error_tolerances%excess_flows_norm))
      
      END IF LagrangianGradient
      
      UsePCG:IF(dual_system%method==dual_network_pcg)THEN
      
      
      IF(solver%max_iterations<0)solver%max_iterations=n_nodes
      solver%min_iterations=1
      
      
      IF(dual_problem%method==SSCNO_TDN)THEN
      
      IF(iteration==1)solver%errors%residual_norm=dual_problem%error_tolerances%flows_infeasibility_norm
      solver%errors=CG_Errors(1.0_r_wp/REAL(iteration,r_wp)**dual_problem%forcing_exponent,&
      MIN(eta*dual_problem%error_tolerances%flows_infeasibility_norm,0.5_r_wp*solver%errors%residual_norm),-1.0_r_wp)
      ELSE
      
      IF(iteration==1)solver%errors%residual_norm=dual_problem%error_tolerances%excess_flows_norm
      solver%errors=CG_Errors(1.0_r_wp/REAL(iteration,r_wp)**dual_problem%forcing_exponent,&
      MIN(eta*dual_problem%error_tolerances%excess_flows_norm,0.5_r_wp*solver%errors%residual_norm),-1.0_r_wp)
      END IF
      WRITE(message_log_unit,*)"Requested tolerances from CG:",solver%errors
      
      CALL StartTimer(dual_problem%timers%solver_timer)
      arrays%nodes_excess_potentials=0.0_r_wp
      CALL SolveDualSystem_PCG(dual_system=dual_system)
      CALL StopTimer(dual_problem%timers%solver_timer)
      
      ELSE UsePCG
      
      CALL StartTimer(dual_problem%timers%solver_timer)
      CALL SolveDualSystem_LLt(dual_system=dual_system)
      CALL StopTimer(dual_problem%timers%solver_timer)
      
      END IF UsePCG
      
      
      NormalizePotentials:IF((dual_problem%n_fixed_potentials>0).AND.(.NOT.dual_system%ground_dummy))THEN
      
      CALL StartTimer(dual_problem%timers%solver_timer)
      CALL VectorShift(vector=arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      shift=-arrays%nodes_excess_potentials(0))
      arrays%nodes_excess_potentials(-n_special_nodes:-(dual_problem%n_free_special_nodes+1))=0.0_r_wp
      CALL StopTimer(dual_problem%timers%solver_timer)
      ELSE IF((dual_problem%n_fixed_potentials<=0).AND.(.NOT.dual_system%ground_dummy))THEN
      CALL StartTimer(dual_problem%timers%solver_timer)
      arrays%nodes_excess_potentials(0)=0.0_r_wp
      CALL VectorShift(vector=arrays%nodes_potentials,&
      shift=-SUM(arrays%nodes_potentials)/REAL(n_nodes+n_special_nodes,r_wp))
      arrays%nodes_excess_potentials(0)=0.0_r_wp
      CALL StopTimer(dual_problem%timers%solver_timer)
      END IF NormalizePotentials
      
      LineSearchDirection:IF(dual_problem%method==SSCNO_TDN)THEN
      
      
      IF(dual_problem%monitor_indicators(3))THEN
      dual_problem%error_tolerances%excess_potentials_norm=&
      SQRT(DOT_PRODUCT(arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes)))
      IF(error_tolerances(3)>0.0_r_wp)&
      convergence_indicators(3)=(dual_problem%error_tolerances%excess_potentials_norm<error_tolerances(3))
      END IF
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL ArcsVoltages(heads_tails=graph%heads_tails,node_offset=n_special_nodes,&
      nodes_potentials=arrays%nodes_excess_potentials,arcs_voltages=arrays%arcs_excess_voltages)
      arrays%arcs_excess_voltages(0)=0.0_r_wp
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      line_search%excess_power=DOT_PRODUCT(arrays%supplies_demands(-dual_problem%n_free_special_nodes:n_nodes),&
      arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes))
      
      ELSE LineSearchDirection
      
      
      
      
      IF(dual_problem%monitor_indicators(3))THEN
      dual_problem%error_tolerances%excess_potentials_norm=&
      SQRT(DOT_PRODUCT(arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes)))
      IF(error_tolerances(3)>0.0_r_wp)&
      convergence_indicators(3)=(dual_problem%error_tolerances%excess_potentials_norm<error_tolerances(3))
      END IF
      
      
      CALL VectorAddition(first=arrays%nodes_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      second=arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes))
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL ArcsVoltages(heads_tails=subgraph%heads_tails,node_offset=subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-subgraph%n_special_nodes:subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(-subgraph%n_special_arcs:subgraph%n_arcs))
      CALL ArcsVoltages(heads_tails=special_subgraph%heads_tails,node_offset=special_subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-special_subgraph%n_special_nodes:special_subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(special_subgraph%n_special_arcs:(-special_subgraph%n_arcs):(-1)))
      
      arrays%arcs_excess_voltages(0)=0.0_r_wp
      CALL VectorSubtraction(from=arrays%arcs_excess_voltages,what=arrays%arcs_voltages,reverse=.FALSE.)
      
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      dual_problem%error_tolerances%voltages_infeasibility_norm=&
      SQRT(DOT_PRODUCT(arrays%arcs_excess_voltages,arrays%arcs_excess_voltages))
      IF(iteration==1)initial_voltages_inf_norm=dual_problem%error_tolerances%voltages_infeasibility_norm
      IF(error_tolerances(1)>0.0_r_wp)&
      convergence_indicators(1)=(dual_problem%error_tolerances%voltages_infeasibility_norm<error_tolerances(1))
      
      dual_problem%error_tolerances%relative_error=dual_problem%error_tolerances%voltages_infeasibility_norm/initial_voltages_inf_&
     &norm
      relative_convergence=(dual_problem%error_tolerances%voltages_infeasibility_norm<=relative_tolerance*initial_voltages_inf_nor&
     &m)
      IF(ALL(convergence_indicators).AND.relative_convergence)THEN
      
      CYCLE TCGN
      END IF
      WRITE(message_log_unit,"(A,I5,A,F10.5,A,E10.3)")"TSQP iteration ",iteration,&
      " with relative error ",dual_problem%error_tolerances%relative_error," and forcing sequence",eta
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL VectorMultiplication(first=arrays%arcs_excess_voltages,second=arrays%arcs_conductances,&
      product=arrays%arcs_excess_flows)
      
      dual_problem%error_tolerances%excess_flows_norm=&
      SQRT(DOT_PRODUCT(arrays%arcs_excess_flows,arrays%arcs_excess_flows))
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      
      CALL StartTimer(dual_problem%timers%array_update_timer)
      CALL ArcsVoltages(heads_tails=subgraph%heads_tails,node_offset=subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-subgraph%n_special_nodes:subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(-subgraph%n_special_arcs:subgraph%n_arcs))
      CALL ArcsVoltages(heads_tails=special_subgraph%heads_tails,node_offset=special_subgraph%n_special_nodes,&
      nodes_potentials=arrays%nodes_potentials(-special_subgraph%n_special_nodes:special_subgraph%n_nodes),&
      arcs_voltages=arrays%arcs_excess_voltages(special_subgraph%n_special_arcs:(-special_subgraph%n_arcs):(-1)))
      
      arrays%arcs_excess_voltages(0)=0.0_r_wp
      
      line_search%excess_power=DOT_PRODUCT(arrays%arcs_excess_voltages,arrays%arcs_excess_flows)
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      END IF LineSearchDirection
      
      PerformSearch:IF(dual_problem%perform_line_search)THEN
      
      IF(line_search%root_finder%max_iterations<0)line_search%root_finder%max_iterations=10
      
      line_search%tolerance=MAX(dual_problem%ls_tolerance_prefactor/REAL(iteration,r_wp)**dual_problem%ls_tolerance_exponent,&
100   .0_r_wp*EPSILON(1.0_r_wp))
      
      LineSearch:IF(dual_problem%method==SSCNO_TDN)THEN
      
      CALL StartTimer(dual_problem%timers%line_search_timer)
      CALL PerformDualSearch_Root(line_search=line_search,ElementalCosts=ElementalCosts,&
      initialize=.FALSE.)
      CALL StopTimer(dual_problem%timers%line_search_timer)
      CALL StartTimer(dual_problem%timers%array_update_timer)
      
      CALL VectorAXPBY(y=arrays%nodes_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      x=arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes),alpha=line_search%step_size)
      CALL StopTimer(dual_problem%timers%array_update_timer)
      
      IF(dual_problem%monitor_indicators(3))dual_problem%error_tolerances%excess_potentials_norm=line_search%step_size*dual_proble&
     &m%error_tolerances%excess_potentials_norm
      
      ELSE LineSearch
      
      IF(dual_problem%force_feasible.AND.(iteration==1))THEN
      
      line_search%step_size=1.0_r_wp
      CALL VectorAddition(first=arrays%arcs_flows,second=arrays%arcs_excess_flows)
      ELSE
      CALL StartTimer(dual_problem%timers%line_search_timer)
      CALL PerformDualSearch_Root(line_search=line_search,ElementalCosts=ElementalCosts,&
      initialize=.FALSE.)
      CALL StopTimer(dual_problem%timers%line_search_timer)
      CALL StartTimer(dual_problem%timers%array_update_timer)
      
      CALL VectorAXPBY(y=arrays%arcs_flows,x=arrays%arcs_excess_flows,alpha=line_search%step_size)
      CALL StopTimer(dual_problem%timers%array_update_timer)
      END IF
      
      dual_problem%error_tolerances%excess_flows_norm=line_search%step_size*dual_problem%error_tolerances%excess_flows_norm
      IF(error_tolerances(4)>0.0_r_wp)&
      convergence_indicators(4)=(dual_problem%error_tolerances%excess_flows_norm<error_tolerances(4))
      
      END IF LineSearch
      
      ELSE
      
      line_search%step_size=1.0_r_wp
      IF(dual_problem%method==SSCNO_TDN)THEN
      CALL VectorAddition(first=arrays%nodes_potentials(-dual_problem%n_free_special_nodes:n_nodes),&
      second=arrays%nodes_excess_potentials(-dual_problem%n_free_special_nodes:n_nodes))
      ELSE
      CALL VectorAddition(first=arrays%arcs_flows,second=arrays%arcs_excess_flows)
      END IF
      
      END IF PerformSearch
      
      END DO TCGN
      dual_problem%n_iterations=iteration
      
      
      ShiftPotentials:IF(dual_problem%n_fixed_potentials<=0)THEN
      
      CALL StartTimer(dual_problem%timers%solver_timer)
      arrays%nodes_potentials(0)=0.0_r_wp
      CALL VectorShift(vector=arrays%nodes_potentials,&
      shift=-SUM(arrays%nodes_potentials)/REAL(n_nodes+n_special_nodes,r_wp))
      arrays%nodes_potentials(0)=0.0_r_wp
      CALL StopTimer(dual_problem%timers%solver_timer)
      END IF ShiftPotentials
      
      
      IF(dual_problem%n_fixed_potentials>0)THEN
      CALL NodesExcessFlows(heads_tails=special_subgraph%heads_tails,node_offset=special_subgraph%n_special_nodes,&
      arcs_flows=arrays%arcs_flows(special_subgraph%n_special_arcs:(-special_subgraph%n_arcs):(-1)),&
      excess_flows=arrays%nodes_excess_flows(-special_subgraph%n_special_nodes:special_subgraph%n_nodes))
      CALL VectorCopy(source=arrays%nodes_excess_flows(-n_special_nodes:-(dual_problem%n_free_special_nodes+1)),&
      target=arrays%supplies_demands(-n_special_nodes:-(dual_problem%n_free_special_nodes+1)))
      WRITE(message_log_unit,*)"The final excess flows at the fixed-potential nodes are:",&
      arrays%supplies_demands(-n_special_nodes:-(dual_problem%n_free_special_nodes+1))
      END IF
      
      IF(dual_problem%log_SSCNO)THEN
      WRITE(message_print_unit,*)
      WRITE(message_print_unit,"(A)")&
      "_______________________________________End of TCGN statistics"
      WRITE(message_print_unit,*)
      END IF
      
      IF(dual_problem%convergence==SSCNO_converged)THEN
      WRITE(message_log_unit,*)"TCG Newton's method successfully converged after: ",&
      iteration," iterations."
      ELSE
      WRITE(message_log_unit,*)"TCG Newton's method did not convege."
      END IF
      
      END SUBROUTINE SolveDualSSCNO_TCGN
      
      
      
      END MODULE Dual_Newton_SSCNO