MODULE Support_Trees USE Precision USE Error_Handling USE System_Monitors USE Network_Matrix_Operations USE Graph_Algorithms IMPLICIT NONE PUBLIC::ST_PropagateNodesPotentials,ST_PropagateArcsFlows,& ST_PropagateArcsConductances,ST_BasicDualNewtonSystem,& ST_DualNewtonPreconditioner,ST_ArcsMask,ST_MatchingNodesMapping PRIVATE CONTAINS SUBROUTINE ST_PropagateNodesPotentials(height,degree,& arcs_voltages,nodes_potentials,potentials_and_voltages) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::height,degree REAL(KIND=r_wp),DIMENSION(0:),INTENT(INOUT),OPTIONAL::nodes_potentials,potentials_and_voltages REAL(KIND=r_wp),DIMENSION(0:),INTENT(IN),OPTIONAL::arcs_voltages INTEGER(KIND=i_wp)::n_nodes,n_tree_leafs INTEGER(KIND=i_wp)::node,parent_node n_tree_leafs=degree**height n_nodes=((n_tree_leafs-1)*degree)/(degree-1) IF(PRESENT(arcs_voltages).AND.PRESENT(nodes_potentials))THEN node=1 DO parent_node=0,n_nodes-n_tree_leafs nodes_potentials(node:(node+degree-1))=nodes_potentials(parent_node)+& arcs_voltages(node:(node+degree-1)) node=node+degree END DO ELSE IF(PRESENT(potentials_and_voltages))THEN node=1 DO parent_node=0,n_nodes-n_tree_leafs potentials_and_voltages(node:(node+degree-1))=potentials_and_voltages(parent_node)+& potentials_and_voltages(node:(node+degree-1)) node=node+degree END DO END IF END SUBROUTINE ST_PropagateNodesPotentials SUBROUTINE ST_PropagateArcsFlows(height,degree,& supplies_demands,arcs_flows,supplies_and_flows) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::height,degree REAL(KIND=r_wp),DIMENSION(0:),INTENT(IN),OPTIONAL::supplies_demands REAL(KIND=r_wp),DIMENSION(0:),INTENT(OUT),OPTIONAL::arcs_flows REAL(KIND=r_wp),DIMENSION(0:),INTENT(INOUT),OPTIONAL::supplies_and_flows INTEGER(KIND=i_wp)::n_nodes,n_tree_leafs INTEGER(KIND=i_wp)::node,parent_node,child_node n_tree_leafs=degree**height n_nodes=((n_tree_leafs-1)*degree)/(degree-1) IF(PRESENT(arcs_flows).AND.PRESENT(supplies_demands))THEN arcs_flows=supplies_demands node=n_nodes PropagateFlows:DO parent_node=n_nodes-n_tree_leafs,0,-1 DO child_node=node,node-degree+1,-1 arcs_flows(parent_node)=arcs_flows(parent_node)+arcs_flows(child_node) END DO node=node-degree END DO PropagateFlows ELSE IF(PRESENT(supplies_and_flows))THEN node=n_nodes DO parent_node=n_nodes-n_tree_leafs,0,-1 DO child_node=node,node-degree+1,-1 supplies_and_flows(parent_node)=supplies_and_flows(parent_node)+supplies_and_flows(child_node) END DO node=node-degree END DO END IF END SUBROUTINE ST_PropagateArcsFlows SUBROUTINE ST_PropagateArcsConductances(height,degree,arc_offset,node_offset,& heads_tails,nodes_mapping,arcs_conductances,& nodes_conductances,nodes_resistances,& ST_arcs_conductances,ST_arcs_resistances,ST_nodes_temp) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::height,degree INTEGER(KIND=i_wp),INTENT(IN)::node_offset,arc_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_mapping REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::& nodes_conductances,nodes_resistances REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_conductances REAL(KIND=r_wp),DIMENSION(0:),INTENT(OUT),OPTIONAL::& ST_arcs_conductances,ST_arcs_resistances REAL(KIND=r_wp),DIMENSION(0:),INTENT(INOUT)::ST_nodes_temp INTEGER(KIND=i_wp)::n_tree_nodes,n_tree_supernodes,n_tree_leafs INTEGER(KIND=i_wp)::tree_node,tree_supernode INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs INTEGER(KIND=i_wp)::node,arc,head,tail,joint,first_climber,second_climber n_special_nodes=-INT(LBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp) n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp) n_tree_leafs=degree**height n_tree_nodes=((n_tree_leafs-1)*degree)/(degree-1) IF(PRESENT(nodes_conductances))nodes_conductances=0.0_r_wp IF(PRESENT(nodes_resistances))nodes_resistances=0.0_r_wp IF(PRESENT(ST_arcs_conductances).AND.PRESENT(nodes_conductances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_CC:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) nodes_conductances(head)=nodes_conductances(head)+arcs_conductances(arc) nodes_conductances(tail)=nodes_conductances(tail)+arcs_conductances(arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_CC:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_CC END IF END DO FindJoint_CC ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_CC CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_conductances) ELSE IF(PRESENT(ST_arcs_resistances).AND.PRESENT(nodes_conductances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_RC:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) nodes_conductances(head)=nodes_conductances(head)+arcs_conductances(arc) nodes_conductances(tail)=nodes_conductances(tail)+arcs_conductances(arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_RC:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_RC END IF END DO FindJoint_RC ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_RC CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_resistances) ST_arcs_resistances=1.0_r_wp/(ST_arcs_resistances+EPSILON(1.0_r_wp)) ELSE IF(PRESENT(ST_arcs_conductances).AND.PRESENT(nodes_resistances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_CR:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) nodes_resistances(head)=nodes_resistances(head)+arcs_conductances(arc) nodes_resistances(tail)=nodes_resistances(tail)+arcs_conductances(arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_CR:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_CR END IF END DO FindJoint_CR ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_CR CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_conductances) nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) ELSE IF(PRESENT(ST_arcs_resistances).AND.PRESENT(nodes_resistances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_RR:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) nodes_resistances(head)=nodes_resistances(head)+arcs_conductances(arc) nodes_resistances(tail)=nodes_resistances(tail)+arcs_conductances(arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_RR:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_RR END IF END DO FindJoint_RR ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_RR CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_resistances) ST_arcs_resistances=1.0_r_wp/(ST_arcs_resistances+EPSILON(1.0_r_wp)) nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) ELSE IF(PRESENT(ST_arcs_conductances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_C:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_C:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_C END IF END DO FindJoint_C ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_C CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_conductances) ELSE IF(PRESENT(ST_arcs_resistances))THEN ST_nodes_temp=0.0_r_wp BoundaryArcs_R:DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) head=nodes_mapping(head) tail=nodes_mapping(tail) IF(head/=tail)THEN first_climber=head second_climber=tail FindJoint_R:DO IF(first_climber>second_climber)THEN first_climber=(first_climber-1)/degree ELSE IF(second_climber>first_climber)THEN second_climber=(second_climber-1)/degree ELSE joint=first_climber EXIT FindJoint_R END IF END DO FindJoint_R ST_nodes_temp(head)=ST_nodes_temp(head)+arcs_conductances(arc) ST_nodes_temp(tail)=ST_nodes_temp(tail)+arcs_conductances(arc) ST_nodes_temp(joint)=ST_nodes_temp(joint)-arcs_conductances(arc)-arcs_conductances(arc) END IF END DO BoundaryArcs_R CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=ST_nodes_temp,arcs_flows=ST_arcs_resistances) ST_arcs_resistances=1.0_r_wp/(ST_arcs_resistances+EPSILON(1.0_r_wp)) END IF END IF END SUBROUTINE ST_PropagateArcsConductances SUBROUTINE ST_BasicDualNewtonSystem(height,degree,& nodes_excess_potentials,nodes_excess_flows,excess_potentials_and_flows,& arcs_conductances,arcs_resistances) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::height,degree REAL(KIND=r_wp),DIMENSION(0:),INTENT(INOUT),OPTIONAL::& nodes_excess_potentials,excess_potentials_and_flows REAL(KIND=r_wp),DIMENSION(0:),INTENT(IN),OPTIONAL::nodes_excess_flows REAL(KIND=r_wp),DIMENSION(0:),INTENT(IN),OPTIONAL::& arcs_conductances,arcs_resistances IF(PRESENT(nodes_excess_potentials).AND.PRESENT(nodes_excess_flows))THEN CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_demands=nodes_excess_flows,arcs_flows=nodes_excess_potentials) IF(PRESENT(arcs_conductances))THEN nodes_excess_potentials=nodes_excess_potentials/arcs_conductances END IF IF(PRESENT(arcs_resistances))THEN nodes_excess_potentials=nodes_excess_potentials*arcs_resistances END IF nodes_excess_potentials(0)=0.0_r_wp CALL ST_PropagateNodesPotentials(height=height,degree=degree,& potentials_and_voltages=nodes_excess_potentials) ELSE IF(PRESENT(excess_potentials_and_flows))THEN CALL ST_PropagateArcsFlows(height=height,degree=degree,& supplies_and_flows=excess_potentials_and_flows) IF(PRESENT(arcs_conductances))THEN excess_potentials_and_flows=excess_potentials_and_flows/arcs_conductances END IF IF(PRESENT(arcs_resistances))THEN excess_potentials_and_flows=excess_potentials_and_flows*arcs_resistances END IF excess_potentials_and_flows(0)=0.0_r_wp CALL ST_PropagateNodesPotentials(height=height,degree=degree,& potentials_and_voltages=excess_potentials_and_flows) END IF END SUBROUTINE ST_BasicDualNewtonSystem SUBROUTINE ST_DualNewtonPreconditioner(height,degree,& node_offset,nodes_mapping,& nodes_excess_flows,nodes_excess_potentials,nodes_resistances,& ST_excess_flows,ST_excess_potentials,& ST_potentials_and_flows,ST_arcs_resistances) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::height,degree INTEGER(KIND=i_wp),INTENT(IN)::node_offset INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_mapping REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_excess_flows REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::nodes_resistances REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::nodes_excess_potentials REAL(KIND=r_wp),DIMENSION(0:),INTENT(INOUT),OPTIONAL::ST_potentials_and_flows REAL(KIND=r_wp),DIMENSION(0:),INTENT(OUT),OPTIONAL::ST_excess_flows,ST_excess_potentials REAL(KIND=r_wp),DIMENSION(0:),INTENT(IN)::ST_arcs_resistances INTEGER(KIND=i_wp)::n_nodes,n_special_nodes INTEGER(KIND=i_wp)::node,special_node n_special_nodes=-INT(LBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(nodes_mapping,DIM=1),KIND=i_wp) TemporaryChoice:IF(PRESENT(ST_excess_potentials).AND.PRESENT(ST_excess_flows))THEN DO node=-n_special_nodes,n_nodes ST_excess_flows(nodes_mapping(node))=ST_excess_flows(nodes_mapping(node))+& nodes_excess_flows(node) END DO CALL ST_BasicDualNewtonSystem(height=height,degree=degree,& nodes_excess_potentials=ST_excess_potentials,nodes_excess_flows=ST_excess_flows,& arcs_resistances=ST_arcs_resistances) DO node=-n_special_nodes,n_nodes ST_excess_flows(nodes_mapping(node))=ST_excess_flows(nodes_mapping(node))-& nodes_excess_flows(node) END DO NodesPotentials:IF(PRESENT(nodes_excess_potentials))THEN IF(PRESENT(nodes_resistances))THEN DO node=-n_special_nodes,n_nodes nodes_excess_potentials(node)=ST_excess_potentials(nodes_mapping(node))+& nodes_excess_flows(node)*nodes_resistances(node) END DO ELSE DO node=-n_special_nodes,n_nodes nodes_excess_potentials(node)=ST_excess_potentials(nodes_mapping(node)) END DO END IF END IF NodesPotentials ELSE IF(PRESENT(ST_potentials_and_flows))THEN ST_potentials_and_flows=0.0_r_wp DO node=-n_special_nodes,n_nodes ST_potentials_and_flows(nodes_mapping(node))=ST_potentials_and_flows(nodes_mapping(node))+& nodes_excess_flows(node) END DO CALL ST_BasicDualNewtonSystem(height=height,degree=degree,& excess_potentials_and_flows=ST_potentials_and_flows,& arcs_resistances=ST_arcs_resistances) NodesPotentials_PF:IF(PRESENT(nodes_excess_potentials))THEN IF(PRESENT(nodes_resistances))THEN DO node=-n_special_nodes,n_nodes nodes_excess_potentials(node)=ST_potentials_and_flows(nodes_mapping(node))+& nodes_excess_flows(node)*nodes_resistances(node) END DO ELSE DO node=-n_special_nodes,n_nodes nodes_excess_potentials(node)=ST_potentials_and_flows(nodes_mapping(node)) END DO END IF END IF NodesPotentials_PF END IF TemporaryChoice END SUBROUTINE ST_DualNewtonPreconditioner SUBROUTINE ST_ArcsMask(node_offset,arc_offset,heads_tails,& arcs_mask,nodes_mapping,mask_internal_arcs) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::node_offset,arc_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(OUT)::arcs_mask INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_mapping LOGICAL,INTENT(IN),OPTIONAL::mask_internal_arcs INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs INTEGER(KIND=i_wp)::node,arc,head,tail LOGICAL::internal_arcs n_special_nodes=-INT(LBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp) n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp) internal_arcs=.TRUE. IF(PRESENT(mask_internal_arcs))internal_arcs=mask_internal_arcs DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) IF((nodes_mapping(head)==nodes_mapping(tail)).AND.internal_arcs)THEN arcs_mask(arc)=.TRUE._l_wp ELSE arcs_mask(arc)=.FALSE._l_wp END IF END DO END SUBROUTINE ST_ArcsMask SUBROUTINE ST_MatchingNodesMapping(node_offset,arc_offset,heads_tails,& arcs_conductances,nodes_mapping,n_random_matchings,& n_components,mapping_offset,heavy_matching) USE Sorting_Ranking IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::node_offset,arc_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::arcs_conductances INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(OUT)::nodes_mapping INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::n_components,mapping_offset INTEGER,INTENT(IN),OPTIONAL::n_random_matchings LOGICAL,INTENT(IN),OPTIONAL::heavy_matching REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::arcs_weights INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::supernodes,heights,& arcs_sequence,arcs_reverse_sequence,end_nodes,first_arcs,next_arcs INTEGER(KIND=i_wp),DIMENSION(:,:),ALLOCATABLE::heads_tails_ INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE::nodes_mask INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs INTEGER(KIND=i_wp)::node,arc,head,tail,next_node,search_node,n_alone,& index,previous_arc,next_arc,first_arc,previous_search_arc,& n_unmatched_nodes,n_unmatched_arcs,counter,first_reverse_arc,& search_arc,adjacent_node,next_search_arc,n_supernodes,& head_supernode,tail_supernode,climber,climber_parent REAL::compression REAL(KIND=r_wp)::max_weight INTEGER::alloc_status,n_contractions,n_random_contractions,debug_unit INTEGER,PARAMETER::n_max_idle_contractions=10 LOGICAL::order_arcs,preorder_arcs,delete_arc,insert_arc,odd_cycle order_arcs=.FALSE. IF(PRESENT(heavy_matching))order_arcs=heavy_matching debug_unit=10 IF(debug_unit/=6)THEN OPEN(UNIT=debug_unit,FILE="Matching.log",ACTION="WRITE",STATUS="UNKNOWN") END IF n_special_nodes=-INT(LBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(nodes_mapping,DIM=1),KIND=i_wp) n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp) n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp) IF(.NOT.ALLOCATED(supernodes))THEN ALLOCATE(supernodes(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(heights))THEN ALLOCATE(heights(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(end_nodes))THEN ALLOCATE(end_nodes(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(nodes_mask))THEN ALLOCATE(nodes_mask(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_byte,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_sequence))THEN ALLOCATE(arcs_sequence(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF ALLOCATE(heads_tails_(1:2,-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=INT(SIZE(heads_tails_),KIND=i_wp),mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) heads_tails_=heads_tails n_random_contractions=0 IF(PRESENT(n_random_matchings))n_random_contractions=n_random_matchings preorder_arcs=.FALSE. IF(order_arcs)THEN IF(.NOT.ALLOCATED(first_arcs))THEN ALLOCATE(first_arcs(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(next_arcs))THEN ALLOCATE(next_arcs(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_reverse_sequence))THEN ALLOCATE(arcs_reverse_sequence(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_weights))THEN ALLOCATE(arcs_weights(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_i_wp,& caller="ST_MatchingNodesMapping",alloc_status=alloc_status) END IF IF(PRESENT(arcs_conductances))THEN preorder_arcs=.TRUE. arcs_weights=arcs_conductances ELSE arcs_weights=1.0_r_wp END IF END IF PreorderArcs:IF(preorder_arcs)THEN IF((n_special_arcs>=0).AND.(n_arcs>=0))arcs_weights(0)=-HUGE(1.0_r_wp) CALL QuickRank(array=arcs_weights,permutation=arcs_reverse_sequence,& partially_ranked=.FALSE.,pivot_selection='R') first_arc=arcs_reverse_sequence(n_arcs) previous_arc=first_arc DO index=n_arcs,-n_special_arcs,-1 arc=arcs_reverse_sequence(index)-n_special_arcs-1 arcs_sequence(previous_arc)=arc previous_arc=arc END DO arcs_sequence(arc)=0 ELSE first_arc=n_arcs DO arc=n_arcs,-n_special_arcs+1,-1 arcs_sequence(arc)=arc-1 END DO arcs_sequence(-n_special_arcs)=0 IF((n_special_arcs>0).AND.(n_arcs>0))arcs_sequence(1)=-1 END IF PreorderArcs DO node=-n_special_nodes,n_nodes supernodes(node)=node nodes_mapping(node)=node end_nodes(node)=node END DO heights=0 nodes_mask=matching_node n_unmatched_nodes=n_special_nodes+n_nodes+1 n_unmatched_arcs=n_special_arcs+n_arcs+1 n_contractions=0 odd_cycle=.TRUE. CoarsenGraph:DO WHILE(first_arc/=0) n_contractions=n_contractions+1 WRITE(debug_unit,*)"At contraction #",n_contractions,& " still unmatched arcs=",n_unmatched_arcs,"nodes=",n_unmatched_nodes CALL ContractMatchedArcs(node_offset=n_special_nodes,arc_offset=n_special_arcs,& heads_tails=heads_tails_,& arcs_sequence=arcs_sequence,first_arc=first_arc,& supernodes=supernodes,heights=heights,& nodes_ordering=nodes_mapping,end_nodes=end_nodes,& nodes_mask=nodes_mask,reuse_matching=.TRUE.,& update_heads_tails=.TRUE.,odd_cycle=odd_cycle) odd_cycle=.NOT.odd_cycle IF(order_arcs.AND.(n_contractions>n_random_contractions))THEN previous_arc=0 arc=first_arc n_unmatched_arcs=0 InitializeIndicators:DO IF(arc==0)EXIT InitializeIndicators n_unmatched_arcs=n_unmatched_arcs+1 head=heads_tails_(1,arc) tail=heads_tails_(2,arc) first_arcs(head)=0 first_arcs(tail)=0 arcs_reverse_sequence(arc)=previous_arc previous_arc=arc arc=arcs_sequence(arc) END DO InitializeIndicators first_reverse_arc=previous_arc arc=first_arc n_unmatched_nodes=0 CompressArcs:DO IF(arc==0)EXIT CompressArcs delete_arc=.FALSE. head=heads_tails_(1,arc) tail=heads_tails_(2,arc) IF(headtail)THEN node=tail adjacent_node=head ELSE delete_arc=.TRUE. END IF SearchList:IF(.NOT.delete_arc)THEN previous_search_arc=0 search_arc=first_arcs(node) IF(search_arc==0)n_unmatched_nodes=n_unmatched_nodes+1 CheckDuplicate:DO IF(search_arc==0)EXIT CheckDuplicate search_node=MAX(heads_tails_(1,search_arc),& heads_tails_(2,search_arc)) DuplicateArc:IF(search_node==adjacent_node)THEN delete_arc=.TRUE. arcs_weights(search_arc)=arcs_weights(search_arc)+arcs_weights(arc) previous_search_arc=arcs_reverse_sequence(search_arc) next_search_arc=arcs_sequence(search_arc) IF(previous_search_arc==0.AND.next_search_arc==0)THEN first_arc=0 first_reverse_arc=0 ELSE IF(previous_search_arc==0)THEN first_arc=next_search_arc arcs_reverse_sequence(first_arc)=0 ELSE IF(next_search_arc==0)THEN first_reverse_arc=previous_search_arc arcs_sequence(first_reverse_arc)=0 ELSE arcs_sequence(previous_search_arc)=next_search_arc arcs_reverse_sequence(next_search_arc)=previous_search_arc END IF ReorderArc:DO IF(previous_search_arc==0)EXIT ReorderArc IF(arcs_weights(search_arc)<=arcs_weights(previous_search_arc))& EXIT ReorderArc previous_search_arc=arcs_reverse_sequence(previous_search_arc) END DO ReorderArc IF(previous_search_arc==0)THEN next_search_arc=first_arc arcs_reverse_sequence(next_search_arc)=search_arc first_arc=search_arc ELSE IF(previous_search_arc==first_reverse_arc)THEN next_search_arc=0 arcs_sequence(previous_search_arc)=search_arc first_reverse_arc=search_arc ELSE next_search_arc=arcs_sequence(previous_search_arc) arcs_sequence(previous_search_arc)=search_arc arcs_reverse_sequence(next_search_arc)=search_arc END IF arcs_sequence(search_arc)=next_search_arc arcs_reverse_sequence(search_arc)=previous_search_arc EXIT CheckDuplicate ELSE IF(search_node>adjacent_node)THEN EXIT CheckDuplicate END IF DuplicateArc previous_search_arc=search_arc search_arc=next_arcs(search_arc) END DO CheckDuplicate END IF SearchList DeleteArc:IF(delete_arc)THEN previous_arc=arcs_reverse_sequence(arc) next_arc=arcs_sequence(arc) IF(previous_arc==0.AND.next_arc==0)THEN first_arc=0 first_reverse_arc=0 ELSE IF(previous_arc==0)THEN first_arc=next_arc arcs_reverse_sequence(first_arc)=0 ELSE IF(next_arc==0)THEN first_reverse_arc=previous_arc arcs_sequence(first_reverse_arc)=0 ELSE arcs_sequence(previous_arc)=next_arc arcs_reverse_sequence(next_arc)=previous_arc END IF ELSE IF(previous_search_arc/=0)THEN next_arc=next_arcs(previous_search_arc) next_arcs(previous_search_arc)=arc ELSE next_arc=first_arcs(node) first_arcs(node)=arc END IF next_arcs(arc)=next_arc END IF DeleteArc arc=arcs_sequence(arc) END DO CompressArcs n_alone=0 arc=first_arc max_weight=arcs_weights(first_arc) UnmatchedNodes:DO IF(arc==0)EXIT UnmatchedNodes node=MIN(heads_tails_(1,arc),heads_tails_(2,arc)) IF(nodes_mask(node)>=& MIN(n_contractions,1+n_max_idle_contractions/n_contractions)+matching_node)THEN n_alone=n_alone+1 arcs_weights(arc)=arcs_weights(arc)+max_weight previous_arc=arcs_reverse_sequence(arc) next_arc=arcs_sequence(arc) IF(previous_arc==0.AND.next_arc==0)THEN first_arc=0 first_reverse_arc=0 ELSE IF(previous_arc==0)THEN first_arc=next_arc arcs_reverse_sequence(first_arc)=0 ELSE IF(next_arc==0)THEN first_reverse_arc=previous_arc arcs_sequence(first_reverse_arc)=0 ELSE arcs_sequence(previous_arc)=next_arc arcs_reverse_sequence(next_arc)=previous_arc END IF search_arc=next_arc next_arc=first_arc DO IF(next_arc==0)EXIT IF(arcs_weights(arc)>=arcs_weights(next_arc))EXIT next_arc=arcs_sequence(next_arc) END DO IF(next_arc==first_arc)THEN previous_arc=0 arcs_reverse_sequence(first_arc)=arc first_arc=arc ELSE IF(next_arc==0)THEN previous_arc=first_reverse_arc arcs_sequence(first_reverse_arc)=arc first_reverse_arc=arc ELSE previous_arc=arcs_reverse_sequence(next_arc) arcs_reverse_sequence(next_arc)=arc arcs_sequence(previous_arc)=arc END IF arcs_sequence(arc)=next_arc arcs_reverse_sequence(arc)=previous_arc arc=search_arc CYCLE UnmatchedNodes END IF arc=arcs_sequence(arc) END DO UnmatchedNodes IF(n_alone>0)WRITE(debug_unit,*)n_alone,& " nodes have still not been matched at contraction #",n_contractions END IF END DO CoarsenGraph IF(PRESENT(mapping_offset))THEN counter=mapping_offset ELSE counter=0 END IF n_supernodes=0 DO node=-n_special_nodes,n_nodes IF(supernodes(node)==node)THEN n_supernodes=n_supernodes+1 search_node=node TraverseList:DO next_node=nodes_mapping(search_node) nodes_mapping(search_node)=counter counter=counter+1 IF(next_node==search_node)EXIT TraverseList search_node=next_node END DO TraverseList WRITE(debug_unit,*)"Supernode #",n_supernodes," [",node,"] ",& " childrens list stored up to index ",counter END IF END DO IF(PRESENT(n_components))THEN compression=REAL(n_components)/REAL(n_nodes+n_special_nodes) nodes_mapping=INT(compression*REAL(nodes_mapping),i_wp) END IF IF(ALLOCATED(supernodes))THEN CALL RecordAllocation(n_elements=-INT(SIZE(supernodes),KIND=i_wp),mold=1_i_wp) DEALLOCATE(supernodes,STAT=alloc_status) END IF IF(ALLOCATED(heights))THEN CALL RecordAllocation(n_elements=-INT(SIZE(heights),KIND=i_wp),mold=1_i_wp) DEALLOCATE(heights,STAT=alloc_status) END IF IF(ALLOCATED(end_nodes))THEN CALL RecordAllocation(n_elements=-INT(SIZE(end_nodes),KIND=i_wp),mold=1_i_wp) DEALLOCATE(end_nodes,STAT=alloc_status) END IF IF(ALLOCATED(nodes_mask))THEN CALL RecordAllocation(n_elements=-INT(SIZE(nodes_mask),KIND=i_wp),mold=1_i_byte) DEALLOCATE(nodes_mask,STAT=alloc_status) END IF IF(ALLOCATED(arcs_sequence))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_sequence),KIND=i_wp),mold=1_i_wp) DEALLOCATE(arcs_sequence,STAT=alloc_status) END IF IF(ALLOCATED(arcs_weights))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_weights),KIND=i_wp),mold=1_i_wp) DEALLOCATE(arcs_weights,STAT=alloc_status) END IF CALL RecordAllocation(n_elements=-INT(SIZE(heads_tails_),KIND=i_wp),mold=1_i_wp,& caller="ST_MatchingNodesMapping") DEALLOCATE(heads_tails_,STAT=alloc_status) IF(order_arcs)THEN IF(ALLOCATED(first_arcs))THEN CALL RecordAllocation(n_elements=-INT(SIZE(first_arcs),KIND=i_wp),mold=1_i_wp) DEALLOCATE(first_arcs,STAT=alloc_status) END IF IF(ALLOCATED(next_arcs))THEN CALL RecordAllocation(n_elements=-INT(SIZE(next_arcs),KIND=i_wp),mold=1_i_wp) DEALLOCATE(next_arcs,STAT=alloc_status) END IF IF(ALLOCATED(arcs_reverse_sequence))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_reverse_sequence),KIND=i_wp),mold=1_i_wp) DEALLOCATE(arcs_reverse_sequence,STAT=alloc_status) END IF END IF END SUBROUTINE ST_MatchingNodesMapping END MODULE Support_Trees