MODULE Graph_Algorithms USE Precision USE Error_Handling USE System_Monitors USE Sorting_Ranking USE Network_Matrix_Operations IMPLICIT NONE PUBLIC::ConnectedComponents,BackboneCycles,ContractMatchedArcs,& BuildSpanningTree,BuildTreeStructures,& PropagateNodesPotentials,PropagateArcsFlows,& CalculatePathLabels,ReBuildSpanningTree,& BuildNetworkMST,BasicDualNewtonSystem,EffectiveCycleResistances,& FactorizeIncompleteQR,FactorizeIncompleteLDLt,& CreateAdjacencyArrays,CalculateCongestionsDilations PRIVATE LOGICAL,PUBLIC,SAVE::debug_graph_algs=.FALSE. INTEGER(KIND=i_byte),PARAMETER,PUBLIC::no_parent=0,& head_is_parent=1,tail_is_parent=2 INTEGER(KIND=i_byte),PARAMETER,PUBLIC::any_tree=1,random_tree=2,& min_cost_tree=3,max_cost_tree=4 INTEGER(KIND=i_byte),PARAMETER,PUBLIC::matching_node=2,dead_node=0 CONTAINS SUBROUTINE ConnectedComponents(heads_tails,node_offset,labels,use_labels,& largest_component_label,largest_component_size) IMPLICIT NONE INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),INTENT(IN)::node_offset INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::labels LOGICAL,INTENT(IN),OPTIONAL::use_labels INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::largest_component_label,largest_component_size INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::heights INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs INTEGER(KIND=i_wp)::node,arc,head,tail,head_root,tail_root,& root,diff_heights,child,parent,& largest_component(1),climber,climber_parent INTEGER::alloc_status LOGICAL::reuse_labels n_special_nodes=-INT(LBOUND(labels,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(labels,DIM=1),KIND=i_wp) n_arcs=INT(SIZE(heads_tails,DIM=2),KIND=i_wp) reuse_labels=.FALSE. IF(PRESENT(use_labels))reuse_labels=use_labels IF(.NOT.reuse_labels)THEN DO node=-n_special_nodes,n_nodes labels(node)=node END DO END IF 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="ConnectedComponents",alloc_status=alloc_status) heights=0 ForAllArcs:DO arc=1,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) climber=head DO WHILE(labels(climber)/=climber) climber=labels(climber) END DO head_root=climber climber=tail DO WHILE(labels(climber)/=climber) climber=labels(climber) END DO tail_root=climber climber=head DO WHILE(labels(climber)/=head_root) climber_parent=labels(climber) labels(climber)=head_root climber=climber_parent END DO climber=tail DO WHILE(labels(climber)/=tail_root) climber_parent=labels(climber) labels(climber)=tail_root climber=climber_parent END DO TreeUnion:IF(head_root/=tail_root)THEN diff_heights=(heights(head_root)-heights(tail_root)) IF(diff_heights>0)THEN parent=head_root child=tail_root ELSE parent=tail_root child=head_root END IF labels(child)=parent IF(diff_heights==0)heights(parent)=heights(parent)+1 END IF TreeUnion END DO ForAllArcs ForAllNodes:DO node=-n_special_nodes,n_nodes climber=node DO WHILE(labels(climber)/=climber) climber=labels(climber) END DO root=climber climber=node DO WHILE(labels(climber)/=root) climber_parent=labels(climber) labels(climber)=root climber=climber_parent END DO END DO ForAllNodes heights=0 DO node=-n_special_nodes,n_nodes heights(labels(node))=heights(labels(node))+1 END DO IF(PRESENT(largest_component_label).OR.PRESENT(largest_component_size))THEN largest_component=INT(MAXLOC(heights),KIND=i_wp) largest_component=largest_component-1-node_offset IF(PRESENT(largest_component_label))& largest_component_label=labels(largest_component(1)) IF(PRESENT(largest_component_size))& largest_component_size=heights(largest_component(1)) END IF CALL RecordAllocation(n_elements=-INT(SIZE(heights),KIND=i_wp),mold=1_i_wp) DEALLOCATE(heights) END SUBROUTINE ConnectedComponents SUBROUTINE BackboneCycles(node_offset,heads_tails,supernodes,& use_supernodes,largest_supernode,largest_supernode_size) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::node_offset INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::supernodes LOGICAL,INTENT(IN),OPTIONAL::use_supernodes INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::largest_supernode,largest_supernode_size INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::cycles,heights INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs INTEGER(KIND=i_wp)::node,arc,head,tail,counter,max_height,& supernode,last_supernode,previous_supernode,next_supernode,& child,parent,head_length,tail_length,head_supernode,tail_supernode,& largest_location(1),climber,climber_parent,supersupernode INTEGER::alloc_status LOGICAL::reuse_supernodes,cycle_closed,increase_height n_special_nodes=-INT(LBOUND(supernodes,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(supernodes,DIM=1),KIND=i_wp) n_arcs=INT(SIZE(heads_tails,DIM=2),KIND=i_wp) reuse_supernodes=.FALSE. IF(PRESENT(use_supernodes))reuse_supernodes=use_supernodes IF(.NOT.reuse_supernodes)THEN DO node=-n_special_nodes,n_nodes supernodes(node)=node END DO END IF ALLOCATE(cycles(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,& caller="ConnectedComponents",alloc_status=alloc_status) DO node=-n_special_nodes,n_nodes cycles(node)=node END DO 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="ConnectedComponents",alloc_status=alloc_status) heights=0 ForAllArcs:DO arc=1,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) head_length=0 max_height=-HUGE(1_i_wp) climber=head DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO supernode=climber climber=head DO WHILE(supernodes(climber)/=supernode) climber_parent=supernodes(climber) supernodes(climber)=supernode climber=climber_parent END DO increase_height=.FALSE. previous_supernode=supernode head_supernode=supernode Uncover_head:DO IF(heights(supernode)>max_height)THEN max_height=heights(supernode) & supersupernode=supernode ELSE IF(heights(supernode)==max_height)THEN increase_height=.TRUE. END IF head=cycles(supernode) cycles(supernode)=previous_supernode; previous_supernode=supernode climber=head DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO supernode=climber climber=head DO WHILE(supernodes(climber)/=supernode) climber_parent=supernodes(climber) supernodes(climber)=supernode climber=climber_parent END DO IF(supernode==previous_supernode)THEN last_supernode=supernode EXIT Uncover_head END IF head_length=head_length+1 END DO Uncover_head tail_length=0 max_height=-HUGE(1_i_wp) climber=tail DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO supernode=climber climber=tail DO WHILE(supernodes(climber)/=supernode) climber_parent=supernodes(climber) supernodes(climber)=supernode climber=climber_parent END DO increase_height=.FALSE. previous_supernode=supernode tail_supernode=supernode Uncover_tail:DO IF(heights(supernode)>max_height)THEN max_height=heights(supernode) & supersupernode=supernode ELSE IF(heights(supernode)==max_height)THEN increase_height=.TRUE. END IF tail=cycles(supernode) cycles(supernode)=previous_supernode; previous_supernode=supernode climber=tail DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO supernode=climber climber=tail DO WHILE(supernodes(climber)/=supernode) climber_parent=supernodes(climber) supernodes(climber)=supernode climber=climber_parent END DO IF(supernode==previous_supernode)THEN last_supernode=supernode EXIT Uncover_tail END IF tail_length=tail_length+1 END DO Uncover_tail AddArc:IF(last_supernode==head_supernode)THEN supernode=head_supernode CondenseCycle:DO supernodes(supernode)=supersupernode IF(supernode==tail_supernode)EXIT CondenseCycle supernode=cycles(supernode) END DO CondenseCycle IF(increase_height)heights(supersupernode)=heights(supersupernode)+1 ELSE IF(head_length>=tail_length)THEN cycles(head_supernode)=tail_supernode ELSE cycles(tail_supernode)=head_supernode END IF END IF AddArc END DO ForAllArcs CALL RecordAllocation(n_elements=-INT(SIZE(cycles),KIND=i_wp),mold=1_i_wp) DEALLOCATE(cycles) ForAllNodes:DO node=-n_special_nodes,n_nodes climber=node DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO supernode=climber climber=node DO WHILE(supernodes(climber)/=supernode) climber_parent=supernodes(climber) supernodes(climber)=supernode climber=climber_parent END DO END DO ForAllNodes heights=0 DO node=-n_special_nodes,n_nodes heights(supernodes(node))=heights(supernodes(node))+1 END DO IF(PRESENT(largest_supernode).OR.PRESENT(largest_supernode_size))THEN largest_location=INT(MAXLOC(heights),KIND=i_wp) largest_location=largest_location-1-node_offset IF(PRESENT(largest_supernode))& largest_supernode=supernodes(largest_location(1)) IF(PRESENT(largest_supernode_size))& largest_supernode_size=heights(largest_location(1)) END IF CALL RecordAllocation(n_elements=-INT(SIZE(heights),KIND=i_wp),mold=1_i_wp) DEALLOCATE(heights) END SUBROUTINE BackboneCycles SUBROUTINE ContractMatchedArcs(node_offset,arc_offset,heads_tails,& arcs_sequence,first_arc,& supernodes,heights,nodes_ordering,end_nodes,& nodes_mask,reuse_matching,update_heads_tails,odd_cycle) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::node_offset,arc_offset INTEGER(KIND=i_wp),INTENT(INOUT)::first_arc INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(INOUT)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-arc_offset:),INTENT(INOUT)::arcs_sequence INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::supernodes,heights INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT),OPTIONAL::& nodes_ordering,end_nodes INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(INOUT)::nodes_mask LOGICAL,INTENT(IN),OPTIONAL::reuse_matching,update_heads_tails LOGICAL,INTENT(IN),OPTIONAL::odd_cycle INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs INTEGER(KIND=i_wp)::node,arc,head,tail,head_root,tail_root,& root,diff_heights,child,parent,index,& climber,climber_parent,end_node,previous_arc INTEGER::alloc_status INTEGER(KIND=i_byte)::matched_node LOGICAL::reuse_supernodes,contract_heads_tails,delete_arc,sequence_nodes n_special_nodes=-INT(LBOUND(supernodes,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(supernodes,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) matched_node=1_i_byte IF(PRESENT(odd_cycle))THEN IF(odd_cycle)matched_node=-matched_node END IF sequence_nodes=(PRESENT(end_nodes).AND.PRESENT(nodes_ordering)) contract_heads_tails=.FALSE. IF(PRESENT(update_heads_tails))contract_heads_tails=update_heads_tails reuse_supernodes=.FALSE. IF(PRESENT(reuse_matching))reuse_supernodes=reuse_matching IF(.NOT.reuse_supernodes)THEN DO node=-n_special_nodes,n_nodes supernodes(node)=node END DO IF(sequence_nodes)THEN DO node=-n_special_nodes,n_nodes nodes_ordering(node)=node end_nodes(node)=node END DO END IF heights=0 END IF previous_arc=0 arc=first_arc ContractArcs:DO IF(arc==0)EXIT ContractArcs head=heads_tails(1,arc) tail=heads_tails(2,arc) delete_arc=.FALSE. climber=head DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO head_root=climber climber=tail DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO tail_root=climber climber=head DO WHILE(supernodes(climber)/=head_root) climber_parent=supernodes(climber) supernodes(climber)=head_root climber=climber_parent END DO climber=tail DO WHILE(supernodes(climber)/=tail_root) climber_parent=supernodes(climber) supernodes(climber)=tail_root climber=climber_parent END DO TreeUnion:IF((head_root==tail_root).OR.& ((nodes_mask(head_root)==dead_node).OR.(nodes_mask(tail_root)==dead_node)))THEN delete_arc=.TRUE. ELSE IF(nodes_mask(head_root)==-matched_node)THEN nodes_mask(head_root)=matching_node ELSE IF(nodes_mask(head_root)>=matching_node)THEN nodes_mask(head_root)=-(MIN(nodes_mask(head_root),HUGE(1_i_byte)-1_i_byte)+1_i_byte) END IF IF(nodes_mask(tail_root)==-matched_node)THEN nodes_mask(tail_root)=matching_node ELSE IF(nodes_mask(tail_root)>=matching_node)THEN nodes_mask(tail_root)=-(MIN(nodes_mask(tail_root),HUGE(1_i_byte)-1_i_byte)+1_i_byte) END IF ContractArc:IF((nodes_mask(head_root)/=matched_node).AND.& nodes_mask(tail_root)/=matched_node)THEN delete_arc=.TRUE. diff_heights=(heights(head_root)-heights(tail_root)) IF(diff_heights>0)THEN parent=head_root child=tail_root ELSE parent=tail_root child=head_root END IF nodes_mask(parent)=matched_node nodes_mask(child)=dead_node supernodes(child)=parent IF(diff_heights==0)heights(parent)=heights(parent)+1 IF(sequence_nodes)THEN end_node=end_nodes(parent) nodes_ordering(end_node)=child end_nodes(parent)=end_nodes(child) END IF END IF ContractArc END IF TreeUnion IF(delete_arc)THEN IF(arc==first_arc)THEN first_arc=arcs_sequence(arc) ELSE arcs_sequence(previous_arc)=arcs_sequence(arc) END IF ELSE previous_arc=arc END IF arc=arcs_sequence(arc) END DO ContractArcs arc=first_arc UpdateNodes:DO IF(arc==0)EXIT UpdateNodes head=heads_tails(1,arc) tail=heads_tails(2,arc) climber=head DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO head_root=climber climber=tail DO WHILE(supernodes(climber)/=climber) climber=supernodes(climber) END DO tail_root=climber climber=head DO WHILE(supernodes(climber)/=head_root) climber_parent=supernodes(climber) supernodes(climber)=head_root climber=climber_parent END DO climber=tail DO WHILE(supernodes(climber)/=tail_root) climber_parent=supernodes(climber) supernodes(climber)=tail_root climber=climber_parent END DO IF(contract_heads_tails)THEN heads_tails(1,arc)=head_root heads_tails(2,arc)=tail_root END IF IF(nodes_mask(head_root)<-matching_node)& nodes_mask(head_root)=ABS(nodes_mask(head_root)) IF(nodes_mask(tail_root)<-matching_node)& nodes_mask(tail_root)=ABS(nodes_mask(tail_root)) arc=arcs_sequence(arc) END DO UpdateNodes END SUBROUTINE ContractMatchedArcs SUBROUTINE BuildSpanningTree(arc_offset,node_offset,arcs_list,heads_tails,& orientations,parents,use_parenthood) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::arcs_list INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(INOUT)::orientations LOGICAL,INTENT(IN),OPTIONAL::use_parenthood INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::heights,forest,cardinalities INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,head,tail,head_root,tail_root,& root,diff_heights,child,parent,parent_node,& climber,climber_parent,diff_cardinalities INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation,climber_orientation LOGICAL::initialize_forest n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) initialize_forest=.TRUE. IF(PRESENT(use_parenthood))initialize_forest=.NOT.use_parenthood 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="BuildSpanningTree",alloc_status=alloc_status) ALLOCATE(forest(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,& caller="BuildSpanningTree",alloc_status=alloc_status) ALLOCATE(cardinalities(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,& caller="BuildSpanningTree",alloc_status=alloc_status) IF(initialize_forest)THEN DO node=-n_special_nodes,n_nodes forest(node)=node END DO orientations=no_parent parents=0 cardinalities=1 ELSE InitializeTrees:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation/=no_parent)THEN forest(node)=heads_tails(orientation,parents(node)) ELSE forest(node)=node END IF END DO InitializeTrees cardinalities=0 ShortenTrees:DO node=-n_special_nodes,n_nodes climber=node DO WHILE(forest(climber)/=climber) climber=forest(climber) END DO root=climber climber=node DO WHILE(forest(climber)/=root) climber_parent=forest(climber) forest(climber)=root climber=climber_parent END DO cardinalities(root)=cardinalities(root)+1 END DO ShortenTrees END IF heights=0 ForAllArcs:DO index=1,INT(SIZE(arcs_list),KIND=i_wp) arc=arcs_list(index) head=heads_tails(1,arc) tail=heads_tails(2,arc) climber=head DO WHILE(forest(climber)/=climber) climber=forest(climber) END DO head_root=climber climber=tail DO WHILE(forest(climber)/=climber) climber=forest(climber) END DO tail_root=climber climber=head DO WHILE(forest(climber)/=head_root) climber_parent=forest(climber) forest(climber)=head_root climber=climber_parent END DO climber=tail DO WHILE(forest(climber)/=tail_root) climber_parent=forest(climber) forest(climber)=tail_root climber=climber_parent END DO AddArc:IF(head_root/=tail_root)THEN diff_heights=(heights(head_root)-heights(tail_root)) IF(diff_heights>0)THEN parent=head_root child=tail_root ELSE parent=tail_root child=head_root END IF forest(child)=parent IF(diff_heights==0)heights(parent)=heights(parent)+1 cardinalities(parent)=cardinalities(parent)+cardinalities(child) diff_cardinalities=(cardinalities(head_root)-cardinalities(tail_root)) IF(diff_cardinalities>0)THEN parent=arc orientation=tail_is_parent climber=tail ELSE parent=arc orientation=head_is_parent climber=head END IF ReverseParenthood:DO climber_orientation=orientations(climber) climber_parent=parents(climber) parents(climber)=parent IF(orientation==head_is_parent)THEN orientations(climber)=tail_is_parent ELSE IF(orientation==tail_is_parent)THEN orientations(climber)=head_is_parent END IF IF(climber_orientation==no_parent)EXIT ReverseParenthood climber=heads_tails(climber_orientation,climber_parent) parent=climber_parent orientation=climber_orientation END DO ReverseParenthood END IF AddArc END DO ForAllArcs CALL RecordAllocation(n_elements=-INT(SIZE(forest),KIND=i_wp),mold=1_i_wp) DEALLOCATE(forest) CALL RecordAllocation(n_elements=-INT(SIZE(heights),KIND=i_wp),mold=1_i_wp) DEALLOCATE(heights) CALL RecordAllocation(n_elements=-INT(SIZE(cardinalities),KIND=i_wp),mold=1_i_wp) DEALLOCATE(cardinalities) END SUBROUTINE BuildSpanningTree SUBROUTINE BuildTreeStructures(arc_offset,node_offset,heads_tails,tree_mask,& orientations,parents,cardinalities,level_ordering,rebuild_cardinalities) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(OUT),OPTIONAL::tree_mask INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::cardinalities INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(OUT)::level_ordering LOGICAL,INTENT(IN),OPTIONAL::rebuild_cardinalities INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::counters INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,head,tail,parent,parent_arc,parent_node,& level_start,level_end,climbers(2),cycle_root INTEGER::side,alloc_status INTEGER(KIND=i_byte)::orientation LOGICAL::new_cardinalities n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) new_cardinalities=.TRUE. IF(PRESENT(rebuild_cardinalities))new_cardinalities=rebuild_cardinalities IF(new_cardinalities)THEN cardinalities=1 ChildDegrees:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation/=no_parent)THEN parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) cardinalities(parent_node)=cardinalities(parent_node)+1 END IF END DO ChildDegrees END IF ALLOCATE(counters(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_nodes+n_special_nodes+1,mold=1_i_wp,& caller="BuildSpanningTree",alloc_status=alloc_status) level_start=INT(UBOUND(level_ordering,DIM=1),KIND=i_wp)+1 level_end=level_start DO node=-n_special_nodes,n_nodes IF(cardinalities(node)==1)THEN level_end=level_end-1 level_ordering(level_end)=node END IF END DO counters=1 BuildLevels:DO WHILE(level_end>=1) level_start=level_start-1 node=level_ordering(level_start) orientation=orientations(node) IF(orientation==no_parent)CYCLE BuildLevels parent_node=heads_tails(orientation,parents(node)) IF(new_cardinalities)THEN counters(parent_node)=counters(parent_node)+1 ELSE counters(parent_node)=counters(parent_node)+counters(node) END IF IF(counters(parent_node)==cardinalities(parent_node))THEN level_end=level_end-1 level_ordering(level_end)=parent_node END IF END DO BuildLevels CALL RecordAllocation(n_elements=-INT(SIZE(counters),KIND=i_wp),mold=1_i_wp) DEALLOCATE(counters) IF(new_cardinalities)THEN cardinalities=1 BuildCardinalities:DO index=n_nodes+n_special_nodes,0,-1 node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)CYCLE BuildCardinalities parent_node=heads_tails(orientation,parents(node)) cardinalities(parent_node)=cardinalities(parent_node)+cardinalities(node) END DO BuildCardinalities END IF IF(PRESENT(tree_mask))THEN tree_mask=.FALSE._l_wp DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation/=no_parent)tree_mask(parents(node))=.TRUE._l_wp END DO END IF END SUBROUTINE BuildTreeStructures SUBROUTINE PropagateNodesPotentials(arc_offset,node_offset,heads_tails,& orientations,parents,level_ordering,& arcs_voltages,tree_arcs_voltages,nodes_potentials,& nodes_multipliers) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(INOUT)::nodes_potentials REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::arcs_voltages REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::tree_arcs_voltages REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::nodes_multipliers INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node INTEGER(KIND=i_byte)::orientation LOGICAL::non_pure n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) non_pure=PRESENT(nodes_multipliers) ArcOrNodeBased:IF(PRESENT(tree_arcs_voltages))THEN PropagateTreePotentials:DO index=0,n_nodes+n_special_nodes node=level_ordering(index) orientation=orientations(node) RootNode:IF(orientation==no_parent)THEN IF(non_pure)THEN nodes_potentials(node)=tree_arcs_voltages(node) ELSE nodes_potentials(node)=nodes_potentials(node)+tree_arcs_voltages(node) END IF ELSE parent_node=heads_tails(orientation,parents(node)) IF(non_pure)THEN nodes_potentials(node)=nodes_multipliers(node)*nodes_potentials(parent_node)+& tree_arcs_voltages(node) ELSE nodes_potentials(node)=nodes_potentials(parent_node)+tree_arcs_voltages(node) END IF END IF RootNode END DO PropagateTreePotentials ELSE IF(PRESENT(arcs_voltages))THEN PropagatePotentials:DO index=0,n_nodes+n_special_nodes node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN IF(non_pure)nodes_potentials(node)=0.0_r_wp CYCLE PropagatePotentials END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) Multipliers:IF(non_pure)THEN IF(orientation==tail_is_parent)THEN nodes_potentials(node)=nodes_multipliers(node)*nodes_potentials(parent_node)+& arcs_voltages(parent_arc) ELSE nodes_potentials(node)=nodes_multipliers(node)*nodes_potentials(parent_node)-& arcs_voltages(parent_arc) END IF ELSE IF(orientation==tail_is_parent)THEN nodes_potentials(node)=nodes_potentials(parent_node)+arcs_voltages(parent_arc) ELSE nodes_potentials(node)=nodes_potentials(parent_node)-arcs_voltages(parent_arc) END IF END IF Multipliers END DO PropagatePotentials END IF ArcOrNodeBased END SUBROUTINE PropagateNodesPotentials SUBROUTINE PropagateArcsFlows(arc_offset,node_offset,heads_tails,& orientations,parents,level_ordering,& supplies_demands,arcs_flows,tree_arcs_flows,& flow_imbalance,nodes_multipliers) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN)::supplies_demands REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(OUT),OPTIONAL::arcs_flows REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::tree_arcs_flows REAL(KIND=r_wp),INTENT(OUT),OPTIONAL::flow_imbalance REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::nodes_multipliers INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,grandparent_arc INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation,parent_orientation LOGICAL::non_pure n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) non_pure=PRESENT(nodes_multipliers) IF(PRESENT(flow_imbalance))flow_imbalance=0.0_r_wp NodeOrArcBased:IF(PRESENT(tree_arcs_flows))THEN tree_arcs_flows=supplies_demands PropagateTreeFlows:DO index=n_nodes+n_special_nodes,0,-1 node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN IF(PRESENT(flow_imbalance))flow_imbalance=flow_imbalance+tree_arcs_flows(node) CYCLE PropagateTreeFlows END IF parent_node=heads_tails(orientation,parents(node)) IF(non_pure)THEN tree_arcs_flows(parent_node)=tree_arcs_flows(parent_node)+& nodes_multipliers(node)*tree_arcs_flows(node) ELSE tree_arcs_flows(parent_node)=tree_arcs_flows(parent_node)+tree_arcs_flows(node) END IF END DO PropagateTreeFlows ELSE IF(PRESENT(arcs_flows))THEN InitialFlows:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)CYCLE InitialFlows parent_arc=parents(node) IF(orientation==tail_is_parent)THEN arcs_flows(parent_arc)=supplies_demands(node) ELSE arcs_flows(parent_arc)=-supplies_demands(node) END IF END DO InitialFlows PropagateFlows:DO index=n_nodes+n_special_nodes,0,-1 node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN IF(PRESENT(flow_imbalance))THEN flow_imbalance=flow_imbalance+supplies_demands(node) END IF CYCLE PropagateFlows END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) parent_orientation=orientations(parent_node) IF(parent_orientation==no_parent)THEN IF(PRESENT(flow_imbalance))THEN flow_imbalance=flow_imbalance+arcs_flows(parent_arc) END IF CYCLE PropagateFlows END IF grandparent_arc=parents(parent_node) Multipliers:IF(non_pure)THEN IF(parent_orientation==orientation)THEN arcs_flows(grandparent_arc)=arcs_flows(grandparent_arc)+& nodes_multipliers(node)*arcs_flows(parent_arc) ELSE arcs_flows(grandparent_arc)=arcs_flows(grandparent_arc)-& nodes_multipliers(node)*arcs_flows(parent_arc) END IF ELSE IF(parent_orientation==orientation)THEN arcs_flows(grandparent_arc)=arcs_flows(grandparent_arc)+arcs_flows(parent_arc) ELSE arcs_flows(grandparent_arc)=arcs_flows(grandparent_arc)-arcs_flows(parent_arc) END IF END IF Multipliers END DO PropagateFlows END IF NodeOrArcBased IF(PRESENT(tree_arcs_flows).AND.PRESENT(arcs_flows))THEN CopyFlows:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)CYCLE CopyFlows parent_arc=parents(node) IF(orientation==tail_is_parent)THEN arcs_flows(parent_arc)=tree_arcs_flows(node) ELSE arcs_flows(parent_arc)=-tree_arcs_flows(node) END IF END DO CopyFlows END IF IF(PRESENT(flow_imbalance))flow_imbalance=ABS(flow_imbalance) END SUBROUTINE PropagateArcsFlows SUBROUTINE CalculatePathLabels(arc_offset,node_offset,heads_tails,level_ordering,& orientations,parents,path_labels,arcs_weights,tree_type) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::path_labels REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_weights INTEGER(KIND=i_byte),INTENT(IN)::tree_type INTEGER(KIND=i_wp)::n_nodes,n_special_nodes INTEGER(KIND=i_wp)::index,node,arc,head,tail,parent_arc,parent_node INTEGER(KIND=i_byte)::orientation INTEGER::alloc_status REAL(KIND=r_wp)::initial_value n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,DIM=1),KIND=i_wp) IF(tree_type==min_cost_tree)THEN initial_value=-HUGE(1.0_r_wp) ELSE initial_value=HUGE(1.0_r_wp) END IF PathLabels:DO index=0,n_nodes+n_special_nodes node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN path_labels(node)=initial_value CYCLE PathLabels END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) IF(tree_type==min_cost_tree)THEN path_labels(node)=MAX(path_labels(parent_node),arcs_weights(parent_arc)) ELSE path_labels(node)=MIN(path_labels(parent_node),arcs_weights(parent_arc)) END IF END DO PathLabels END SUBROUTINE CalculatePathLabels SUBROUTINE ReBuildSpanningTree(arc_offset,node_offset,heads_tails,arcs_list,& arcs_weights,orientations,parents,cardinalities,& percent_pivoted,tree_type) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::arcs_list REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_weights INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::& parents,cardinalities INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(INOUT)::orientations REAL,INTENT(OUT),OPTIONAL::percent_pivoted INTEGER(KIND=i_byte),INTENT(IN)::tree_type INTEGER(KIND=i_dp)::n_traces,n_pivots INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::node,arc,arc_index,entering_arc,leaving_arc,climbers(2),head,tail,& climber,parent,climber_parent,climber_cardinality,cardinality,& cycle_root,subtree_root,hooking_node,unhooking_node,subtree_size INTEGER(KIND=i_byte)::orientation,climber_orientation,leaving_orientation INTEGER::alloc_status,side,leaving_side LOGICAL::weights_test,tree_arc,climbers_at_root(2) REAL(KIND=r_wp)::max_tree_weight n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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_traces=0 n_pivots=0 TestPivot:DO arc_index=INT(LBOUND(arcs_list,DIM=1),KIND=i_wp),INT(UBOUND(arcs_list,DIM=1),KIND=i_wp) entering_arc=arcs_list(arc_index) head=heads_tails(1,entering_arc) tail=heads_tails(2,entering_arc) tree_arc=((orientations(head)==tail_is_parent).AND.(parents(head)==entering_arc)).OR.& ((orientations(tail)==head_is_parent).AND.(parents(tail)==entering_arc)) IF(tree_arc.OR.(entering_arc==0))CYCLE TestPivot n_traces=n_traces+1 leaving_arc=entering_arc climbers=(/head,tail/) climbers_at_root=.FALSE. TraceCycle:DO IF(climbers(1)==climbers(2))THEN cycle_root=climbers(1) EXIT TraceCycle END IF ChooseSide:IF(climbers_at_root(1).AND.climbers_at_root(2))THEN CYCLE TestPivot ELSE IF(climbers_at_root(1))THEN side=2 ELSE IF(climbers_at_root(2))THEN side=1 ELSE IF(cardinalities(climbers(1))<=cardinalities(climbers(2)))THEN side=1 ELSE side=2 END IF END IF ChooseSide parent=parents(climbers(side)) orientation=orientations(climbers(side)) ClimbUp:IF(orientation==no_parent)THEN climbers_at_root(side)=.TRUE. ELSE IF(tree_type==min_cost_tree)THEN weights_test=(arcs_weights(leaving_arc)arcs_weights(parent)) END IF IF(weights_test)THEN leaving_arc=parent leaving_orientation=orientation leaving_side=side END IF climbers(side)=heads_tails(orientation,parent) END IF ClimbUp END DO TraceCycle PivotArc:IF(leaving_arc==entering_arc)THEN CYCLE TestPivot ELSE n_pivots=n_pivots+1 IF(leaving_orientation==head_is_parent)THEN unhooking_node=heads_tails(1,leaving_arc) subtree_root=heads_tails(2,leaving_arc) ELSE unhooking_node=heads_tails(2,leaving_arc) subtree_root=heads_tails(1,leaving_arc) END IF subtree_size=cardinalities(subtree_root) IF(leaving_side==1)THEN hooking_node=heads_tails(2,entering_arc) climber=heads_tails(1,entering_arc) orientation=head_is_parent ELSE hooking_node=heads_tails(1,entering_arc) climber=heads_tails(2,entering_arc) orientation=tail_is_parent END IF parent=entering_arc cardinality=0 ReverseParenthood:DO climber_parent=parents(climber) climber_orientation=orientations(climber) climber_cardinality=cardinalities(climber) parents(climber)=parent IF(orientation==head_is_parent)THEN orientations(climber)=tail_is_parent ELSE IF(orientation==tail_is_parent)THEN orientations(climber)=head_is_parent END IF cardinalities(climber)=subtree_size-cardinality IF(climber==subtree_root)EXIT ReverseParenthood climber=heads_tails(climber_orientation,climber_parent) parent=climber_parent orientation=climber_orientation cardinality=climber_cardinality END DO ReverseParenthood climber=unhooking_node DO WHILE(climber/=cycle_root) cardinalities(climber)=cardinalities(climber)-subtree_size climber=heads_tails(orientations(climber),parents(climber)) END DO climber=hooking_node DO WHILE(climber/=cycle_root) cardinalities(climber)=cardinalities(climber)+subtree_size climber=heads_tails(orientations(climber),parents(climber)) END DO END IF PivotArc END DO TestPivot IF(debug_graph_algs)THEN WRITE(*,*)"In ReBuildSpanningTree, n_traces=",n_traces,& ", n_pivots=",n_pivots,", percent=",REAL(n_pivots)/REAL(n_traces) END IF IF(PRESENT(percent_pivoted))percent_pivoted=REAL(n_pivots)/REAL(n_traces) END SUBROUTINE ReBuildSpanningTree SUBROUTINE BuildNetworkMST(arc_offset,node_offset,heads_tails,arcs_mask,arcs_weights,& level_ordering,orientations,parents,cardinalities,path_labels,tree_mask,& from_scratch,defragment_ordering,tree_type) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::arcs_mask REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::arcs_weights INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(OUT)::level_ordering INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::parents,cardinalities INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(INOUT)::orientations REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::path_labels LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(INOUT)::tree_mask LOGICAL,INTENT(IN),OPTIONAL::from_scratch LOGICAL,INTENT(IN),OPTIONAL::defragment_ordering INTEGER(KIND=i_byte),INTENT(IN)::tree_type INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::arcs_ordering,candidate_arcs,& nodes_roots,nodes_ordering INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::node,arc,arc_index,index,parent_arc,parent_node,& head,tail,n_candidate_arcs,offset,root_node INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation LOGICAL::initialize_tree,mask_arcs,rebuild_cardinalities,reorder_nodes,& tree_arc_test,path_labels_test n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) initialize_tree=.FALSE. IF(PRESENT(from_scratch))initialize_tree=from_scratch mask_arcs=PRESENT(arcs_mask) reorder_nodes=.FALSE. IF(PRESENT(defragment_ordering))reorder_nodes=defragment_ordering BuildMST:IF((tree_type==max_cost_tree).OR.(tree_type==min_cost_tree))THEN OrderArcs:IF(initialize_tree)THEN rebuild_cardinalities=.TRUE. ALLOCATE(arcs_ordering(1:n_special_arcs+n_arcs+1),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1_i_wp,& caller="BuildNetworkMST",alloc_status=alloc_status) CALL QuickRank(array=arcs_weights,permutation=arcs_ordering,& pivot_selection='R',partially_ranked=.FALSE.) arcs_ordering=arcs_ordering-n_special_arcs-1 IF(mask_arcs)THEN n_candidate_arcs=0 DO arc_index=1,n_special_arcs+n_arcs+1 arc=arcs_ordering(arc_index) IF(arcs_mask(arc))THEN n_candidate_arcs=n_candidate_arcs+1 arcs_ordering(n_candidate_arcs)=arc END IF END DO ELSE n_candidate_arcs=n_special_arcs+n_arcs+1 END IF IF(tree_type==max_cost_tree)THEN CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,& arcs_list=arcs_ordering(n_candidate_arcs:1:(-1)),heads_tails=heads_tails,& orientations=orientations,parents=parents,use_parenthood=.FALSE.) ELSE CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,& arcs_list=arcs_ordering(1:n_candidate_arcs),heads_tails=heads_tails,& orientations=orientations,parents=parents,use_parenthood=.FALSE.) END IF CALL RecordAllocation(n_elements=-INT(SIZE(arcs_ordering),KIND=i_wp),mold=1_i_wp) DEALLOCATE(arcs_ordering,STAT=alloc_status) ELSE OrderArcs rebuild_cardinalities=.FALSE. ALLOCATE(candidate_arcs(1:n_special_arcs+n_arcs+1),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1_i_wp,& caller="UpdateSpanningTree",alloc_status=alloc_status) CALL CalculatePathLabels(arc_offset=n_special_arcs,node_offset=n_special_nodes,& heads_tails=heads_tails,orientations=orientations,parents=parents,& level_ordering=level_ordering,path_labels=path_labels,& arcs_weights=arcs_weights,tree_type=tree_type) n_candidate_arcs=0 DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) tree_arc_test=.NOT.tree_mask(arc) IF(mask_arcs)tree_arc_test=tree_arc_test.AND.arcs_mask(arc) IF(tree_arc_test)THEN IF(tree_type==min_cost_tree)THEN path_labels_test=& (arcs_weights(arc)<=MAX(path_labels(head),path_labels(tail))) ELSE path_labels_test=& (arcs_weights(arc)>=MIN(path_labels(head),path_labels(tail))) END IF IF(path_labels_test)THEN n_candidate_arcs=n_candidate_arcs+1 candidate_arcs(n_candidate_arcs)=arc END IF END IF END DO CALL ReBuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,& arcs_list=candidate_arcs(1:n_candidate_arcs),heads_tails=heads_tails,& arcs_weights=arcs_weights,orientations=orientations,parents=parents,& cardinalities=cardinalities,tree_type=tree_type) CALL RecordAllocation(n_elements=-INT(SIZE(candidate_arcs),KIND=i_wp),mold=1_i_wp) DEALLOCATE(candidate_arcs,STAT=alloc_status) END IF OrderArcs ELSE BuildMST rebuild_cardinalities=.TRUE. ALLOCATE(candidate_arcs(1:n_special_arcs+n_arcs+1),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1_i_wp,& caller="UpdateSpanningTree",alloc_status=alloc_status) IF(mask_arcs)THEN n_candidate_arcs=0 DO arc=-n_special_arcs,n_arcs IF(arcs_mask(arc))THEN n_candidate_arcs=n_candidate_arcs+1 candidate_arcs(n_candidate_arcs)=arc END IF END DO ELSE DO arc=-n_special_arcs,n_arcs candidate_arcs(arc+n_special_arcs+1)=arc END DO n_candidate_arcs=n_special_arcs+n_arcs+1 END IF CALL BuildSpanningTree(arc_offset=n_special_arcs,node_offset=n_special_nodes,& arcs_list=candidate_arcs(1:n_candidate_arcs),heads_tails=heads_tails,& orientations=orientations,parents=parents,& use_parenthood=.NOT.initialize_tree) CALL RecordAllocation(n_elements=-INT(SIZE(candidate_arcs),KIND=i_wp),mold=1_i_wp) DEALLOCATE(candidate_arcs,STAT=alloc_status) END IF BuildMST ReorderNodes:IF(reorder_nodes)THEN ALLOCATE(nodes_ordering(0:n_special_nodes+n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1_i_wp,& caller="BuildNetworkMST",alloc_status=alloc_status) ALLOCATE(nodes_roots(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1_i_wp,& caller="BuildNetworkMST",alloc_status=alloc_status) CALL BuildTreeStructures(arc_offset=n_special_arcs,node_offset=n_special_nodes,& heads_tails=heads_tails,orientations=orientations,parents=parents,& cardinalities=cardinalities,level_ordering=nodes_ordering,tree_mask=tree_mask,& rebuild_cardinalities=rebuild_cardinalities) offset=0 DefragmentOrdering:DO index=0,n_nodes+n_special_nodes node=nodes_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN level_ordering(offset)=node parents(node)=offset nodes_roots(node)=node offset=offset+cardinalities(node) CYCLE DefragmentOrdering END IF parent_node=heads_tails(orientation,parents(node)) root_node=nodes_roots(parent_node) nodes_roots(node)=root_node parents(root_node)=parents(root_node)+1 level_ordering(parents(root_node))=node END DO DefragmentOrdering CALL RecordAllocation(n_elements=-INT(SIZE(nodes_ordering),KIND=i_wp),mold=1_i_wp) DEALLOCATE(nodes_ordering,STAT=alloc_status) CALL RecordAllocation(n_elements=-INT(SIZE(nodes_roots),KIND=i_wp),mold=1_i_wp) DEALLOCATE(nodes_roots,STAT=alloc_status) ELSE ReorderNodes CALL BuildTreeStructures(arc_offset=n_special_arcs,node_offset=n_special_nodes,& heads_tails=heads_tails,orientations=orientations,parents=parents,& cardinalities=cardinalities,level_ordering=level_ordering,tree_mask=tree_mask,& rebuild_cardinalities=rebuild_cardinalities) END IF ReorderNodes END SUBROUTINE BuildNetworkMST SUBROUTINE BasicDualNewtonSystem(node_offset,arc_offset,heads_tails,& orientations,parents,level_ordering,tree_mask,& nodes_excess_potentials,nodes_excess_flows,& arcs_conductances,arcs_resistances,& arcs_excess_voltages,tree_arcs_excess_voltages,& nodes_conductances,nodes_resistances,nodes_multipliers) IMPLICIT NONE INTEGER,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)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN)::tree_mask REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(INOUT)::nodes_excess_potentials REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_excess_flows REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::& arcs_conductances,arcs_resistances REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(OUT),OPTIONAL::arcs_excess_voltages REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::tree_arcs_excess_voltages REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::& nodes_conductances,nodes_resistances,nodes_multipliers INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,head,tail INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation INTEGER::temp n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) ArcOrNodeBased:IF(PRESENT(arcs_excess_voltages))THEN IF(PRESENT(nodes_multipliers))THEN CALL PropagateArcsFlows(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& supplies_demands=nodes_excess_flows,& arcs_flows=arcs_excess_voltages,nodes_multipliers=nodes_multipliers) ELSE CALL PropagateArcsFlows(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& supplies_demands=nodes_excess_flows,& arcs_flows=arcs_excess_voltages) END IF IF(PRESENT(arcs_conductances))THEN WHERE(tree_mask)arcs_excess_voltages=arcs_excess_voltages/arcs_conductances END IF IF(PRESENT(arcs_resistances))THEN WHERE(tree_mask)arcs_excess_voltages=arcs_excess_voltages*arcs_resistances END IF IF(PRESENT(nodes_conductances))THEN Conductances:DO node=-n_special_nodes,n_nodes IF(orientations(node)==no_parent)CYCLE Conductances parent_arc=parents(node) arcs_excess_voltages(parent_arc)=arcs_excess_voltages(parent_arc)/& nodes_conductances(node) END DO Conductances END IF IF(PRESENT(nodes_resistances))THEN Resistances:DO node=-n_special_nodes,n_nodes IF(orientations(node)==no_parent)CYCLE Resistances parent_arc=parents(node) arcs_excess_voltages(parent_arc)=arcs_excess_voltages(parent_arc)*& nodes_resistances(node) END DO Resistances END IF IF(PRESENT(nodes_multipliers))THEN CALL PropagateNodesPotentials(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& arcs_voltages=arcs_excess_voltages,nodes_potentials=nodes_excess_potentials,& nodes_multipliers=nodes_multipliers) ELSE CALL PropagateNodesPotentials(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& arcs_voltages=arcs_excess_voltages,nodes_potentials=nodes_excess_potentials) END IF ELSE IF(PRESENT(tree_arcs_excess_voltages))THEN IF(PRESENT(nodes_multipliers))THEN CALL PropagateArcsFlows(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& supplies_demands=nodes_excess_flows,& tree_arcs_flows=tree_arcs_excess_voltages,nodes_multipliers=nodes_multipliers) ELSE CALL PropagateArcsFlows(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& supplies_demands=nodes_excess_flows,& tree_arcs_flows=tree_arcs_excess_voltages) END IF IF(PRESENT(arcs_conductances))THEN DO node=-n_special_nodes,n_nodes IF(orientations(node)==no_parent)THEN tree_arcs_excess_voltages(node)=0.0_r_wp ELSE tree_arcs_excess_voltages(node)=tree_arcs_excess_voltages(node)/arcs_conductances(parents(node)) END IF END DO END IF IF(PRESENT(arcs_resistances))THEN DO node=-n_special_nodes,n_nodes IF(orientations(node)==no_parent)THEN tree_arcs_excess_voltages(node)=0.0_r_wp ELSE tree_arcs_excess_voltages(node)=tree_arcs_excess_voltages(node)*arcs_resistances(parents(node)) END IF END DO END IF IF(PRESENT(nodes_conductances))THEN tree_arcs_excess_voltages=tree_arcs_excess_voltages/nodes_conductances END IF IF(PRESENT(nodes_resistances))THEN tree_arcs_excess_voltages=tree_arcs_excess_voltages*nodes_resistances END IF IF(PRESENT(nodes_multipliers))THEN CALL PropagateNodesPotentials(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& tree_arcs_voltages=tree_arcs_excess_voltages,nodes_potentials=nodes_excess_potentials,& nodes_multipliers=nodes_multipliers) ELSE CALL PropagateNodesPotentials(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,& orientations=orientations,parents=parents,& level_ordering=level_ordering,& tree_arcs_voltages=tree_arcs_excess_voltages,nodes_potentials=nodes_excess_potentials) END IF END IF ArcOrNodeBased END SUBROUTINE BasicDualNewtonSystem SUBROUTINE EffectiveCycleResistances(arc_offset,node_offset,heads_tails,& parents,cardinalities,orientations,level_ordering,& tree_mask,arcs_resistances,cycles_resistances) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN)::tree_mask INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents,cardinalities INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_resistances REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(OUT)::cycles_resistances REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::nodes_resistances INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,& climbers(2),head,tail,parent,cycle_root INTEGER(KIND=i_byte)::orientation INTEGER::side,alloc_status n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) ALLOCATE(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="EffectiveCycleResistances",alloc_status=alloc_status) PathResistances:DO index=0,n_nodes+n_special_nodes node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN nodes_resistances(node)=0.0_r_wp CYCLE PathResistances END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) nodes_resistances(node)=nodes_resistances(parent_node)+arcs_resistances(parent_arc) END DO PathResistances CycleResistances:DO arc=-n_special_arcs,n_arcs IF(tree_mask(arc))THEN cycles_resistances(arc)=0.0_r_wp ELSE head=heads_tails(1,arc) tail=heads_tails(2,arc) climbers=(/head,tail/) TraceCycle:DO IF(climbers(1)==climbers(2))THEN cycle_root=climbers(1) EXIT TraceCycle END IF IF(cardinalities(climbers(1))<=cardinalities(climbers(2)))THEN side=1 ELSE side=2 END IF parent=parents(climbers(side)) orientation=orientations(climbers(side)) climbers(side)=heads_tails(orientation,parent) END DO TraceCycle cycles_resistances(arc)=nodes_resistances(head)-nodes_resistances(cycle_root)+& nodes_resistances(tail)-nodes_resistances(cycle_root) END IF END DO CycleResistances END SUBROUTINE EffectiveCycleResistances SUBROUTINE FactorizeIncompleteQR(arc_offset,node_offset,heads_tails,& orientations,parents,level_ordering,tree_mask,& arcs_conductances,nodes_conductances,nodes_resistances,& nodes_multipliers) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN)::tree_mask REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_conductances REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::& nodes_conductances,nodes_resistances REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::nodes_multipliers INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,head,tail INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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(PRESENT(nodes_conductances))THEN nodes_conductances=0.0_r_wp Conductances_C:DO arc=-n_special_arcs,n_arcs IF(.NOT.tree_mask(arc))THEN 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) END IF END DO Conductances_C Multipliers_C:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)THEN nodes_multipliers(node)=1.0_r_wp CYCLE Multipliers_C END IF parent_arc=parents(node) nodes_conductances(node)=nodes_conductances(node)+arcs_conductances(parent_arc) nodes_multipliers(node)=arcs_conductances(parent_arc)/nodes_conductances(node) END DO Multipliers_C ELSE IF(PRESENT(nodes_resistances))THEN nodes_resistances=0.0_r_wp Conductances_R:DO arc=-n_special_arcs,n_arcs IF(.NOT.tree_mask(arc))THEN 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) END IF END DO Conductances_R Multipliers_R:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)THEN nodes_multipliers(node)=1.0_r_wp CYCLE Multipliers_R END IF parent_arc=parents(node) nodes_resistances(node)=nodes_resistances(node)+arcs_conductances(parent_arc) nodes_multipliers(node)=arcs_conductances(parent_arc)/nodes_resistances(node) END DO Multipliers_R nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) END IF END SUBROUTINE FactorizeIncompleteQR SUBROUTINE FactorizeIncompleteLDLt(arc_offset,node_offset,heads_tails,& arcs_mask,orientations,parents,level_ordering,tree_mask,& arcs_conductances,nodes_conductances,nodes_resistances,& nodes_multipliers,diagonal_nodes_conductances) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::arc_offset,node_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::arcs_mask INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN)::parents INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN),OPTIONAL::tree_mask REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_conductances REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::& nodes_conductances,nodes_resistances REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::nodes_multipliers REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::& diagonal_nodes_conductances INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,head,tail INTEGER::alloc_status INTEGER(KIND=i_byte)::orientation n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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(PRESENT(nodes_conductances))THEN InitializeConductances_C:IF(PRESENT(diagonal_nodes_conductances))THEN nodes_conductances=diagonal_nodes_conductances AddDiagonal_C:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)CYCLE AddDiagonal_C nodes_conductances(node)=nodes_conductances(node)+arcs_conductances(parents(node)) END DO AddDiagonal_C ELSE nodes_conductances=0.0_r_wp IF(PRESENT(arcs_mask))THEN MaskedConductances_C:DO arc=-n_special_arcs,n_arcs IF(arcs_mask(arc))THEN 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) END IF END DO MaskedConductances_C ELSE Conductances_C: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) END DO Conductances_C END IF END IF InitializeConductances_C Multipliers_C:DO index=n_nodes+n_special_nodes,0,-1 node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN nodes_multipliers(node)=1.0_r_wp CYCLE Multipliers_C END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) nodes_multipliers(node)=arcs_conductances(parent_arc)/nodes_conductances(node) nodes_conductances(parent_node)=nodes_conductances(parent_node)-& nodes_multipliers(node)*arcs_conductances(parent_arc) END DO Multipliers_C ELSE IF(PRESENT(nodes_resistances))THEN InitializeConductances_R:IF(PRESENT(diagonal_nodes_conductances))THEN nodes_resistances=diagonal_nodes_conductances AddDiagonal_R:DO node=-n_special_nodes,n_nodes orientation=orientations(node) IF(orientation==no_parent)CYCLE AddDiagonal_R nodes_resistances(node)=nodes_resistances(node)+arcs_conductances(parents(node)) END DO AddDiagonal_R ELSE nodes_resistances=0.0_r_wp IF(PRESENT(arcs_mask))THEN MaskedConductances_R:DO arc=-n_special_arcs,n_arcs IF(arcs_mask(arc))THEN 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) END IF END DO MaskedConductances_R ELSE Conductances_R: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) END DO Conductances_R END IF END IF InitializeConductances_R Multipliers_R:DO index=n_nodes+n_special_nodes,0,-1 node=level_ordering(index) orientation=orientations(node) IF(orientation==no_parent)THEN nodes_multipliers(node)=1.0_r_wp CYCLE Multipliers_R END IF parent_arc=parents(node) parent_node=heads_tails(orientation,parent_arc) nodes_multipliers(node)=arcs_conductances(parent_arc)/nodes_resistances(node) nodes_resistances(parent_node)=nodes_resistances(parent_node)-& nodes_multipliers(node)*arcs_conductances(parent_arc) END DO Multipliers_R nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) END IF END SUBROUTINE FactorizeIncompleteLDLt SUBROUTINE CreateAdjacencyArrays(arc_offset,node_offset,& heads_tails,neighbours,my_neighbours,incident_arcs,& n_edges,include_loops) IMPLICIT NONE INTEGER,INTENT(IN)::node_offset,arc_offset INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(OUT)::neighbours,my_neighbours,incident_arcs INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::n_edges LOGICAL,INTENT(IN),OPTIONAL::include_loops INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::nodes_degrees INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs,& n_vertices INTEGER(KIND=i_wp)::node,arc,head,tail,special_node,special_arc,& degree,head_neighbour,tail_neighbour,& node_shift,arc_shift,arc_counter,node_counter INTEGER::alloc_status,status LOGICAL::graph_OK,add_loops n_special_nodes=node_offset n_nodes=INT(UBOUND(my_neighbours,DIM=1),KIND=i_wp)-2-n_special_nodes n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp) n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp) n_vertices=n_nodes+n_special_nodes+1 node_shift=n_special_nodes+1 arc_shift=n_special_arcs+1 add_loops=.TRUE. IF(PRESENT(include_loops))add_loops=include_loops IF(.NOT.ALLOCATED(nodes_degrees))THEN ALLOCATE(nodes_degrees(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,& caller="CreateAdjacencyArrays",alloc_status=alloc_status) END IF nodes_degrees=0 DO arc=-n_special_arcs,n_arcs head=heads_tails(1,arc) tail=heads_tails(2,arc) IF((head/=tail).OR.add_loops)THEN nodes_degrees(head)=nodes_degrees(head)+1 nodes_degrees(tail)=nodes_degrees(tail)+1 END IF END DO arc_counter=1 DO node=-n_special_nodes,n_nodes my_neighbours(node+node_shift)=arc_counter arc_counter=arc_counter+nodes_degrees(node) END DO my_neighbours(0)=0 my_neighbours(n_vertices+1)=arc_counter IF(PRESENT(n_edges))n_edges=arc_counter-1 DO arc=n_arcs,-n_special_arcs,-1 head=heads_tails(1,arc) tail=heads_tails(2,arc) IF((head/=tail).OR.add_loops)THEN nodes_degrees(head)=nodes_degrees(head)-1 head_neighbour=my_neighbours(head+node_shift)+nodes_degrees(head) nodes_degrees(tail)=nodes_degrees(tail)-1 tail_neighbour=my_neighbours(tail+node_shift)+nodes_degrees(tail) neighbours(head_neighbour)=tail+node_shift neighbours(tail_neighbour)=head+node_shift incident_arcs(head_neighbour)=arc incident_arcs(tail_neighbour)=arc END IF END DO neighbours(0)=0 incident_arcs(0)=0 IF(ALLOCATED(nodes_degrees))THEN CALL RecordAllocation(n_elements=-INT(SIZE(nodes_degrees),KIND=i_wp),mold=1_i_wp) DEALLOCATE(nodes_degrees,STAT=alloc_status) END IF END SUBROUTINE CreateAdjacencyArrays SUBROUTINE CalculateCongestionsDilations(node_offset,arc_offset,heads_tails,& orientations,parents,cardinalities,level_ordering,tree_mask,joints_and_roots,& arcs_conductances,arcs_congestions_dilations,nodes_conductances,normalize) IMPLICIT NONE INTEGER,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)::parents,cardinalities INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(IN)::orientations INTEGER(KIND=i_wp),DIMENSION(0:),INTENT(IN)::level_ordering LOGICAL(KIND=l_wp),DIMENSION(-arc_offset:),INTENT(IN)::tree_mask INTEGER(KIND=i_wp),DIMENSION(-arc_offset:),INTENT(OUT)::joints_and_roots REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_conductances REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(OUT)::arcs_congestions_dilations REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::nodes_conductances LOGICAL,INTENT(IN),OPTIONAL::normalize INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_special_arcs,n_arcs INTEGER(KIND=i_wp)::index,node,arc,parent_arc,parent_node,head,tail,& climbers(2),cycle_root,parent INTEGER::side,alloc_status INTEGER(KIND=i_byte)::orientation INTEGER::temp n_special_nodes=-INT(LBOUND(parents,DIM=1),KIND=i_wp) n_nodes=INT(UBOUND(parents,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) BEPJoints:DO arc=-n_special_arcs,n_arcs IF(tree_mask(arc))THEN joints_and_roots(arc)=0 ELSE head=heads_tails(1,arc) tail=heads_tails(2,arc) climbers=(/head,tail/) TraceCycle:DO IF(climbers(1)==climbers(2))THEN cycle_root=climbers(1) EXIT TraceCycle END IF IF(cardinalities(climbers(1))<=cardinalities(climbers(2)))THEN side=1 ELSE side=2 END IF parent=parents(climbers(side)) orientation=orientations(climbers(side)) IF(orientation==no_parent)THEN cycle_root=0 EXIT TraceCycle END IF climbers(side)=heads_tails(orientation,parent) END DO TraceCycle joints_and_roots(arc)=cycle_root END IF END DO BEPJoints nodes_conductances=0.0_r_wp DO arc=-n_special_arcs,n_arcs IF(.NOT.tree_mask(arc))THEN nodes_conductances(heads_tails(1:2,arc))=nodes_conductances(heads_tails(1:2,arc))+& arcs_conductances(arc) nodes_conductances(joints_and_roots(arc))=nodes_conductances(joints_and_roots(arc))-& arcs_conductances(arc)-arcs_conductances(arc) END IF END DO CALL PropagateArcsFlows(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,orientations=orientations,& parents=parents,level_ordering=level_ordering,& supplies_demands=nodes_conductances,arcs_flows=arcs_congestions_dilations) IF(PRESENT(normalize))THEN IF(normalize)THEN WHERE(tree_mask) arcs_congestions_dilations=& arcs_congestions_dilations/(arcs_conductances+EPSILON(1.0_r_wp)) END WHERE END IF END IF nodes_conductances(level_ordering(0:1))=0.0_r_wp CALL PropagateNodesPotentials(arc_offset=arc_offset,node_offset=node_offset,& heads_tails=heads_tails,orientations=orientations,& parents=parents,level_ordering=level_ordering,& arcs_voltages=arcs_congestions_dilations,nodes_potentials=nodes_conductances) DO arc=-n_special_arcs,n_arcs IF(.NOT.tree_mask(arc))THEN arcs_congestions_dilations(arc)=& nodes_conductances(heads_tails(1,arc))-nodes_conductances(joints_and_roots(arc))+& nodes_conductances(heads_tails(2,arc))-nodes_conductances(joints_and_roots(arc)) END IF END DO WHERE(tree_mask)arcs_congestions_dilations=ABS(arcs_congestions_dilations) END SUBROUTINE CalculateCongestionsDilations END MODULE Graph_Algorithms