Chapter 3
Implementing the Premultiplier Algorithm

In this chapter we will explain how we implement the three versions of the premultiplier algorithm introduced in Chapter 2.

3.1 Nonscaling Premultiplier Algorithm

3.1.1 Pseudo Code

Here we break up the procedure simplex-p_ivot into 4 steps:

1. find_in_arc will find the entering arc in_arc; in_arc=-1 if we could not find any entering arc.
2. update_tree1 will change the spanning tree structure so that the tail of in_arc becomes the root
3. find_out_arc will add the in_arc into the tree; construct the basic cycle W; identify the leaving arc out_arc.
4. update_tree2 will construct the new spanning tree so that the tail of out_arc becomes the root.

Algorithm nonscaling premultiplier
begin

initialization;
in_arc = find_in_arc;
while flag=1 do
begin
while in_arc -1 do
begin

update_tree1;
out_arc = find_out_arc;
update_tree2;
in_arc = find_in_arc;
end
flag = update_pi;
if flag=1 then in_arc = find_in_arc;
end

Initialization

read network files, get all nodes & arcs data;
compute C'=10 (|C|+1) where C=max{c_ij | (i,j) belongs to A};
add a supernode as root, with =C'+1; for each supply and transshipment node i, construct an artificial arc to the supernode with flow=b(i); for each demand node j, construct an artificial arc from the supernode to it with flow=-b(j); all the artificial arcs have cost=C' and capacity=M, where M is a sufficiently big number; (999999999 is sufficiently big in our tests);
set p_i=2C'+1 for supply nodes; p_i=1 for demand nodes;
set CurrentArc(i)=FirstArc(i) for all nodes;
mark all nodes except supernode to be ineligible;
get initial tree data structure;
turn on flag to be 1;

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.

Finding the Entering Arc

let E denote the set of eligible nodes; initialize E as {root};
set D=M where M is a sufficient large number;
scan nontree arc adjacent to root starting from CurrentArc(root);
if found a flowable eligible arc then return it; leave this procedure;
starting from the root, perform a depth-first-search along the spanning tree T(v);
while there is still some eligible node which has not been exhausted do
begin
if found an exhausted eligible node then put it in E, skip it, keep searching;
else if we find an inexhausted eligible node i then
begin
scan each nontree arc adjacent to node i starting from CurrentArc(i);
if found a flowable eligible arc then return it; leave this procedure;
else put node i in E, mark it to be exhausted, and keep searching;
end
else
we find an ineligible node j;
set D=min{ D, absolute value of the reduced cost of pred_arc(j) };
end
end
return -1; leave this procedure;

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.

Update Tree Structure: The First Part

Procedure update_tree1
begin

if the tail i of in_arc is not the root then
begin
make i the new root;
update the tree structure;
end
set the predecessor of the new root to be the head of in_arc;
set FirstChild(head of in_arc) to be the new root;
construct the basic cycle W;

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.

Find the Leaving Arc

Procedure find_out_arc
begin

set v=pred(root), u=root, d=M, where M is a sufficiently large number;
while v is not the root do
begin

if r_uv = 0 then return (u,v), and leave this procedure;
else
let d = min{ d, r_uv }; record (u,v) if old d>r_uv;
set v=pred(v);
end
send d along W;
begin
if the direction of the tree arc is the same as W then
decrease its residual capacity by d;
else increase its residual capacity by d;
compute Sum(c^p_ijd), along each (u,v) on W;

end
return (u,v);

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

say, the in_arc is (u,v) where u is the current root; out_arc is (p, q) where p will become the new root; (note that pred(u)=v, and pred(p)=q now)
if u=q and v=p then
begin
if (u,v) and (p,q) correspond to the same arc in the original network then
the tree structure remains the same, reset pred(p)=-1, pred_arc(p)=-1, FirstChld(q)=RightSibling(p); leave this procedure;
else
make p the new root, reset pred(p)=-1, pred_arc(p)=-1;
if LeftSibling(p)=-1 then FirstChild(q)=RightSibling(p), RightSibling(p)=-1; leave this procedure;
else
RightSibling(LeftSibling(p))=RightSibling(p), RightSibling(p)=-1; leave this procedure;
end
else
begin
make p the new root, reset pred(p)=-1, pred_arc(p)=-1;
if LeftSibling(p)=-1 then FirstChild(q)=RightSibling(p), RightSibling(p)=-1; leave this procedure;
else
RightSibling(LeftSibling(p))=RightSibling(p), RightSibling(p)=-1; leave this procedure;
end

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.

Modify the Premultipliers

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

if D = M then return flag=0, leave this procedure;
else
add D to p_i ,for each i belongs to E;
reset CurrentArc(i)=FirstArc(i),for each i belongs to E;
return flag=1, leave this procedure;

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.

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.

Claim : If node i is not exhausted yet, that is, CurrentArc(i)<>-1, there will be no entering arc between FirstArc(i) and CurrentArc(i). If CurrentArc(i) = -1, we can skip i, and go to scan the next eligible node.
Proof:
We first consider the case that CurrentArc(i)<>-1. Assume there are arcs between FirstArc(i) and CurrentArc(i). Since CurrentArc(i) was the entering arc when we scanned node i last time, which means, those arcs between FirstArc(i) and CurrentArc(i) can be classified into the following three groups:
1. nontree arcs with c^p_ij>= 0, u_ij>= r_ij>= 0;.
2. nontree arcs with c^p_ij < 0, r_ij=0
3. tree arcs, either (a) point upward with c^p_ij =0 or (b) point downward with c^p_ij>= 0
Now we discuss each group of arcs respectively for the later iterations when we scan node i again.
The arcs of group 1 will fall into group 1 or 3(b) when we scan node i:
Since c^p_ij will not decrease, (i,j) will never become eligible when we scan node i. However, (j,i) may be eligible when we scan node j if node j is eligible, c^p_ji<0 and u_ji>=0, r_ji>0. In such case, (j,i) will become a tree arc, and belongs to group 3(b). If node j is not eligible or even if it is eligible but (j,i) is not p_ivoted in in later iteration, arc (i,j) will remain in group1. Therefore, arcs of group 1 won't become eligible.
The arcs of group 2 will fall into group 1 or 2 when we scan node i:
Since c^p_ij will not decrease, (i,j) will remain in group 2 if p_j does not increase too much. If p_j increases much enough, c^p_ij will become nonnegative and (i,j) will be in group 1. According to our discussion for the arcs in group 1, the arcs in group 2 may also fall into group 3(b) later.
The arcs of group 3 will fall into group 1 or 3 when we scan node i:
For some tree arc (i,j) in group 3, it will remain in group 3 if it is still a tree arc. If (i,j) is p_ivoted out to be a nontree arc, c^p_ij will be nonnegative and thus (i,j) may fall into group 1. Similarly, (i,j) may fall into group 3(b) again according to our discussion for the arcs in group 1.
Therefore, any arc belonging to those 3 groups will never become an entering arc during our discussion.
Now consider the case that CurrentArc(i) = -1. In this case, all the arcs adjacent to node i will not be eligible by using the same discussion above, thus we can skip this node and go to scan the next eligible node.

Tips to speed up the algorithm

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.

3.2 Scaling Premultiplier Algorithm

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.

---

3.2.1 Pseudo Code for the Generic Scaling Version

Algorithm generic scaling premultiplier begin

initialization
compute e=max{c^p_ij | (i, j) belongs to A };
in_arc = find_in_arc;
mark all nodes to be awake; N*=set of all nodes;
while e < 1/n do
begin
while N* is not empty do
begin

while in_arc -1 do
begin
update_tree1;
out_arc = find_out_arc;
update_tree2;
set method=1;
in_arc = find_in_arc;
end
flag = update_pi;
if flag=1 then
in_arc = find_in_arc;
if (in_arc=-1 and method=2) then set mode=1;
else set mode=0;
end
goto next scaling phase, get new e, reset N*=set of all nodes, mark all nodes awake, reset CurrentArc(i)=FirstArc(i)for each i belongs to N;
in_arc = find_in_arc;
set mode=0;
end

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_ivoting 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.

Finding the Entering Arc

Procedure find_in_arc
begin

if method=2 then
begin
while there is any eligible and awakened node i do
begin
scan each nontree arc adjacent to the root starting from CurrentArc(i);
if found a flowable eligible arc then return it; leave this procedure;
end
return -1; leave this procedure;
end
else method=1 then
begin
let E be the set of eligible nodes; empty E; put root into E; set D₁= M where M is a very big number;
if root is awake then
begin
scan each nontree arc adjacent to the root starting from CurrentArc(root);
if found a flowable eligible arc then return it; leave this procedure;
end
starting from the root, perform a depth-first-search along the spanning tree T(v);
while there is still some eligible and awakened node do
begin
if found an eligible and asleep node then put it in E, skip it, keep searching;
else if we find an un exhausted eligible and awaken node i then
begin
scan nontree arc adjacent to node i starting from CurrentArc(i);
if found a flowable eligible arc then return it; leave this procedure;
else put it in E, mark it to be asleep, take it out of N* if its premultiplier has not been updated during this scaling phase, keep searching;
end
else we find an ineligible node j then
begin
we do not search further for that subtree;
set D₁=min{ D₁, absolute value of the reduced cost of pred_arc(j) };
end
end
return -1; leave this procedure;
end

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.

Modify the Premultipliers

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.

if D₁ = M then return flag=0, leave this procedure;
if mode=0 then
begin
set D₂(i)=e/4 - p_i (mod e/4)for each i belongs to E;
let D₂=min{ D₂(i),for each i belongs to E };
end
else mode=1 then
begin
scan all eligible nodes,
decrease D₂(i) by previous D ,for each i belongs to E;
reset D₂(i)=e/4 if D₂(i)=0,for each i belongs to E;
let D₂=min{ D₂(i),for each i belongs to E };
decrease D₁ by previous D;
end
if D₁ > D₂ then set D = D₂ , method=2;
else set D = D₁ , method=1;
for each j belongs to E, increase p_j by D , if p_j becomes multiple of e/4, we put node j into R, mark j to be awake;
return flag=1;

Now we will explain the controllers "method" and "mode". Consider the case that in some iteration after we update premultipliers by D₂ , 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₂(i) by subtracting the old D₂(i) with the previous D (which was the previous D₂). We use mode=1 to denote this situation. If mode=0, we have to compute all the D₂(i) by (e/4 - p_i (mod e/4)).

Note that if D=D₂, 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₂. 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].

Tips to speed up the algorithm

From our computational testing results, we found this scaling algorithm often uses D₂ to update premultipliers, and may scan the same eligible nodes set often without finding any flowable admissible arc until D₁ << D₂(i) for all eligible nodes i. This is unprofitable since we will waste too much time in updating D₂ but finally we still need D₁. Although we use method=2 to avoid some computations in obtaining the new D₂ , 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₂.

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

if D₁ = M then return flag=0, leave this procedure;
set D₂(i)=e/4 - p_i (mod e/4)for each i belongs to E ;
set tc(i)=min{c^p_ij | (i,j) does not belong to T(v), r_ij>0 } or min{-c^p_ji | (j,i) does not belong to T(v), r_ji>0 },for each i belongs to E;
reset D₂(i)=old D₂(i) + kie/4 such that ki is the minimum nonnegative integer making tc(i) - D₂<-e/4, i belongs to E;
get D₂=min{D₂(i),for each i belongs to E};
if D₁ > D₂ then set D = D₂ , method=2;
else set D = D₁ , method=1;

for each j belongs to E, increase p_j by D , if p_j becomes multiple of e/4, we put node j into R, mark j to be awake;
return flag=1;

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⁵, 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₂(i) that is bigger than tc(i) for each eligible node i. And choose the smallest D₂(i) to be D₂. Then, we can guarantee to find an entering arc (i,j) when we scan node i if D₂<= D₁.

Therefore, after update_pi, we will either find an admissible arc if D₂<= D₁ or increase the number of eligible nodes if D₁ < D₂. Basically, we can say this modified version combines those unprofitable iterations of updating D₂ in the generic version into one iteration.

Tips to speed up the algorithm

We use the same tips used in the generic version here. Also, for node i, we can stop scanning for D₂(i) whenever we find D₂(i) < D₁.

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.

Previous: Chapter 1 & 2
Next: Go to Chapter 4
Go to Premultiplier C Codes