In this chapter we will explain how we implement the three versions of the premultiplier algorithm introduced in Chapter 2.
Here we break up the procedure simplex-p_{i}vot into 4 steps:
Algorithm nonscaling premultiplier
begin
As we described in Chapter 2, we add a supernode as the initial root node, and construct artificial arcs between all nodes and the supernode. Each artificial arc has a big cost C', and big capacity M. We want to set C' sufficiently large so that no optimal solution will have an artificial arc. In theory, C' should be nC+1; however, such a large value of C' may cause overflow, and we determined that C'=10(C+1) was satisfactory for our tests.
Basically, we perform a depth-first-search in this procedure. Referring to the definition 2 in Chapter 2, we can find an eligible node by searching along the tree T(v) starting from the root if all the arcs in the path from that node to the root have zero reduced cost. Standing at an eligible node, we scan its adjacent nontree arcs to identify eligible arcs. Note that we do not scan all the adjacent nontree arcs of a node every time; instead, we only start scanning from CurrentArc(i). The reasons will be explained in subsection 3.1.2. If an eligible node is exhausted, that is, CurrentArc(i)=-1, we will skip scanning this node and go to scan the next eligible node. By using CurrentArc(i), we avoid many redundant scans and speed up the algorithm.
Starting from any eligible node, we keep searching further in its subtree. If we find an ineligible node, we will not search further in its subtree. We record the absolute value of the reduced cost of the arc connecting that ineligible node with its predecessor. If we have searched all the eligible nodes but have not found any eligible arc, we will set D to be the minimum absolute value of the reduced cost of those arcs connecting those ineligible nodes and their eligible predecessors. This D will be used for updating the premultipliers of eligible nodes in the procedure update_pi.
We define in_arc to be the flowable eligible arc in the residual network pointing outward from the eligible node we are scanning. This procedure will either return the eligible arc if there is any, or return -1 if no eligible arc is found. More detailed procedures about this pseudo code are included in Appendix A.1.
Procedure update_tree1
begin
Since we use the FirstChild and the RightSibling tree structure for this algorithm, this part will be a little complicated. To make node i the root, the arcs in the path from i to the root have to befor eachreversed". An example is illustrated in Figure 3.1.
Recall that the tail of in_arc is the eligible node we are scanning. In the end of this procedure, we set the predecessor of the root to be the head of in_arc so that we can easily traverse around the basic cycle W later. If we can push flow along in_arc, the cost will decrease. The orientation of the basic cycle W is the same as in_arc. The next procedure is to push flow along W, and identify the leaving arc, out_arc.
Procedure find_out_arc
begin
Now we have identified the leaving arc, out_arc. We set the tail of out_arc to be the new root. Since we have defined the predecessor of the current root to be the head of in_arc, it is fairly easy to update the new tree structure in which the tail of out_arc ia the new root.
Procedure update_tree2
begin
Here we define a new node pointer, LeftSibling(i), which points to the sibling of i immediately to the left. We can find LeftSibling(i) by traverse from FirstChild(pred(i)) to the right. We have to be careful about the case that in_arc and out_arc have the same node pairs. Figure 3.2 illustrates the tree before and after procedure update_tree2. Note that the subtrees of node 4 and 6 are easily maintained.
After we have searched all the eligible nodes but still did not find any eligible arc, we have to modify the premultipliers for all the eligible nodes.
Procedure update_pi
begin
The premultipliers will only be updated after we could not find any eligible arcs in the procedure find_in_arc. Therefore, the vector of reduced costs will also only be changed after the premultipliers have been changed. We also need to reset the CurrentArc(i) to be FirstArc(i) for every eligible node i. This is because we can not guarantee if there is no eligible arc between FirstArc(i) and CurrentArc(i) after we changed p_{i} . The reason will be explained later in subsection 3.1.2.
Remember that we initialize D to be M in the beginning of the procedure find_in_arc. D will remain M if we have scanned all nodes but still could not find any eligible arc. In such case, we have reached the optimality condition since all the tree arcs have zero reduced cost, and all the nontree arcs are either ineligible or eligible with zero residual capacity.
Here we will explain more about some tips used in our nonscaling code which speed up this algorithm.
Why we can skip scanning arcs between FirstArc(i) and CurrentArc(i)In this algorithm, we use CurrentArc(i) to be the starting pointer to scan for eligible arcs instead of FirstArc(i). This is because there will be no eligible arcs between FirstArc(i) and CurrentArc(i) as long as p_{i} has not been updated. That is, during the two iterations of updating p_{i} , we can skip scanning those arcs between FirstArc(i) and CurrentArc(i). Furthermore, we can even skip scanning this node and go to scan the next eligible node in E if CurrentArc(i)=-1.
Since we reset CurrentArc(i) to be FirstArc(i) after we update p_{i} , we will focus our following discussion on any iteration between the two iterations of updating p_{i}. Also, (i,j) and (j,i) correspond to the same arc in the original network in our discussion. However, we will only use the outgoing arc (i,j) to represent that arc when we scan node i.
Recall that FirstArc(i) is the first arc adjacent to node i. CurrentArc(i) points to the arc which was the entering arc when we scanned node i last time.
Since c^{p}_{ij}=c_{ij}-p_{i}+p_{j} and p_{i} remains the same during our discussion, c^{p}_{ij} will either increase or remain the same due to the nondecreasing property of premultipliers.
We use two ideas to speed up this algorithm: (1) using CurrentArc(i) (2) computing reduced cost only when necessary.
As described earlier in this subsection, we can skip scanning those arcs between FirstArc(i) and CurrentArc(i). This did save much running time in the procedure find_in_arc and still gives us the correct result.
We implemented two versions of scaling premultiplier algorithm. We call them thefor eachgeneric scaling version"(scaling 1) and thefor eachmodified scaling version"(scaling 2). They differ in the procedure update_pi. Both of them have the same procedures initialization, update_tree1, find_out_arc and update_tree2 as the nonscaling version described in section 3.1.
Algorithm generic scaling premultiplier begin
The main loop is for the scaling phase. Each main loop corresponds to a scaling phase with a new value of e which decreases by at least a factor of 2 in each loop. The algorithm terminates when e < 1/n. During each scaling phase, it performs similar procedures as the nonscaling algorithm but with different criteria in identifying eligible arcs. We also use some controllers such as method and mode to switch the p_{i}voting and updating policies so that we can save some CPU running time.
Since this algorithm uses the same procedures of initialization, update_tree1, find_out_arc, and update_tree2 as the nonscaling version did, here we will only state the procedure find_in_arc and update_pi.
Procedure find_in_arc
begin
This procedure is very similar to its nonscaling version. In fact, they are exactly the same except the way of identifying an entering arc. In the nonscaling version, we scan an inexhausted eligible node for a flowable nontree arc with reduced cost < 0. Similarly, here in this scaling version, we scan an eligible and awakened node for a flowable admissible arc whose reduced cost < -e/4 during the e scaling phase.
We will discuss about the controllerfor eachmethod" after we explain the procedure update_pi.
After we have searched all the eligible and awakened nodes but still did not find any flowable admissible arc, we have to modify the premultipliers for all the eligible nodes.
Procedure update_piNow we will explain the controllers "method" and "mode". Consider the case that in some iteration after we update premultipliers by D_{2} , we still can not find any admissible arc, and we have to call the procedure update_pi again. At this moment, all the eligible nodes remain the same, so, we can compute the new D_{2}(i) by subtracting the old D_{2}(i) with the previous D (which was the previous D_{2}). We use mode=1 to denote this situation. If mode=0, we have to compute all the D_{2}(i) by (e/4 - p_{i} (mod e/4)).
Note that if D=D_{2}, we only need to scan those reawakened nodes in the procedure find_in_arc. We set method=2 for such case so that we do not need to perform the depth first search on the tree. Otherwise, we set method=1 and perform the depth first search along the tree as did the nonscaling version.
Now we can explain why we set all the initial premultipliers to be positive(1, C'+1, and 2C'+1). It is because we need to compute p_{i} (mod e/4) for D_{2}. It is easier for us to compute the modulus and make the premultipliers multiple of e/4 if p_{i} is positive.
We use CurrentArc(i) as in the nonscaling algorithm. It can be shown that when N* is empty, c^{p}_{ij}>= -e/2 for each (i, j) belongs to G [2].
From our computational testing results, we found this scaling algorithm often uses D_{2} to update premultipliers, and may scan the same eligible nodes set often without finding any flowable admissible arc until D_{1} << D_{2}(i) for all eligible nodes i. This is unprofitable since we will waste too much time in updating D_{2} but finally we still need D_{1}. Although we use method=2 to avoid some computations in obtaining the new D_{2} , this does not save much time. Therefore, we change this updating policy and do another slightly different scaling version which will focus on how to find a flowable admissible arc earlier and avoid some unprofitable loops of updating D_{2}.
This modified scaling version only differs with the generic one in the procedure update_pi. The following is the detailed pseudo code for the procedure update_pi:
Procedure update_pi
begin
This procedure is a little more complicated than its generic version. For each eligible node i, the algorithm scans each of its adjacent flowable nontree arcs and records the one with smallest reduced cost^{5}, say, (i,j). For node i, (i,j) is the best candidate to be admissible after we update p_{i} later. In other words, if we increase p_{i} by tc(i), (i,j) will be the first arc whose reduced cost becomes zero. However, in order to maintain the polynomial running time, we need to reawaken node i. Thus we try to get the smallest D_{2}(i) that is bigger than tc(i) for each eligible node i. And choose the smallest D_{2}(i) to be D_{2}. Then, we can guarantee to find an entering arc (i,j) when we scan node i if D_{2}<= D_{1}.
Therefore, after update_pi, we will either find an admissible arc if D_{2}<= D_{1} or increase the number of eligible nodes if D_{1} < D_{2}. Basically, we can say this modified version combines those unprofitable iterations of updating D_{2} in the generic version into one iteration.
We use the same tips used in the generic version here. Also, for node i, we can stop scanning for D_{2}(i) whenever we find D_{2}(i) < D_{1}.
As for how to compute the reduced cost efficiently, it is not easy to analyze which way is better: (1) to compute the reduced cost for the arcs adjacent to eligible nodes after update_pi or (2) compute the reduced cost of an arc when it is necessary as the generic version did. It is possible that we may compute the reduced cost for the same arc many times if we use (2) . However, we still use (2) to compute the reduced cost in our code.