\section{Parallel Cost-scaling Algorithm} \label{scale} \subsection{Description of the Algorithm} \label{desc_cost_scale} Recently, we have designed Assignment algorithms for SIMD implementation, building upon theoretical advances in algorithms based on the PRAM model \cite{br88}. The technique is an adaptation of the sublinear time parallel PRAM algorithm of \cite{gpv88}, which uses a two-phase cost scaling approach. We obtain an $O(n^{5/3} \log^2 n C)$ time algorithm, where weights are integral and $C$ represents the weight of maximum magnitude (one extra factor of $\log n C$ is incurred by the cost scaling iterations). The best previous SIMD algorithms appear to have been $O(n^2 \log n C)$ on a bit-serial $n \times n$ mesh. This algorithm is fairly memory efficient, using only a small constant number of $\log n$-bit registers per processor. To adapt these ideas to the strict SIMD model, we employ a novel technique in which the data is permuted in each pass to make optimal use of the local communication. A step-by-step description of the algorithm is given below. The outer loop of the algorithm performs cost scaling, in which a solution is obtained through iterative improvement. Steps 1 through 4 in ``Match()'' (Figure~\ref{match_procedure}) begin with a solution which is 2 $ \epsilon$-optimal and improve it to an $ \epsilon$-optimal solution (see \cite{be86} for a discussion of the cost scaling technique). The algorithm has two major phases. The first is a very efficient "Push and Relabel" strategy which takes one step on many augmenting paths in parallel. The second is a more powerful but slower (in constant factor) routine which searches for complete disjoint augmenting paths. Each pass of the second phase finds all of the one-step paths of the first phase, some multiple-step paths, and at least one complete augmenting path. Thus, while each pass of the second phase takes longer, it may require fewer total passes. Theoretically, the first phase is more efficient at the beginning when there are many free vertices, and later the second phase becomes more efficient because it is guaranteed to increase the matching size by at least one during every pass. There is, therefore, an optimal ``cutoff'' point at which one should switch from phase 1 to 2. \begin{figure} \begin{tabbing} \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \kill Match( ) \\ \{ \\ \> $\forall i \forall j : c(i, j) \leftarrow n \times w(i, j)$ \\ \> $\epsilon \leftarrow C$ \\ \> {\bf while} ($ \epsilon > 1$){\bf do } \\ \> \> Step 1: \> \> \> Find the minimum unmatched edge cost $C_j$ in each column $j$. \\ \> \> Step 2: \> \> \> Subtract $C_j$ from all $c(i, j)$ in each column $j$. \\ \> \> Step 3: \> \> \> Unmatch arcs $(i, j)$ with $c(i, j) \geq \epsilon$. \\ \> \> Step 4: \> \> \> {\bf while} there are unmatched vertices (number % of matches $< n$) {\bf do} \\ \> \> \> \> \> \> {\bf if} ( number of matches $<$ Crossover ) \\ \> \> \> \> \> \> \> CALL Push\_Relabel\_Matching( ); \\ \> \> \> \> \> \> {\bf else} \\ \> \> \> \> \> \> \> CALL Shortest\_Path\_Matching( ); \\ \> \> \> \> \> \> {\bf endif} \\ \> \> \> \> \> {\bf endwhile} \\ \> \> $\epsilon \leftarrow \epsilon/2$; \\ \> {\bf endwhile} \\ \} \\ \end{tabbing} \caption{Procedure Match()} \label{match_procedure} \end{figure} Let us describe the computation required in each phase. Phase 1 is described below as ``Push\_Relabel\_Matching'' (Figure~\ref{push_relabel_matching}). Phase 1 begins by marking eligible edges from free $S$-vertices, using simple row and column operations. Next, a maximal matching of eligible edges is computed, by sweeping the matrix diagonally from the top left to the bottom right. (This subroutine is not much more time consuming than the row/column operations.) Finally, after selecting new matched edges and unmatching the old matched edges with which they conflict, we update the dual variables. This consists of increasing the costs in each column until the column's matched edge has cost $= \epsilon$. We then decrease the costs in each row until the cost of some unmatched edge in the row is $= - \epsilon$. Again, these require only simple row and column operations. Each step is $O(n)$, times the bit length of the operation, with very small constants. \begin{figure} \begin{tabbing} \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \kill Push\_Relabel\_Matching( ) \\ \{ \\ \> Step 4.1: \> \> \> $matched\_cost(i, j) \leftarrow$ column broadcast of the cost of the matched edge. \\ \> \> \> \> (If the column is unmatched, broadcast $\epsilon$.) \\ \> Step 4.2: \> \> \> For all unmatched edges, mark eligible if % $c(i, j) < matched\_cost(i, j)$. \\ \> Step 4.3: \> \> \> Find a maximal matching of the eligible edges. \\ \> Step 4.4: \> \> \> Column broadcast the newly matched edges $new\_M$. \\ \> Step 4.5: \> \> \> Unmatch the old edges from newly matched columns. \\ \> \> \> \> Match the new edges in these columns. \\ \> Step 4.6: \> \> \> Modify the column costs (increase them as high as % possible without making the \\ \> \> \> \> matched edge in the column have $c(i, j) > \epsilon$). Column broadcast the cost of the \\ \> \> \> \> newly matched edges. \\ \> Step 4.7: \> \> \> Modify the row costs (decrease them as low as possible % without making the \\ \> \> \> \> unmatched edge in the row have $c(i, j) <\epsilon$).\\ \} \\ \end{tabbing} \caption{Procedure Push\_Relabel\_Matching( )} \label{push_relabel_matching} \end{figure} The second phase, called ``Shortest\_Path\_Matching'' (Figure~\ref{shortest_path_matching}), works much harder at the dual variable update stage, finding all feasible augmenting paths. In place of the maximal matching step, this allows us to take as many steps as possible along augmenting paths from the free $S$-vertices. Furthermore, by reordering the rows and columns, we can find these augmenting paths in one diagonal sweep of the matrix, in time not much greater than a simple maximal matching. The time-consuming portion of phase 2 is in finding all feasible augmenting paths. The heart of the computation is the calculation of the shortest path from every $S$-vertex to some free $T$-vertex. Theoretically, by carefully pipelining the operations, such a distance computation can be done in $O(n \log n C)$ steps on a bit-serial $n \times n$ array \cite{va80}, and our implementation achieves this bound. Therefore, phase 2 has the same asymptotic complexity as phase 1. However, the constant for the distance computation is nearly two orders of magnitude greater than that of the row/column operations, which limits its competitiveness. \begin{figure} \begin{tabbing} \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \hspace*{.25in} \= \kill Shortest\_Path\_Matching( ) \\ \{ \\ \> Step 4.1: \> \> \> Compute the distance labels (shortest paths to free $T$-vertices). \\ \> Step 4.2: \> \> \> Permute rows and columns to order $S$-vertices and $T$-vertices by \\ \> \> \> \> distance label, with the largest distance labels on top/left. \\ \> Step 4.3: \> \> \> Mark eligible edges (those which appear in shortest augmenting paths). \\ \> Step 4.4: \> \> \> Set reduced costs $cb(i, j)$ as: \\ \> \> \> \> $cb(i, j) \leftarrow c(i, j) - d_s(i) + d_t(j)$ where arcs $(i, j)$ are unmatched. \\ \> \> \> \> $cb(i, j) \leftarrow c(i, j) - d_t(j) + d_s(i)$ where arcs $(j, i)$ are matched. \\ \> Step 4.5: \> \> \> Find and augment along maximal $S$-paths, starting at the largest distance labels. \\ \> Step 4.6: \> \> \> Set $c(i, j) \leftarrow cb(i, j).$ \\ \} \\ \end{tabbing} \caption{Procedure Shortest\_Path\_Matching( )} \label{shortest_path_matching} \end{figure} \subsection{Performance Results of the Parallel Cost-Scaling Algorithm} \label{results_scale} We have implemented the above algorithm on the Wavetracer, and report some results (Tables~\ref{fig_results_scale}, ~\ref{fig_distance_scale}, ~\ref{fig_combine_scale} and ~\ref{fig_xxx_scale}) for each of the two major subroutines (in msec) and for the entire algorithm (in sec) below. These results are for random data in the range $0 - 1000$. \begin{table} \begin{center} \begin{tabular}{|c|c|} \hline \hline {\em Size } & {\em Time (in sec)} \\ \hline $8 \times 8$ & 0.02 \\ \hline $16 \times 16$ & 0.07 \\ \hline $32 \times 32$ & 0.31 \\ \hline $64 \times 64$ & 1.11 \\ \hline $128 \times 128$& 6.37 \\ \hline \hline \end{tabular} \end{center} \caption{Push Relabel Matching Algorithm ($w(i, j) = i \times j$)} \label{fig_results_scale} \end{table} \begin{table} \begin{center} \begin{tabular}{|c|c|} \hline \hline {\em Size } & {\em Time (in sec)} \\ \hline $8 \times 8$ & 0.06 \\ \hline $16 \times 16$ & 0.28 \\ \hline $32 \times 32$ & 1.06 \\ \hline $64 \times 64$ & 4.34 \\ \hline $128 \times 128$& 27.8 \\ \hline \hline \end{tabular} \end{center} \caption{Shortest Path Matching Algorithm ($w(i, j) = i \times j$)} \label{fig_distance_scale} \end{table} \begin{table} \begin{center} \begin{tabular}{|c|c|} \hline \hline $128 \times 128$ & {\em Time (in sec)} \\ \hline Crossover at 115 & 7.3 \\ \hline Crossover at 121 & 6.1 \\ \hline Crossover at 125 & 3.9 \\ \hline Crossover at 126 & 3.7 \\ \hline Crossover at 127 & 3.4 \\ \hline No Crossover & 6.4 \\ \hline \hline \end{tabular} \end{center} \caption{Combination of Shortest Path and Push Relabel Matching Algorithm ($w(i, j) = i \times j$)} \label{fig_combine_scale} \end{table} \begin{table} \begin{center} \begin{tabular}{|c|c|c|c|} \hline \hline {\em Size } & {\em Push\_Relabel cycles} & {\em Distance cycles} & {\em Algorithm Cycles} \\ \hline - & (in msec) & (in msec) & (in sec) \\ \hline $256 \times 256$ & 1.9 & 93 & 1.12 \\ \hline $512 \times 512$ & 3.5 & 19 & 3.33 \\ \hline $1024 \times 1024$ & 6.9 & 404 & 13.1 \\ \hline $2048 \times 2048$ & 13.8 & 843 & 34.9 \\ \hline \hline \end{tabular} \end{center} \caption{Average Performance of the Push Relabel and Shortest Path passes, and of the entire algorithm for large input size (random 11-bit data)} \label{fig_xxx_scale} \end{table} Note that the second phase (distance) is about 50 to 60 times slower than the first. For the random data tested, the average time was faster when using only phase 1 to completion, rather than invoking phase 2. It appears that for some examples it may be advantageous to use the second phase, only for the very last matched edge (i.e., for one pass, to avoid a potentially very large number of passes of phase 1). One particularly difficult problem is obtained by setting each cost $c(i,j)$ to $i \times j$. For this problem, we find that a 256-node instance requires 333 sec using ``Push\_Relabel\_Matching'' alone, and only 156 sec when crossing over to full distance calculation when all but one pair of nodes are matched. Crossing over sooner degrades this improvement, however. The speed of these ``theoretically good'' algorithms does not appear competitive with the Munkres implementation of the previous section. However, it must be noted that we no longer have the use of row and column busses, providing row and column OR in one cycle. The Push Relabel pass would be much improved by row/column broadcast, since row/column minimum and maximum are heavily used. For the 16K processor, machine, this would mean a factor of 128 improvement in row and column operations. The shortest path matching algorithm could probably also be improved by using the busses. The majority of the time is spent in the distance computation, which, as written, uses only local communication. However, with broadcast we could write a much simpler non-pipelined algorithm with perhaps an order of magnitude improvement in the constant. A second source slowdown is in the cost scaling. To simplify implementation, we used only a scale factor of 2. Therefore, a great many cost scaling passes were required, often only to bring a few new edges into the lower epsilon range. For instance, consider a 1K by 1K problem with data in the range 0 to 1000. The number of passes with scale factor 2 would be about 20, while increasing the scale to just 10 would result in 6 passes. Finally, one critical difficulty with these algorithms is that they require more storage per processor than the Munkres algorithm. This limits the size of problem which can be efficienty run on the parallel array, and is a reason why we have added the Auction algorithm to our implementations, as discussed in the next section.