Combinatorial Algorithms for Computational Scienceand Engineering
Erik G. Boman1 , Doruk Bozdag2 , Umit V. Catalyurek2 ,Karen D. Devine1 , Assefaw H. Gebremedhin3 , Paul D. Hovland4 ,and Alex Pothen3
1Discrete Algorithms and Math Department, Sandia National Laboratories
Biomedical Informatics, and Electrical and Computer Engineering, Ohio State University
3
Department of Computer Sciences and Computing Research Institute, Purdue University
4
Mathematics and Computer Science Division, Argonne National Laboratory
2Abstract. We discuss recent activities in the SciDAC Applied Math Institute for
Combinatorial Scientific Computing and Petascale Simulations (CSCAPES).

1. Introduction
CSCAPES was established in 2006 to help
contribute to the SciDAC mission by developing infrastructural algorithmic andsoftware technologies for supporting highperformance computing. Its scope encompasses four areas: load-balancing and related tasks in parallel processing (such asdata migration and task scheduling); autoScope of CSCAPES
matic differentiation (AD), a technology fortransforming computer source code for computing a function into code for computing its Figure 1. Key research areas in CSCAPES andderivatives; basic ingredients in linear solver their relationship to a typical SciDAC application.
technology; and runtime data and iteration An arrow from A to B indicates that A in somereordering to improve performance in irreg- sense uses B.
ular computation.
The focus in CSCAPES in each of these areas is on the formulation (modeling) and solution(algorithms) of the underlying fundamental combinatorial problems, which are often phrased
using graphs or hypergraphs. Load-balancing, the task of distributing data and work in a largescale computation among the processors of a machine to minimize total execution time, can bephrased as a partitioning problem on a graph, or with more accuracy, on a hypergraph. Whencomputational dependency among subtasks is modeled using a graph, a distance-1 coloringof the graph can be used to determine a scheduling of the tasks for concurrent computation.
Distance-2 coloring and several other specialized variants are used to model matrix partitioning
problems that arise in the efficient computation of sparse Jacobians and Hessians using AD.
In AD, exploitation of associativity and commutativity in the chain rule of differential calculus
Computational Science Application

Automatic Differentiation

Data Migration

Elimination

Optimization

Task Scheduling

Coloring

Scientific Computing Tool

Linear Solver

Load Balancing

Matching

HPC Task

Eigen Solver

Performance Improvement

Partitioning

Combinatorial Problem

Orderingleads to a variety of vertex or edge elimination problems in the directed acyclic graph used torepresent the computation of a function and its derivative. Various types of matchings in graphs
are used in the solution of sparse linear systems (numerical preprocessing and block triangulardecomposition) and in the coarsening phase of multilevel methods for graph partitioning. Workand storage reduction in direct solvers for large sparse systems calls for vertex ordering problems
in graphs. Graph and hypergraph based models are used to capture the needs in reordering thenodes and elements within an unstructured mesh to improve runtime performance.
CSCAPES researchers are developing sequential and parallel algorithms and software forthe combinatorial problems outlined in the previous paragraph. The load-balancing andparallelization capabilities are being deployed through the Sandia-housed Zoltan software toolkit.
Similarly, the AD capabilities are being deployed via the tools being developed at Argonne within
the OpenAD framework. In addition, serial as well as parallel codes (for coloring, matching, andre-ordering) are being separately developed at the academic institutions within CSCAPES, withthe ultimate goal of being integrated with the major software outlets, Zoltan and OpenAD.
Both Zoltan and the tools in OpenAD pre-date CSCAPES and have been under continuous
development for several years with sustained support from DOE; similarly, the coloring workin CSCAPES was previously supported by NSF. Zoltan’s internal interface is currently being
redesigned to improve support for integration of third party libraries. For example, it has
recently acquired an interface to the parallel ordering functionality available in PT-Scotch.
The load-balancing and partitioning facilities in Zoltan have been used to enable awide variety of applications in structural mechanics, chemical engineering, groundwater flow,
biological systems, electronic circuit simulations, and molecular dynamics. Likewise, the tools
in OpenAD have been applied to a broad range of applications: modeling breast cancer, climate,weather, semiconductor devices, power networks, and groundwater; atmospheric chemistry;computational fluid dynamics; the network-enabled optimization server (NEOS); water reservoir
simulation; and chemical kinetics.
In the remainder of this short paper we provide a few “nuggets” of our recent works, togive a flavor of the ongoing research and development activities within CSCAPES. For more
information, consult the references listed in this paper and the website http://www.cscapes.org.
2. Dynamic Hypergraph Repartitioning and 2D Matrix Partitioning
In parallel adaptive computations, processor work loads may become unbalanced over time due
to changes in the computations. In adaptive mesh refinement, for example, the refinement andcoalescing of elements may cause significant load imbalance. This kind of application requires
dynamic load balancing. Dynamic load balancing has several goals with complicated trade-offs:
load balance in the new data distribution; low communication cost within the application; lowcost to migrate data from the old distribution to a new one; and short repartitioning time.
CSCAPES researchers have recently developed a hypergraph-based repartitioning method [5].
In this method, new vertices are added to represent each part. Then, each of these vertices isconnected to the objects in the same part with new (weighted) edges to model the data migration
cost. The resulting hypergraph accurately accounts for both the communication volume and themigration cost to move the data, and can be partitioned using existing algorithms and software.
In Figure 2, results from applying this hypergraph repartitioning method (HG-Repart) to an
adaptive mesh refinement problem from the ALEGRA shock physics code [2] is shown. HGRepart is compared with other connectivity-based approaches—“scratch-remap” hypergraphpartitioning (HG-Scratch) and graph-based partitioning (G-Scratch) and repartitioning (GRepart)—and the coordinate-based methods Recursive Coordinate Bisection (RCB) and Hilbert
Space-Filling Curve partitioning (HSFC). Scratch-remap methods compute partitions without
accounting for existing part assignments, and then remap parts to try to minimize migration
costs. While the coordinate-based methods do not model communication and migration

300
10250
HG−repart mig
HG−repart com
HG−scratch mig
HG−scratch com
G−repart mig
G−repart com
G−scratch mig
G−scratch com
RCB mig
RCB com
HSFC mig
HSFC com

200150100500

p=16

p=32

p=64

Number of processors

Partitioning time

Communication Volume (MBytes)

1001

HG−repart
HG−scratch
G−repart
G−scratch
RCB
HSFC

0.1

0.01

p=16

p=32

p=64

Number of processors

Figure 2. (Top) Sample adaptive finite element meshes at time-steps 0, 54, and 108; the smallest
mesh has 132,209 nodes, and the largest has 1,380,266 nodes. (Bottom) A comparsion of communication
volume and execution time over 109 time-steps of an adaptive finite element simulation on 16, 32 and64 processors of a Linux cluster (dual 2.2GHz AMD Opterons with 4 GB RAM and a Myrinet network).
costs explicitly, they are implicitly incremental, making them viable for repartitioning meshbased problems. The repartitioning hypergraph method HG-repart produced lower total cost(application communication volume plus data migration communication volume) than all othermethods. The coordinate-based methods were much faster than graph- and hypergraph-based
methods, but their total cost was high. Execution time for HG-repart was greater than G-repart,
indicating the need for faster heuristics in the hypergraph implementation for applications withrelatively low and homogeneous connectivity.
The multiplication of a sparse matrix by a vector is an important kernel in scientificcomputing. CSCAPES researchers have long studied how to optimize the performance of
this operation in parallel by reducing communication volume [6]. Recently, new partitioning
algorithms have been developed for novel two-dimensional data distributions that can becomputed efficiently using current partitioning software [1, 11].
3. Parallel Coloring
CSCAPES researchers have recently developed a framework for parallelizing greedy distance-1
coloring algorithms on distributed-memory computers [4]. Formulated in a generic manner, thetechniques employed in the framework include: careful exploitation of features of the initial data
distribution; speculation—maximizing concurrency by tentatively tolerating inconsistencies andthen detecting and resolving them; randomization; and infrequent, coarse-grain communication
among processors as opposed to frequent, fine-grained communication. In [4], several specializedalgorithms designed using the framework have been presented and evaluated; experiments carriedout on moderate size machines demonstrated good scalability.
In a work that is currently being finalized [3], the distance-1 coloring framework has been
extended for the distance-2 coloring problem, an archetypal model for partitioning problems
that arise in the efficient computation of sparse Jacobian and Hessian matrices using AD (seeFigure 3). To address the challange of accessing color information of distance-2 neighbors,

Figure 3.

c1c5

a11 a12 0 0 a15
a21 a22 0 0 0
a31 0 0 a34 00 0 a43 a44 a45

c3

r1
r2
r3
r4

c4
A

Ga

|V |
94,893
179,860
217,918
503,712952,203|E|
3,260,967
4,966,618
5,708,253
18,156,315
22,785,136

Degree
max
299
125
179
842
76

colors
(seq)
303
140
180
843
112

Gb
60

50

ldoor
pwtk
shipsec5
inline_1
pkustk13

Normalized number of colors

Agraphnamepkustk13
shipsec5
pwtk
inline 1ldoorc2

Speedup

a11 0 0 0 a15
0 a22 a23 0 00 a32 a33 a34 a35
0 0 a43 a44 a45
a51 0 a53 a54 a55

Structurally orthogonal column partition of
c1 a derivative matrix A and its
c2 representation as a distance-2
c3 coloring in the adjacency graph
c4 Ga (A) when A is symmetc5 ric, and bipartite graph Gb (A)when A is nonsymmetric.

403020100020406080

Number of processors

1001.15
1.1
1.05
10.95
0.90.85ldoor
pwtk
shipsec5
inline_1
pkustk130.8020406080100

Number of processors

Figure 4. Performance of the parallel distance-2 coloring algorithm on select application graphs. Themaximum degree in a graph is a lower bound on the optimal number of colors needed for distance-2coloring. A greedy sequential algorithm gave near optimal (or perhaps optimal) solution for these testgraphs. The parallel algorithm used in the worst case only 12% more colors.
efficient means of information exchange between processors hosting a pair of vertices that are
two edges away from each other are introduced. Experiments run on upto a hundred processors
have shown that the resulting parallel distance-2 coloring algorithm scales well and gives asolution very close to the sequential variant, which in turn often gives near optimal solution (seeFigure 4). MPI implementations of the distance-2 coloring algorithm as well as the distance-1
coloring algorithm discussed earlier have been incorporated in Zoltan.
In the case of Hessian computation, the distance-2 coloring model does not exploit theavailable symmetry. Models that exploit symmetry to various degrees exist. The most accurate
model for a direct method is star coloring and that for a substitution method is acyclic coloring.
CSCAPES investigators have in a recent work exploited the structures of two-colored induced
subgraphs in star and acyclic coloring to design novel sequential algorithms that have been shown
to be superior to previously known approaches [9]. The serial software they have developed has
been incorporated into the operator-overloading-based AD tool ADOL-C. Experiments using
ADOL-C on Jacobian computation in a Simulated Moving Bed process, a method used inchromatographic separation in chemical engineering, and Hessian computation on an electric
power grid problem have shown that the use of the coloring techniques in large-scale problems
reduces overall run time by several orders of magnitude [7, 8].
4. Automatic Differentiation Q
Suppose the derivatives of y = ni=1 xi with respect to all xi need to be computed. One
algorithm to compute this function is y1 = x1 ; yi = xi ∗ yi−1 , i = 2, . . . , n. One can thencompute the derivatives using y n = 1; y i = xi+1 ∗ y i+1 , xi+1 = yi ∗ y i+1 , i = n − 1, . . . , 1.
Then, ∂y/∂xi = xi . This algorithm to compute the derivatives corresponds to the reverse
mode of automatic differentiation. The reverse mode is optimal in this case andQwithin afactor two of optimal for all scalar functions. Now, consider the case of y = wT ni=1 Mi v,
where w and v are vectors and the Mi are dense, rectangular matrices (y remains a scalar).
Then one might compute y using b0 = wT ; bi = bi−1 Mi , i = 1, . . . , n; y = bn v. Then, thederivatives with respect to the elements of the Mi can be computed optimally (or nearly so)using bn = v; bi−1 = Mi bi , M i = bi bi−1 , = n, . . . , 1. However, if one computes the moreP
generalTTT
sequence a1 = w M1 v1 , a2 = w M1 M2 v2 , . . ., an = w M1 M2 . . . Mn vn , where y = ni=1 ai ,and then computes the derivatives of each a0 using the algorithm above (summing to produce
the derivatives of y), the result is far from optimal. The problem is that although the chained
vector-matrix products from the left can be reused from one ai to the next, the chained matrixvector products from the right cannot be reused because the vi vary. To illustrate the problem,
suppose all of the Mi are square matrices of size n × n. Then, the cost of computing the functionis O(n3 ) while the cost of computing the derivative matrices is O(n4 ). If, however, one recasts
the computation as b0 = wT ; bi = bi−1 Mi , ai = bi vi , i = 1, . . . , n; then, one can simply recur
backward on the bk sequence, with an extra term thrown in to account for the inner productwith xk . This brings the cost back down to O(n3 ). Thus, in developing methods to computethe derivatives of a function, it is important to consider the full algorithm used to compute thefunction, not just subcomponents. We identified a way to use a similar technique to reduce thecost of computing the gradient of a function arising in computational chemistry from O(n5 ) toO(n4 ) [10]. This reduction will have a major impact on the size of problems that can be solved.
Acknowledgments
Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin
Company, for the U. S. Department of Energy’s National Nuclear Security Administration undercontract DE-AC04-94AL85000. The work at Argonne was supported by the Mathematical,Information, and Computational Sciences Division subprogram of the Office of Advanced
Scientific Computing Research, U.S. Department of Energy under Contract W-31-109-Eng-38.
The CSCAPES Institute is supported by the U.S. Department of Energy’s Office of Sciencethrough grant DE-FC-0206-ER-25774, as part of its SciDAC program.
References[1] E. G. Boman. A nested dissection approach to sparse matrix partitioning. In Proc. of Appl. Math. Mech.,
to appear, 2008. (Presented at ICIAM07, Zurich, July 2007).
[2] E. A. Boucheron, K. H. Brown, K. G. Budge, S. P. Burns, D. E. Carroll, S. K. Carroll, M. A. Christon, R. R.
Drake, C. G. Garasi, T. A. Haill, J. S. Peery, S. V. Petney, J. Robbins, A. C. Robinson, R. Summers, T. E.
Voth, and M. K. Wong. ALEGRA: User Input and Physics Descriptions Version 4.2. Sandia National
Laboratories, Albuquerque, NM, 2002. Tech. Report SAND2002-2775.
[3] D. Bozdag, U. Catalyurek, A. Gebremedhin, F. Manne, E. Boman, and F. Ozguner. Distributed-memory
parallel graph coloring algorithms for Jacobian and Hessian computation. To be submitted to SIAM J.
Sci. Comput., 2008.
[4] D. Bozdag, A. Gebremedhin, F. Manne, E. Boman, and U. Catalyurek. A framework for scalable greedy
coloring on distributed memory parallel computers. Journal of Parallel and Distributed Computing,68(4):515–535, 2008.
[5] U. Catalyurek, E. Boman, K. Devine, D. Bozdag, R. Heaphy, and L. Riesen. Hypergraph-based dynamic
load balancing for adaptive scientific computations. In Proc. of 21th International Parallel and DistributedProcessing Symposium (IPDPS’07). IEEE, 2007. 11 pages.
[6] U. V. Çatalyürek and C. Aykanat. Hypergraph-partitioning based decomposition for parallel sparse-matrix
vector multiplication. IEEE Transactions on Parallel and Distributed Systems, 10(7):673–693, 1999.
[7] A. Gebremedhin, A. Pothen, A. Tarafdar, and A. Walther. Efficient computation of sparse Hessians usingcoloring and automatic differentiation. INFORMS Journal on Computing (to appear), 2008. 28 pages.
[8] A. Gebremedhin, A. Pothen, and A. Walther. Exploiting sparsity in Jacobian computation via coloring andautomatic differentiation: a case study in a simulated moving bed process. In The Fifth International
Conference on Automatic Differentiation, August 2008, Bonn, Germany (to appear), 2008. 10 pages.
[9] A. Gebremedhin, A. Tarafdar, F. Manne, and A. Pothen. New acyclic and star coloring algorithms withapplication to computing Hessians. SIAM Journal on Scientific Computing, 29(3):1042–1072, 2007.
[10] Personal communication with Ron Shepherd.
[11] M. Wolf, E. Boman, and B. Hendrickson. Optimizing sparse matrix-vector multiplication by cornerpartitioning. In Proc. of PARA08, to appear, 2008.