MODULE Dual_Network_Solvers
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Graph_Algorithms
      USE Network_Matrix_Operations
      USE Network_Data_Types
      USE SCOTCH_Interface
      USE CHACO_Interface
      USE CHACO_Constants
      USE TAUCS_Interface
      USE Support_Trees
      USE Conjugate_Gradient
      USE Vector_Operations
      USE Network_Graphics
      IMPLICIT NONE
      PUBLIC::InitializeDualSolver,InitializeDualPreconditioner,&
      DestroyDualSolver,DestroyDualPreconditioner,&
      SolveDualSystem_PCG,SolveDualSystem_LLt
      PUBLIC::Dual_Network_Preconditioner,Dual_Network_Factorization,Dual_System_Timers,&
      Dual_Network_System,Dual_System_Handle
      
      
      PRIVATE
      
      INTEGER,PARAMETER,PUBLIC::no_precond=0,diagonal_precond=1,&
      MST_BHBt_precond=2,MST_QR_precond=3,MST_LDLt_precond=4,&
      ST_precond=5,ST_MST_precond=6,&
      TAUCS_MWB_precond=7,TAUCS_LLt_precond=8
      
      
      PUBLIC::Support_Tree_Arrays,Support_Tree_Preconditioner
      INTEGER,PARAMETER,PUBLIC::Default_ordering=0,Matching_ordering=1,&
      SCOTCH_ordering=2,CHACO_ordering=3
      
      TYPE Support_Tree_Arrays
      INTEGER(KIND=i_wp),DIMENSION(:),POINTER::nodes_mapping=>NULL()
      
      LOGICAL(KIND=l_wp),DIMENSION(:),POINTER::crossing_arcs_mask=>NULL()
      
      REAL(KIND=r_wp),DIMENSION(:),POINTER::&
      ST_arcs_resistances=>NULL(),ST_potentials_and_flows=>NULL()
      ENDTYPE
      TYPE Support_Tree_Preconditioner
      TYPE(SCOTCH_Mapping_Ordering),POINTER::scotch_mapping=>NULL()
      TYPE(CHACO_Mapping_Partitioning),POINTER::chaco_mapping=>NULL()
      TYPE(Directed_Graph),POINTER::graph=>NULL()
      TYPE(Support_Tree_Mapping)::mapping
      TYPE(Support_Tree_Arrays),POINTER::arrays=>NULL()
      INTEGER(KIND=i_wp)::height=0,max_height=0,degree=2
      INTEGER(KIND=i_wp)::n_leafs=0,n_nodes=0
      INTEGER::mapping_method=Default_ordering
      LOGICAL::CHACO_sequence=.FALSE.
      LOGICAL::support_special_nodes=.TRUE.,allocated_ST_arrays=.FALSE.
      ENDTYPE
      
      TYPE Dual_Network_Preconditioner
      INTEGER::method=diagonal_precond
      TYPE(Spanning_Tree),POINTER::tree=>NULL()
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree=>NULL()
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_preconditioner=>NULL()
      LOGICAL::use_spanning_tree=.FALSE.,use_support_tree=.FALSE.,reoptimize_ST_MST=.TRUE.,&
      allocated_taucs=.FALSE.,allocated_tree=.FALSE.,allocated_support_tree=.FALSE.
      REAL(KIND=r_wp),DIMENSION(:),POINTER::&
      nodes_resistances=>NULL(),tree_arcs_excess_voltages=>NULL(),&
      nodes_conductances=>NULL(),nodes_multipliers=>NULL()
      
      ENDTYPE
      
      INTEGER,PARAMETER,PUBLIC::TAUCS_fact_llt=TAUCS_incomplete_fact,&
      TAUCS_fact_ll=TAUCS_complete_fact_ll,TAUCS_fact_mf=TAUCS_complete_fact_mf
      
      INTEGER,PARAMETER,PUBLIC::no_fact_ordering=0,TAUCS_fact_ordering=1,SCOTCH_fact_ordering=2
      
      TYPE Dual_Network_Factorization
      INTEGER::factorization=TAUCS_fact_mf,ordering=SCOTCH_fact_ordering
      TYPE(Fill_Reducing_Ordering),POINTER::graph_ordering
      TYPE(SCOTCH_Mapping_Ordering),POINTER::scotch_ordering=>NULL()
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_factorization=>NULL()
      LOGICAL::reorder_nodes=.TRUE.,allocated_ordering=.FALSE.,&
      allocated_scotch_ordering=.FALSE.,allocated_factorization=.FALSE.
      INTEGER(KIND=i_wp),DIMENSION(:),POINTER::nodes_reordering=>NULL(),&
      nodes_renumbering=>NULL()
      
      ENDTYPE
      
      TYPE Dual_System_Timers
      TYPE(CG_Timers)::pcg_timers
      INTEGER::precond_initialization_timer=-1,precond_reoptimization_timer=-1,&
      precond_factorization_timer=-1,pcg_timer=-1,&
      llt_ordering_timer=-1,llt_factorization_timer=-1,llt_solver_timer=-1
      
      ENDTYPE
      
      INTEGER,PARAMETER,PUBLIC::dual_network_LLt=1,dual_network_PCG=2
      
      TYPE Dual_Network_System
      INTEGER::PIN=1
      INTEGER::method=dual_network_LLt
      TYPE(CG_System)::system
      TYPE(Directed_Graph),POINTER::graph=>NULL()
      TYPE(Network_Arrays),POINTER::arrays=>NULL()
      TYPE(Dual_Network_Preconditioner)::preconditioner
      TYPE(Dual_Network_Factorization)::factorization
      TYPE(Dual_System_Timers)::timers
      REAL(KIND=r_wp)::grounding_conductance=1.0_r_wp
      LOGICAL::rooted_graph=.FALSE.,ground_dummy=.FALSE.,&
      allocated_solver=.FALSE.
      ENDTYPE
      
      TYPE Dual_System_Handle
      TYPE(Dual_Network_System),POINTER::system=>NULL()
      ENDTYPE
      TYPE(Dual_System_Handle),DIMENSION(:),ALLOCATABLE,PUBLIC::&
      dual_systems
      
      CONTAINS
      
      
      SUBROUTINE InitializeDualSolver(dual_system)
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(CG_Solver),POINTER::solver
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(SCOTCH_Mapping_Ordering),POINTER::scotch_ordering
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_factorization
      
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs
      LOGICAL::valid_problem
      INTEGER::alloc_status
      
      graph=>dual_system%graph
      arrays=>dual_system%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_excess_flows).AND.&
      ASSOCIATED(arrays%nodes_excess_potentials)))valid_problem=.FALSE.
      END IF
      IF(.NOT.valid_problem)THEN
      CALL CriticalError(message="Network system not fully specified",&
      caller="InitializeDualSolver")
      RETURN
      END IF
      
      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((n_special_nodes>=0).AND.(n_nodes>=0))THEN
      dual_system%rooted_graph=.TRUE.
      ELSE
      dual_system%ground_dummy=.FALSE.
      END IF
      
      UsePCG:IF(dual_system%method==dual_network_pcg)THEN
      
      
      IF(.NOT.ASSOCIATED(dual_system%system%solver))THEN
      dual_system%allocated_solver=.TRUE.
      ALLOCATE(dual_system%system%solver)
      END IF
      solver=>dual_system%system%solver
      solver%timers=dual_system%timers%pcg_timers
      solver%lb=-n_special_nodes
      solver%ub=n_nodes
      solver%cg_x=>arrays%nodes_excess_potentials
      dual_system%system%cg_b=>arrays%nodes_excess_flows
      
      
      
      IF(.NOT.ASSOCIATED(solver%cg_Vy))THEN
      ALLOCATE(solver%cg_Vy(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(solver%cg_Vx))THEN
      ALLOCATE(solver%cg_Vx(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(solver%cg_residuals))THEN
      ALLOCATE(solver%cg_residuals(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(solver%cg_Vz))THEN
      ALLOCATE(solver%cg_Vz(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(dual_system%system%matrix))THEN
      ALLOCATE(dual_system%system%matrix)
      END IF
      
      
      IF(.NOT.ALLOCATED(dual_systems))THEN
      ALLOCATE(dual_systems(1:MAX(10,dual_system%PIN)))
      END IF
      
      ELSE UsePCG
      
      IF(.NOT.ASSOCIATED(dual_system%factorization%taucs_factorization))THEN
      dual_system%factorization%allocated_factorization=.TRUE.
      ALLOCATE(dual_system%factorization%taucs_factorization,STAT=alloc_status)
      END IF
      taucs_factorization=>dual_system%factorization%taucs_factorization
      
      taucs_factorization%graph=>graph
      taucs_factorization%conductances=>arrays%arcs_conductances
      taucs_factorization%excess_potentials=>arrays%nodes_excess_potentials
      taucs_factorization%excess_flows=>arrays%nodes_excess_flows
      taucs_factorization%precond_method=TAUCS_LLt
      taucs_factorization%factorization=dual_system%factorization%factorization
      
      
      IF(.NOT.ASSOCIATED(taucs_factorization%eliminated_nodes))THEN
      
      ALLOCATE(taucs_factorization%eliminated_nodes(1:1),STAT=alloc_status)
      GroundGraph:IF(dual_system%rooted_graph)THEN
      taucs_factorization%eliminated_nodes=(/0/)
      ELSE
      taucs_factorization%eliminated_nodes=(/-n_special_nodes/)
      END IF GroundGraph
      END IF
      taucs_factorization%diagonal_conductance=dual_system%grounding_conductance
      
      
      IF((dual_system%factorization%ordering==SCOTCH_fact_ordering).OR.&
      (dual_system%factorization%ordering==no_fact_ordering))&
      taucs_factorization%ordering_method="none"
      
      dual_system%factorization%reorder_nodes=.FALSE.
      IF(dual_system%factorization%ordering==SCOTCH_fact_ordering)THEN
      
      dual_system%factorization%reorder_nodes=.TRUE.
      IF(.NOT.ASSOCIATED(dual_system%factorization%graph_ordering))THEN
      dual_system%factorization%allocated_ordering=.TRUE.
      ALLOCATE(dual_system%factorization%graph_ordering,STAT=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(dual_system%factorization%nodes_renumbering))THEN
      ALLOCATE(dual_system%factorization%nodes_renumbering(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(dual_system%factorization%nodes_reordering))THEN
      ALLOCATE(dual_system%factorization%nodes_reordering(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="InitializeDualSolver",alloc_status=alloc_status)
      END IF
      dual_system%factorization%graph_ordering%nodes_renumbering=>dual_system%factorization%nodes_renumbering
      dual_system%factorization%graph_ordering%nodes_reordering=>dual_system%factorization%nodes_reordering
      
      taucs_factorization%graph_ordering=>dual_system%factorization%graph_ordering
      
      END IF
      
      UseSCOTCH:IF(dual_system%factorization%ordering==SCOTCH_fact_ordering)THEN
      
      
      IF(.NOT.ASSOCIATED(dual_system%factorization%scotch_ordering))THEN
      dual_system%factorization%allocated_scotch_ordering=.TRUE.
      ALLOCATE(dual_system%factorization%scotch_ordering,STAT=alloc_status)
      END IF
      scotch_ordering=>dual_system%factorization%scotch_ordering
      
      scotch_ordering%graph=>graph
      scotch_ordering%map_graph=.FALSE.
      scotch_ordering%order_graph=.TRUE.
      scotch_ordering%graph_ordering=>dual_system%factorization%graph_ordering
      
      IF(scotch_ordering%scotch_ordering%ordering_string=="")THEN
      scotch_ordering%scotch_ordering%ordering_string="m"
      END IF
      WRITE(message_log_unit,*)&
      "Using SCOTCH fill-reducing ordering ",TRIM(scotch_ordering%scotch_ordering%ordering_string)
      CALL StartTimer(dual_system%timers%llt_ordering_timer)
      CALL SCOTCH_InitializeMappingOrdering(mapping_ordering=scotch_ordering)
      CALL SCOTCH_ComputeMappingOrdering(mapping_ordering=scotch_ordering)
      CALL StopTimer(dual_system%timers%llt_ordering_timer)
      
      CALL SCOTCH_DestroyMappingOrdering(mapping_ordering=scotch_ordering)
      
      IF(dual_system%factorization%allocated_scotch_ordering)THEN
      DEALLOCATE(dual_system%factorization%scotch_ordering,STAT=alloc_status)
      END IF
      
      
      END IF UseSCOTCH
      
      WRITE(message_log_unit,*)&
      "Initializing a complete factorization structure in TAUCS."
      CALL TAUCS_InitializePreconditioner(taucs_preconditioner=taucs_factorization)
      
      END IF UsePCG
      
      END SUBROUTINE InitializeDualSolver
      
      
      
      SUBROUTINE InitializeDualPreconditioner(dual_system)
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      
      INTEGER::alloc_status
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(CG_Solver),POINTER::solver
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_preconditioner
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs,node,arc
      
      IF(dual_system%method/=dual_network_pcg)THEN
      CALL Warning("Preconditioner initialization attempted without PCG",&
      caller="InitializeDualPreconditioner")
      RETURN
      END IF
      
      graph=>dual_system%graph
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      arrays=>dual_system%arrays
      solver=>dual_system%system%solver
      
      CALL StartTimer(dual_system%timers%precond_initialization_timer)
      
      IF((dual_system%preconditioner%method==MST_BHBt_precond).OR.(dual_system%preconditioner%method==ST_MST_precond).OR.&
      (dual_system%preconditioner%method==MST_QR_precond).OR.(dual_system%preconditioner%method==MST_LDLt_precond))THEN
      
      dual_system%preconditioner%use_spanning_tree=.TRUE.
      END IF
      IF((dual_system%preconditioner%method==ST_precond).OR.(dual_system%preconditioner%method==ST_MST_precond))THEN
      
      dual_system%preconditioner%use_support_tree=.TRUE.
      END IF
      
      
      IF((dual_system%preconditioner%method==MST_QR_precond).OR.(dual_system%preconditioner%method==MST_LDLt_precond).OR.&
      (dual_system%preconditioner%method==diagonal_precond).OR.(dual_system%preconditioner%method==ST_precond).OR.&
      (dual_system%preconditioner%method==ST_MST_precond))THEN
      
      
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%nodes_resistances))THEN
      ALLOCATE(dual_system%preconditioner%nodes_resistances(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      IF((dual_system%preconditioner%method==MST_QR_precond).OR.(dual_system%preconditioner%method==MST_LDLt_precond).OR.&
      (dual_system%preconditioner%method==ST_MST_precond))THEN
      
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%nodes_multipliers))THEN
      ALLOCATE(dual_system%preconditioner%nodes_multipliers(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      END IF
      
      END IF
      
      BuildST:IF(dual_system%preconditioner%use_support_tree)THEN
      WRITE(message_log_unit,*)"Initializing a support tree for the network preconditioner"
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%support_tree))THEN
      dual_system%preconditioner%allocated_support_tree=.TRUE.
      ALLOCATE(dual_system%preconditioner%support_tree,STAT=alloc_status)
      END IF
      support_tree=>dual_system%preconditioner%support_tree
      
      IF(.NOT.ASSOCIATED(support_tree%arrays))THEN
      support_tree%allocated_ST_arrays=.TRUE.
      ALLOCATE(support_tree%arrays,STAT=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(support_tree%arrays%nodes_mapping))THEN
      ALLOCATE(support_tree%arrays%nodes_mapping(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(support_tree%arrays%crossing_arcs_mask))THEN
      ALLOCATE(support_tree%arrays%crossing_arcs_mask(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=.TRUE._l_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(support_tree%mapping_method==CHACO_ordering)THEN
      IF(.NOT.ASSOCIATED(support_tree%chaco_mapping))THEN
      ALLOCATE(support_tree%chaco_mapping,STAT=alloc_status)
      END IF
      IF(ASSOCIATED(arrays%nodes_coordinates))THEN
      support_tree%chaco_mapping%coordinates=>arrays%nodes_coordinates(:,1:n_nodes)
      END IF
      ELSE IF(support_tree%mapping_method==SCOTCH_ordering)THEN
      IF(.NOT.ASSOCIATED(support_tree%scotch_mapping))THEN
      ALLOCATE(support_tree%scotch_mapping,STAT=alloc_status)
      END IF
      END IF
      
      InitializeST:IF((support_tree%mapping_method==SCOTCH_ordering).OR.&
      (support_tree%mapping_method==CHACO_ordering).OR.&
      (support_tree%mapping_method==Matching_ordering))THEN
      IF(support_tree%support_special_nodes.OR.(.NOT.dual_system%rooted_graph))THEN
      
      support_tree%support_special_nodes=.TRUE.
      support_tree%graph=>graph
      support_tree%mapping%nodes_mapping=>support_tree%arrays%nodes_mapping
      support_tree%mapping%arcs_weights=>dual_system%arrays%arcs_conductances
      
      ELSE
      
      ALLOCATE(support_tree%graph,STAT=alloc_status)
      support_tree%graph%heads_tails=>graph%heads_tails(1:2,1:n_arcs)
      support_tree%graph%n_special_nodes=-1
      support_tree%graph%n_special_arcs=-1
      support_tree%graph%n_nodes=n_nodes
      support_tree%graph%n_arcs=n_arcs
      support_tree%mapping%nodes_mapping=>support_tree%arrays%nodes_mapping(1:n_nodes)
      support_tree%mapping%arcs_weights=>dual_system%arrays%arcs_conductances(1:n_arcs)
      
      support_tree%arrays%nodes_mapping(-n_special_nodes:0)=0
      support_tree%arrays%crossing_arcs_mask(-n_special_arcs:0)=.FALSE._l_wp
      
      END IF
      END IF InitializeST
      
      IF(support_tree%mapping_method==SCOTCH_ordering)THEN
      support_tree%scotch_mapping%graph=>support_tree%graph
      support_tree%scotch_mapping%ST_mapping=>support_tree%mapping
      ELSE IF(support_tree%mapping_method==CHACO_ordering)THEN
      support_tree%chaco_mapping%graph=>support_tree%graph
      support_tree%chaco_mapping%ST_mapping=>support_tree%mapping
      END IF
      
      IF(support_tree%height<0)THEN
      support_tree%height=INT(LOG(REAL(n_nodes))/LOG(REAL(support_tree%degree)),i_wp)
      IF(support_tree%max_height>0)THEN
      support_tree%height=MAX(1_i_wp,MIN(support_tree%max_height,support_tree%height))
      ELSE
      support_tree%height=MAX(1_i_wp,support_tree%height)
      END IF
      END IF
      support_tree%n_leafs=support_tree%degree**support_tree%height
      support_tree%n_nodes=((support_tree%n_leafs-1)*support_tree%degree)/(support_tree%degree-1)
      
      
      IF(.NOT.ASSOCIATED(support_tree%arrays%ST_arcs_resistances))THEN
      ALLOCATE(support_tree%arrays%ST_arcs_resistances(0:support_tree%n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(support_tree%n_nodes)-(0)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(support_tree%arrays%ST_potentials_and_flows))THEN
      ALLOCATE(support_tree%arrays%ST_potentials_and_flows(0:support_tree%n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(support_tree%n_nodes)-(0)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      SELECT CASE(support_tree%mapping_method)
      CASE(SCOTCH_ordering)
      
      support_tree%scotch_mapping%map_graph=.TRUE.
      support_tree%scotch_mapping%order_graph=.FALSE.
      IF(support_tree%scotch_mapping%scotch_mapping%architecture_string=="")THEN
      WRITE(support_tree%scotch_mapping%scotch_mapping%architecture_string,*)"leaf ",support_tree%height,support_tree%height," 1"
      END IF
      IF(support_tree%scotch_mapping%scotch_mapping%mapping_string=="")THEN
      support_tree%scotch_mapping%scotch_mapping%mapping_string="b{job=t,map=t,poli=S,"//&
      "strat=m{low=(g x),asc=f{type=b} x,rat=0.95,type=h,vert=100}}"
      END IF
      WRITE(message_log_unit,*)&
      "Using SCOTCH partitioning with architecture ",TRIM(support_tree%scotch_mapping%scotch_mapping%architecture_string),&
      " and strategy ",TRIM(support_tree%scotch_mapping%scotch_mapping%mapping_string)
      CALL SCOTCH_InitializeMappingOrdering(mapping_ordering=support_tree%scotch_mapping)
      
      CASE(CHACO_ordering)
      
      CHACOSequence:IF(support_tree%CHACO_sequence)THEN
      
      support_tree%chaco_mapping%strategy%sequence=1
      support_tree%chaco_mapping%strategy%rqi_flag=CHACO_Lanczos
      support_tree%chaco_mapping%strategy%vmax=MAX(1_i_wp,n_nodes/10_i_wp)
      support_tree%chaco_mapping%architecture%ndims_tot=support_tree%n_leafs
      WRITE(message_log_unit,*)"Using spectral CHACO to sequence the nodes into",&
      support_tree%n_leafs," slices."
      
      ELSE
      
      support_tree%chaco_mapping%strategy%sequence=0
      support_tree%chaco_mapping%architecture%architecture=CHACO_mesh1D
      support_tree%chaco_mapping%architecture%mesh_dims(1)=support_tree%n_leafs
      IF(support_tree%chaco_mapping%strategy%global_method==CHACO_dummy)THEN
      IF(ASSOCIATED(arrays%nodes_coordinates))THEN
      support_tree%chaco_mapping%strategy%global_method=CHACO_inertial
      ELSE
      support_tree%chaco_mapping%strategy%global_method=CHACO_linear
      END IF
      END IF
      IF(support_tree%chaco_mapping%strategy%local_method==CHACO_dummy)THEN
      support_tree%chaco_mapping%strategy%local_method=CHACO_local_KL
      END IF
      support_tree%chaco_mapping%strategy%ndims=1
      IF(support_tree%chaco_mapping%strategy%vmax<=1)&
      support_tree%chaco_mapping%strategy%vmax=10
      WRITE(message_log_unit,*)&
      "Using CHACO partitioning with 1D mesh of length ",support_tree%n_leafs
      
      END IF CHACOSequence
      CALL CHACO_InitializeMapping(mapping_partitioning=support_tree%chaco_mapping)
      
      CASE(Matching_ordering)
      
      
      
      CASE DEFAULT
      
      IF(support_tree%support_special_nodes.OR.(.NOT.dual_system%rooted_graph))THEN
      
      WRITE(message_log_unit,"(A,F5.1,A)")" Using default (SFC) nodes ordering with",&
      REAL(n_nodes+n_special_nodes)/REAL(support_tree%n_leafs)," nodes per partition"
      DO node=-n_special_nodes,n_nodes
      support_tree%arrays%nodes_mapping(node)=INT(REAL(node+n_special_nodes)*&
      REAL(support_tree%n_leafs)/REAL(n_nodes+n_special_nodes+1))+&
      (support_tree%n_nodes-support_tree%n_leafs)+1
      END DO
      IF(dual_system%ground_dummy.AND.dual_system%rooted_graph)support_tree%arrays%nodes_mapping(0)=0
      CALL ST_ArcsMask(node_offset=n_special_nodes,arc_offset=n_special_arcs,&
      heads_tails=graph%heads_tails,arcs_mask=support_tree%arrays%crossing_arcs_mask,&
      nodes_mapping=support_tree%arrays%nodes_mapping,mask_internal_arcs=.TRUE.)
      
      ELSE
      
      WRITE(message_log_unit,"(A,F5.1,A)")" Using default (SFC) nodes ordering with",&
      REAL(n_nodes)/REAL(support_tree%n_leafs)," nodes per partition"
      DO node=1,n_nodes
      support_tree%arrays%nodes_mapping(node)=INT(REAL(node-1)*REAL(support_tree%n_leafs)/REAL(n_nodes+1))+&
      (support_tree%n_nodes-support_tree%n_leafs)+1
      END DO
      CALL ST_ArcsMask(node_offset=-1,arc_offset=-1,&
      heads_tails=graph%heads_tails(:,1:),arcs_mask=support_tree%arrays%crossing_arcs_mask(1:),&
      nodes_mapping=support_tree%arrays%nodes_mapping(1:),mask_internal_arcs=.TRUE.)
      
      support_tree%arrays%nodes_mapping(-n_special_nodes:0)=0
      support_tree%arrays%crossing_arcs_mask(-n_special_arcs:0)=.FALSE._l_wp
      
      END IF
      
      ENDSELECT
      
      END IF BuildST
      
      BuildMST:IF(dual_system%preconditioner%use_spanning_tree)THEN
      WRITE(message_log_unit,*)&
      "Initializing a spanning forest for the network preconditioner"
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%tree))THEN
      dual_system%preconditioner%allocated_tree=.TRUE.
      ALLOCATE(dual_system%preconditioner%tree)
      END IF
      tree=>dual_system%preconditioner%tree
      
      tree%graph=>graph
      tree%PIN=graph%PIN
      
      
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%tree_arcs_excess_voltages))THEN
      ALLOCATE(dual_system%preconditioner%tree_arcs_excess_voltages(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(tree%parents))THEN
      ALLOCATE(tree%parents(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(tree%cardinalities))THEN
      ALLOCATE(tree%cardinalities(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(tree%level_ordering))THEN
      ALLOCATE(tree%level_ordering(0:n_nodes+n_special_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes+n_special_nodes)-(0)+1,mold=1_i_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(tree%orientations))THEN
      ALLOCATE(tree%orientations(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_byte,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(tree%path_labels))THEN
      ALLOCATE(tree%path_labels(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(tree%mask))THEN
      ALLOCATE(tree%mask(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_l_wp,&
      caller="InitializeDualPreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(dual_system%preconditioner%method==ST_MST_precond)THEN
      
      IF(dual_system%preconditioner%reoptimize_ST_MST)THEN
      DO node=-n_special_nodes,n_nodes
      tree%parents(node)=node
      END DO
      tree%orientations=no_parent
      END IF
      ELSE
      CALL BuildNetworkMST(arc_offset=tree%graph%n_special_arcs,node_offset=tree%graph%n_special_nodes,&
      heads_tails=dual_system%preconditioner%tree%graph%heads_tails,&
      orientations=tree%orientations,level_ordering=tree%level_ordering,&
      parents=tree%parents,cardinalities=tree%cardinalities,&
      path_labels=tree%path_labels,tree_mask=tree%mask,&
      from_scratch=.TRUE.,tree_type=any_tree)
      END IF
      
      END IF BuildMST
      
      InitializeTAUCS:IF((dual_system%preconditioner%method==TAUCS_LLt_precond).OR.&
      (dual_system%preconditioner%method==TAUCS_MWB_precond))THEN
      
      IF(.NOT.ASSOCIATED(dual_system%preconditioner%taucs_preconditioner))THEN
      dual_system%preconditioner%allocated_taucs=.TRUE.
      ALLOCATE(dual_system%preconditioner%taucs_preconditioner,STAT=alloc_status)
      END IF
      taucs_preconditioner=>dual_system%preconditioner%taucs_preconditioner
      
      taucs_preconditioner%graph=>graph
      taucs_preconditioner%conductances=>arrays%arcs_conductances
      taucs_preconditioner%excess_potentials=>solver%cg_Vz
      taucs_preconditioner%excess_flows=>solver%cg_residuals
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%eliminated_nodes))THEN
      ALLOCATE(taucs_preconditioner%eliminated_nodes(1:1),STAT=alloc_status)
      GroundGraph:IF(dual_system%rooted_graph)THEN
      taucs_preconditioner%eliminated_nodes=(/0/)
      ELSE
      taucs_preconditioner%eliminated_nodes=(/-n_special_nodes/)
      END IF GroundGraph
      END IF
      taucs_preconditioner%diagonal_conductance=dual_system%grounding_conductance
      
      IF(dual_system%preconditioner%method==TAUCS_LLt_precond)THEN
      WRITE(message_log_unit,*)&
      "Initializing a TAUCS incomplete LLt factorization preconditioner."
      taucs_preconditioner%precond_method=TAUCS_LLt
      taucs_preconditioner%factorization=TAUCS_incomplete_fact
      ELSE
      taucs_preconditioner%precond_method=TAUCS_MWB
      taucs_preconditioner%factorization=TAUCS_complete_fact_mf
      WRITE(message_log_unit,*)&
      "Initializing a TAUCS MWB support-graph preconditioner."
      END IF
      
      CALL TAUCS_InitializePreconditioner(taucs_preconditioner=taucs_preconditioner)
      
      END IF InitializeTAUCS
      
      CALL StopTimer(dual_system%timers%precond_initialization_timer)
      
      END SUBROUTINE InitializeDualPreconditioner
      
      
      
      SUBROUTINE DestroyDualSolver(dual_system)
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      INTEGER::alloc_status
      TYPE(CG_Solver),POINTER::solver
      TYPE(Network_Arrays),POINTER::arrays
      
      solver=>dual_system%system%solver
      arrays=>dual_system%arrays
      
      UsePCG:IF(dual_system%method==dual_network_pcg)THEN
      
      IF(dual_system%allocated_solver)THEN
      solver=>dual_system%system%solver
      
      IF(ASSOCIATED(solver%cg_Vy))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(solver%cg_Vy),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(solver%cg_Vy,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(solver%cg_Vx))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(solver%cg_Vx),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(solver%cg_Vx,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(solver%cg_residuals))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(solver%cg_residuals),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(solver%cg_residuals,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(solver%cg_Vz))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(solver%cg_Vz),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(solver%cg_Vz,STAT=alloc_status)
      END IF
      DEALLOCATE(dual_system%system%solver,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(dual_system%system%matrix))THEN
      DEALLOCATE(dual_system%system%matrix,STAT=alloc_status)
      END IF
      
      ELSE UsePCG
      
      IF(dual_system%factorization%allocated_ordering)THEN
      
      IF(ASSOCIATED(dual_system%factorization%nodes_renumbering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(dual_system%factorization%nodes_renumbering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(dual_system%factorization%nodes_renumbering,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(dual_system%factorization%nodes_reordering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(dual_system%factorization%nodes_reordering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(dual_system%factorization%nodes_reordering,STAT=alloc_status)
      END IF
      DEALLOCATE(dual_system%factorization%graph_ordering,STAT=alloc_status)
      END IF
      
      CALL TAUCS_DestroyPreconditioner(taucs_preconditioner=dual_system%factorization%taucs_factorization)
      IF(dual_system%factorization%allocated_factorization)THEN
      DEALLOCATE(dual_system%factorization%taucs_factorization,STAT=alloc_status)
      END IF
      
      END IF UsePCG
      
      END SUBROUTINE DestroyDualSolver
      
      
      
      SUBROUTINE DestroyDualPreconditioner(dual_system)
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_preconditioner
      INTEGER::alloc_status
      
      
      IF(ASSOCIATED(dual_system%preconditioner%nodes_resistances))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(dual_system%preconditioner%nodes_resistances),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(dual_system%preconditioner%nodes_resistances,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(dual_system%preconditioner%nodes_multipliers))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(dual_system%preconditioner%nodes_multipliers),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(dual_system%preconditioner%nodes_multipliers,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(dual_system%preconditioner%tree_arcs_excess_voltages))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(dual_system%preconditioner%tree_arcs_excess_voltages),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(dual_system%preconditioner%tree_arcs_excess_voltages,STAT=alloc_status)
      END IF
      
      tree=>dual_system%preconditioner%tree
      IF(dual_system%preconditioner%allocated_tree)THEN
      
      IF(ASSOCIATED(tree%parents))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%parents),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(tree%parents,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(tree%cardinalities))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%cardinalities),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(tree%cardinalities,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(tree%orientations))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%orientations),KIND=i_wp),mold=1_i_byte)
      DEALLOCATE(tree%orientations,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(tree%level_ordering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%level_ordering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(tree%level_ordering,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(tree%path_labels))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%path_labels),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(tree%path_labels,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(tree%mask))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tree%mask),KIND=i_wp),mold=1_l_wp)
      DEALLOCATE(tree%mask,STAT=alloc_status)
      END IF
      DEALLOCATE(dual_system%preconditioner%tree,STAT=alloc_status)
      END IF
      
      support_tree=>dual_system%preconditioner%support_tree
      IF(dual_system%preconditioner%use_support_tree)THEN
      SELECT CASE(support_tree%mapping_method)
      CASE(SCOTCH_ordering)
      CALL SCOTCH_DestroyMappingOrdering(mapping_ordering=support_tree%scotch_mapping)
      CASE(CHACO_ordering)
      CALL CHACO_DestroyMapping(mapping_partitioning=support_tree%chaco_mapping)
      ENDSELECT
      IF(support_tree%allocated_ST_arrays)THEN
      
      IF(ASSOCIATED(support_tree%arrays%nodes_mapping))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(support_tree%arrays%nodes_mapping),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(support_tree%arrays%nodes_mapping,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(support_tree%arrays%crossing_arcs_mask))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(support_tree%arrays%crossing_arcs_mask),KIND=i_wp),mold=.TRUE._l_wp)
      DEALLOCATE(support_tree%arrays%crossing_arcs_mask,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(support_tree%arrays%ST_arcs_resistances))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(support_tree%arrays%ST_arcs_resistances),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(support_tree%arrays%ST_arcs_resistances,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(support_tree%arrays%ST_potentials_and_flows))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(support_tree%arrays%ST_potentials_and_flows),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(support_tree%arrays%ST_potentials_and_flows,STAT=alloc_status)
      END IF
      DEALLOCATE(support_tree%arrays,STAT=alloc_status)
      END IF
      IF(.NOT.support_tree%support_special_nodes)DEALLOCATE(support_tree%graph,STAT=alloc_status)
      END IF
      IF(dual_system%preconditioner%allocated_support_tree)THEN
      DEALLOCATE(dual_system%preconditioner%support_tree,STAT=alloc_status)
      END IF
      
      taucs_preconditioner=>dual_system%preconditioner%taucs_preconditioner
      IF(ASSOCIATED(taucs_preconditioner))THEN
      DEALLOCATE(taucs_preconditioner%eliminated_nodes)
      CALL TAUCS_DestroyPreconditioner(taucs_preconditioner=taucs_preconditioner)
      IF(.NOT.ASSOCIATED(taucs_preconditioner%graph,dual_system%graph))DEALLOCATE(taucs_preconditioner%graph,STAT=alloc_status)
      END IF
      IF(dual_system%preconditioner%allocated_taucs)DEALLOCATE(dual_system%preconditioner%taucs_preconditioner,STAT=alloc_status)
      
      END SUBROUTINE DestroyDualPreconditioner
      
      
      
      SUBROUTINE SolveDualSystem_LLt(dual_system)
      IMPLICIT NONE
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_factorization
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs,node,arc
      
      graph=>dual_system%graph
      arrays=>dual_system%arrays
      taucs_factorization=>dual_system%factorization%taucs_factorization
      
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      
      CALL StartTimer(dual_system%timers%llt_factorization_timer)
      CALL TAUCS_CreatePreconditioner(taucs_preconditioner=taucs_factorization)
      CALL StopTimer(dual_system%timers%llt_factorization_timer)
      CALL StartTimer(dual_system%timers%llt_solver_timer)
      CALL TAUCS_ApplyPreconditioner(taucs_preconditioner=taucs_factorization)
      CALL StopTimer(dual_system%timers%llt_solver_timer)
      CALL TAUCS_FreePreconditioner(taucs_preconditioner=taucs_factorization)
      
      END SUBROUTINE SolveDualSystem_LLt
      
      
      
      SUBROUTINE SolveDualSystem_PCG(dual_system)
      IMPLICIT NONE
      TYPE(Dual_Network_System),INTENT(INOUT),TARGET::dual_system
      
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree
      TYPE(CG_System),POINTER::system
      TYPE(CG_Solver),POINTER::solver
      TYPE(Network_Arrays),POINTER::arrays
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs,node,arc
      
      IF((dual_system%PIN<INT(LBOUND(dual_systems,DIM=1),KIND=i_wp)).OR.&
      (dual_system%PIN>INT(UBOUND(dual_systems,DIM=1),KIND=i_wp)))THEN
      CALL CriticalError(message="Not enough dual system handles available or PIN out of range",&
      caller="InitializeDualSolver")
      RETURN
      ELSE
      IF(ASSOCIATED(dual_systems(dual_system%PIN)%system))THEN
      CALL CriticalError(message="Dual system handle already taken",&
      caller="InitializeDualSolver")
      RETURN
      ELSE
      dual_systems(dual_system%PIN)%system=>dual_system
      END IF
      END IF
      dual_system%system%PIN=dual_system%PIN
      dual_system%system%matrix%PIN=dual_system%PIN
      
      graph=>dual_system%graph
      tree=>dual_system%preconditioner%tree
      support_tree=>dual_system%preconditioner%support_tree
      system=>dual_system%system
      solver=>system%solver
      arrays=>dual_system%arrays
      
      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(dual_system%preconditioner%method==diagonal_precond)THEN
      
      CALL StartTimer(dual_system%timers%precond_factorization_timer)
      CALL Diagonal_ADAt(heads_tails=graph%heads_tails,node_offset=n_special_nodes,&
      arcs_conductances=arrays%arcs_conductances,nodes_resistances=dual_system%preconditioner%nodes_resistances)
      CALL StopTimer(dual_system%timers%precond_factorization_timer)
      END IF
      
      UpdateST:IF(dual_system%preconditioner%use_support_tree)THEN
      IF((support_tree%mapping_method==SCOTCH_ordering).OR.(support_tree%mapping_method==CHACO_ordering).OR.&
      (support_tree%mapping_method==Matching_ordering))THEN
      CALL StartTimer(dual_system%timers%precond_reoptimization_timer)
      SELECT CASE(support_tree%mapping_method)
      CASE(SCOTCH_ordering)
      CALL SCOTCH_ComputeMappingOrdering(mapping_ordering=support_tree%scotch_mapping,&
      mapping_offset=(support_tree%n_nodes-support_tree%n_leafs)+1)
      CASE(CHACO_ordering)
      CALL CHACO_ComputeMapping(mapping_partitioning=support_tree%chaco_mapping,&
      mapping_offset=(support_tree%n_nodes-support_tree%n_leafs)+1)
      CASE(Matching_ordering)
      CALL ST_MatchingNodesMapping(node_offset=support_tree%graph%n_special_nodes,&
      arc_offset=support_tree%graph%n_special_arcs,&
      heads_tails=support_tree%graph%heads_tails,&
      arcs_conductances=support_tree%mapping%arcs_weights,nodes_mapping=support_tree%mapping%nodes_mapping,&
      n_components=support_tree%n_leafs,mapping_offset=(support_tree%n_nodes-support_tree%n_leafs)+1,&
      heavy_matching=.TRUE.)
      ENDSELECT
      
      IF(dual_system%ground_dummy.AND.dual_system%rooted_graph)support_tree%arrays%nodes_mapping(0)=0
      CALL ST_ArcsMask(node_offset=support_tree%graph%n_special_nodes,&
      arc_offset=support_tree%graph%n_special_arcs,&
      heads_tails=support_tree%graph%heads_tails,&
      arcs_mask=support_tree%arrays%crossing_arcs_mask(-support_tree%graph%n_special_arcs:support_tree%graph%n_arcs),&
      nodes_mapping=support_tree%mapping%nodes_mapping,mask_internal_arcs=.TRUE.)
      CALL StopTimer(dual_system%timers%precond_reoptimization_timer)
      END IF
      
      
      CALL StartTimer(dual_system%timers%precond_factorization_timer)
      IF(dual_system%preconditioner%method==ST_MST_precond)THEN
      CALL ST_PropagateArcsConductances(height=support_tree%height,degree=support_tree%degree,&
      arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=graph%heads_tails,arcs_conductances=arrays%arcs_conductances,&
      nodes_mapping=support_tree%arrays%nodes_mapping,&
      ST_arcs_resistances=support_tree%arrays%ST_arcs_resistances,ST_nodes_temp=support_tree%arrays%ST_potentials_and_flows)
      ELSE
      CALL ST_PropagateArcsConductances(height=support_tree%height,degree=support_tree%degree,&
      arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=graph%heads_tails,arcs_conductances=arrays%arcs_conductances,&
      nodes_mapping=support_tree%arrays%nodes_mapping,nodes_resistances=dual_system%preconditioner%nodes_resistances,&
      ST_arcs_resistances=support_tree%arrays%ST_arcs_resistances,ST_nodes_temp=support_tree%arrays%ST_potentials_and_flows)
      END IF
      CALL StopTimer(dual_system%timers%precond_factorization_timer)
      
      
      END IF UpdateST
      
      UpdateMST:IF(dual_system%preconditioner%use_spanning_tree)THEN
      IF(dual_system%preconditioner%method==ST_MST_precond)THEN
      
      
      CALL StartTimer(dual_system%timers%precond_reoptimization_timer)
      ReoptimizeMST:IF(dual_system%preconditioner%reoptimize_ST_MST)THEN
      
      DO node=-n_special_nodes,n_nodes
      IF(.NOT.support_tree%arrays%crossing_arcs_mask(tree%parents(node)))THEN
      tree%orientations(node)=no_parent
      
      END IF
      END DO
      END IF ReoptimizeMST
      
      CALL BuildNetworkMST(arc_offset=tree%graph%n_special_arcs,node_offset=tree%graph%n_special_nodes,&
      heads_tails=dual_system%preconditioner%tree%graph%heads_tails,&
      arcs_mask=support_tree%arrays%crossing_arcs_mask(-tree%graph%n_special_arcs:tree%graph%n_arcs),&
      orientations=tree%orientations,level_ordering=tree%level_ordering,&
      parents=tree%parents,cardinalities=tree%cardinalities,&
      path_labels=tree%path_labels,tree_mask=tree%mask,&
      from_scratch=.NOT.dual_system%preconditioner%reoptimize_ST_MST,defragment_ordering=.FALSE.,&
      tree_type=any_tree)
      CALL BuildNetworkMST(arc_offset=tree%graph%n_special_arcs,node_offset=tree%graph%n_special_nodes,&
      heads_tails=dual_system%preconditioner%tree%graph%heads_tails,&
      arcs_weights=arrays%arcs_conductances(-tree%graph%n_special_arcs:tree%graph%n_arcs),&
      arcs_mask=support_tree%arrays%crossing_arcs_mask(-tree%graph%n_special_arcs:tree%graph%n_arcs),&
      orientations=tree%orientations,level_ordering=tree%level_ordering,&
      parents=tree%parents,cardinalities=tree%cardinalities,&
      path_labels=tree%path_labels,tree_mask=tree%mask,&
      from_scratch=.FALSE.,defragment_ordering=.TRUE.,&
      tree_type=max_cost_tree)
      CALL StopTimer(dual_system%timers%precond_reoptimization_timer)
      
      ELSE
      
      CALL StartTimer(dual_system%timers%precond_reoptimization_timer)
      
      CALL BuildNetworkMST(arc_offset=tree%graph%n_special_arcs,node_offset=tree%graph%n_special_nodes,&
      heads_tails=dual_system%preconditioner%tree%graph%heads_tails,&
      arcs_weights=arrays%arcs_conductances(-tree%graph%n_special_arcs:tree%graph%n_arcs),&
      orientations=tree%orientations,level_ordering=tree%level_ordering,&
      parents=tree%parents,cardinalities=tree%cardinalities,&
      path_labels=tree%path_labels,tree_mask=tree%mask,&
      from_scratch=.FALSE.,tree_type=max_cost_tree)
      
      CALL StopTimer(dual_system%timers%precond_reoptimization_timer)
      
      END IF
      
      CALL StartTimer(dual_system%timers%precond_factorization_timer)
      IF(dual_system%preconditioner%method==MST_QR_precond)THEN
      
      CALL FactorizeIncompleteQR(node_offset=tree%graph%n_special_nodes,&
      arc_offset=tree%graph%n_special_arcs,&
      heads_tails=tree%graph%heads_tails,orientations=tree%orientations,&
      parents=tree%parents,level_ordering=tree%level_ordering,&
      tree_mask=tree%mask,arcs_conductances=arrays%arcs_conductances,&
      nodes_resistances=dual_system%preconditioner%nodes_resistances,&
      nodes_multipliers=dual_system%preconditioner%nodes_multipliers)
      ELSE IF((dual_system%preconditioner%method==MST_LDLt_precond).OR.(dual_system%preconditioner%method==ST_MST_precond))THEN
      
      CALL FactorizeIncompleteLDLt(node_offset=tree%graph%n_special_nodes,&
      arc_offset=tree%graph%n_special_arcs,&
      heads_tails=tree%graph%heads_tails,orientations=tree%orientations,&
      parents=tree%parents,level_ordering=tree%level_ordering,&
      tree_mask=tree%mask,arcs_conductances=arrays%arcs_conductances,&
      nodes_resistances=dual_system%preconditioner%nodes_resistances,&
      nodes_multipliers=dual_system%preconditioner%nodes_multipliers)
      END IF
      CALL StopTimer(dual_system%timers%precond_factorization_timer)
      END IF UpdateMST
      
      CALL StartTimer(dual_system%timers%precond_factorization_timer)
      UpdateTAUCS:IF((dual_system%preconditioner%method==TAUCS_LLt_precond).OR.&
      (dual_system%preconditioner%method==TAUCS_MWB_precond))THEN
      
      IF(dual_system%preconditioner%taucs_preconditioner%droptol<=0.0_r_wp)THEN
      dual_system%preconditioner%taucs_preconditioner%factorization=TAUCS_complete_fact_mf
      ELSE
      dual_system%preconditioner%taucs_preconditioner%factorization=TAUCS_incomplete_fact
      END IF
      CALL TAUCS_CreatePreconditioner(taucs_preconditioner=dual_system%preconditioner%taucs_preconditioner)
      END IF UpdateTAUCS
      CALL StopTimer(dual_system%timers%precond_factorization_timer)
      
      
      CALL StartTimer(dual_system%timers%pcg_timer)
      SELECT CASE(dual_system%preconditioner%method)
      CASE(diagonal_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualDiagonalPrecond)
      CASE(MST_BHBt_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualMSTPrecond)
      CASE(MST_QR_precond,MST_LDLt_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualFactorizationPrecond)
      CASE(ST_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualSupportTreePrecond)
      CASE(TAUCS_LLt_precond,TAUCS_MWB_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualTAUCSPrecond)
      CASE(ST_MST_precond)
      
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply,&
      Preconditioner=DualSupportedMSTPrecond)
      CASE DEFAULT
      CALL PreconditionedConjugateGradient(system=system,&
      MatrixVectorMultiplication=DualMatrixMultiply)
      ENDSELECT
      CALL StopTimer(dual_system%timers%pcg_timer)
      
      IF(solver%convergence==CG_converged)THEN
      WRITE(message_log_unit,*)"--> CG converged after ",solver%n_iterations,&
      " iterations with errors ",solver%errors
      ELSE
      WRITE(message_log_unit,*)"--> CG did NOT converge.",&
      " Last recorded errors are ",solver%errors
      END IF
      
      IF((dual_system%preconditioner%method==TAUCS_LLt_precond).OR.&
      (dual_system%preconditioner%method==TAUCS_MWB_precond))THEN
      CALL TAUCS_FreePreconditioner(taucs_preconditioner=dual_system%preconditioner%taucs_preconditioner)
      END IF
      
      NULLIFY(dual_systems(dual_system%PIN)%system)
      
      END SUBROUTINE SolveDualSystem_PCG
      
      
      SUBROUTINE DualMatrixMultiply(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      TYPE(Dual_Network_System),POINTER::dual_system
      
      dual_system=>dual_systems(system%PIN)%system
      
      CALL Multiply_ADAt(heads_tails=dual_system%graph%heads_tails,&
      node_offset=dual_system%graph%n_special_nodes,&
      arcs_conductances=dual_system%arrays%arcs_conductances,&
      nodes_flows=dual_system%system%solver%cg_Vy,&
      nodes_potentials=dual_system%system%solver%cg_Vx)
      
      
      IF(dual_system%ground_dummy)&
      dual_system%system%solver%cg_Vy(0)=dual_system%system%solver%cg_Vy(0)+dual_system%grounding_conductance*dual_system%system%s&
     &olver%cg_Vx(0)
      
      END SUBROUTINE DualMatrixMultiply
      
      
      SUBROUTINE DualDiagonalPrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      
      dual_system=>dual_systems(system%PIN)%system
      
      CALL VectorMultiplication(first=dual_system%preconditioner%nodes_resistances,&
      second=dual_system%system%solver%cg_residuals,&
      product=dual_system%system%solver%cg_Vz)
      
      
      
      IF(dual_system%ground_dummy)THEN
      CALL VectorShift(vector=dual_system%system%solver%cg_Vz,&
      shift=-dual_system%system%solver%cg_Vz(0))
      END IF
      
      
      END SUBROUTINE DualDiagonalPrecond
      
      
      SUBROUTINE DualMSTPrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(CG_Solver),POINTER::solver
      TYPE(Network_Arrays),POINTER::arrays
      INTEGER(KIND=i_wp)::lb,ub
      
      dual_system=>dual_systems(system%PIN)%system
      tree=>dual_system%preconditioner%tree
      solver=>dual_system%system%solver
      arrays=>dual_system%arrays
      
      
      lb=-tree%graph%n_special_nodes
      ub=tree%graph%n_nodes
      solver%cg_Vz(tree%level_ordering(0))=0.0_r_wp
      CALL BasicDualNewtonSystem(node_offset=tree%graph%n_special_nodes,&
      arc_offset=tree%graph%n_special_arcs,heads_tails=tree%graph%heads_tails,&
      orientations=tree%orientations,tree_mask=tree%mask,&
      parents=tree%parents,level_ordering=tree%level_ordering,&
      tree_arcs_excess_voltages=dual_system%preconditioner%tree_arcs_excess_voltages,&
      nodes_excess_potentials=solver%cg_Vz(lb:ub),&
      nodes_excess_flows=solver%cg_residuals(lb:ub),&
      arcs_conductances=arrays%arcs_conductances(-tree%graph%n_special_arcs:tree%graph%n_arcs))
      
      IF(dual_system%ground_dummy)THEN
      CALL VectorShift(vector=solver%cg_Vz,&
      shift=-solver%cg_Vz(0))
      END IF
      
      
      END SUBROUTINE DualMSTPrecond
      
      
      SUBROUTINE DualFactorizationPrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(IN)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(CG_Solver),POINTER::solver
      TYPE(Network_Arrays),POINTER::arrays
      INTEGER(KIND=i_wp)::lb,ub
      
      dual_system=>dual_systems(system%PIN)%system
      tree=>dual_system%preconditioner%tree
      solver=>dual_system%system%solver
      arrays=>dual_system%arrays
      
      
      lb=-tree%graph%n_special_nodes
      ub=tree%graph%n_nodes
      CALL BasicDualNewtonSystem(node_offset=tree%graph%n_special_nodes,&
      arc_offset=tree%graph%n_special_arcs,tree_mask=tree%mask,&
      heads_tails=tree%graph%heads_tails,orientations=tree%orientations,&
      parents=tree%parents,level_ordering=tree%level_ordering,&
      tree_arcs_excess_voltages=dual_system%preconditioner%tree_arcs_excess_voltages,&
      nodes_excess_potentials=solver%cg_Vz(lb:ub),&
      nodes_excess_flows=solver%cg_residuals(lb:ub),&
      nodes_resistances=dual_system%preconditioner%nodes_resistances(lb:ub),&
      nodes_multipliers=dual_system%preconditioner%nodes_multipliers(lb:ub))
      
      IF(dual_system%ground_dummy)THEN
      CALL VectorShift(vector=solver%cg_Vz,&
      shift=-solver%cg_Vz(0))
      END IF
      
      
      END SUBROUTINE DualFactorizationPrecond
      
      
      
      SUBROUTINE DualSupportTreePrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree
      TYPE(CG_Solver),POINTER::solver
      
      dual_system=>dual_systems(system%PIN)%system
      solver=>dual_system%system%solver
      support_tree=>dual_system%preconditioner%support_tree
      
      CALL ST_DualNewtonPreconditioner(height=support_tree%height,degree=support_tree%degree,&
      node_offset=dual_system%graph%n_special_nodes,nodes_mapping=support_tree%arrays%nodes_mapping,&
      nodes_excess_flows=solver%cg_residuals,&
      nodes_excess_potentials=solver%cg_Vz,&
      nodes_resistances=dual_system%preconditioner%nodes_resistances,&
      ST_potentials_and_flows=support_tree%arrays%ST_potentials_and_flows,&
      ST_arcs_resistances=support_tree%arrays%ST_arcs_resistances)
      
      IF(dual_system%ground_dummy)THEN
      CALL VectorShift(vector=solver%cg_Vz,&
      shift=-solver%cg_Vz(0))
      END IF
      
      
      END SUBROUTINE DualSupportTreePrecond
      
      
      SUBROUTINE DualSupportedMSTPrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(Network_Arrays),POINTER::arrays
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Spanning_Tree),POINTER::tree
      TYPE(Support_Tree_Preconditioner),POINTER::support_tree
      TYPE(CG_Solver),POINTER::solver
      INTEGER::node,arc
      
      dual_system=>dual_systems(system%PIN)%system
      arrays=>dual_system%arrays
      solver=>dual_system%system%solver
      graph=>dual_system%graph
      support_tree=>dual_system%preconditioner%support_tree
      tree=>dual_system%preconditioner%tree
      
      CALL ST_DualNewtonPreconditioner(height=support_tree%height,degree=support_tree%degree,&
      node_offset=graph%n_special_nodes,nodes_mapping=support_tree%arrays%nodes_mapping,&
      nodes_excess_flows=solver%cg_residuals,&
      ST_potentials_and_flows=support_tree%arrays%ST_potentials_and_flows,&
      ST_arcs_resistances=support_tree%arrays%ST_arcs_resistances)
      
      CALL BasicDualNewtonSystem(node_offset=tree%graph%n_special_nodes,&
      arc_offset=tree%graph%n_special_arcs,&
      heads_tails=tree%graph%heads_tails,orientations=tree%orientations,&
      parents=tree%parents,level_ordering=tree%level_ordering,&
      tree_mask=tree%mask,&
      nodes_excess_potentials=solver%cg_Vz,&
      nodes_excess_flows=solver%cg_residuals,&
      tree_arcs_excess_voltages=dual_system%preconditioner%tree_arcs_excess_voltages,&
      nodes_resistances=dual_system%preconditioner%nodes_resistances,&
      nodes_multipliers=dual_system%preconditioner%nodes_multipliers)
      
      solver%cg_Vz=solver%cg_Vz+support_tree%arrays%ST_potentials_and_flows(support_tree%arrays%nodes_mapping)
      
      IF(dual_system%ground_dummy)THEN
      CALL VectorShift(vector=solver%cg_Vz,&
      shift=-solver%cg_Vz(0))
      END IF
      
      
      END SUBROUTINE DualSupportedMSTPrecond
      
      
      SUBROUTINE DualTAUCSPrecond(system)
      IMPLICIT NONE
      TYPE(CG_System),INTENT(INOUT)::system
      
      
      TYPE(Dual_Network_System),POINTER::dual_system
      TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_preconditioner
      
      dual_system=>dual_systems(system%PIN)%system
      taucs_preconditioner=>dual_system%preconditioner%taucs_preconditioner
      
      CALL TAUCS_ApplyPreconditioner(taucs_preconditioner=taucs_preconditioner)
      
      IF(dual_system%ground_dummy)THEN
      
      END IF
      
      
      END SUBROUTINE DualTAUCSPrecond
      
      
      
      END MODULE Dual_Network_Solvers