## 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-pivot 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.
flag is a variable used to check if the current tree satisfies the optimality conditions. We will explain this more in the procedure update_pi.

The following is pseudo code for the nonscaling premultiplier algorithm:

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
end

Now we will explain the details of each procedure.

### Initialization

Procedure Initialization
begin
read network files, get all nodes & arcs data;
compute C'=10 (|C|+1) where C=max{cij | (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 pi=2C'+1 for supply nodes; pi=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;
end

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

This is the procedure to find an eligible arc as the entering arc.
Procedure find_in_arc
begin
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;
end

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

In this procedure, we let the tail of the entering arc be the root.

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;
end

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.

Figure 3.1: An example for the procedure update_tree1

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

In this procedure, we walk along W starting from the root and identify the blocking arc. The original network simplex algorithm also uses the same procedure to identify 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 ruv = 0 then return (u,v), and leave this procedure;
else
let d = min{ d, ruv }; record (u,v) if old d>ruv;
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(cpijd), along each (u,v) on W;
end
return (u,v);
end

Update Tree Structure: The Second Part

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

Figure 3.2: An example for the procedure update_tree2

### 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 pi ,for each i belongs to E;
reset CurrentArc(i)=FirstArc(i),for each i belongs to E;
return flag=1, leave this procedure;
end

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

### 3.1.2 More about our implementing

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 pi has not been updated. That is, during the two iterations of updating pi , 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 pi , we will focus our following discussion on any iteration between the two iterations of updating pi. 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 cpij=cij-pi+pj and pi remains the same during our discussion, cpij 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 cpij>= 0, uij>= rij>= 0;.
2. nontree arcs with cpij < 0, rij=0
3. tree arcs, either (a) point upward with cpij =0 or (b) point downward with cpij>= 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 cpij 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, cpji<0 and uji>=0, rji>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 pivoted 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 cpij will not decrease, (i,j) will remain in group 2 if pj does not increase too much. If pj increases much enough, cpij 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 pivoted out to be a nontree arc, cpij 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

The following is pseudo code of our code:

Algorithm generic scaling premultiplier begin

initialization
compute e=max{cpij | (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
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 pivoting 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

This is the procedure to find an eligible arc as 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 D1= 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 D1=min{ D1, absolute value of the reduced cost of pred_arc(j) };
end
end
return -1; leave this procedure;
end
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.

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

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

Note that if D=D2, 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 pi (mod e/4) for D2. It is easier for us to compute the modulus and make the premultipliers multiple of e/4 if pi is positive.

We use CurrentArc(i) as in the nonscaling algorithm. It can be shown that when N* is empty, cpij>= -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 D2 to update premultipliers, and may scan the same eligible nodes set often without finding any flowable admissible arc until D1 << D2(i) for all eligible nodes i. This is unprofitable since we will waste too much time in updating D2 but finally we still need D1. Although we use method=2 to avoid some computations in obtaining the new D2 , 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 D2.

### 3.2.2 Pseudo Code for the Modified Scaling Version

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 D1 = M then return flag=0, leave this procedure;
set D2(i)=e/4 - pi (mod e/4)for each i belongs to E ;
set tc(i)=min{cpij | (i,j) does not belong to T(v), rij>0 } or min{-cpji | (j,i) does not belong to T(v), rji>0 },for each i belongs to E;
reset D2(i)=old D2(i) + kie/4 such that ki is the minimum nonnegative integer making tc(i) - D2<-e/4, i belongs to E;
get D2=min{D2(i),for each i belongs to E};
if D1 > D2 then set D = D2 , method=2;
else set D = D1 , method=1;
for each j belongs to E, increase pj by D , if pj becomes multiple of e/4, we put node j into R, mark j to be awake;
return flag=1;
end

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 cost5, say, (i,j). For node i, (i,j) is the best candidate to be admissible after we update pi later. In other words, if we increase pi 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 D2(i) that is bigger than tc(i) for each eligible node i. And choose the smallest D2(i) to be D2. Then, we can guarantee to find an entering arc (i,j) when we scan node i if D2<= D1.

Therefore, after update_pi, we will either find an admissible arc if D2<= D1 or increase the number of eligible nodes if D1 < D2. Basically, we can say this modified version combines those unprofitable iterations of updating D2 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 D2(i) whenever we find D2(i) < D1.

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.

5.without loss of generality, here we only consider outgoing arcs (i,j) for node i.