A Separator-Based Framework for Automated Partitioning and Mapping of Parallel Algorithms in Scientific Computing

Guy Blelloch^{1}, Anja Feldmann^{1}, Omar Ghattas^{2}, John Gilbert^{3}, Gary Miller^{1}, David O’Hallaron^{1},

Eric J. Schwabe^{1}, Jonathan Shewchuk^{1}, Shang-Hua Teng^{3}

1School of Computer Science ^{∗} 2Department of Civil Engineering

Carnegie Mellon University Pittsburgh, PA 15213

3Xerox PARC 3333 Coyote Hill Road Palo Alto, CA 94304

Abstract

This paper is a report on ongoing work in developing automated systems for the partitioning, placement, and routing of data that is necessary for the efficient parallel solution of large unstructured problems in the area of scientific computing, specifically the numerical solution of partial differential equations. Many of these problems have as an iterated inner loop the multiplication of a vector of values by a large sparse matrix. This problem of sparse matrix-vector multiplication has an underlying combinatorial graph structure that can be exploited. Using geometric information from the original problem, we can partition this combinatorial graph using provably good two- or threedimensional graph separators (depending on the dimension of the problem). The resulting partitions into subproblems have good load balancing properties and a relatively small amount of communication between subproblems. In order to develop effective heuristics for the placement of these subproblems on the available processors and the routing of messages between them, we must also carefully consider the characteristics of the target architectures. The first parallel machine we are considering is the iWarp system. The novel communication mechanism of the iWarp system allows us to draw an analogy between our placement and routing problem and certain area minimization problems in the field of VLSI circuit layout, giving us an additional collection of insights and heuristics which can be brought to bear on our problem.

1 Introduction

There is widespread interest in the engineering and scientific communities in using parallel computers to solve larger problems than were previously feasible. Many problems in computational science and engineering can be modeled by partial differential equations. Numerical approximation of these equations over the spatial domain of interest by finite difference, finite element, or finite volume methods gives rise to algebraic systems with large, sparse coefficient matrices. Often, however, the distribution of nonzeroes is unstructured, as in problems with geometric complexity or differing scales of resolution (typical of problems in solid and fluid mechanics). In these cases, achieving maximum performance on distributed memory parallel computers (which generally are limited in both memory and communication resources) requires careful partitioning of problems along with efficient placement of and routing between the resulting subproblems. Such domain-based problems have an underlying combinatorial graph with useful properties derived from its geometric structure, among them the existence of small separators which can be used to partition the problem. The placement of subproblems on the processors of a parallel computer and the routing between them, however, require more detailed knowledge of the target machine if maximum performance is to be obtained.

∗Acknowledge everyone’s grants and contracts here. This research was supported in part by the Defense Advanced Research Projects Agency (DOD) monitored by DARPA/CMO under Contract MDA972-90-C-0035, the Air Force Office of Scientific Research under Contract F49620-92-J-0131, and the National Science Foundation under Grant CCR– 8713489.

Multiplying a vector by a sparse matrix is a frequently-used inner loop in many PDE solvers, such as Krylov subspace methods (e.g. conjugate gradients) and other iterative methods. There is a natural combinatorial graph representation of this problem which reduces the communication requirements of this multiplication to exchanging data across all edges of the graph. The unstructured meshes associated with finite element calculations have underlying graphs with embeddings in two or three dimensions already defined, which will in general have small separators that may be used recursively to partition these graphs into as many subgraphs as desired. The separator properties of the underlying graph in two or three dimensions will give bounds on the load balancing and communication requirements of the resulting partition.

As was stated earlier, effective performance in placement and routing requires careful attention to the properties of the target machines. Our first target is the iWarp system at Carnegie Mellon [1, 2], whose novel communication mechanism adds some interesting elements to the problem of placement and routing. Its communication mechanism, which uses circuit-switched logical pathways multiplexed over each physical channel, changes the measures we must consider when choosing a set of routing paths. In particular, the maximum congestion of a set of routing paths is of far greater importance than its maximum dilation. This allows us to draw a strong analogy between the problem of placement and routing and established work in the area of circuit layout [4, 5], and bring insights and techniques from this area to bear on our problem. In addition, geometric information from the original problem is helpful in the placement of subproblems. Preliminary experiments suggest that even relatively simple heuristic approaches to placement and routing for these graphs that take advantage of their geometric structure perform nearly twice as well as the current placement and routing packages on iWarp [7] in terms of the quantity of routing resources required to route a given underlying graph.

2 A Combinatorial Graph Representation of Sparse Matrix-Vector Multiplication

The problem of multiplying a vector of elements by a large sparse matrix is ubiquitous in scientific computation. Together with inner products and element-wise vector operations (e.g., addition, multiplication by scalars), it is one of the basic primitives used in iterative numerical methods, and is among the most communicationintensive. For example, the basic iterative step of the conjugate-gradient method for finite-element calculations consists of several inner-product calculations followed by one multiplication of a sparse matrix by a vector (and some simpler element-wise operations). It is important to note at this point that although the nonzero entries of the matrix may change from iteration to iteration, the nonzero structure of the matrix is fixed.

We now consider the communication requirements for calculating the product y of a sparse symmetric n × n matrix A and an n-element vector x. Assume that there are n processors, and each processor i stores a single element x_{i} of the vector. In order to calculate y_{i}, processor i must receive the values x_{j} from each processor j such that A_{ij} ≠ 0. The value x_{j} must be multiplied by A_{ij}, which can be thought of as an edge weight, at either the source or the destination. Thus the set of routing paths which we must implement on the network of processors includes a bidirectional edge between processors i and j whenever A_{ij} ≠ 0. This is just the undirected graph on n nodes whose adjacency matrix is A. This reduces the problem of calculating the product of a sparse matrix and a vector to that of exchanging data across all edges of a graph.

Given this graph, the problem of performing an inner product is relatively simple, assuming that the processor containing each x_{i} also contains the ith element of the other factor in the inner product. Then all that is required is for the ith processor to calculate x_{i}y_{i}, and then for the resulting values to be accumulated over the n processors using a spanning tree in the graph just described. Element-wise operations on the vectors require no communication, but distribution of scalars for scalar-vector multiplication will require a broadcast, which can use the same spanning tree. These two primitives, data exchange and accumulation/broadcast, are sufficient to implement most iterative methods. Since the data exchange graph contains the spanning tree for the accumulation/broadcast operation, we need only consider the implementation of the communication needed for the sparse matrixvector multiplication.

3 Graph Partitioning Using Separators

3.1 Partitioning  of  Sparse

Matrix-Vector Multiplication

In general, the sparse matrix-vector multiplications used in scientific computing can have tens of thousands or even hundreds of thousands of variables, which is more than the number of processors in current parallel machines. Therefore we must partition the problem and map the subproblems to the available processors, while allowing for the necessary communication between subproblems.

The most natural way to partition a graph is to remove some subset of its edges, leaving some connected components of vertices remaining. If this method of partitioning is used, and the connected components are assigned to processors, then the local computation in each processor consists of for each node in the component accumulating the values of its neighbors that are in the same component. To complete the computation, each node with a neighbor in another component must send its value to the processor containing that neighbor. After these two phases, the multiplication can be completed.

Current finite-element implementations suggest a different partitioning approach. Since the weight along each edge, which corresponds to a nonzero entry in the matrix, depends on all corners of the element to which that edge belongs (e.g., in two dimensions, all nodes on the face containing that edge), cutting edges so that the nodes of a given face are divided among several processors will mean additional communication is necessary to calculate the nonzero entries of the matrix. It is better to partition the graph by cutting between elements and putting copies of the nodes and edges along the cut in both bordering blocks of the partition, so as to keep entire elements on single processors. This will allow the computation of matrix entries to be performed locally, and although it will require some extra memory, the borders of the partition blocks will be shown to be small so this is not a problem. In this case, the multiple copies of a node will each be able to collect some partial sum of their entry in the product, and only these partial sums need to be communicated. Therefore for this partitioning scheme, the communication is between the multiple copies of each node, rather than between neighboring nodes.

In order to get an efficient implementation there are two measures we must keep in mind. First, we must keep the sizes of the partition blocks relatively balanced — note that for this problem, the amount of local work done on a processor is the number of edges contained in its partition block, not the number of nodes. Second, we want to keep the boundaries of the partition blocks as small as possible, as this corresponds to the amount of interprocessor communication we must perform. In the next subsection we will discuss how graph separators allow us to satisfy both of these concerns.

If we have partitioned the original graph into P blocks, where P is the number of available processors, then we can define the P -node minor derived from the original graph and the partition. The minor is a weighted graph with a supernode for each block of the partition, and an edge between two supernodes if their corresponding blocks need to communicate with each other (i.e., they either have a cut edge between them or contain copies of the same node, depending on the partitioning approach). The weight on a edge is the number of values which must be communicated between those two subproblems in order to perform a multiplication. This weighted minor on the P supernodes is what we must place and route on the target machine.

3.2 Implementation of Separators in Two and Three Dimensions

Miller, Teng, Thurston and Vavasis [6] give provably good separators in two or three dimensions for the sort of combinatorial graphs generated by these problems. These separators can be used recursively to obtain a partition of the graph with both good balance among block sizes (corresponding to processor load) and a relatively small number of cut edges or duplicated nodes (corresponding to data which must be routed), depending on the method of partitioning. This approach to graph partitioning has been implemented in MATLAB.

(John and Shang-Hua, could you put in a couple of paragraphs here with a little more information about the separator algorithm and implementation?)

4 Placement and Routing

4.1 Background

Once we have partitioned the graph into subgraphs that correspond to subproblems to be solved locally, and identified which subproblems will have to communicate and how much information they will have to send to each other, we must embed the resulting P -node weighted graph into the underlying graph of the interconnection network of the target machine. To implement the exchange of data across each edge of the original graph needed to perform the sparse matrix-vector multiplication, we must communicate between each pair of adjacent nodes in the P -node graph an amount of data proportional to the weight of the edge between them.

For our purposes, an embedding of a weighted graph G into a graph H (as distinguished from an embedding of a graph in two- or three-dimensional space) is a mapping of the nodes of G to the nodes of H, corresponding to an assignment of subproblems to processors, together with a mapping of the edges of G to paths in H, corresponding to routing paths between processors containing subproblems which need to communicate. We are only considering one-to-one node mappings because our first target machine does not support multi-tasking and our partitioning approach already guarantees us load balancing. In general, the two most important measures of the efficiency of such an embedding are its dilation, which is the maximum length of the paths in H corresponding to edges in G, and its congestion, which is the maximum over all edges in H of the total weight of the edges in G whose corresponding paths contain that edge. (Since we are only considering one-to-one mappings, the load of our embeddings will always be one.)

The goal of keeping dilation and congestion low is to minimize the total time required for communication, so the relative importance of these measures for a particular machine is dependent on its communication mechanism and routing method (e.g., store-and-forward versus circuit-switched routing). When we consider different machines, we also expect to discover additional measures of the efficiency (or even feasibility) of embeddings that depend on specific features of their architectures.

4.2 The iWarp System

Our first target machine is the iWarp system. The iWarp system is based on a single chip VLSI processor [1, 2] developed jointly by Carnegie Mellon and Intel Corporation. The iWarp component contains a computation agent, and a communication agent for transferring data to and from other iWarp processors. Computing rates up to 20 Megaflops and 20 MIPS per cell are possible in the computation agent. The communication agent can sustain a bandwidth of 320 Megabytes/second per cell with a latency of 0.2 milliseconds per hop. iWarp systems are 2-dimensional tori of iWarp processors, ranging in size from 4 to 1024 processors. All implementations are done on a prototype 64-processor system at Carnegie Mellon that runs at half-speed.

iWarp processors communicate with neighboring processors over structures called logical channels, which can be chained together to form pathways [2]. Although the physical connection pattern of the iWarp is a 2-dimensional torus, logical channels and pathways allow for logical connection patterns such as rings, trees, or even hypercubes [8]. A network of logical channels can be created statically at compile time or created and destroyed dynamically at runtime. Each iWarp processor can support at most 20 logical channels incident to it at any point in time. (Note that a logical channel is a unidirectional communication medium.) Since pathways are handled completely by the communication agent the delay at each processor along such a pathway is very small, and data communicated thusly does not interrupt the computation agent.

From these characteristics of iWarp systems we can now discuss some measurements of the efficiency of embeddings. Dilation is of minor importance because of the pathways. Congestion is a very important measurement because all pathways sharing a certain physical connection are multiplexed over that link, and the total available bandwidth is limited. In addition to these traditional measurements we have to face the fact that each iWarp processor can only support up to 20 logical channels at any point in time. Therefore if we cannot limit the number of logical channels needed to be at most 20, then we have to use dynamic setup of connections (or a message passing system), which is slower than the static setup of connections. This makes the node congestion of our embeddings, defined for each node in the host graph H as the total number of paths in H that are incident to that node, of paramount importance since it determines which style of connection setup to use. (We are currently investigating whether there are more efficient ways to implement precomputed routing patterns which require more than 20 logical channels per node.)

Code for the direct solution of a threedimensional finite-element problem is currently running on the iWarp, but it was designed for one specific (and very regular) three-dimensional grid. The next step is to modify this code to handle more general and useful irregular grids in two and three dimensions.

4.3 Approaches to Placement and Routing

Because of the communication mechanism of the iWarp system, the congestion of our embeddings are more important than their dilation to the amount of time it will take to communicate the needed messages. This is analogous to some of the issues considered in Gate Array layout (see [5]), where the placement problem reduces to an assignment problem. In the Gate Array layout problem, each gate in the netlist of the circuit has to be assigned to a set of cells on the master that will implement the gate. In our placement problem, each subtask needs to be assigned to a processor. The routing problem between gates in the Gate Array layout problem is characterized by the fact that the routing channels on the master have a fixed width (often the same over the whole chip). This is in contrast to standard-cell and full-custom layout, where the channel widths can be chosen by the layout algorithm. Routing a connection in the layout problem is analogous to embedding a graph edge in the processor network. The limited number of logical channels is analogous to the fixed uniform channel width in the Gate Array layout problem. We use the analogy between these problems to bring some ideas from layout algorithms to our heuristics for embedding the P -node weighted graph of subproblems into the graph representing the interconnection network of the target machine. Placement corresponds to the mapping of nodes to nodes, and routing to the mapping of edges to paths.

For placement, it is clearly to our advantage to take into account as much of the geometric information given in the original problem as possible. After all, the first step in layout is often to come up with an embedding of the graph in the plane, and our original graph is given to us already embedded in two- or three-dimensional space. One simple approach to placing two-dimensional graphs in the grid is to simply halve the set of nodes repeatedly with alternating vertical and horizontal cuts, and mapping the 64 nodes to the torus in the natural way. On the other hand, we can take advantage of the structure of the graph in a deeper way by using the recursion tree from the partitioning algorithm to map the P nodes of the graph of partition blocks to the nodes of the torus. There are many other heuristics for this placement problem, including simulated annealing and distance minimization, which we continue to investigate.

There are also a variety of approaches which can be taken for routing the needed connections. The simplest are one- and two-bend routing schemes with or without randomization, but algorithms using weighted shortest paths to minimize congestion seem to work considerably better. By rerouting paths which pass through areas of high congestion and varying the node weights and cost functions in appropriate ways, routings with very good balance of congestion across the network can be obtained. Local optimization heuristics yield further improvements. We are currently experimenting extensively with these methods.

Preliminary experiments have shown that it is not likely that we will be able to embed even the weighted graphs which are derived from partitions of two-dimensional problems within the current numerical limits that would allow us to implement our communication as a single set of static paths.

This does not bode well at all for three-dimensional problems, where the average degree of the weighted graphs nearly doubles. However, in two dimensions it seems that in general we can divide up the communication into two sets of statically routable paths (as it turns out, simply putting one orientation of each bidirectional edge into each set will work); two sets may suffice for three dimensions as well, but this requires more study. iWarp researchers are currently studying the question of how efficiently one can switch between precomputed sets of routes such as these when the communication phases are tightly synchronized. Recent work on barrier synchronization on iWarp systems [3] may be helpful here.

4.4 The Connection Machine CM-5

The next target machine for this approach is the Connection Machine CM-5. The routing methods used on this machine make the choice of routing paths no longer a concern, but partitioning (and placement?) are still important.

(Guy and Anja, this section could use a paragraph or two about the CM-5 and where the implementation stands and where it is going. The last part might belong in Section 6 instead.)

5 A Step-by-Step Example

On the following page we illustrate the partitioning, placement, and routing of a unstructured finite-element mesh with roughly 5000 nodes on an 8 × 8 torus. The mesh was designed to calculate stresses in a cracked plate, and is hundreds of times more dense at the center (at the tip of the crack) than near the edges. The purpose of this small example is to show the behavior of our partitioning method on an unstructured mesh with a large variation in density.

The first figure illustrates the original unstructured mesh for the crack problem. Note that the mesh is so dense near the center that the figure is completely black. To its right is the same mesh after having been partitioned into 64 blocks be the removal. Due to density variations in the mesh, many of the partition blocks are indistinguishable. The third figure illustrates the 64-node derived from collapsing all the nodes in each block of the partition to a single supernode. Finally, the last figure shows a placement and routing of this graph on an 8 × 8 torus. Edges extending past the boundary of the grid wrap around to the other side.

6 Future Directions

Incorporate into finite-element program generator

C and NESL implementations of separator code

Other routing approaches on iWarp

Implementing only a statically routable subgraph, and embedding the full communication graph into it.

References

[1] S. Borkar et al. “iWarp: an Integrated Solution to High-Speed Parallel Computing”. In Proceedings of Supercomputing ’88.

[2] S. Borkar et al. “Supporting Systolic and Memory Communication in iWarp”. In Proceedings of the 17th Annual ACM International Symposium on Computer Architecture, 1990.

[3] A. Feldmann, T. Gross, D. O’Hallaron, T. Stricker. “Subset Barrier Synchronization on a Private Memory Parallel System.” To appear in Proceedngs of the 4th Annual ACM Symposium on Parallel Algorithms and Architectures, 1992.

[4] F. T. Leighton. Complexity Issues in VLSI. MIT Press, 1983.

[5] T. Lengauer. Combinatorial Algorithms for Circuit Layout. Publisher, 19??.

[6] G. Miller, S.-H. Teng, W. Thurston, S. Vavasis. “A Unified Geometric Approach to Graph Separators.” In Proceedings of the 32nd IEEE FOCS, 1991.

[7] D. R. O’Hallaron. “The ASSIGN Parallel Program Generator”. Technical Report CMU-CS-91-141.

[8] T. Stricker. “Supporting the Hypercube Programming Model on Mesh Architectures (A Fast Sorter for iWarp Tori).” To appear in Proceedngs of the 4th Annual ACM Symposium on Parallel Algorithms and Architectures, 1992.