Implementation of Partial Separability in a
Source-to-Source Transformation AD Tool
Sri Hari Krishna Narayanan1 , Boyana Norris1 , Paul Hovland1 , and Assefaw
Gebremedhin2

Abstract A significant number of large optimization problems exhibit structure
known as partial separability, for example, least squares problems, where elemental
functions are gathered into groups that are then squared. The sparsity of the Jacobian of a partially separable function can be exploited by computing the smaller
Jacobians of the elemental functions and then assembling them into the full Jacobian. We implemented partial separability support in ADIC2 by using pragmas to
identify partially separable function values, applying source transformations to subdivide the elemental gradient computations, and using the ColPack coloring toolkit
to compress the sparse elemental Jacobians. We present experimental results for anelastic-plastic torsion optimization problem from the MINPACK-2 test suite.
Key words: Forward mode, partial separability, ADIC2, ColPack

1 Introduction
As introduced by Griewank and Toint [13], a function f : Rn 7→ R is consideredpartially separable if f can be represented in the formmf (x) = ∑ fi (x),

(1)

i=1

where fi depends on pi ≪ n variables. Bouaricha and Moré [5] and Bischof and ElKhadiri [3] , among others, have explored different ways to exploit the sparsity of
the Jacobians and Hessians of partially separable functions. To compute the (usuallydense) gradient ∇ f of f , one can first compute the much smaller (and possibly
1 Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL, USA,[snarayan,norris,hovland]@mcs.anl.gov.
2 Purdue University, West Lafayette, IN , USA, agebreme@purdue.edu

12Narayanan et al.

sparse) gradients of the fi elementals, then assemble the full gradient of f . Thisapproach can significantly reduce the memory footprint and number of floatingpoint operations for overall gradient computation compared with computing dense
gradients.
This paper describes the following new capabilities of the ADIC2 source transformation tool.
• Pragma-guided source transformations to perform scalar expansion of the elemental components of partially separable scalar-valued functions.
• Exposing of the sparsity present in the elemental Jacobians.
• Calculation of the compressed elemental Jacobians using ColPack.
• Combining of the elementals into the scalar valued result.

1.1 ADIC2
ADIC2 is a source transformation tool for
ADIC2
automatic differentiation of C and C++
EDG/ROSE
Open64
codes, with support for both the forward
front-ends Sage III
whirl
and reverse modes of AD [15]. ADIC2 uses
the ROSE compiler framework [17], whichOpen
relies on the EDG C/C++ parsers [10].
whirlToXAIF
SageToXAIF
Open64
Analysis
ADIC2 is part of the OpenAD frameworkat Argonne, whose general structure is ilXAIF
lustrated in Figure 1. Briefly, the process oftransforming the original source code into
OpenAD
Angel
code computing the derivatives consists of
xaifBooster
boost
several steps: (1) canonicalization (seman(AD source transformation)
xercestically equivalent transformations for removing features that hamper analysis orsubsequent AD transformations); (2) pro- Fig. 1 OpenAD component structure andgram analysis (e.g., control flow graph con- source transformation workflow.
struction, def-use chains); (3) generation of
the language independent XAIF intermediate representation; (4) AD transformationof the XAIF representation; (5) conversion of the resulting AD XAIF code back to
the ROSE intermediate representation; and (6) generation of C/C++ derivative code.
The general differentiation process as implemented by ADIC2 is discussed in detail in [15]. To exploit the sparsity of the gradients of partially separable functions,
we have implemented several extensions of the AD process, which are described in
Section 3.

Partial Separability in ADIC2

3

1.2 ColPack
When a Jacobian (or a Hessian) matrix is sparse, the runtime and memory efficiency
of its computation can be improved through compression by avoiding storing andcomputing with zeros. Curtis, Powell, and Reid demonstrated that when two or morecolumns of a Jacobian are structurally orthogonal, they can be approximated simultaneously using finite differences by perturbing the corresponding independent variables simultaneously [9]. Two columns are structurally orthogonal if there is no row
in which both columns have a nonzero. Coleman and Moré showed that the problemof partitioning the columns of a Jacobian into the fewest groups, each consisting ofstructurally orthogonal columns, can be modeled as a graph coloring problem [7].
The methods developed for finite-difference approximations are readily adapted toautomatic differentiation with appropriate initialization of the seed matrix [2].
ColPack is a software package containing algorithms for various kinds of graphcoloring and related problems arising in compression-based computation of sparse
Jacobians and Hessians [12]. In ColPack, the Jacobian is represented using a bipartite graph. Thus, a partitioning of the columns of a Jacobian into groups of structurally orthogonal columns is obtained using a distance-2 coloring of the column
vertices of the bipartite graph.
The coloring algorithms in ColPack are fast, yet effective, greedy heuristics. They
are greedy in the sense that vertices are colored sequentially one at a time and thecolor assigned to a vertex is never changed. The number of colors used by the heuristic depends on the order in which vertices are processed. Hence, ColPack contains
implementations of various effective ordering techniques for each of the coloring
problems it supports.
The rest of this paper is organized as follows. Section 2 contains a brief overview
of related work. Section 3 describes our implementation approach. We show experimental results for an optimization application use case in Section 4, and we conclude
Section 5 with a brief summary and discussion of future work.

2 Related Work
Bischof and El-Khadiri [3] describe the approach they took in implementing partial separability support in ADIFOR. Our approach, though similar in spirit, has anumber of significant differences. The ADIFOR approach assumed that the elemental functions were encoded in separate loops, while our approach does not rely on
this assumption and supports partial separability when multiple elemental functions
are computed in the same loop nest. To determine the sparsity pattern automaticallywhen the Jacobian structure is unknown, both ADIFOR and ADIC2 use runtime
detection through different versions of the SparsLinC library; in addition however
ADIC2 also relies on ColPack to compute a coloring, which is used to initialize theseed matrix for computing a compressed Jacobian (or Hessian) using the forward
mode of AD.

4

Narayanan et al.

Järvi [14] describes an object-oriented model for parameter estimation of a partially separable function. Conforti et al. [8] describe a master-slave approach to effect a parallel implementation of AD. The master process on the basis of the partially separable structure of the function, splits the workload among the slaves, andcollects the results of the distributed computation as soon as they are available.
Gay [11] describes the detection of partially separable structure automatically by
walking expression graphs. The structure is then used to facilitate explicit Hessiancomputations.
To exploit the sparsity in Jacobians of the elemental functions, we perform scalar
expansion, which is the conversion of a scalar value into a temporary array. For
example, scalar expansion can convert a scalar variable with a value of 1 to a vector
or matrix where all the elements are 1 (See Section 3). Typically scalar expansion
is used in compiler optimizations to remove scalar data dependences across loop
iterations to enable vectorization or automated parallelization. This transformation
is usually limited to counter-controlled for loops without control flow or function
calls. The size of the temporary arrays is typically determined through polyhedral
analysis of the iteration space of the loops containing the scalar operations thatare candidates for expansion. Polyhedral analysis is implemented in a number ofcompilers and analysis tools, including Omega [16], CHiLL [6], and PLuTo [4]. Ourcurrent approach to the implementation of scalar expansion does not use polyhedral
analysis. We describe our approach in more detail in Section 3.

3 Implementation
The changes required to support partial separability were implemented in our
ADIC2 source-to-source transformation infrastructure introduced in Section 1.1.
While this source translation can be performed in a standalone manner, it was convenient to implement it before the canonicalization step in ADIC2. The ROSE compiler framework, on which ADIC2 is based, parses the input code and generates
an abstract syntax tree (AST). ROSE provides methods to traverse and modify the
AST through the addition and deletion of AST nodes representing data structuresand program statements.
The translation is implemented as the following three traversals of the nodes
representing the function definitions within the AST.
1. Identification of partial separability (traversal T1 )2. Promotion of the elementals within loops and creation of a summation function(traversal T2 )3. Generation of elemental initialization loop (traversal T3 )In the first traversal, (T1 ), the statements within each function definition are examined. If the pragma $adic partiallyseparable is found, then the statement immediately following the pragma is an assignment statement whose left-hand side is
the dependent variable and the scalar-valued result of a partially-separable func-

Partial Separability in ADIC2

5

tion computation. The right-hand side of the assignment statement is an expression
involving the results of the elemental function computations. The names of the variables representing the elementals are specified in the pragma.
The second traversal, (T2 ), is initiated if the pragma $adic partiallyseparable isfound in T1 . The second traversal visits the nodes of the AST in two phases. The
first phase, called TopDown, visits the nodes starting at the node representing thefunction definition and proceeds down the tree to the leaf nodes. The second phase,called BottomU p, visits the leaf nodes first and proceeds up the tree all the way to
the node representing the function definition. Therefore T2 visits each node twice. In
both the phases information gained by visiting nodes can be passed in the directionin which nodes are visited. The following transformations occur in T2 .
1. Scalar expansion of elementals. In this transformation, the declaration of each of
the scalar elementals is changed into a dynamically allocated array. Next, eachreference to the scalar-variable elemental is modified into a reference to the arrayelemental.
To allocate memory for the array, its size is determined to be the maximum number of updates to the value inside any loop nest within the function body andis calculated by the BottomU p phase. In the BottomU p phase, when an innermost loop is visited, we create a parametrized expression whose value will be thenumber of times that loop executes (based only on its own loop bounds). This
expression is passed to its parent. Each parent loop multiplies its own local expression by the maximum of the expressions received from its children. Whenthe BottomU p phase is concluded, the function definition node will contain the
maximum number of updates to the elementals. This value is used to allocate
memory for the type-promoted elementals.
To modify the references of the scalar-variable elemental into a reference to thearray elemental, the index of each array reference is determined by the bounds ofthe surrounding loops. For example,

double elemental ;for ( j = lb0 ; j < ub0; j++) {for ( i = lb1 ; i < ub1; i++) {elemental = ...
( omitted )}}is transformed to

double ∗ elemental ;ADIC SPARSE Create1DimArray(&elemental, (ub1−lb1) ∗ (ub0−lb0));for ( j = lb0 ; j < ub0; j++) {for ( i = lb1 ; i < ub1; i++) {temp0 = j ∗( ub0 − lb0) + i ;elemental [temp0] = ...
( omitted )}}6

Narayanan et al.

This transformation is made possible using both phases. The expression that calculates the value of the index is created by passing the bounds of outer loops in
the TopDown phase to the inner loops. Then, in the BottomU p phase, if a reference is to an elemental is encountered, it is converted into an array reference andan appropriate assignment of the index-calculation expression to the array index
variable is inserted at the beginning of the inner loop body.
2. Creation of the result vector. In this transformation, the assignment statement
immediately following the pragma $adic partiallyseparable is modified. This
assignment statement is not affected by the previous transformation. A for loop
is created to replace the assignment statement. The assignment statement itself
is inserted into the body of the for loop. The for loop iterates as many times as
the maximum number of updates to the elementals, which is determined in the
BottomU p phase. For example,#pragma $adic partiallyseparable , elemental∗ scalar = temp ∗ ( elemental );is transformed to#pragma $adic partiallyseparable , elemental
for (k = 0; k< (ub1−lb1) ∗ (ub0−lb0); k++) {∗ scalar = temp ∗ ( elemental );}The elemental variable references within the assignment undergo scalar expansion, and the LHS of the assignment statement is replaced by an array reference.
The dimensions of this array are the maximum number of updates to the elementals, which was determined in the previous transformation, and the loop becomes#pragma $adic partiallyseparable , elemental
for (k = 0; k< (ub1−lb1) ∗ (ub0−lb0); k++) {temp vector [k] = temp∗(elemental )[ k ];}3. Summation of the result vector. Last, a call to a summation function is added
to the code. The arguments to the summation function are the scalar dependent
variable and the temporary array reference that forms the left-hand side of themodified assignment statement. For example,

ADIC SPARSE Summation(scalar, temp vector);
is a call that can result from this transformation.
In the third traversal (T3 ), the statements within each function definition are examined again. If the pragma $adic partialelemental is found, then the statement
immediately following the pragma is an assignment statement whose left-hand side
is an elemental and whose right-hand side is an initialization value. Such an assignment statement is not modified by any earlier transformation. Similar to the creation

Partial Separability in ADIC2

7

of a result vector, a for loop is created that iterates as many times as the maximum
number of updates to the elementals, which was determined in the previous transformation (T2 ). The assignment statement itself is inserted into the body of the forloop. Finally the loop replaces the annotated assignment statements. For example,#pragma $adic partialelemental
elemental = 0.0;
is transformed to#pragma $adic partiallyseparable , elemental
for (k = 0; k< (ub1−lb1) ∗ (ub0−lb0); k++) {elemental [k] = 0.0;}

4 Experimental Results
We evaluated the performance of the partial separability support in ADIC2 by using
a two-dimensional elastic-plastic torsion model from the MINPACK-2 test problem
collection [1]. This model uses a finite-element discretization to compute the stressfield on an infinitely long cylindrical bar to which a fixed angle of twist per unit
length has been applied. The resulting unconstrained minimization problem can be
expressed as min( f (u)), where f : ℜn → ℜ, f (u) is the quadraticZ1{ k∇u(x)k2 − cu(x)}dx,(2)2
c is a constant, and D is a bounded domain with a smooth boundary.
In our experiments, we applied ADIC2 to the C version of the function implementation after the original Fortran code was manually translated into C. As described in more detail in Section 3, we insert two types of pragmas to define (i)the elementals of the partially separable function and (ii) the initialization of theelementals.
f (u) =

D...
#pragma $adic partialelemental
fquad = zero ;#pragma $adic partialelemental
flin = zero ;...
/∗ computation of the elementals fquad and fln of thefunction f and their sparse Jacobians ∗/...
#pragma $adic partiallyseparable , fquad, flin
∗ f = area ∗(p5∗fquad+flin );8

Narayanan et al.

...
In the initialization portion, the value of the double constant zero is 0.0. Afterthe transformation (described in Section 3), the last portion of the code above (thefunction computation) is transformed to...
#pragma $adic partiallyseparable , fquad, flin
ADIC SPARSE Create1DimArray(& adic temp f, ad var max);for ( adic temp0 = 0; adic temp0 < ad var max; adic temp0++) {adic temp f [ adic temp0 ] = ( area ∗ (( p5 ∗ fquad[ adic temp0 ])+ flin [ adic temp0 ]));}ADIC SPARSE Summation(f, adic temp f,ad var max);...
where ad var max is the size of the gradient vector array for the full Jacobian.
We validated the correctness of the sparse computation by comparing it with thevalues produced by the analytical version. For example, for an input array size 100,
the error estimated by the norm of the difference was approximately 1.5e-16, or near
the limit of machine precision for floating-point computations.
We measured the execution time of the gradient computation on an Intel Xeon
workstation with dual quad-core E5462 Xeon processors (8 cores total) running
at 2.8 GHz (1600 MHz FSB) with 32 KB L1 cache, 12 MB of L2 cache (6 MB
shared per core pair), and 16 GB of DDR2 FBDIMM RAM, running Linux kernel
version 2.6.35 (x86-64). All codes were compiled with gcc version 4.4.5 with -O2
optimizations enabled.

103
102

Time (secs)

101

105

Analytical
Dense
PSS (Total)- Sparsity Detection- Seed Generation- Gradient Computation

Dense
PSS

104

Gradient Array Size

104

100

10-1

103

10-2
10-3

102

10-4
10-5
10-6 210103

Input Array Size

104

105

101 2
10

103

Input Array Size

104

105

Fig. 2 Left: comparison between the runtimes of the three gradient versions: hand-coded analyticalgradients, dense AD, and partially separable sparse AD (PSS). A breakdown of the steps in the PSS
version is also shown: sparsity detection, seed generation, and gradient computation and Jacobianrecovery. Right: gradient array sizes for the dense and PSS versions.

Partial Separability in ADIC2

9

Figure 2 (left) shows the execution times for computing the Jacobian of the function in Eq. 2 using three approaches: (1) manually implemented analytical derivative computation; (2) dense, forward mode AD (using ADIC2); and (3) sparse, forward mode AD (using ADIC2) with additional sparsity detection (using SparsLinC)
and Jacobian compression (using ColPack). Additionally, the graph shows the time
taken to perform the steps involved in approach (3): sparsity detection, seed generation, and gradient computation and Jacobian recovery. The analytic version performs
best, as expected. The forward-mode dense AD version is between 500 and 3,400
times slower than the manually optimized analytical derivatives computation, while
the partially separable sparse AD version achieves performance within a factor of
6 of the analytic version for small array sizes, and less than a factor of 2 slower
for larger array sizes. The right side of Fig. 2 shows the reduction in memory requirements (ranging from 500- to 4000-fold) for storing the gradients using the PSS
approach compared with dense gradients. The size of the gradient array in the PSS
approach is the number of colors used by ColPack to compresses the sparse Jacobian.The implication of this result is that for arrays sizes beyond a certain size, thedense approach is not feasible as it’s memory requirement will exceed the available
machine capacity. However, the approach exploiting partial becomes infeasible at amuch larger input size.

5 Conclusion
We presented an approach to exploiting sparsity in the computation of gradients ofpartially separable functions, which are common in large-scale optimization. We
identify partially separable computations by using pragmas, which guide our source
transformation system to perform scalar expansion and generate efficient forward
mode AD code for computing the gradients of the elemental functions. In addition, we exploit sparsity in these gradients by using the SparsLinC library and the
ColPack coloring toolkit to enable efficient forward mode AD by using staticallyallocated compressed dense vectors for computing the gradients of the elemental
functions. We evaluated the performance of our implementation using a case studyof the elastic-plastic torsion problem in the MINPACK-2 test suite, demonstratingthat (1) exploiting partial separability and sparsity significantly reduces the memory
requirements of the generated code, enabling the solution of larger problems thanpossible with dense forward mode, and (2) the performance of the best AD version
compares favorably with that of the hand-coded gradients. In future work we will
extend ADIC2 to remove certain restrictions, for example, the assumption that thegradient vectors of different elemental functions are of the same size. We also plan
to integrate the polyhedral analysis currently being developed in ROSE and add
support for exploiting partial separability when using the reverse mode in ADIC2.
Acknowledgments This work was supported by the U.S. Dept. of Energy
Office of Science Applied Mathematics Program under Contract No. DE-AC0206CH11357.

10

Narayanan et al.

References
1. Averick, B.M., Carter, R.G., Moré, J.J., Xue, G.L.: The MINPACK-2 test problem collection.
Tech. Rep. Preprint MCS-P152-0694, Mathematics and Computer Science Division, Argonne
National Laboratory (1992)
2. Averick, B.M., Moré, J.J., Bischof, C.H., Carle, A., Griewank, A.: Computing large sparse
Jacobian matrices using automatic differentiation. SIAM J. Sci. Comput. 15(2), 285–294(1994). DOI 10.1137/0915020. URL http://link.aip.org/link/?SCE/15/285/
1
3. Bischof, C., El-Khadiri, M.: On exploiting partial separability and extending the compile-timereverse mode in ADIFOR. Technical Memorandum ANL/MCS–TM–163, Mathematics andComputer Science Division, Argonne National Laboratory, Argonne, Ill. (1992). ADIFOR
Working Note # 7.
4. Bondhugula, U., Hartono, A., Ramanujam, J., Sadayappan, P.: PLuTo: A practical and fully
automatic polyhedral program optimization system. In: Proc. ACM SIGPLAN 2008 Conference on Programming Language Design and Implementation (PLDI) (2008)
5. Bouaricha, A., Moré, J.J.: Impact of partial separability on large-scale optimization. Comp.
Optim. Appl. 7(1), 27–40 (1997)
6. Chen, C.: Model-guided empirical optimization for memory hierarchy. Ph.D. thesis (2007)
7. Coleman, T.F., Moré, J.J.: Estimation of sparse Jacobian matrices and graph coloring problems. SIAM Journal on Numerical Analysis 20(1), 187–209 (1983)
8. Conforti, D., Luca, L.D., Grandinetti, L., Musmanno, R.: A parallel implementation of automatic differentiation for partially separable functions using pvm. Parallel Computing
22(5), 643 – 656 (1996). DOI 10.1016/0167-8191(96)00014-2. URL http://www.
sciencedirect.com/science/article/pii/0167819196000142
9. Curtis, A.R., Powell, M.J.D., Reid, J.K.: On the estimation of sparse Jacobian matrices. Journal of the Institute of Mathematics and Applications 13, 117–119 (1974)
10. Edison Design Group C++ Front End.
http://www.edg.com/index.php?
location=c_frontend
11. Gay, D.M.: More AD of nonlinear AMPL models: Computing Hessian information and exploiting partial separability. In: M. Berz, C. Bischof, G. Corliss, A. Griewank (eds.) Computational Differentiation: Techniques, Applications, and Tools, pp. 173–184. SIAM
12. Gebremedhin, A.H., Nguyen, D., Patwary, M., Pothen, A.: ColPack: Software for graph coloring and related problems in scientific computing. Tech. rep., Purdue University (2011)
13. Griewank, A., Toint, P.L.: On the unconstrained optimization of partially separable functions.
In: M.J.D. Powell (ed.) Nonlinear Optimization 1981, pp. 301–312. Academic Press, New
York, NY (1982)
14. Järvi, J.: Object-oriented model for partially separable functions in parameter estimation. Acta
Cybernetica 14(2), 285–302 (1999)
15. Narayanan, S.H.K., Norris, B., Winnicka, B.: ADIC2: Development of a component
source transformation system for differentiating C and C++. Procedia Computer Science
1(1), 1845–1853 (2010). DOI DOI:10.1016/j.procs.2010.04.206. URL http://www.
sciencedirect.com/science/article/pii/S1877050910002073.
ICCS
2010
16. Pugh, W.: The Omega Project web page. http://www.cs.umd.edu/projects/
omega/
17. Quinlan, D.: ROSE: Compiler support for object-oriented frameworks. Tech. Rep. UCRLID-136515, Lawrence Livermore National Laboratory (1999). URL http://www.osti.
gov/bridge/servlets/purl/793936-hdq2WX/native/793936.PDF