@z This file was created by Aleksandar Donev as part of the Network Optimization project. Feel free to use any portion of it and contact me at donev@pa.msu.edu @x \Title{Graph and Network Algorithms} \author{Aleksandar Donev} \date{February 2000} @*0 Module |Graph_Algorithms|. This module provides several graph algorithms to be used with network optimization codes. For example, a connected-components labeling routine is provided for extracting a single connected component from a graph. Also provided are very efficient novel routines for manipulating spanning trees related to network optimization. Finally, some routines related to updating minimal spanning trees under small perturbations to the weights are provided, as well as routines for conversion from heads-and-tails format into the more commonly used adjacency-list-based data structures. A spanning tree in the sense used here is a subset of the arcs that contains no cycles, and has a chosen {\em root node}. Once the root is known, parenthood relations are determined uniquely. In this library, parenthood relations are usually stored in the array called |parents(-n_special_nodes:n_nodes)|, where the parent of each node is in fact a pointer (the index) of the {\em parent arc}, i.e. the arc that connects the node to its parent. Since our graph is stored as usual in a heads-tails array, we need to know whether the head is the parent, or the tail of the parent arc is the parent. Therefore each parent arc has its orientation, usually stored in the array |orientations(-n_special_nodes:n_nodes)|, which can be one of the predefined public parameters |head_is_parent|, |tail_is_parent|, or |no_parent| (for example, for rooting arcs). Root nodes should not have a parent, and this is especially imprortant if dealing with spanning forests, which this library can handle in most routines. One would usually use the sign of the entry in |parents| to denote the orientation, but since I allow for arcs with negative numbers (special arcs), I do not do this in this code. The parenthood relations described above uniquely determine a tree and a way to go from a node up to the root of the tree. A final important data structure in manipulating the tree is a {\em thread-like ordering of the nodes}, which is a permutation of the nodes, usually stored in the array |level_ordering(0:n_nodes+n_special_nodes)| in these codes. The requirment here is that {\em every node comes after its parent} in |level_ordering|... In other words, by traversing the nodes in forward order, one would always visit the parent {\em before} visiting the child. Vice versa, by traversing the nodes in the reverse-thread order, one would always visit {\em all} children of a node before visiting the node itself. These are crucial operations in traversing a tree when calculating flows or potentials in network flow algorithms. Please note that this array is not the usual thread-order defined by computer scientists. My definition is a very loose one, and so the ordering is not unique. I have developed a way to create this ordering in a two-pass procedure over the nodes of the tree, as given in the routine |BuildTreeStructures|, suited for a serial computer. Finally, in maintaining and working with the trees, two more auxilliary node labels are used. One is the {\em cardinality} label for the nodes, which gives the number of nodes in the subtree rooted at the given node, {\em counting the node itself}, and is usually stored in |cardinalities(-n_special_nodes:n_nodes)|. When talking about minimal or maximal weight spanning trees, denoted by MST from now on (I support both min and max ones, not just the usual min-weight trees), we have weights on the arcs, and the {\em path node labels} array |path_labels(-n_special_nodes:n_nodes)| stores either the minimal (for MaxST) or maximal (for MinST) weight on the path from the node to the root. These paths do not have to be exact, only exact upper or lower bounds. To tell the codes what kind of tree one wants to build or manipulate, the predefined parameters |any_tree| (meaning an arbitrary spanning tree), |random_tree| (which is like |any_tree|, only a different tree will be generated each time because an expansive random permutation of the arcs is used), |min_cost_tree| (for MinST) or |max_cost_tree| (for MaxST). These are further used in the front-end module |Network_Spanning_Trees|. Finally, the logical flag |debug_graph_algs| tells whether some useful statistics should be gathered (not complete yet). The list of routines in this module along with a short description (please look below for more details) is: \begin{description} \item[|ConnectedComponents|] labels each node with a number that is a chosen unique node from the connected component it belongs to. \item[|BuildSpanningTree|] makes a spanning tree or forest by using the union-find algorithm. It accepts a list of arcs to consider in entering the tree, and checks whether each arc forms a cycle, and includes it in the tree if it does. This allows great flexibility in using this routine. For example, on can start with half-made trees and add a collection of arcs. Or, one can order the arcs according to weight and then build the MST, which is Kruskal's algorithm for building MST's. But remember to {\em always initialize the tree parenthood relations to whatever it is that you want to start from}. This routine creates just the parenthood relations for the tree. \item[|BuildTreeStructures|] completes a spanning tree with the cardinality labels and creates a thread like level ordering of the nodes. The cardinalities may be provided as input if they are known, and this speeds the routine (for example, the routine |UpdateSpanningTree| maintains cardinality labels exactly, so these do not need to be rebuilt). \item[|CalculatePathLabels|] calculates the path labels for the nodes in the weighted spanning tree. \item[|UpdateSpanningTree|] takes a tree and arcs weights and a list of arcs to consider for entering either the minimal or maximal spanning tree. It can be used to update MST's after small perturbation of the weights, or to build an MST from scratch by starting with a random tree. It uses the cardinality and parenthood relations to traverse Basic Equivalent Path (BEP) cycles and performs arc pivots. \item[|PropagateNodesPotentials|] propagates potentials from the root of a tree down to the leaf nodes, given potentials for the roots and potential drops accross each of the arcs. In terms of network optimization, this is a multiplication of an arc vector with the matrix $B^{-T}$, where $B$ is the basis for the node-arc incidence matrix. This routine can deal with forests as well. \item[|PropagateArcsFlows|] propagates flows from the leaf nodes to the root by using the given supplies-demands on each node. This is in fact multiplication of a nodal vector with $B^{-1}$. \end{description} @f index _index @ @a 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. @; // Whether to print useful statistics INTEGER(KIND=i_byte), PARAMETER, PUBLIC :: no_parent=0, & head_is_parent=1, tail_is_parent=2 @; // Possible orientations for parenthood relations INTEGER(KIND=i_byte), PARAMETER, PUBLIC :: any_tree=1, random_tree=2, & min_cost_tree=3, max_cost_tree=4 @; // Types of spanning trees INTEGER(KIND=i_byte), PARAMETER, PUBLIC :: matching_node=2, dead_node=0 @; // Statuses used in matching code CONTAINS @; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; @@; END MODULE Graph_Algorithms @; @*0 Graph Connectivity Labeling. The routine |ConnectedComponents| takes a graph and just labels the nodes according to the connected cluster they belong to, where each cluster get's its own label. The somewhat more expensive routine |BackboneCycles| is a bit more detailed and labels each cluster of cycles (blob in some terminology) with a unique label, so that dangling ends and the percolating backbone can be identified in a given graph. The first routine uses plain union-find algorithm as its basis, while backbone indentification requires both a union-find and a matching algorithm. @*1 Labeling Connected Clusters. This procedure will label the connected components of a graph with integers (not controllable at this point) and possibly return the label for the largest connected-component. Therefore, all nodes that are connected to eachother will be given the same integer label in the array |labels| on exit, and this integer will be returned for the largest component in the optional argument |largest_component_label| along with (optionally) the size of this component in |largest_component_size|. The algorithm can start with already partially labeled graphs--forests (for example, when the edges are processed in batches, as in a parallel system where inter-processor edges must be handled separately), in which case |use_labels| should be set to |.TRUE.| (default is |.FALSE.|). Otherwise, it starts from a forest of all single-node trees and builds from there. But notice that here |labels| uses the root node for each connected-component tree as the label, and this can not be changed. So not just any integers can be used as labels. This was done in order to avoid using two arrays, one for the labels, the other for the parenthood relations, but it can easily be changed. @ The following two macros perform two basic operation in union-find algorithms and basically and tree algorithm--they find the root of a given node |node| using the parenthood relations in the array |parents|. The second macro performs path compression along the path to the root: @m _FindRoot(parents, node, root, climber) @; climber = node @; // Start at the node DO WHILE (parents[climber] != climber) @; // Until you reach a node that's its own parent climber = parents[climber] @; // Climb one level up the tree END DO @; root = climber @; @m _CompressPath(parents, node, root, climber, climber_parent) @; climber = node @; // Start at the node DO WHILE (parents[climber] != root) @; // Until you reach a node that's its own parent @%% counter = counter + 1 @; // This counts the {\bf nodes at depth 2 or higher} climber_parent = parents[climber] @; // Climb one level up the tree parents[climber] = root @; // Jump pointer to root climber = climber_parent @; END DO @; @ Here is the routine |ConnectedComponents|. The array labels is of dimension |(-node_offset:)| and should be addressed as such in the given heads-and-tails array |heads_tails|. This is not to be confused or the code will crash! @= 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 @; // Heads-tails array for $G$ INTEGER(KIND=i_wp), INTENT(IN) :: node_offset @; // Lower bound for nodal arrays INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: labels @; // With displaced counting LOGICAL, INTENT(IN), OPTIONAL :: use_labels @; INTEGER(KIND=i_wp), INTENT(OUT), OPTIONAL :: largest_component_label, largest_component_size @; // Largest component, if needed 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=-_LBOUND(labels, i_wp, DIM=1) @; n_nodes=_UBOUND(labels, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_arcs=_SIZE(heads_tails, i_wp,DIM=2) @; // $m$ 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 @; // Start with forest of all nodes 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 @; // Initialize heights ForAllArcs: DO arc = 1, n_arcs @; // Not masked for now head = heads_tails [1, arc] @; tail = heads_tails [2, arc] @; _FindRoot(labels, head, head_root, climber) @; _FindRoot(labels, tail, tail_root, climber) @; _CompressPath(labels, head, head_root, climber, climber_parent) @; _CompressPath(labels, tail, tail_root, climber, climber_parent) @; TreeUnion: IF (head_root != tail_root) THEN @; // Perform union operation diff_heights = (heights [head_root] - heights [tail_root]) @; IF (diff_heights > 0) THEN @; // Choose which tree merges with which parent = head_root @; child = tail_root @; ELSE @; parent = tail_root @; child = head_root @; END IF @; labels (child) = parent @; // Hook the two trees--union IF (diff_heights == 0) heights[parent] = heights[parent]+1 @; // Otherwise heights are not-changed by union operation END IF TreeUnion @; END DO ForAllArcs @; ForAllNodes: DO node = -n_special_nodes, n_nodes @; // Final compression of all paths _FindRoot(labels, node, root, climber) @; _CompressPath(labels, node, root, climber, climber_parent) @; END DO ForAllNodes @; /* We sometimes need to extract the largest component, so here I find this component and its label. This will make who ever uses this routine much happier in many cases, as sorting would be needed to determine this in the general case. */ heights=0 @; // Tree sizes will be stored here DO node = -n_special_nodes, n_nodes @; // Count tree sizes heights(labels[node])=heights(labels[node])+1 @; END DO @; IF(PRESENT(largest_component_label).OR.PRESENT(largest_component_size)) THEN @; largest_component=_MAXLOC(heights, i_wp) @; // Assumes numbering starting at 1! largest_component=largest_component-1-node_offset @; // Correct for numberin 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=-_SIZE(heights, i_wp), mold=1_i_wp) @; DEALLOCATE (heights) @; END SUBROUTINE ConnectedComponents @; @*1 Condensing Cycles. The structure of this routine is almost the same as |ConnectedComponents|, only now each backbone gets it's own label, as opposed to each connected cluster. These labels are returned in the array supernodes, and one can establish preexisting supernodes by passing |use_supernodes=.TRUE.|. The algorithm that is used is a relatively novel approach and it is based on previous work in rigidity percolation. It is essentially based on matching, and is an algorithm in which unique paths are established by making a forest in the network by adding the arcs one-by-one. Before adding an arc we must check whether there is a cycle formed in this way in the {\whole} network. This is done by reversing the parenthood relations (stored here in the array |cycles|) on the path from the head to the root of its forest tree (called a supernode here) and then try to do the same with the tail. If there is a cycle through the network, then this will identify it, if not, the edge can be added to the cycles. Once a cycle is identified, the nodes in the cycle are merged together into a supernode using the famed union-find algorithm (set relations here are stored in the array |supernodes|). This process is extremely efficient and soo clever, but there is little documentation of it here or elsewhere. Contact me if it sounds interesing. Note that this particular version does some redundant operations for beauty and simplicity's sake, and is not fully optimized. @ @= SUBROUTINE BackboneCycles(node_offset, heads_tails, supernodes, & use_supernodes, largest_supernode, largest_supernode_size) @; IMPLICIT NONE @; INTEGER(KIND=i_wp), INTENT(IN) :: node_offset @; // Lower bound for nodal arrays INTEGER(KIND=i_wp), DIMENSION(:,:), INTENT(IN) :: heads_tails @; // Heads-tails array INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: supernodes @; // Cycle supernodes LOGICAL, INTENT(IN), OPTIONAL :: use_supernodes @; // Are there already some supernodes? INTEGER(KIND=i_wp), INTENT(OUT), OPTIONAL :: largest_supernode, largest_supernode_size @; INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: cycles, heights @; // Path labels for tracing cycles and the heights of the union-find trees INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_arcs @; // Counters 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=-_LBOUND(supernodes, i_wp, DIM=1) @; n_nodes=_UBOUND(supernodes, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_arcs=_SIZE(heads_tails, i_wp,DIM=2) @; // $m$ 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 @; // Start with forest of all nodes 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 @; // Start with no cycles 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 @; // Initialize heights ForAllArcs: DO arc = 1, n_arcs @; // Not masked for now! head = heads_tails [1, arc] @; tail = heads_tails [2, arc] @; _UncoverNode(head, supernodes, cycles, head, head_supernode, last_supernode, head_length, & supersupernode, increase_height) @; // Uncover the head first and see if a cycle was formed _UncoverNode(tail, supernodes, cycles, tail, tail_supernode, last_supernode, tail_length, & supersupernode, increase_height) @; // Now also uncover the tail @%% head_supernode=supernodes[head] @; // The supernode for the head @%% tail_supernode=supernodes[tail] @; // The supernode for the tail /* We know that the supertrees on the cycle have already been compressed, so that we only need to climb at most one step up to find the supernode: */ AddArc: IF(last_supernode==head_supernode) THEN @; // Condense the identified cycle @%% WRITE(*,*) "Condensing arc: ", arc, supernodes(heads_tails[1:2,arc]) @; supernode=head_supernode @; CondenseCycle: DO @; // Trace back the cycle supernodes[supernode]=supersupernode @; // Label all with the shortest supernode IF(supernode==tail_supernode) @~ EXIT CondenseCycle @; // Cycle was traced supernode=cycles[supernode] @; // Move along the cycle (already compressed) @%% WRITE(*,*) supernodes(cycles[supernode]), cycles[supernode] @; END DO CondenseCycle @; IF(increase_height) @~ heights[supersupernode]=heights[supersupernode]+1 @; ELSE @; IF(head_length>=tail_length) THEN @; // Try to make as short of cycles as possible cycles[head_supernode]=tail_supernode @; // Add the arc to the cycles ELSE @; cycles[tail_supernode]=head_supernode @; // Add this arc to the cycles END IF @; END IF AddArc @; END DO ForAllArcs @; CALL RecordAllocation(n_elements=-_SIZE(cycles, i_wp), mold=1_i_wp) @; DEALLOCATE(cycles) @; ForAllNodes: DO node = -n_special_nodes, n_nodes @; // Final compression of all supertrees _FindRoot(supernodes, node, supernode, climber) @; _CompressPath(supernodes, node, supernode, climber, climber_parent) @; END DO ForAllNodes @; /* We sometimes need to extract the largest component, so here I find this component and its supernode label. This will make who ever uses this routine much happier in many cases, as sorting would be needed to determine this in the general case. */ heights=0 @; // Tree sizes will be stored here DO node = -n_special_nodes, n_nodes @; // Count tree sizes heights(supernodes[node])=heights(supernodes[node])+1 @; END DO @; IF(PRESENT(largest_supernode).OR.PRESENT(largest_supernode_size)) THEN @; largest_location=_MAXLOC(heights, i_wp) @; // Assumes numbering starting at 1! largest_location=largest_location-1-node_offset @; // Correct for numbering 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=-_SIZE(heights, i_wp), mold=1_i_wp) @; DEALLOCATE (heights) @; END SUBROUTINE BackboneCycles @; @ The following macro is a crucial operation in traversing the potential cycles formed by additions of arcs. It starts at |node| and it traverses its cycle until it reaches an uncovered node, i.e. a parent supernode. Notice that each node in the cycle has its supernode, and all operations must be performed on the supernodes to avoid updating the nodes all the time. Also notice that full path-compression is performed here, and in my experience this is a good optimization. I also decided to use height labels for the supernode union-sets, even though these do not bring any improvement and demand some extra memeory. The nice thing about using heights in union-find is that they guarantee a logarithmic time bound, and may be useful in pathological graphs. @m _UncoverNode(ID, supernodes, cycles, node, first_supernode, last_supernode, path_length, supersupernode, equal_heights) @; path_length=0 @; max_height=-HUGE(1_i_wp) @; _FindRoot(supernodes, node, supernode, climber) @; // Find the first cluster supernode _CompressPath(supernodes, node, supernode, climber, climber_parent) @; // Shorten supertree increase_height=.FALSE. @; previous_supernode=supernode @; // Start from uncovered node first_supernode=supernode @; Uncover_@e@&ID: DO @; // Reverse the cycle path of this node @%% supersupernode=MIN(supersupernode, supernode) @; // Choose a common supernode label IF(heights[supernode]>max_height) THEN @; // Find the tallest supertree and merge all to it max_height=heights[supernode] @; supersupernode=supernode @; ELSE IF(heights[supernode]==max_height) THEN @; // There is a subtree of the same height equal_heights=.TRUE. @; END IF @; node=cycles[supernode] @; // Find the next node in the``half'' cycle cycles[supernode]=previous_supernode; // Reverse the parenthood relation previous_supernode=supernode @; // Save for later _FindRoot(supernodes, node, supernode, climber) @; // Find the cluster supernode _CompressPath(supernodes, node, supernode, climber, climber_parent) @; // Shorten supertrees IF(supernode==previous_supernode) THEN @; // We reached the end of a path or cycle last_supernode=supernode @; EXIT Uncover_@e@&ID @; END IF @; path_length=path_length+1 @; // The length of this piece of a cycle END DO Uncover_@e@&ID @; @*0 Graph Contraction Via Matching. The following routine is useful for ordering the nodes or arcs in a graph based on their connectivity proximity, as used in the module |Support_Trees| in the routine |ST_MatchingNodesMapping|. The routine |ContractMatchedArcs| is a rather general routine though, and it performs matching of the arcs contained in the linked list of arcs |arcs_sequence|. An arc in the sequence will be matched only if its two endpoints are not matched already. It is meant to be called inside a loop from another routine that updates the list of arcs to try to match, so all of the arrays needed are passed to it as arguments. Also, it uses some clever schemes to mask which nodes to match, and this is done based on the parity of the number of graph contractions performed so far. I made these matching routines myself and it would take a while to fully document them. They are not really complete and optimized yet either, so I leave this for some future time. One thing to notice is that I use arc number 0 explicitly as a sentinel end-of-list value, so it can not be included in the matching even if you want to (for some silly reason). @ @= 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 @; // First arc in the list 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(KIND=l_wp), DIMENSION(-arc_offset:), INTENT(INOUT) :: arcs_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=-_LBOUND(supernodes, i_wp, DIM=1) @; n_nodes=_UBOUND(supernodes, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ matched_node=1_i_byte @; IF(PRESENT(odd_cycle)) THEN @; IF(odd_cycle) @~ matched_node=-matched_node @; // On even cycles END IF @; @%% WRITE(*,*) "Matched nodes mask in this loop:", matched_node @; 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 @; // Initialize the arrays DO node = -n_special_nodes, n_nodes @; supernodes[node] = node @; // Start with forest of all nodes END DO @; IF(sequence_nodes) THEN @; DO node = -n_special_nodes, n_nodes @; nodes_ordering[node] = node @; // Every node by itself end_nodes[node]=node @; END DO @; END IF @; heights=0 @; END IF @; /* Now we start matching the supernodes using the arcs in the arcs sequence (linked list): */ previous_arc=0 @; arc=first_arc @; ContractArcs: DO @; // Not masked IF(arc==0) @~ EXIT ContractArcs @; // We exhausted all arcs @; head = heads_tails [1, arc] @; tail = heads_tails [2, arc] @; delete_arc=.FALSE. @; _FindRoot(supernodes, head, head_root, climber) @; _FindRoot(supernodes, tail, tail_root, climber) @; _CompressPath(supernodes, head, head_root, climber, climber_parent) @; _CompressPath(supernodes, tail, tail_root, climber, climber_parent) @; TreeUnion: IF ((head_root==tail_root).OR. & ((nodes_mask[head_root]==dead_node).OR.(nodes_mask[tail_root]==dead_node))) THEN @; // The arc is redundant delete_arc=.TRUE. @; ELSE @; // Perform union operation 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) @; // Mark as unmatched in this round 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) @; // Mark as unmatched in this round END IF @; ContractArc: IF((nodes_mask[head_root]!=matched_node).AND. & nodes_mask[tail_root]!=matched_node) THEN @; // Both endpoints are unmatched and available @%% WRITE(*,*) "Matching arc:", arc, " (", head, tail," )", & @%% " Supernodes:", head_root, tail_root, & @%% " Mask: ", nodes_mask[head_root], nodes_mask[tail_root] @; delete_arc=.TRUE. @; // Skip this arc in the list diff_heights = (heights [head_root] - heights [tail_root]) @; IF (diff_heights > 0) THEN @; // Choose which tree merges with which parent = head_root @; child = tail_root @; ELSE @; parent = tail_root @; child = head_root @; END IF @; nodes_mask[parent]=matched_node @; // Mark the node as a matched supernode nodes_mask[child]=dead_node @; // We are done with this node now supernodes[child]=parent @; // Hook the two trees--union IF (diff_heights == 0) heights[parent] = heights[parent]+1 @; IF(sequence_nodes) THEN @; // Keep a time record of the matching end_node=end_nodes[parent] @; // The end of one of the ordering lists nodes_ordering[end_node]=child @; // Merge the two lists for the matching end_nodes[parent]=end_nodes[child] @; END IF @; END IF ContractArc @; END IF TreeUnion @; IF(delete_arc) THEN @; // Arc is to be skipped IF(arc==first_arc) THEN @; // We must move the pointer to the first arc first_arc=arcs_sequence[arc] @; // The next arc becomes first ELSE @; // Shortcut this arc in the list arcs_sequence[previous_arc]=arcs_sequence[arc] @; // Skip this arc in the list END IF @; ELSE @; previous_arc=arc @; END IF @; arc=arcs_sequence[arc] @; // Move to the next arc in the list END DO ContractArcs @; arc=first_arc @; UpdateNodes: DO @; IF(arc==0) @~ EXIT UpdateNodes @; // We exhausted all arcs @; head=heads_tails[1,arc] @; tail=heads_tails[2,arc] @; _FindRoot(supernodes, head, head_root, climber) @; _FindRoot(supernodes, tail, tail_root, climber) @; /* The following two may be redundant work (?): */ _CompressPath(supernodes, head, head_root, climber, climber_parent) @; _CompressPath(supernodes, tail, tail_root, climber, climber_parent) @; IF(contract_heads_tails) THEN @; heads_tails[1,arc]=head_root @; // Update HT 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] @; // Move to next arc END DO UpdateNodes @; END SUBROUTINE ContractMatchedArcs @; @*0 Bulding Spanning Trees. This section deals with routines for building spanning trees of the kind needed in network flow algorithms. Then I will give algorithms for maintaning these spanning trees under changes of the weights of the arcs. @*1 Bulding Parenthood Relations. This routine will build a spanning tree by examining the arcs in |arcs_list| one by one and adding all arc that do not cause a cycle. It also builds parenthood relations for the tree. It also works for forests as well, as do all spanning-tree related routines in this module. The logical flag |use_parenthood| (defaulting to |.FALSE.|) should be used to tell whether the routine should use the already provided parenthood relations (in |parents| and |orientations|), or initialize these to a forest of all single nodes and start building from there (the default). As in all routines I give, the nodes are numbered starting from |-n_special_nodes|, and the arcs are numbered from |-n_special_arcs|, which are called |node_offset| and |arc_offset| in the argument lists from now on. These do not default to 1, but have to always be given (Fortran restrictions). @ @= 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 @; // Lower bound for arc and nodal arrays INTEGER(KIND=i_wp), DIMENSION(:), INTENT(IN) :: arcs_list @; // Order in which to examine arcs INTEGER(KIND=i_wp), DIMENSION(:,-arc_offset:), INTENT(IN) :: heads_tails @; // Heads-tails array for $G$ INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: parents @; // With displaced counting INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(INOUT) :: orientations @; // Orientations for basic arcs in the tree, one of |no_parent|, |tail_is_parent| or |head_is_parent| LOGICAL, INTENT(IN), OPTIONAL :: use_parenthood @; INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: heights, forest, cardinalities @; // The heights of the union-find trees and the union sets in the spanning forest, and approximate cardinalities of the tree nodes INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_special_arcs, n_arcs @; // Counters 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 @; // Temporaries INTEGER :: alloc_status @; INTEGER(KIND=i_byte) :: orientation, climber_orientation @; // Temporaries LOGICAL :: initialize_forest @; n_special_nodes=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ initialize_forest=.TRUE. @; // Start from a single-node forest 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 @; // Start with forest of all nodes separately END DO @; orientations=no_parent @; // No parenthood relations at start parents=0 @; // Dummy arc is parent of all at start cardinalities=1 @; // Each node is by itself ELSE @; // Find the root of each node to initalize the forest InitializeTrees: DO node = -n_special_nodes, n_nodes @; // Initialize the forest orientation=orientations[node] @; IF(orientation!=no_parent) THEN @; forest[node]=heads_tails(orientation, parents[node]) @; // Find the parent ELSE @; forest[node] = node @; // Node is by itself in the forest END IF @; END DO InitializeTrees @; cardinalities=0 @; // This may be an expensive operation if we start with an almost complete spanning tree: ShortenTrees: DO node = -n_special_nodes, n_nodes @; // Compress all paths _FindRoot(forest, node, root, climber) @; _CompressPath(forest, node, root, climber, climber_parent) @; cardinalities[root]=cardinalities[root]+1 @; END DO ShortenTrees @; END IF @; heights=0 @; // All trees are of height 0 (or 1, as you please...) ForAllArcs: DO index=1, _SIZE(arcs_list, i_wp) @; // Try to add each arc in |arcs_list| to the tree arc=arcs_list[index] @; head = heads_tails [1, arc] @; tail = heads_tails [2, arc] @; /* To decide whether to add this arc to the tree I use the standard union-find algorithm, as in the connected-components labeling algorithm. This requires extra memory, but speeds the decision process considerably: */ _FindRoot(forest, head, head_root, climber) @; _FindRoot(forest, tail, tail_root, climber) @; _CompressPath(forest, head, head_root, climber, climber_parent) @; _CompressPath(forest, tail, tail_root, climber, climber_parent) @; AddArc: IF (head_root != tail_root) THEN @; // Add this arc to the tree diff_heights = (heights [head_root] - heights [tail_root]) @; // Which set tree is higher IF (diff_heights > 0) THEN @; // Choose which tree merges with which parent = head_root @; child = tail_root @; ELSE @; parent = tail_root @; child = head_root @; END IF @; forest [child] = parent @; // Hook the two trees--union IF (diff_heights == 0) heights[parent] = heights[parent]+1 @; // Otherwise heights are not-changed by union operation cardinalities[parent]=cardinalities[parent]+cardinalities[child] @; // Update cardinalities /* Now the difficult and expensive part is maintaining the parenthood relations after this arc was added to the tree--this requires reversing parenthood orientations on the path from either the head or tail to its true parent root. Here I assume that a smaller tree (smaller cardinality) implies a shorter path to the true parent root, which is usually true for trees with no special structure: */ diff_cardinalities = (cardinalities[head_root] - cardinalities[tail_root]) @; IF (diff_cardinalities > 0) THEN @; // Choose which tree needs parents updating parent = arc @; orientation=tail_is_parent @; climber = tail @; ELSE @; parent = arc @; orientation=head_is_parent @; climber = head @; END IF @; ReverseParenthood: DO @; // Reverse the parenthood relations on one of the root paths climber_orientation=orientations[climber] @; // Save this climber_parent=parents[climber] @; // Save parents[climber]=parent @; // Reverse the parenthood relations IF(orientation==head_is_parent) THEN @; // Reverse orientation of parent-relation 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 @; // Reached the root climber=heads_tails[climber_orientation, climber_parent] @; // Climb up the tree parent=climber_parent @; // Go one step up to the parent root orientation=climber_orientation @; END DO ReverseParenthood @; END IF AddArc @; END DO ForAllArcs @; CALL RecordAllocation(n_elements=-_SIZE(forest, i_wp), mold=1_i_wp) @; DEALLOCATE(forest) @; CALL RecordAllocation(n_elements=-_SIZE(heights, i_wp), mold=1_i_wp) @; DEALLOCATE(heights) @; CALL RecordAllocation(n_elements=-_SIZE(cardinalities, i_wp), mold=1_i_wp) @; DEALLOCATE(cardinalities) @; END SUBROUTINE BuildSpanningTree @; @*1 Building Additional Tree Data Structures. Although the parenthood relations built by |BuildSpanningTree| uniquely determine the spanning tree, a few other data-structures, such as the cardinalities of the nodes or a thread-like root-to-leaf ordering of the nodes is needed in most tree algorithms. These are built in the routine |BuildTreeStructures|. If the cardinalities, are already known, |rebuild_cardinalities| should be set to |.FALSE.|, and the routine will be speeded up. Please note that the dummy node at position 0 is not treated as a regular node but is left over as the first node in the ordering always. The way that the node ordering |level_ordering| is built is a novel simple procedure which starts from the leafs of the trees (identified through their cardinality being 1--first pass through the nodes) and adds successive levels of parent nodes to the ordering until the root is reached (second pass through nodes). To determine when a node is ready to be added, I make sure that all the children of the node have been already added (simply via counting them!), and then add the node. Other schemes are possible too, such as a breath-first search from the root node, but we do not have the data structures for that. This procedure is actually very fast and simple and works well. @ @= 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 @; // Lower bound for arc and nodal arrays INTEGER(KIND=i_wp), DIMENSION(:,-arc_offset:), INTENT(IN) :: heads_tails @; // Heads-tails array LOGICAL(KIND=l_wp), DIMENSION(-arc_offset:), INTENT(OUT), OPTIONAL :: tree_mask @; // A tree mask INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(IN) :: parents @; // With displaced counting INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(IN) :: orientations @; INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: cardinalities @; // Cardinalities of tree nodes INTEGER(KIND=i_wp), DIMENSION(0:), INTENT(OUT) :: level_ordering @; // A thread-like ordering of the tree nodes by levels--{\em starting from the root and going to the leaf nodes} LOGICAL, INTENT(IN), OPTIONAL :: rebuild_cardinalities @; // Build cardinalities from scratch or use old ones?--default is build from scratch INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: counters @; // A temporary array 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 @; // An indicator n_special_nodes=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ new_cardinalities=.TRUE. @; IF(PRESENT(rebuild_cardinalities)) @~ new_cardinalities=rebuild_cardinalities @; /* We need either the cardinalities or the nodal child degrees in order to initialize the level ordering of the nodes. If there are no cardinalities provided, it is easier to build child degrees (the number of children of each node counting the node itself): */ IF(new_cardinalities) THEN @; // Initialize the cardinalities temporarily with nodal child degrees cardinalities=1 @; @%% cardinalities[0]=0 @; // Dummy node is special ChildDegrees: DO node=-n_special_nodes, n_nodes @; // Count the number of children of each node orientation=orientations[node] @; IF(orientation!=no_parent) THEN @; // Add 1 to the number of children of the parent 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) @; // Temporary counters CALL RecordAllocation( n_elements=n_nodes+n_special_nodes+1, mold=1_i_wp, & caller="BuildSpanningTree", alloc_status=alloc_status) @; /* This way of building an ordering of the nodes into levels may not be the best if we are dealing with a forest of many smaller trees, since it puts all the leafs first, then all the one-level-higher nodes, and so on. A better scheme in that case would be to put all the nodes of one tree in the forest next to one another--this will speed cache performance most likely. This can be done easily since we have cardinalities and we can calculate offsets. But in general using the ordering computed here the nodes can be reordered trivially later on, so I leave this task up to an external user routine if the user thinks this will be beneficial. */ level_start=_UBOUND(level_ordering, i_wp)+1 @; // Initialize to end of array level_end=level_start @; // Sweeping pointer as nodes are added to the ordering DO node=-n_special_nodes, n_nodes @; // Find the leaf nodes to start the level structure @%% IF(orientations[node]==no_parent) @~ WRITE(*,*) "Root node:", node, cardinalities[node] @; IF(cardinalities[node]==1) THEN @; // This is a leaf non-dummy node @%% WRITE(*,*) "Leaf node:", node @; level_end=level_end-1 @; level_ordering[level_end]=node @; // Add to first level END IF @; END DO @; @%% WRITE(*,*) "Level ordering at present:", level_ordering[level_end:] counters=1 @; @%% counters[0]=0 @; @%% level_ordering[0]=0 @; // Dummy node always goes in front BuildLevels: DO WHILE(level_end>=1) @; level_start=level_start-1 @; // Add the parent of this node to the ordering node=level_ordering[level_start] @; orientation=orientations[node] @; IF(orientation==no_parent) @~ CYCLE BuildLevels @; // Root node parent_node=heads_tails(orientation, parents[node]) @; // Try to add this node to the list @%% WRITE(*,*) "level_bounds:", level_start, level_end, " node:", node @; IF(new_cardinalities) THEN @; // |counters| is child-degree counters[parent_node]=counters[parent_node]+1 @; ELSE @; // |counters| is in fact the cardinality counters[parent_node]=counters[parent_node]+counters[node] @; END IF @; IF(counters[parent_node]==cardinalities[parent_node]) THEN @; // Add the node to the level list level_end=level_end-1 @; level_ordering[level_end]=parent_node @; END IF @; END DO BuildLevels @; @%% CALL QuickRank(array=level_ordering, permutation=counters) @; @%% WRITE(*,*) "Level order:", level_ordering @; @%% WRITE(*,*) "Ordered :", level_ordering(counters-1) @; CALL RecordAllocation(n_elements=-_SIZE(counters, i_wp), mold=1_i_wp) @; DEALLOCATE(counters) @; IF(new_cardinalities) THEN @; // Now calculate the true cardinalities cardinalities=1 @; @%% cardinalities[0]=0 @; // Dummy node 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 @; // Build a mask that selects the tree arcs from the non-tree arcs 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 @; @*0 Linear Algebra with $B^{-1}$ and $B^{-T}$. After a tree is built, it determines a basis $B$ for the network. The routine |PropagateNodesPotentials| multiplies a tension vector $t$ with $B^{-T}$ to give a potential vector $\lambda=B^{-T} \, t$. The routine |PropagateArcsFlows| multiplies a supply-demand vector $b$ with $B^{-1}$ to give a flow vector $x=B^{-1} \, b$. These are essentially routines that traverse the tree from the root down to leaf nodes or vice versa and propagate either potentials or flows along the way. These are essential steps in algorithms based on trees. Some of the routines support a non-pure basis $F$, which has the same structure as a basis $B$, but the off-diagonal entries are not $1$ or $-1$ but some real multiplier (associated here with the node whose parent-arc we are looking at via the array called |nodes_multipliers|). These are meant to be used with some of the preconditioning/factorization routines in the next section. Note that the routines in this module are not meant to perform atomistic operations, such as changing the potential of one node, or pivoting one arc, as used in simplex codes. Rather, these are streamlined operations acting on the whole network meant to be used in conjugate-gradient solvers to perform matrix-vector multiplications. @*1 Propagating Potentials. This routine will essentially mutiply the arc vector |arcs_voltages| with $B^{-T}$ to give the nodal vector |nodes_potentials|. It traverses the tree from the root to the leafs using |level_ordering|. The most important thing to watch for is that this routine expects the potentials of the root nodes (i.e. nodes which do not have real parents) to be pre-set before calling this routine. The usual way of setting these to 0 is not a good design choice in this case, especially since this rouitine can deal with forests, which have multiple root. Also, it is important to note that although the routine accepts the potential drops on all arcs as input, only the values for the tree arcs will actually be used, the rest will simply be ignored. An optional argument |nodes_multipliers| can be used when the basis $B$ is not a pure basis, but some tree matrix $F$ with non-unit off-diagonal entries. Here each multiplier is associated to the corresponding tree arc via the node whose parent arc it is (to save on wasted storage for any multipliers of non-tree arcs). {\bf Important note on efficiency:} Although I made efforts to have the code be efficient, clarity and versatility were my focus. In particular, some of the |IF| statements inside the loops can be eliminated if we are certain that we are dealing with a spanning tree, and not a spanning forest. Such an assumption is made in the next section, but not here because these routines make sense ona forest as well. Also, statements such as: @ @= 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 @; @ could be replaced with code that avoids the conditional and instead uses extra flops: |nodes_potentials[node]=nodes_potentials[parent_node]+(3-2*orientation)*arcs_voltages[parent_arc]| which would be more efficient on a lot of machines. However, on the Pentium I benchmarked this code on this was not the case and conditionals seemed to be handled well inside loops, so I chose to go with the elegant version. But this may need to be changed on some other architectures and compilers! @ @= 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 @; // Lower bound for arc and nodal arrays INTEGER(KIND=i_wp), DIMENSION(:,-arc_offset:), INTENT(IN) :: heads_tails @; // Heads-tails array for $G$ INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(IN) :: parents @; // With displaced counting INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(IN) :: orientations @; INTEGER(KIND=i_wp), DIMENSION(0:), INTENT(IN) :: level_ordering @; // A thread-like ordering of the tree nodes by levels--{\em starting from the root} REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: nodes_potentials @; // The potentials of the nodes--the roots should have their potentials pre-set! REAL(KIND=r_wp), DIMENSION(-arc_offset:), INTENT(IN), OPTIONAL :: arcs_voltages @; // The potential drops accross the arcs--only tree arcs will be used REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(IN), OPTIONAL :: tree_arcs_voltages @; // Only tree arcs REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(IN), OPTIONAL :: nodes_multipliers @; // Arc multipliers, used when $B$ is not a pure basis but has arc gains 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ non_pure=PRESENT(nodes_multipliers) @; ArcOrNodeBased: IF(PRESENT(tree_arcs_voltages)) THEN @; // Faster and better PropagateTreePotentials: DO index=0, n_nodes+n_special_nodes @; // Start from root node and propagate the potentials down the tree to the leaf nodes /* Here I assume that |nodes_potentials[node]| for root nodes is the potential of the rooting virtual node connected to (the root of) this tree in the forest. Please note that this is handled differently below, where I assume that |nodes_potentials[node]| for root nodes is the potential of the root node itself, since in that case I do not explicitly have virtual rooting nodes (no storage array for that). These are both equivalent if we are dealing with one tree only and not a forest: WRONG: THE POTENTIALS OF THE ROOTS ARE OVERWRITTEN IN THE CASE WITH MULTIPLIERS */ node=level_ordering[index] @; orientation=orientations[node] @; RootNode: IF(orientation==no_parent) THEN @; IF(non_pure) THEN @; // Use a virtual parent (the node itself) nodes_potentials[node]=tree_arcs_voltages[node] @; // Overwrite the input ELSE @; nodes_potentials[node]=nodes_potentials[node]+tree_arcs_voltages[node] @; END IF @; ELSE @; // Use the parent potential parent_node=heads_tails[orientation, parents(node)] @; // This is the parent 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 @; // Slower! PropagatePotentials: DO index=0, n_nodes+n_special_nodes @; // Start from root node and propagate the potentials down the tree to the leaf nodes node=level_ordering[index] @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; IF(non_pure) @~ nodes_potentials[node]=0.0_r_wp @; // We must ground it CYCLE PropagatePotentials @; END IF @; // Leave root potentials alone! parent_arc=parents[node] @; parent_node=heads_tails[orientation, parent_arc] @; // This is the parent node 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 @; @%% WRITE(*,*) nodes_potentials[parent_node], nodes_potentials[node] @; END IF Multipliers @; END DO PropagatePotentials @; END IF ArcOrNodeBased @; END SUBROUTINE PropagateNodesPotentials @; @*1 Propagating Flows. This routine will essentially mutiply the nodal vector |supplies_demands| with $B^{-1}$ to give the arc vector |tree_flows|. It traverses the tree from the leafs to the root using |level_ordering|. This routine will leave the rooting arcs, i.e. the ficticious arc |parents[node]| for which |orientations[node]==no_parent| alone. If the supplies-demands vector adds up to a total in/outflow of zero, the flow that should go through the rooting arc will also be zero. To allow for this not to be the case and give the user some useful info, the actual flow that is left over after all the pushing is done is returned in the optional |flow_imbalance|. Again, non-pure bases are supported via the optional |nodes_multipliers|. No array temporary is needed here for storing the excess flows at each node, even though such a temporary would speed the loops somewhat. It is also important to note that I leave the roots and their rooting arcs alone in this collection of libraries. The user will probably have to deal with them separately before or after calling these routines. @ @= 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 // The feasible flows in the tree arcs. The non-tree arcs are {\em not} assigned flows! REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(OUT), OPTIONAL :: tree_arcs_flows @; // Only tree arcs are of interest REAL(KIND=r_wp), INTENT(OUT), OPTIONAL :: flow_imbalance @; // Total excess flow in the forest REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(IN), OPTIONAL :: nodes_multipliers @; // Arc multipliers, used when $B$ is not a pure basis but has arc gains 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ non_pure=PRESENT(nodes_multipliers) @; IF(PRESENT(flow_imbalance)) @~ flow_imbalance=0.0_r_wp @; NodeOrArcBased: IF(PRESENT(tree_arcs_flows)) THEN @; // Faster and simpler--try it first // We assume |tail_is_parent| for all tree arcs!!! tree_arcs_flows=supplies_demands @; // Initialize PropagateTreeFlows: DO index=n_nodes+n_special_nodes, 0, -1 @; // Start from leaf nodes and propagate the flows up to the root, without intermediate storage 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 @; // This is a tree root END IF @; parent_node=heads_tails[orientation, parents(node)] @; // This is the parent 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 @; // More complicated and slower 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 @; // Start from leaf nodes and propagate the flows up to the root, without intermediate storage 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] @; // This is the parent node parent_orientation=orientations[parent_node] @; IF(parent_orientation==no_parent) THEN @; // The parent is a root IF(PRESENT(flow_imbalance)) THEN @; flow_imbalance=flow_imbalance+arcs_flows[parent_arc] @; END IF @; CYCLE PropagateFlows @; // Leave root potentials alone! 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 @; // Copy the result 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) @; // Total flow excess END SUBROUTINE PropagateArcsFlows @; @*0 Maintaining Minimal/Maximal Spanning Trees. The following set of routines provides algorithms for efficient maintainance of minimal spanning trees under changes in the edge weights. @*1 Initializing Node Path Labels. The routine |CalculatePathLabels| calculates path labels for the nodes in a spanning forest. It does this with an efficient single-pass loop through the nodes from the root to the leafs (using |level_ordering|). Recall that for a minimal-weight spanning trees path labels are defined as the maximal weight (or an upper bound on it, to be precise) on the path from the node to the root, while for maximal-weight STs it is a lower bound on the minimal weight on the path to the root. These labels are not really used in this module, but they can be efficiently used to reject a lot of candidate arcs for entry into a spanning tree just by looking at the path labels, as is done in the module |Network_Spanning_Trees|: @ @= 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 @; // Heads-tails INTEGER(KIND=i_wp), DIMENSION(0:), INTENT(IN) :: level_ordering @; // Ordering from root to nodes 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 @; // Max or Min REAL(KIND=r_wp), DIMENSION(-arc_offset:), INTENT(IN) :: arcs_weights @; // Weights on arcs INTEGER(KIND=i_byte), INTENT(IN) :: tree_type @; // |min_cost_tree| or |max_cost_tree| 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; 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 @; // Start from root node and propagate the potentials down the tree to the leaf nodes node=level_ordering[index] @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; path_labels[node]=initial_value @; // Initialize the root path labels CYCLE PathLabels @; // Leave root potentials alone! END IF @; parent_arc=parents[node] @; parent_node=heads_tails[orientation, parent_arc] @; // This is the parent node 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 @; @*1 Rebuilding a Spanning Tree. The following routine will update a spanning tree by checking whether the arcs in |arcs_list| offer advantegous exchanges with some arc in their basic equivalent path (BEP) cycle. If there is such an advantegous change, the examined arc enters the spanning tree, while the arc with the maximal or minimal weight in the BEP leaves the tree. This is a pivoting operation, and requires updating the parenthood and cardinality relations. Cardinalities are used to trace the BEP without having to go all the way to the root, as described in other people's work on simplex-like codes and similar {\bf non-greedy} MST algorithms. The user can choose whether he wants a minimal or maximal weight ST via |tree_type|, and a helpful statistics on the number of actual pivots performed, divided by the number of arcs (cycles) examined, is returned in the optional |percent_pivoted|. If this is small the starting tree can be judged near-optimal with good confidence. @ @= 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 @; // Order in which to examine non-tree arcs for entry into the tree. These do {\em not} have to be non-tree arcs though. REAL(KIND=r_wp), DIMENSION(-arc_offset:), INTENT(IN) :: arcs_weights @; // Weights on arcs INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(INOUT) :: & parents, cardinalities @; // For tree nodes INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(INOUT) :: orientations @; // For tree arcs REAL, INTENT(OUT), OPTIONAL :: percent_pivoted @; // How many traces resulted in pivots INTEGER(KIND=i_byte), INTENT(IN) :: tree_type @; // |min_cost_tree| or |max_cost_tree| /* Local variables: */ INTEGER(KIND=i_dp) :: n_traces, n_pivots @; // Operation counters INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_special_arcs, n_arcs @; // Counters 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 @; @%% REAL :: dice @; // A random dice for random pivoting operations n_special_nodes=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ n_traces=0 @; // Counter for the number of BEP cycle traces n_pivots=0 @; // Counter for the number of pivots @%% CALL ResetTimer(50) @; @%% CALL ResetTimer(51) @; @%% CALL StartTimer(40) @; TestPivot: DO arc_index=_LBOUND(arcs_list, i_wp), _UBOUND(arcs_list, i_wp) @; entering_arc=arcs_list[arc_index] @; head=heads_tails[1,entering_arc] @; tail=heads_tails[2,entering_arc] @; // We test to make sure this is not a tree arc in order to make the code more robust: 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 @; // This is a tree arc! @%% WRITE(*,*) "Considering arc: ", entering_arc @; /* Test BEP cycle for an exchange of arcs: */ n_traces=n_traces+1 @; leaving_arc=entering_arc @; climbers=(/head, tail/) @; // Start tracing the cycle formed by |entering_arc| climbers_at_root=.FALSE. @; // In case we are dealing with a forest! TraceCycle: DO @; IF(climbers[1]==climbers[2]) THEN @; // Found a joint in the tree cycle_root=climbers[1] @; @%% WRITE(*,*) "Found joint at: ", cycle_root @; EXIT TraceCycle @; END IF @; ChooseSide: IF(climbers_at_root[1].AND.climbers_at_root[2]) THEN @; // This arc crossed between trees in the forest--we have different tree roots @%% WRITE(*,*) "Arc crossed between forests" @; CYCLE TestPivot @; ELSE IF(climbers_at_root[1]) THEN @; // One climber finished at a root side=2 @; ELSE IF(climbers_at_root[2]) THEN @; // One climber finished at a root side=1 @; ELSE @; // Use cardinalities to make a smart decision IF(cardinalities[climbers[1]]<=cardinalities[climbers[2]]) THEN @; side=1 @; // Climb on side 1 ELSE @; side=2 @; // Climb on side 2 END IF @; END IF ChooseSide @; parent=parents(climbers[side]) @; // Candidate to leave the tree--an arc in the cycle orientation=orientations(climbers[side]) @; ClimbUp: IF(orientation==no_parent) THEN @; // Reached a root climbers_at_root[side]=.TRUE. @; ELSE @; @%% @e#!_CompareWeights(weights_test, , arcs_weights[parent]) @; 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] @; // Climb up one level on this side END IF ClimbUp @; END DO TraceCycle @; PivotArc: IF(leaving_arc==entering_arc) THEN @; // Arc did not qualify to enter the MST @%% WRITE(*,*) "Arc: ", entering_arc, " did not enter the tree" @; CYCLE TestPivot @; ELSE @; // We need to perform a pivot operation @%% WRITE(*,*) "Arc: ", entering_arc, " replacing ", leaving_arc @; n_pivots=n_pivots+1 @; IF(leaving_orientation==head_is_parent) THEN @; unhooking_node=heads_tails[1, leaving_arc] @; // This is where subtree unhooks from subtree_root=heads_tails[2,leaving_arc] @; ELSE @; unhooking_node=heads_tails[2, leaving_arc] @; // This is where subtree unhooks from 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 @; // Initialize for |ReverseParenthood| cardinality=0 @; // Initialize @%% WRITE(*,*) "Reversing parenthood..." ReverseParenthood: DO @; // Reverse the parenthood relations on part of the cycle climber_parent=parents[climber] @; // Save climber_orientation=orientations[climber] @; // Save this too climber_cardinality=cardinalities[climber] @; // Save parents[climber]=parent @; // Reverse the parenthood relations IF(orientation==head_is_parent) THEN @; // Reverse orientation of parent-relation orientations[climber]=tail_is_parent @; ELSE IF(orientation==tail_is_parent) THEN @; orientations[climber]=head_is_parent @; END IF @; cardinalities[climber]=subtree_size-cardinality @; // Update the cardinality IF(climber==subtree_root) @~ EXIT ReverseParenthood @; // Reached the root climber=heads_tails[climber_orientation, climber_parent] @; // Climb up the tree parent=climber_parent @; // Go one step up to the parent root orientation=climber_orientation @; cardinality=climber_cardinality @; END DO ReverseParenthood @; @%% WRITE(*,*) "Updating cardinalities..." climber=unhooking_node @; DO WHILE (climber!=cycle_root) @; // Update cardinalities cardinalities[climber]=cardinalities[climber]-subtree_size @; climber=heads_tails(orientations[climber], parents[climber]) @; // Move up END DO @; climber=hooking_node @; DO WHILE (climber!=cycle_root) @; // Update cardinalities cardinalities[climber]=cardinalities[climber]+subtree_size @; climber=heads_tails(orientations[climber], parents[climber]) @; // Move up END DO @; @%% EXIT TestPivot @; @%% CALL StopTimer(51) @; END IF PivotArc @; END DO TestPivot @; IF(debug_graph_algs) THEN @; WRITE(*,*) "In ReBuildSpanningTree, n_traces=", n_traces, & ", n_pivots=", n_pivots, ", percent=", @E REAL(n_pivots)/@E REAL(n_traces) @; @%% WRITE(*,*) "Tracing time:", ReadTimer(50)," Pivoting time:", ReadTimer(51) @; END IF @; IF(PRESENT(percent_pivoted)) @~ percent_pivoted=@E REAL(n_pivots)/@E REAL(n_traces) @; END SUBROUTINE ReBuildSpanningTree @; @*0 Preconditioning Laplacian Matrices via Incomplete MST Factorizations. For network optimization problems whose solution has a tree structure, such as linear min-cost problems, several incomplete preconditioners have been developed in the literature and are implemented here. These basically form a preconditioner of the form $M=F D_F F^T$ where $F$ is either a pure basis (called $B$ above), in which case we have the MST preconditioner of Vaidya/Vega, or a tree matrix with non-unit off-diagonal entries, such as one coming from the incomplete zero-fill $Q R$ factorization of $A D_A A^T$ of Resende and Pardalos, or the complete root-free Cholesky factorization of $B D_B B^T + \Lambda$ of Mehrotra and Wang, where $\Lambda$ is diagonal, usually the non-tree contribution to the diagonal of $A D_A A^T$. Here all $D$ matrices are diagonal conductivity matrices, although they could also be supplied as resistivities in some routines to avoid more expensive divisions. @*1 Building the Preconditioner Spanning Tree. This procedure can be used to initialize or rebuild the spanning tree used by the preconditioners. It does not contain much new material, but is mostly a wrapper for the other routines in this module.If the tree is built from scratch, then Quick Ranking is used to sort the weights of the arcs and build the MST, however, if this is not the case, path labels are created and used to select the list of arcs candidates for entry into the MST. By setting |tree_type=any_tree| an initial dummy MST can be created very quickly. @ @= 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 @; // Heads-tails LOGICAL(KIND=l_wp), DIMENSION(-arc_offset:), INTENT(IN), OPTIONAL :: arcs_mask @; // Arc mask REAL(KIND=r_wp), DIMENSION(-arc_offset:), INTENT(IN), OPTIONAL :: arcs_weights @; // Arc weights INTEGER(KIND=i_wp), DIMENSION(0:), INTENT(OUT) :: level_ordering @; // From root to nodes 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 @; // Max or Min LOGICAL(KIND=l_wp), DIMENSION(-arc_offset:), INTENT(INOUT) :: tree_mask @; LOGICAL, INTENT(IN), OPTIONAL :: from_scratch @; // Is this a brand new tree? (|.FALSE.|) LOGICAL, INTENT(IN), OPTIONAL :: defragment_ordering @; // Pack each forest in |level_ordering| or not? (|.FALSE.|) INTEGER(KIND=i_byte), INTENT(IN) :: tree_type @; // |min_cost_tree|, |max_cost_tree| or |any_tree| INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: arcs_ordering, candidate_arcs, & nodes_roots, nodes_ordering @; // Temporary storage arrays INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_special_arcs, n_arcs @; // Counters 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ 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 @; // Building an MST is a bit complicated because arcs need to be sorted OrderArcs: IF(initialize_tree) THEN @; // We need to order the arcs weights in this case rebuild_cardinalities=.TRUE. @; // Cardinalities are not maintained here 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 @; // The correctly numbered permutation IF(mask_arcs) THEN @; @%% WRITE(*,*) "Considering masked arcs only in sorted order" @; // Compress the arc list to contain only masked arcs 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 @; @%% WRITE(*,*) "Considering all arcs in sorted order" @; n_candidate_arcs=n_special_arcs+n_arcs+1 @; // Consider all arcs END IF @; @%% WRITE(*,*) "List of candidate arcs to enter spanning tree:", & @%% arcs_ordering[1:n_candidate_arcs] @; @%% WRITE(*,*) "Masked OK:", arcs_mask(arcs_ordering[1:n_candidate_arcs]) IF(tree_type==max_cost_tree) THEN @; // Consider arcs in reverse order 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 @; // Arcs are in right order--from small to large weights 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=-_SIZE(arcs_ordering, i_wp), mold=1_i_wp) @; DEALLOCATE(arcs_ordering, STAT=alloc_status) @; ELSE OrderArcs @; /* In this case we want to reoptimize a legal spanning forest. Masking the arcs is not supported, since arcs between trees can never enter the tree in this scheme (one must first rebuild a forest with the new mask and then reoptimize this forest in this case!): */ rebuild_cardinalities=.FALSE. @; // Cardinalities are not maintained here 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 @; // How many arcs are eligible to enter the MST DO arc=-n_special_arcs, n_arcs @; // Find eligible arcs for entry into tree head=heads_tails[1,arc] @; tail=heads_tails[2,arc] @; /* An eligible arc is a non-tree arc whose weight compares well with the path labels: */ 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 @; // The arc is eligible n_candidate_arcs=n_candidate_arcs+1 @; candidate_arcs[n_candidate_arcs]=arc @; END IF @; END IF @; END DO @; @%% WRITE(*,*) "List of candidate arcs to re-enter spanning tree:", & @%% candidate_arcs[1:n_candidate_arcs] @; 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=-_SIZE(candidate_arcs, i_wp), mold=1_i_wp) @; DEALLOCATE(candidate_arcs, STAT=alloc_status) @; END IF OrderArcs @; ELSE BuildMST @; // We just need a random spanning forest // I do not support random trees here, only |any_tree|, since sorting is expensive rebuild_cardinalities=.TRUE. @; // Cardinalities are not maintained here 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 @; @%% WRITE(*,*) "Considering masked arcs only" @; // Compress the arc list to contain only masked arcs in sequence 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 @; @%% WRITE(*,*) "Considering all arcs" @; 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 @; // Consider all arcs END IF @; @%% WRITE(*,*) "List of candidate arcs to re-enter spanning tree:", & @%% candidate_arcs[1:n_candidate_arcs] @; 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=-_SIZE(candidate_arcs, i_wp), mold=1_i_wp) @; DEALLOCATE(candidate_arcs, STAT=alloc_status) @; END IF BuildMST @; /* Now we can build the rest of the tree data-structures. In the case of a spanning forest, there may be a better ordering of the nodes onto levels, where each tree in the forest is stored in a contiguous section of |level_ordering|. So it may be wise to recompute the ordering here, for cache performance and such? :*/ ReorderNodes: IF(reorder_nodes) THEN @; // We will reorder the nodes to improve cache performance 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) @; @%% WRITE(*,*) "Initial nodes ordering:", nodes_ordering @; offset=0 @; // Start here DefragmentOrdering: DO index=0, n_nodes+n_special_nodes @; // Offset each root into level_ordering accordingly node=nodes_ordering[index] @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; // This is a root node @%% WRITE(*,*) "Root, cardinality, offset:", node, cardinalities[node], offset @; level_ordering[offset]=node @; parents[node]=offset @; // Save this offset--I am overwriting parents here! nodes_roots[node]=node @; // This is a root offset=offset+cardinalities[node] @; // Next root will go here CYCLE DefragmentOrdering @; // Go to next node END IF @; parent_node=heads_tails(orientation, parents[node]) @; // This is the parent node root_node=nodes_roots[parent_node] @; nodes_roots[node]=root_node @; // Save the root of this forest parents[root_node]=parents[root_node]+1 @; // One more child has been added level_ordering(parents[root_node])=node @; // Add this child to the ordering END DO DefragmentOrdering @; @%% WRITE(*,*) "New nodes ordering:", level_ordering @; @%% WRITE(*,*) "Nodes roots:", nodes_roots @; CALL RecordAllocation(n_elements=-_SIZE(nodes_ordering, i_wp), mold=1_i_wp) @; DEALLOCATE(nodes_ordering, STAT=alloc_status) @; CALL RecordAllocation(n_elements=-_SIZE(nodes_roots, i_wp), mold=1_i_wp) @; DEALLOCATE(nodes_roots, STAT=alloc_status) @; ELSE ReorderNodes @; // Use the default ordering 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 @; @*1 Incomplete $Q R$ preconditioner. This routine computes the Resende/Pardalos incomplete $Q R$ factorization of the matrix $A D_A A^T$ in the form $F D_F F^T$, where $F$ is given via the non-unit off-diagonal |nodes_mutlipliers| and the diagonal $D_F$ is given via |nodes_conductances|. The construction of this preconditioner is very simple, but in practice the incomplete Cholesky factorization of Mehrotra/Wang seems more efficient. This preconditioner approximates a diagonal preconditioner when the MST preconditioner $B D_B B^T$ is very bad, and approximates the MST preconditioner when it is good, so it is sort of an adaptive transition between the two. Similar arguments apply to the the Cholesky factorization, though less formally than the $Q R$ preconditioner. @ @= 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 @; // $y$ 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ IF(PRESENT(nodes_conductances)) THEN @; _FactorizeIncompleteQR(C, nodes_conductances) @; ELSE IF(PRESENT(nodes_resistances)) THEN @; _FactorizeIncompleteQR(R, nodes_resistances) @; nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) @; END IF @; END SUBROUTINE FactorizeIncompleteQR @; @ @m _FactorizeIncompleteQR(ID, nodes_conductances) @; nodes_conductances=0.0_r_wp @; Conductances_@e@&ID: DO arc=-n_special_arcs, n_arcs @; IF(.NOT.tree_mask[arc]) THEN @; // Tree arcs only contribute to the child node! 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_@e@&ID @; Multipliers_@e@&ID: DO node=-n_special_nodes, n_nodes @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; nodes_multipliers[node]=1.0_r_wp @; // Just a sentinel value CYCLE Multipliers_@e@&ID @; 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_@e@&ID @; @*1 Diagonal $L D L^T$ Cholesky preconditioner. This routine computes a root-free full $F D F^T$ Cholesky factorization of $B D_B B^T + \Lambda$ as described by Mehrotra and Wang (here $L=F$). The interesting thing is that this factorization does not require any fill so long as $\Lambda$ is diagonal. Here I take $\Lambda=Diag(N D_N N^T)$ to be the non-basic contribution to the diagonal of $A D_A A^T$, though Mehrotra/Wang propose some other choices as well without sufficient support. @ @= 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 @; // Not really needed, but added to make the same interface as above 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 @; // $y$ REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(OUT):: nodes_multipliers @; REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(IN), OPTIONAL :: & diagonal_nodes_conductances @; // $\Lambda$ if given explicitly 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ IF(PRESENT(nodes_conductances)) THEN @; _FactorizeIncompleteLDLt(C, nodes_conductances) @; ELSE IF(PRESENT(nodes_resistances)) THEN @; _FactorizeIncompleteLDLt(R, nodes_resistances) @; nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp)) @; END IF @; END SUBROUTINE FactorizeIncompleteLDLt @; @ @m _FactorizeIncompleteLDLt(ID, nodes_conductances) @; InitializeConductances_@e@&ID: IF(PRESENT(diagonal_nodes_conductances)) THEN @; // The diagonal term $\Lambda$ is given explicitly nodes_conductances=diagonal_nodes_conductances @; // Start with non-basic contribution AddDiagonal_@e@&ID: DO node=-n_special_nodes, n_nodes @; // Add the basic diagonal contribution orientation=orientations[node] @; IF(orientation==no_parent) @~ CYCLE AddDiagonal_@e@&ID @; nodes_conductances[node]=nodes_conductances[node]+arcs_conductances(parents[node]) @; END DO AddDiagonal_@e@&ID @; ELSE @; nodes_conductances=0.0_r_wp @; IF(PRESENT(arcs_mask)) THEN @; MaskedConductances_@e@&ID: 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_@e@&ID @; ELSE @; Conductances_@e@&ID: 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_@e@&ID @; END IF @; END IF InitializeConductances_@e@&ID @; Multipliers_@e@&ID: DO index=n_nodes+n_special_nodes, 0, -1 @; node=level_ordering[index] @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; @%% WRITE(*,*) node, nodes_resistances[node], nodes_multipliers[node] @; nodes_multipliers[node]=1.0_r_wp @; // Just a sentinel value CYCLE Multipliers_@e@&ID @; END IF @; parent_arc=parents[node] @; parent_node=heads_tails[orientation, parent_arc] @; // This is the parent node 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_@e@&ID @; @*1 Solving a Basic Dual Newton System. The following routine actually solves a system of the form $(F D_F F^T)x=b$, where $F$ is either a pure or non-pure basis of the network. Both |nodes_multipliers| for $F$ and the conductances/resistances in $D_F$ are optional. In fact, $D_F$ can be given as either arc or nodal conductances/resistances, via optional arguments (usually one should supply only 1 argument for $D_F$). The potentials of the root nodes should be pre-specified in |nodes_excess_potentials| because the solution can only be determined up to a constant displacement--the potential of the root node (under the assumption that the supply-demand vctor $b$ has zero sum, which it {\em must be} in order for the system to be feasible!). Here I solve the system by step-wise solving the three different systems with $F$, $D_F$ and $F^T$. I am not aware of a way to solve this system in one pass through the arcs, and I believe such a pass does not exist. @ @= 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 @; // Heads-tails array INTEGER(KIND=i_wp), DIMENSION(-node_offset:), INTENT(IN) :: parents @; // With displaced counting 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 @; // Temporary REAL(KIND=r_wp), DIMENSION(-node_offset:), INTENT(OUT), OPTIONAL :: tree_arcs_excess_voltages @; // Temporary 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ ArcOrNodeBased: IF(PRESENT(arcs_excess_voltages)) THEN @; // Arc-based temporary IF(PRESENT(nodes_multipliers)) THEN @; // Solve $F y =b$ using |arcs_excess_voltages| as a temporary for the tree flows 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) @; // Basic arc flows END IF @; /* Now solve $D_F t =y$ for $t$ using a different method depending on how $D_F$ is specified: */ 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 @; // Finally solve $F^T x=t$ for $x$ 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 @; // Node-based temporary 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 @; @%% WRITE(*,*) "Propagated flows:", tree_arcs_excess_voltages @; /* Now solve $D_F t =y$ for $t$ using a different method depending on how $D_F$ is specified: */ IF(PRESENT(arcs_conductances)) THEN @; DO node=-n_special_nodes, n_nodes @; IF(orientations[node]==no_parent) THEN @; // The virtual rooting arc has infinite conductance 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 @; // The virtual rooting arc has zero resistance 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 @; @%% WRITE(*,*) "Propagated voltages:", tree_arcs_excess_voltages @; IF(PRESENT(nodes_multipliers)) THEN @; // Finally solve $F^T x=t$ for $x$ 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 @; @%% WRITE(*,*) "Propagated potentials:", nodes_excess_potentials @; END IF ArcOrNodeBased @; END SUBROUTINE BasicDualNewtonSystem @; @*1 Calculating Effective Cycle Resistances. The following routine will sum the resistances of the arcs in the BEP cycle of each non-basic arc. This is an operation that can for example be used in preconditioning certain systems of equations: @ @= 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 @; // Heads-tails LOGICAL(KIND=l_wp), DIMENSION(-arc_offset:), INTENT(IN) :: tree_mask @; INTEGER(KIND=i_wp), DIMENSION(0:), INTENT(IN) :: level_ordering @; // Ordering from root to nodes 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 @; // Weights on arcs REAL(KIND=r_wp), DIMENSION(-arc_offset:), INTENT(OUT) :: cycles_resistances @; REAL(KIND=r_wp), DIMENSION(:), ALLOCATABLE :: nodes_resistances @; // Temporary 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ 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) @; @%% nodes_resistances[0]=0.0_r_wp @; PathResistances: DO index=0, n_nodes+n_special_nodes @; // Total resistance of path to root node=level_ordering[index] @; orientation=orientations[node] @; IF(orientation==no_parent) THEN @; nodes_resistances[node]=0.0_r_wp @; // Initialize the root path labels CYCLE PathResistances @; // Leave root potentials alone! END IF @; parent_arc=parents[node] @; parent_node=heads_tails[orientation, parent_arc] @; // This is the parent node 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 @; // For basic arcs there is no cycle cycles_resistances[arc]=0.0_r_wp @; ELSE @; // For non-basic arcs we need to trace the BEP cycle head=heads_tails[1,arc] @; tail=heads_tails[2,arc] @; climbers=(/head, tail/) @; // Start tracing the cycle formed by |arc| 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 to climb on side=1 @; // Climb on side 1 ELSE @; side=2 @; // Climb on side 2 END IF @; parent=parents(climbers[side]) @; // Candidate to leave the tree--an arc in the cycle orientation=orientations(climbers[side]) @; climbers[side]=heads_tails[orientation, parent] @; // Climb up one level on this side 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 @; @*0 Graph Format Conversion Routines. The following routine will take a heads and tails array and convert it into an adjacency-list-like representation, suitable for use with the SCOTCH library: @*1 Adjacency lists. @ @= 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 @; // The actual number of edges (without self-loops) LOGICAL, INTENT(IN), OPTIONAL :: include_loops @; // Should self loops be included or not (|.TRUE.|) INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: nodes_degrees @; // The nodal degrees INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_arcs, n_special_arcs, & n_vertices @; // Counters 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=_UBOUND(my_neighbours, i_wp, DIM=1)-2-n_special_nodes @; // The {\em first and last entry in |my_neighbours| are sentinel values!} n_special_arcs=-_LBOUND(heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND(heads_tails, i_wp, DIM=2) @; n_vertices=n_nodes+n_special_nodes+1 @; // For now I add all the vertices always! node_shift=n_special_nodes+1 @; // Shift in node numbering arc_shift=n_special_arcs+1 @; // Shift in arc numbering @%% WRITE(*,*) n_special_nodes, n_nodes, n_special_arcs, n_arcs, n_vertices, n_edges @; add_loops=.TRUE. @; // By default I include self loops IF(PRESENT(include_loops)) @~ add_loops=include_loops @; _AllocateNodalArray(nodes_degrees, 1_i_wp, .NOT.ALLOCATED, "CreateAdjacencyArrays") @; /* First we need to calculate nodal degrees: */ nodes_degrees=0 @; DO arc=-n_special_arcs, n_arcs @; // Note: Here we need to separately add to the head and tail due to loop 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 @; /* Now we do a |SUM_PREFIX| on the nodal degrees to figure out the offsets in |my_neighbours|: */ 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 @; // Sentinel value my_neighbours[n_vertices+1]=arc_counter @; // Sentinel end-of-array value IF(PRESENT(n_edges)) @~ n_edges=arc_counter-1 @; DO arc=n_arcs, -n_special_arcs, -1 @; // Now we traverse the arcs in reverse order to make the adjacency lists 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 @; _DeallocateArray(nodes_degrees, 1_i_wp, ALLOCATED) @; END SUBROUTINE CreateAdjacencyArrays @; @*1 Tree Partitioning. @ @= 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 @; // Heads-tails array 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 @; // Divide by conductances or not 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=-_LBOUND(parents, i_wp, DIM=1) @; n_nodes=_UBOUND(parents, i_wp, DIM=1) @; // The number of nodes, $n$ @; n_special_arcs=-_LBOUND (heads_tails, i_wp, DIM=2) @; n_arcs=_UBOUND (heads_tails, i_wp, DIM=2) @; // $m$ BEPJoints: DO arc=-n_special_arcs, n_arcs @; IF(tree_mask[arc]) THEN @; // For basic arcs there is no cycle joints_and_roots[arc]=0 @; ELSE @; // For non-basic arcs we need to trace the BEP cycle head=heads_tails[1,arc] @; tail=heads_tails[2,arc] @; climbers=(/head, tail/) @; // Start tracing the cycle formed by |arc| 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 to climb on side=1 @; // Climb on side 1 ELSE @; side=2 @; // Climb on side 2 END IF @; parent=parents(climbers[side]) @; // Candidate to leave the tree--an arc in the cycle orientation=orientations(climbers[side]) @; IF(orientation==no_parent) THEN @; // The arc is between two different trees cycle_root=0 @; EXIT TraceCycle @; END IF @; climbers[side]=heads_tails[orientation, parent] @; // Climb up one level on this side 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) @; // Calculate the congestions for tree arcs arcs_congestions_dilations= & arcs_congestions_dilations/(arcs_conductances+EPSILON(1.0_r_wp)) @; // We must take |ABS| later because orientation does not matter END WHERE @; END IF @; END IF @; /* HOW ABOUT THE ORIENTATIONS OF THE ARCS!?! */ 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 @; @%%