Sparse Matrix Code Dependence Analysis Simplification at Compile Time
Mahdi Soltan Mohammadi Kazem Cheshmi Ganesh Gopalakrishnan Mary Hall Maryam Mehri Dehnavi Anand Venkat Tomofumi Yuki Michelle Mills Strout



         
                                        
                                        

                                        
                                        


                                        
                                        















                                        



  Department of Computer Science              
  University of Arizona            






kingmahdi@cs.arizona.edu          



                                        



  Department of Computer Science              
  University of Toronto            






kazem@cs.toronto.edu          



                                        



  School of Computing              
  University of Utah            






ganesh@cs.utah.edu          



                                        



  School of Computing              
  University of Utah            






mhall@cs.utah.edu          



                                        



  Department of Computer Science              
  University of Toronto            






mmehride@cs.toronto.edu          



                                        




  Intel Corporation            






anand.venkat@intel.com          



                                        



                
  Univ. Rennes, Inria, CNRS, IRISA            






tomofumi.yuki@inria.fr          



                                        



  Department of Computer Science              
  University of Arizona            






mstrout@cs.arizona.edu          












Analyzing array-based computations to determine
data dependences is useful for many applications
including automatic parallelization, race detection, 
computation and communication overlap, verification, and shape analysis.
For sparse matrix codes, array data dependence
analysis is made more difficult by the use of
index arrays that make it possible to store only the nonzero entries of 
the matrix (e.g., in A[B[i]], B is an index array). 
Here, dependence analysis 
is often stymied by such indirect array accesses due to the
values of the index array not being available at compile time. 
Consequently, many dependences cannot be proven unsatisfiable 
or determined until runtime.
Nonetheless, index arrays in sparse matrix codes often have properties 
such as monotonicity of index array elements that can be exploited to
reduce the amount of runtime analysis needed.
In this paper, we contribute a formulation of array data dependence 
analysis that includes encoding  index array properties as 
universally quantified constraints.
This makes it possible to leverage
existing SMT solvers to determine
whether such dependences are unsatisfiable and
significantly reduces the number of dependences that
require runtime analysis in a set of eight sparse matrix kernels.
Another contribution is an algorithm for simplifying 
the remaining satisfiable data dependences by discovering
equalities and/or subset relationships.
These simplifications are essential to make
a runtime-inspection-based approach feasible.






<ccs2012>
<concept>
<concept_id>10011007.10011006.10011008</concept_id>
<concept_desc>Software and its engineering General programming languages</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>
<concept_id>10003456.10003457.10003521.10003525</concept_id>
<concept_desc>Social and professional topics History of programming languages</concept_desc>
<concept_significance>300</concept_significance>
</concept>
</ccs2012>


[500]Software and its engineering General programming languages
[300]Social and professional topics History of programming languages





  







Sparse Matrix Code Dependence Analysis Simplification at Compile Time
    Michelle Mills Strout
    December 30, 2023
=====================================================================








§ INTRODUCTION


Data dependence analysis answers questions about which memory accesses 
in which loop iterations access the same memory location thus creating 
a partial ordering (or dependence) between loop iterations.
Determining this information enables 
iteration space slicing <cit.>,
provides input to race detection, makes automatic parallelization 
and associated optimizations such as tiling or communication/computation 
overlap possible,
and enables more precise data-flow analysis, or abstract interpretation.
A data dependence exists
between two array accesses that (1) access the same array element with
at least one access being a write and (2) that access occurs within the loop
bounds for each of the accesses' statement(s).
These conditions for a data dependence has been posed 
as  a constraint-based decision problem <cit.>, a data-flow
analysis with polyhedral set information <cit.>, and
linear memory access descriptors <cit.>.
However, such approaches require a runtime component when analyzing
codes with indirect memory accesses (e.g., A[B[i]], B being an index array)
such as those that occur in sparse matrix codes.
In this paper, we present an approach to improve the precision of
compile-time dependence analysis for sparse matrix codes
and simplification techniques for decreasing the complexity of
any remaining runtime checks.


Sparse matrix computations occur in many codes, such as graph analysis,
partial differential equation solvers, and molecular dynamics
solvers. Sparse matrices save on storage and computation by only
storing the nonzero values in a matrix.
Figure <ref>
illustrates one example of how the sparse matrix vector multiplication 
(y⃗ = Ax⃗) can be written to use a common sparse matrix format 
called compressed sparse row (CSR). In CSR, the nonzeros are organized by row,
and for each nonzero, the column for that nonzero is explicitly stored.
Since the percentage of nonzeros in a matrix can be less than 1%,
sparse matrix formats are critical for many computations to fit
into memory and execute in a reasonable amount of time.
Unfortunately, the advantage sparse matrices have on saving storage compared
to dense matrices comes with the cost of complicating program analysis. 
Compile-time only approaches resort to conservative 
approximation <cit.>.
Some approaches use runtime dependence tests to complement compile time
analysis and obtain more precise answers <cit.>. 
Runtime dependence information is also used to detect task-graph 
or wavefront parallelism that arises due to the 
sparsity of dependences <cit.>.







The data dependence analysis approach presented here is constraint-based.
Some constraint-based data dependence analysis 
approaches <cit.>
represent the index arrays in sparse matrix computations as uninterpreted functions 
in the data dependence relations. For example, the loop bounds for the k loop in 
Figure <ref> can be represented as rowptr(i) ≤ k < rowptr(i+1).
Previous work by others generates simplified constraints at compile time
that can then be checked at runtime 
with the goal of finding fully parallel loops <cit.>.
We build on the previous work of <cit.>.
In that work, dependences in a couple of sparse matrix codes were
determined unsatisfiable manually, simplified using equalities found through
a partial ordering of uninterpreted function terms, or approximated by
removing enough constraints to ensure a reasonable runtime analysis complexity.
In this paper, we automate the determination of unsatisfactory data dependences, 
 find equalities using the Integer Set Library (ISL) <cit.>, 
have developed ways to detect data dependence subsets for simplifying runtime analysis, 
and perform an evaluation with eight popular sparse kernels.




In this paper, we have two main goals:
(1) prove as many data dependences as possible
to be unsatisfiable, thus reducing the number of dependences 
that require runtime tests; and 
(2) simplify the satisfiable data dependences so that a
runtime inspector for them has complexity less than or equal
to the original code.
Figure <ref> shows an overview of our approach.
We use the ISL library <cit.> like an SMT solver to determine which
data dependences are unsatisfiable. Next, we manipulate any
remaining dependence 
relations using IEGenLib <cit.> and ISL libraries to discover equalities 
that lead to simplification.



Fortunately, much is known about the index arrays that 
represent sparse matrices as well as the assumptions made by 
the numerical algorithms that operate on those matrices. 
For example, in the CSR representation shown in Figure <ref>, 
the rowptr index array is monotonically strictly increasing.
In Section <ref> we explain how such information 
can be used to add more inequality and equality constraints 
to a dependence relation. The new added constraints in some cases
cause conflicts, and hence we can detect that
those relations are unsatisfiable.






The dependences that cannot be shown as unsatisfiable
at compile time, still require runtime tests.
For those dependences, the goal is to simplify the constraints and reduce 
the overhead of any runtime test. Sometimes index array properties 
can be useful for reducing the complexity of runtime inspector 
by finding equalities that are implied by the dependence constraints 
in combination with assertions about the index arrays.
The equalities such as i = col(m') can help us remove
a loop level from the inspector, i in the example. Since some equality
constraint would allow us deduce value of an iterator in the inspector from
another one, e.g, we can deduce i values from m' values using i = col(m').
Another simplification involves determining when data dependence relations
for a code are involved in a subset relationship.  When this occurs,
runtime analysis need only occur for the superset.



This paper makes the following contributions:


  * An automated approach for the determination of unsatisfiable
 dependences in sparse codes.


  * An implementation of an instantiation-based decision procedure that discovers
equality relationships in satisfiable dependences.


  * An approach that discovers
subset relationships in satisfiable dependences thus reducing run-time
analysis complexity further.


  * A description of common properties of index arrays arising in sparse
matrix computations, expressed as universally quantified constraints.


  * Evaluation of the utility of these properties for determining
unsatisfiability or simplifying dependences from a 
suite of real-world sparse matrix codes.












§ BACKGROUND: DATA DEPENDENCE ANALYSIS



Data dependence analysis of a loop nest is a common code analysis that is used
in different applications, such as automatic parallelization <cit.>
and data race detection <cit.>.
This section reviews the data dependence analysis process and how that process
differs when analyzing sparse matrix codes. Then, we review some of the applications
of data dependence analysis including an example of its use for
finding wavefront parallelism in sparse codes.



 §.§ Data Dependence Analysis



 A data dependence occurs between two
iterations of a loop when both of the iterations access the same memory
location and at least one of the accesses is a write.  Data dependence
constraints are of the following form:

    Dep : (∃I⃗,I⃗'⃗)(I⃗≺I⃗'⃗∧ F(I⃗)=G(I⃗'⃗)
    		∧ Bounds(I⃗) ∧ Bounds(I⃗'⃗)),

where I⃗ and I⃗'⃗ are iteration vector instances from the same loop
nest, F and G are array index expressions to the same array
with at least one of the accesses being a write, and Bounds(I⃗) expands
to the loop nest bounds for the I⃗ iteration vectors.
In this paper, the term dependence relation is
used interchangeably with dependence constraints by viewing them as a relation
between  I⃗ and I⃗'⃗.




For example, consider the dense matrix implementation for forward solve
in Figure <ref>. Forward solve solves for the vector
u⃗ in the equation Au⃗=f⃗ assuming that the matrix is lower
triangular (i.e., nonzeros are only on the diagonals or below as shown in the
example in Figure <ref>).  The dense forward solve code
has the following dependences for the outermost loop i:


  * A loop-carried dependence due to the scalar tmp
variable.  However, since tmp is written before being read in
each iteration of the i loop, it is privatizable, 
which means each processor in a potential parallelization of the i 
loop can be given its own private copy of tmp.

  * A loop-carried dependence between the write u[i]
in Statement S2 and the read u[j] in Statement S1
with constraints

    (∃ i,j, i',j')(   i<i' ∧ i=j' ∧ 0 ≤ i,i' <N ∧
    0 ≤ j <i ∧ 0 ≤ j' <i').




The iterators i' and j' are different instances of i and j.
This dependence due to the writes and reads to array u is satisfiable
because the computation for any row i' depends on all previous rows i< i'.
This means that the outer loop in dense forward solve is fully ordered
due to data dependences
and therefore not parallelizable.






 §.§ Sparse Codes and Runtime Parallelism



For sparse codes, after compile time dependence analysis, some remaining 
dependences may involve index arrays as subscript expressions.
The data dependence constraints can use uninterpreted functions to 
represent the index arrays at compile time.
Because the values of the index arrays are unknown until run time, 
proving such dependences are unsatisfiable may require runtime dependence testing.  
However, even when dependences arise at runtime, it still may be possible to implement
a sparse parallelization called wavefront parallelization.
Identifying wavefront parallelizable loops combines
compile time and runtime analyses.  The compiler generates inspector
code to find the data dependence graph at runtime.


We now consider the sparse forward solve with Compressed Sparse Row
CSR format in Figure <ref>. We are interested in detecting
loop-carried dependences of the outermost loop.  There are two pairs of
accesses on array u in S1 and S2 that  can potentially cause
loop-carried dependences: u[col[k]] (read), u[i] (write); and u[i] (write), u[i] (write). The constraints for the two dependence tests
are shown in the following.





Dependences for the write/write u[i] in S2:

    (1)    i = i' ∧ i < i' ∧ 0 ≤  i < N ∧ 0 ≤ i' < N  
    ∧ rowptr(i) ≤ k < rowptr(i+1) 
    ∧ rowptr(i') ≤ k'  < rowptr(k'+1)


    (2)   i = i' ∧ i' < i ∧ 0 ≤  i < N ∧ 0 ≤ i' < N  
    ∧ rowptr(i) ≤ k < rowptr(i+1) 
    ∧ rowptr(i') ≤ k'  < rowptr(k'+1)


Dependences for read u[col[k]] and write u[i] in S1, and S2:

    (3)   i = col(k') ∧ i < i' ∧ 0 ≤  i < N ∧ 0 ≤ i' < N  ∧
     rowptr(i') ≤ k' < rowptr(i'+1)


    (4)   i = col(k') ∧ i' < i ∧ 0 ≤  i < N ∧ 0 ≤ i' < N  ∧
     rowptr(i') ≤ k' < rowptr(i'+1)


These dependences can be tested at runtime when concrete interpretations for
the index arrays (contents of arrays rowptr and col) are available.
The runtime dependence analyzers, called inspectors <cit.>,
may be generated from the dependence constraints <cit.>. 

Suppose the matrix 
in the forward solve code in Figure <ref>
had the nonzero pattern as in Figure <ref>. The runtime check 
would create the dependence graph for this example based on the
four dependences above as shown in Figure <ref>.
Once the dependence graph is constructed
a breadth-first traversal of the dependence
graph can derive sets of iterations that may be safely 
scheduled in parallel without a dependence violations, with each level set
being called a wavefront as shown in Figure <ref>.








 §.§ Applications of the Sparse Data Dependence Analysis


Besides wavefront parallelism, there are many other uses for
sparse data dependence analysis.  Any application of sparse data dependence
analysis would benefit from a reduction in the number of data dependences
that need to be inspected at runtime and from any complexity reduction
of data dependences that do require runtime inspection.
Here we summarize some of those
applications.

Race detection:
Dynamic race detection is an essential prerequisite to the
parallelization of existing sequential codes.
While the front-end static analysis methods employed
in these checkers <cit.> can often suppress race
checks on provably race-free loops, they fail
to do so when presented with non-affine access patterns
that occur in sparse matrix codes.
In addition to significantly increasing runtimes,
the shadow memory cells employed by dynamic race
checkers also increases memory pressure, often by a factor of four.
The techniques presented in this paper can help suppress
race checks when we can prove the independence of loop iterations.

Dynamic program slicing: 
Pugh and Rosser introduced the concept of iteration space slicing where
program slicing is done on a loop iteration basis using 
Presburger representations <cit.>.
Similar dynamic approaches for tiling across loops in sparse codes 
were presented by various groups <cit.>.
All of these techniques would require runtime data dependence
analysis, thus disproving dependences or reducing the complexity
of inspecting dependences at runtime would be applicable.

High-level synthesis:
Optimizations in high-level synthesis (HLS) uses runtime dependence checks.
In HLS, it is important to pipeline the innermost loops to get
efficient hardware.  Alle et al. have proposed using runtime dependence checks
to dynamically check if an iteration is in conflict with those currently in the
pipeline, and add delays only when necessary <cit.>.

Distributed memory parallelization: 
Another possible application of our work can be found in the work by <cit.>.
The authors produce distributed parallel code that uses MPI
for loops where there might be indirect loop bounds, 
and/or array accesses.
The read and write 
sets/data elements of each process are computed via an inspector 
where indirect accesses are involved to determine if each process is 
reading/writing data that is owned by other processes. 
Basumallik and Eigenmann use run-time inspection of data dependences
to determine how to reorder a loop to perform computation and communication
overlap <cit.>.





§ AUTOMATING (UN)SATISFIABILITY ANALYSIS FOR SPARSE DATA DEPENDENCES


For any application of data dependence analysis for sparse codes,
the best outcome is to determine that a potential data dependence
is unsatisfiable.  Any dependence that is unsatisfiable does not
for runtime analysis.
Previous work used domain-specific knowledge about the index
arrays used to represent sparse matrices to guide manual determination
of unsatisfactory data dependences <cit.>.
In this paper, we show how to automate this process 
by specifying the domain-specific knowledge as universally
quantified constraints on uninterpreted functions and then
using instantiation methods similar to those used by 
SMT solvers to produce more constraints that can cause
conflicts.



 §.§ Detecting Unsatisfiable Dependences Using Domain Information



As an example of how domain information can be used to
show dependences are unsatisfiable consider the following 
constraints from a dependence relation:

    i' < i ∧  k = m' ∧  0 ≤ i,i' <n ∧ 
    rowptr(i) ≤ k < rowptr(i+1) ∧  rowptr(i'-1) ≤ m' < rowptr(i').

Some relevant domain information is that 
the rowptr index array 
is strictly monotonically  increasing:

    (∀ x_1,x_2)(x_1 < x_2  rowptr(x_2) < rowptr(x_2)).


Since the dependence relation in question has the constraints
i' < i. Then, using the above strict monotonicity information 
would result in adding rowptr(i') < rowptr(i).
But, considering  the constraint, k = k', rowptr(i) ≤ k, and
m' < rowptr(i') we know that, rowptr(i) < rowptr(i'). 
This leads to a conflict, 

    rowptr(i) < rowptr(i')  ∧  rowptr(i') < rowptr(i).

This conflict would indicate the dependence in question was
unsatisfiable and therefore does not require any runtime analysis.




 §.§ Universally Quantified Assertions about Index Arrays

Even if a formula that includes uninterpreted function calls is satisfiable in
its original form, additional constraints about the uninterpreted functions may
make it unsatisfiable. This has been exploited abundantly in
program verification community to obtain more precise results <cit.>. A common approach to express such additional
constraints is to formulate them as universally quantified assertions.  
For instance, <cit.> use following to indicate 
that array a is sorted within a certain domain:


    (∀ x_1, x_2)(i < x_1 ≤ x_2 < j  a(x_1) ≤ a(x_2) ).


There are several methods that SMT solvers use to reason about quantified
formulas, the most common one being quantifier
instantiation <cit.>.
In quantifier instantiation, instances of universally quantified assertions, 
where the universally quantified variables are replaced with ground terms, 
are added to the original formula. 
Any of the added constraints might contradict constraint(s) in the formula
that would show the original formula is unsatisfiable.
For the general case of quantified first order logic, there is no complete 
instantiations procedure. That means the instantiation can go on forever
not exactly knowing whether the formula is satisfiable or unsatisfiable.
In some limited cases, the quantified assertions can be completely
replaced by a set of quantifier instances to construct an equisatisfiable
quantifier-free formula <cit.>.


Combining the constraints from dependences with arbitrary universally
quantified assertions would create a first order 
logic theory that in general is undecidable. Undecidability would
imply that we cannot implement an algorithm for deciding the formulas
that would always terminate.
Numerous works such as <cit.> present 
different decidable fragments of first order logic. 
The approach that these works use to make decidable fragments is to
put restrictions what type of universally quantified assertions
can be used. 
The restriction are usually on on syntax of the allowed 
assertions <cit.>,
and sometimes specific properties that a specific instantiation
procedure for assertions must have <cit.>.
We perform a terminating instantiation that is sound but incomplete.
In other words, the dependences we determine unsatisfiable
are in fact unsatisfiable, but we may characterize some unsatisfiable
constraints as may satisfiable.




 §.§ Domain Information about Index Arrays


We represent domain information about index arrays as universally
quantified assertions. In this section, we illustrate 
some assertions relevant to numerical benchmarks
and relate the corresponding assertions to the existing
theory fragments.  
Table <ref> in Section <ref> lists all the 
assertions we use in the evaluation.  Below are some example properties.


For the forward solve with compressed sparse
row (CSR) code in
Figure <ref>, we know the following:



  * Monotonic index arrays: 
The row index array values increase monotonically.
This property of index arrays can be expressed with an
assertion about the uninterpreted function symbol
that represents the index array. 
For instance, in the example the rowptr() function is monotonically increasing. 
If we assume that all the sparse matrix rows have at least one nonzero,
then rowptr() is strictly monotonically increasing.
This assertion can be encoded as follows: 
	
    (∀ x_1,x_2)(x_1<x_2  rowptr(x_1)<rowptr(x_2)).



  * Lower Triangular Matrix:
The forward solve algorithm shown in Figure <ref>
operates on lower triangular matrices.
For the CSR format that leads to the following domain-specific
assertion:
        
    (∀ x_1,x_2)(x_1<rowptr(x_2)  col(x_1)<x_2 )

This indicates that nonzeros of rows before row i have columns less than i.


The domain information in Table <ref> in Section <ref> 
can be represented with following 
general forms:


    1.    (∀ x_1, x_2)(x_1+c_1 ≤ x_2  f(x_1)+c_2 ≤ f(x_2)


    2.    (∀ x_1, x_2)(x_1+c_1 ≤ x_2  f(x_1)+c_2 ≤ g(x_2)


    3.    (∀ x_1, x_2)(x_1+c_1 ≤ f(x_2)  g(x_1)+c_2 ≤ x_2)

Where c_1 and c_2 can be 0 or 1. 
The first and second assertions fit the decidable LIA fragment that is
presented by <cit.>.  However, to the best of
our knowledge the third assertion form does not fit any previously presented
decidable fragment, and its decidability remains open. 


Modern SMT solvers are equipped with heuristic-based quantifier instantiations
to reason about quantified formulas. Existing techniques for quantifier
instantiation can construct the set of instantiations for deciding some of our
assertions, e.g., non-strict monotonicity, but not for all of them.  For
unsatisfiable formulas with universal quantifiers where the solver only needs a
small set of relevant instances to find contradicting constraints, the existing
heuristics can work well. For all our examples, both Z3 and CVC4 were able to
identify all unsatisfiable dependences.  The solvers also time out for
satisfiable ones given a small timeout (5 seconds).  This is as expected, since
specific instances of universally quantified formulas usually do not help in
proving that the quantified formula is satisfiable.


Nonetheless, we cannot just use a conventional SMT solver like Z3 in
our context.
The key reason is that we are not just interested in satisfiability of the dependence
constraints.  If unsatisfiability cannot be proven statically, runtime checks
will be generated.  It is important for these runtime checks to be as fast as
possible, and hence we are also interested in using the assertions to decrease
the cost of runtime checks.  
For example, additional equalities means the data dependence inspector
iteration space
has lower dimensionality, thus reducing the algorithmic complexity of
runtime checks.  We illustrate the complexity reduction through instantiation
of assertions with two examples in Section <ref>.




 §.§ Detecting Unsatisfiable Sparse Dependences


Instantiation-based quantifier elimination is a natural choice for our context,
since we seek to either prove unsatisfiability or find additional constraints
that simplifies runtime checks. Unfortunately, our assertions are not fully
covered by decidable fragments <cit.> where equisatisfiable
quantifier-free formulas can always be obtained.  Nonetheless, using
inspiration from the decidable fragments <cit.> we have a
procedure that detects all unsatisfiable examples from our benchmark suite that
represent a wide range of numerical analysis codes. 

Note that we can show our general assertions (1), (2), and (3), presented
in Section <ref> as:


    ∀x⃗, φ_I(x⃗) φ_V(x⃗)


Where x⃗ denotes vector of quantified variables,
φ_I(x⃗) denotes antecedent of the assertion,
and φ_V(x⃗) denotes consequent of the assertion.
Then the following definitions define our procedure to instantiate
quantified variables, 
and potentially use a SMT to detect their unsatisfiability.




Definition 1 (E) 
We define E to be the set of expressions used as arguments to 
all uninterpreted function calls in the original set of constraints.
We use this set to instantiate quantified assertions.



Definition 2 (UNSAT_ψ)

2pt

  * The inference rules for turning the universally
quantified predicates into quantifier-free predicates is as follows:



[forall] 

ψ[ ∀x⃗,  φ_I(x⃗) φ_V(x⃗) ]


ψ[ ⋀_x⃗ ∈E^n    (φ_I(x⃗) φ_V(x⃗)) ]




where E^n is the set of vectors of size n=|x⃗| produced as Cartesian product of E.

  * Solve the quantifier-free formula ψ output of step with an SMT solver that
decide union of quantifier-free theories of uninterpreted functions with equality 
and Presburger Arithmetics.


Correctness:
Although the above procedure is incomplete, we do have soundness.
This means if a dependence is determined unsatisfiable, it in fact 
is not a dependence.  However, if a dependence is determined
satisfiable at compile time, it could be that at runtime the actual values of index
arrays lead to the dependence not being satisfiable.
Since our procedure is conservatively correct, it is sound.


To show that the decidability procedure UNSAT_ψ is sound,
we need to show that if the original formula ψ is 
satisfiable, then so is the unquantified formula ψ',
	
    ψ∈ SAT ψ' ∈ SAT.

This is equivalent to
	
    ψ' ∉ SAT ψ∉ SAT.

Since universal quantification is being replaced with 
specific expression instantiations to create ψ',
ψ' is a potentially weaker set of constraints
than ψ.
This means that ψ' is a conservative approximation
of ψ.  As such, if ψ' is not satisfiable,
then ψ is not satisfiable.





§ SIMPLIFYING THE DEPENDENCES UTILIZING EQUALITIES

The finite instantiation proposed in Section <ref> can
prove many of the dependence relations to be unsatisfiable. However, some of
the relations remain satisfiable, thus requiring runtime checks. It is then
important to minimize the runtime cost by simplifying the dependence relations
as much as possible. In this section, we discuss one of such simplifications
utilizing additional equalities after finite instantiations.






 §.§ Discovering New Equality Constraints and Their Usefulness

Sometimes index array properties 
can help reduce the complexity of runtime inspectors
through introducing equalities to the dependence's constraints. 
The new equalities are discoverable after instantiating the
universally quantified assertions and combining those with other inequality
and equality relationships.
For instance, consider the following set of constraints; 
it is a satisfiable dependence that needs a runtime inspector with 
complexity of O(n^2) to traverse the space of values for i and i':

    (i ≤ i')  ∧  (f(i') ≤ f(i)) ∧ (0 ≤ i,i' <n).

And assume we also know following universally quantified rule 
about the uninterpreted function f (strict monotonicity):

    (∀ x_1,x_2), (x_1 < x_2)  (f(x_1) < f(x_2)).

With any universally quantified implication, 
if the left side of the implication is true, then the right side must 
be true to satisfy the assertion (i.e., p  q). It is also 
the case that the contrapositive is true (i.e., q  p). 
For this example, the negation of the right-hand side of the implication is 
f(x_2) ≤ f(x_1), which matches one of the constraints in the dependence.
Thus the negation of the left-hand side must be true and therefore 
x_2 ≤ x_1. With x_1 matching i and x_2 matching i', 
we find i' ≤ i. Thus an equality has been found:

    (i ≤ i' ∧ i' ≤ i)   i = i'

Using this equality we can iterate over either i or i'
in the inspector and calculate the other by taking advantage 
the equality. The runtime inspection would 
have complexity of only O(n).






 §.§ Finding Equalities Example: Left Cholesky


For a more realistic example from one of the benchmarks used
in the evaluation, consider a maybe satisfiable dependence from 
our Static Left Cholesky shown in Figure <ref>. Following 
dependence is coming from a read in S1 (lValues[l]), and a
write in S2 (lValues[j]):


    { [colNo] ⇒ [colNo']  :  ∃ j, i', l', (j = l') ∧ (colNo < colNo') 
    ∧ (0 ≤ colNo < n) ∧ (0 ≤ colNo' < n) 
    ∧ (lcolptr(pruneSet(i')) ≤ l' < lcolptr(pruneSet(i')+1)) 
    ∧ (prunePtr(colNo') ≤ i' < prunePtr(colNo'+1)) 
    ∧ (lcolptr(colNo) <  j < lcolptr(colNo+1)) }



An inspector for this dependence is 
shown in Figure <ref>.
We do not need loops for j and l'
in the inspector, because they can be projected out.
Note, index array prunePtr points to nonzeros in
the sparse matrix, ranging from 0 to number of nonzeros, nnz,
and n denotes the number of column.
The two loops, colNop and ip, combined are traversing 
all the nonzero values and hence have a combined complexity of nnz, 
followed by the colNo loop traversing over columns, n. 
Consequently, the complexity of this inspector is n(nnz).

The equality colNo = pruneSet(i') is found using the additional knowledge
that lcolptr is strictly monotonically increasing as demonstrated in the
following.

We have the following constraints in the original dependence:

    lcolptr(pruneSet(i')) <= l' < lcolptr(pruneSet(i')+1)
       ∧   j = lp ∧ lcolptr(colNo) < j < lcolptr(colNo+1),

which gives the following through transitivity:

    lcolptr(pruneSet(i')) < lcolptr(colNo+1)  
       ∧   lcolptr(colNo) < lcolptr(pruneSet(i')+1).


We have the following assertion:
	
    (∀ x_1,x_2)(lcolptr(x_1) < lcolptr(x_2)  x_1 < x_2)

where two instances x_1=pruneSet(i'),x_2=colNo+1 and x_1=colNo,x_2=pruneSet(i')+1 give new constraints:
   
    pruneSet(i')<colNo+1 ∧ colNo<pruneSet(i')+1
    ⇒     pruneSet(i')≤ colNo ∧ colNo≤ pruneSet(i')
    ⇒     colNo = pruneSet(i')


The optimized inspector based on new discoveries
is shown in Figure  <ref>.
We do not need loop for colNo, since we can get its values from 
pruneSet(i') based on
colNo = pruneSet(i').
This simplified inspector
has a complexity of (nnz),
 which is significantly better than the original n(nnz).














§ SIMPLIFYING THE DEPENDENCES UTILIZING SUPERSET RELATIONSHIP



Another way to deal with data dependence relations that cause complex
runtime analysis is to remove it from consideration by determining it is a
 subset of a less expensive relation.
Consider two dependence relations R1 and R2, and two iterations of the
outermost loop i and i'. If we can show that for all i and i' that are
dependent according to R2, the same pairs of i and i' are also dependent
according to R1, then it is sufficient to only test R1. We say that R1 is
a superset of R2, written R1 ⊇ R2, in such cases, and remove R2
from runtime check. Note that in the above definition, R1 may have more pairs
of outermost iterators that are dependent than R2.

Taking advantage of this redundancy can result in lower complexity runtime analysis. 
As an example, consider the  Incomple Cholesky code shown in
Figure <ref>. In section, we refer to an array access A at statement S as
A@S for brevity. One of the dependence tests is between the write
val[k]@S3 and the read val[m]@S3. This test is
redundant with the test between the write val[k]@S3 and the read
val[m]@S2. This is because an iteration of the i loop
that make the read from val[m] in S3 is guaranteed to access the same
memory location while executing the loop surrounding S2. Thus, the more
expensive check between accesses in S3 can be removed.



In this section, we describe our approach to identify redundant dependence
relations. The key challenge is to determine superset relations between two
dependence tests involving uninterpreted functions. We present two approaches
that cover important cases, and discuss possible extensions.




 §.§ Trivial Superset Relations

Given a polyhedral dependence relation, it is easy to characterize the pairs of
loop iterations that are dependent. All the indices that do not correspond to
the loop iterators in question can be projected out to obtain the set of
dependent iterations. These sets can be compared to determine if a dependence
test is subsumed by another.  In principle, this is what we do to check if a
dependence relation is redundant with another. However,  dependence relations
from sparse codes
may have variables passed as parameters to uninterpreted functions.
Such variables  cannot be projected out.
Thus, we employ
an approach based on similarities in the constraint systems. The trivial case
is when the dependence relation R1 is expressed with a subset of constraints
in another relation R2. If this is the case, then R1 can be said to be
superset equal to R2.

We illustrate this approach with the earlier example from Incomplete Cholesky.
We take two dependence relations, R1 between val[k]@S3
and val[m]@S2, and R2 between val[k]@S3 and
val[m]@S3. The relations—omitting the obviously common constraints 
for val[k]@S3—are:

    R1 = { [i,m,k,l] → [i',m'] : k = m' ∧ i < i' ∧ 0 ≤ i' < n ∧  col(i')+1 ≤ m' < col(i'+1) }
    
    R2 = { [i,m,k,l] → [i',m',k',l'] : k = m' ∧ i < i' ∧ 0 ≤ i' < n ∧  col(i')+1 ≤ m' < col(i'+1) 
    ∧  col(row(m'))≤ k' < col(row(m)+1) ∧ m'≤ l' < col(i+1)
    ∧  row(l')=row(k') ∧ row(l'+1)≤ row(k') }


Since R1 is expressed with a subset of constraints in R2, 
we may conclude that R1 ⊇ R2. 




 §.§ Superset Relation due to Overlapped Accesses



The trivial check is sufficient for many pairs of relations. However, some
relations require a more involved process. We use a different dependence
relation from Incomplete Cholesky (Figure <ref>) as an example of
such cases.  We consider the dependence relation R3 between val[k]@S3
and val[l]@S3 that is redundant with R1. This can be intuitively
observed from the code structure: the set of memory locations that may be
accessed by the read of val[l] when l=m, i.e., the first iteration of
the l loop, is exactly the same as the reads by val[m]@S2.  This
guarantees that even if the guard on S3 always evaluated to true, the
dependence between iterations of the i loop would be redundant with that
imposed by S2.


The constraints for R3 (omitting those for val[k]@S3) are as follows:

    R3 = { [i,m,k,l] → [i',m',k',l'] : k = l' ∧ i < i' ∧ 0 ≤ i' < n ∧  col(i')+1 ≤ m' < col(i'+1)
    ∧  col(row(m'))≤ k' < col(row(m)+1) ∧ m'≤ l' < col(i+1)
    ∧  row(l')=row(k') ∧ row(l'+1)≤ row(k') }




  * We first identify that k=m' in R1 is not a constraint in R3.

  * The equality k=m' has a similar (one side of the equality is the same) equation k=l' in R3.

  * The bounds on m' and l' are collected from the respective constraints.

  * Because the bound on m' subsumes that of l', and since k=m' was the
only constraint that was not in R3, we may conclude that  R1 ⊇ R3.


It is important to note that the bounds on l'—the set of values accessed in
the subset relation—can be conservative, i.e., may accesses, but the
bounds on m'—the set of values accessed in the superset relation—must be
exact. If both bounds represent may accesses, then the superset relation does
not hold. This is important for situations as illustrated in the example above,
where statements have data-dependent guards.



 §.§ Limitations and Extensions

Although the approach presented above was able to cover all the important cases
we encountered, it is by no means complete. The difficulty of manipulating
integer sets with uninterpreted function symbols have led us to work directly
on the constraints. This may cause our approach to miss some superset
relations, since the same relation can be expressed in many different ways. Adding some
form of normalization to the constraint system will help us avoid such pitfalls.

The overlapped iterator approach to finding a superset in 
Section <ref> was developed specifically for the problematic
data dependence relation R3.  Future work includes developing a more general
simplification approach based on this overlapping iterator concept.

In terms of scaling, there  is potentially a 
problem of selecting the pairs of dependence relations to
test for redundancy. We currently try all possible candidate pairs, which does
not pose a problem since a large number of dependence relations are filtered
out through unsatisfiability test described in
Section <ref>. However, selecting promising pairs to limit
the number of tests would be an useful extension.




§ IMPLEMENTATION

The data dependence analysis and simplification have been automated 
except for the superset simplification.  
This section summarizes the software packages the implementation relies on, 
discusses some important optimization to make our implementation scalable,

and compares the ISL-based implementation with that of an SMT solver.






 §.§ Software Description

The artifact for reproducing the results presented in this article
 can be found at the following public github repository:
https://github.com/CompOpt4Apps/Artifact-datadepsimplify-arXiv-July2018

We use three software packages to automate applying methods described 
in this paper: 
IEGenLib library <cit.>,
ISL library <cit.>,
and CHILL compiler framework <cit.>.  
CHiLL is a source-to-source compiler framework for composing and applying
high level loop transformations to improve the performance of nested loop 
written in C. We utilize CHILL to extracted the dependence relations
from the benchmarks.
ISL is a library for manipulating integer sets and relations that 
only contain affine constraints.
It can act as a constraint solver by testing the
emptiness of integer sets. It is also equipped with a function for detecting
equalities in sets and relations.  ISL does not support uninterpreted
functions, and thus cannot directly represent the dependence constraints in
sparse matrix code. 
IEGenLib is a set manipulation library that can  manipulate 
integer sets/relations that  contain uninterpreted function symbols. 
The IEGenLib library utilizes ISL for some of its fundamental functionalities. 
We have implemented detecting unsatisfiable dependences and finding 
the equalities utilizing the IEGenLib and ISL libraries.

The following briefly describes how the automation works.
First, we extract the dependence relations utilizing CHILL, and
store them in IEGenLib. The user defined index array properties 
are also stored in IEGenLib.
Next, the instantiation procedure is carried out in IEGenLib.
Then inside IEGenLib, the uninterpreted functions are removed by replacing each
call with a fresh variable, and adding constraints that encode functional
consistency <cit.>. 
Next, ISL can be utilized by IEGenLib to
find the empty sets, i.e, unsatisfiable relations.
Additionally, equality detection is available
as one of many operations supported by ISL
The finite instantiations described in
Section <ref> 
are intersections of the assertions with the dependence
relation. 








 §.§ Optimization

A straightforward approach to implementing the procedure in
Section <ref> would be to take the  quantifier-free
formula resulting from instantiation,
replace the uninterpreted functions, and directly pass it to ISL.
However, this approach does not scale to large numbers of instantiations.
An instantiated assertion is encoded as a union of two constraints (p
∨ q).  Given n instantiations, this approach introduces 2^n
disjunctions to the original relation, although many of the clauses may be
empty.  In some of our dependence relations, the value of n may exceed
1000, resulting in a prohibitively high number of disjunctions.  We have
observed that having more than 100 instantiations causes ISL to start
having scalability problems.

We apply an optimization to avoid introducing disjunctions when possible.
Given a set of instantiations, the optimization adds the instantiations to 
the dependence relation in two phases. 
The first phase only instantiates those that do not introduce disjunctions 
to the dependence relation. 
During this phase, we check if the antecedent is already part of 
the dependence constraint, and thus is always true. 
If this is the case, then q can be directly added to the dependence relation. 
We also perform the same for q p 
and add p to the dependence relation if q is always true. 
The second phase adds the remaining instantiations that introduce disjunctions.
This optimization helps reducing the cost of dependence testing in two ways:
(1) if the relation is unsatisfiable after the first phase, disjunctions are completely avoided; 
and 
(2) the second phase only instantiates the remainder, reducing the number of disjunctions.

If the dependence relation remains non-empty after the second phase, then the
relation is checked at runtime. All equalities in a relation is made explicit
before inspector code generation with ISL so that the code generator can take
advantage of the equalities to simplify the generated code.




 §.§ Contrasting SMT with ISL


SMT solvers are specialized for solving satisfiability problems expressed as a
combination of background theories. ISL is a library for manipulating integer
sets, and is specialized for the theory of Presburger arithmetic over integers.

The finite instantiation in Section <ref> is well-suited
for SMT solvers. In fact, SMT solvers are equipped with their own instantiation
algorithms that also work well for our dependence relations. However, SMT
solvers do not provide any equality relationships they might derive while answering
the satisfiability question.
Although it is possible to use SMT solvers to test if an equation is true for a
set of constraints, we cannot search for an equation given the constraints
(unless all candidates are enumerated—but there are infinite candidates in general). 

For our implementation, the choice was between adding finite instantiation to
ISL or adding equality detection to SMT solvers. We have chosen the former
option as it seemed simpler to do, and also because we are more familiar with
ISL.






§ EVALUATION OF UNSATISFIABILITY AND SIMPLIFICATION APPROACHES


In this section, we study the impact of our approach of utilizing domain information
about index arrays on the data dependence analysis of eight sparse kernels.
Our approach may help data dependence analysis in three ways: 
(1) The runtime check can be completely removed
if the dependences are proven unsatisfiable; 
(2) Deriving equalities from instantiated universally quantified
assertions about domain information 
can simplify dependences and reduce respected runtime check complexity; and
(3) Reducing all maybe satisfiable relations of a given code to
a set of dependence relations that encompass all potential dependences.
We do this by finding relations that are superset equal of other relations.
This can discard even more dependence relations that potentially
might need expensive runtime checks. 

We first describe the suite of numerical kernels that we have compiled to
evaluate our approach. Then we evaluate the impact of each step in our
approach, from the relevance of the index property assertions to the simplification
using superset relations. Finally, we report the complexity of inspectors
with and without our proposed simplifications.



 §.§ Numerical Algorithms in Benchmark


We have included some of the most popular sparse kernels 
in a benchmark suite: (1) The Cholesky factorization, Incomplete LU0 and
Incomplete Cholesky0, and the sparse triangular solver, which are commonly used
in direct solvers and as preconditioners in iterative solvers; (2)  sparse
matrix vector multiplication, and
Gauss-Seidel methods, which are often used in iterative solvers.
Table <ref> summarizes the benchmarks indicating which library each
algorithm came from and how the benchmark compares with the implementations in
existing libraries. 







We modified onr of the benchmarks, left Cholesky, to make temporary arrays
privatizable and to remove dynamic index array updates so that the compiler can
analyze the sparse code.

Left Cholesky:  
This code has following changes compared to 
a more common implementation in CSparse <cit.>:
(i) Privatization of temporary arrays: We analyzed dependences between 
reads and writes to temporary arrays to detect privatizable arrays.  This
can be challenging for a compiler to do with  
sparse codes since accesses to these arrays are irregular.
We set the values of these arrays 
to zero at the beginning of each loop so 
a compiler could identify them as privatizable.
(ii) Removal of dynamic index array updates: 
Previous data dependence analysis work 
focuses on cases where index arrays are not updated.
However, in some numerical codes, updating 
index arrays is a common optimization.
We refer to this as dynamic index array updates,
and it usually occurs 
when the nonzero structure of an output matrix
is modified in the sparse code during the computation. 
This would make 
dependence analysis very complicated for the compiler.
We removed dynamic index arrays by partially decoupling symbolic
analysis from the numerical code in these benchmarks.
Symbolic analysis here refers to terminology used in the numerical
computing community. Symbolic analysis uses the nonzero pattern of the matrix
to determine computation patterns and the sparsity of the resulting data.
To remove dynamic index array updates, we decouple symbolic analysis from
the code similar to the approach used by
 <cit.>.


Performance Impact:
The changes made to Left Cholesky do not have a 
noticeable effect on the code performance. Based on our 
experiments using five matrices[Problem1, rdb450l,
wang2, ex29, Chebyshev2] from the Florida Sparse 
Matrix Collection <cit.>  
the performance cost of these modifications is on average less 
than 10% than the original code. 




 §.§ Relevance of Index Array Properties

We have extracted the constraints to test for dependences that are carried by
the outermost loop for the sparse matrix codes in Table <ref>. A
total of 124 data dependences relations were collected from the benchmarks.
Of those 124, only 83 of them were unique, the repetition coming from accesses with 
same access indices in the same statements, or other situations.
Table <ref> summarizes the index array assertions relevant to the
benchmarks.



We are not claiming to have found all the array 
properties that exist either in our example suite nor in general. 
Also, we only consider dependence relations
for outermost loops, however, dependence relations can be extracted for other 
loop levels in a loop nest and can be used 
for vectorization 
and in other applications of dependence analysis.








 §.§ Detecting Unsatisfiability


In this section, we show the impact of using index array properties
to detect unsatisfiability for the relations collected 
from dependences from our benchmark suite.
To not conflate the impact of the index
array properties that we are evaluating with what 
traditional methods are capable of, we first apply functional consistency 
in the theory of Presburger arithmetic combined with 
uninterpreted functions  <cit.>.
This detects 12 dependences as unsatisfiable. 
Nevertheless, we must note that, all of the 12 dependences have inconsistencies in 
their affine parts and functional consistency does not help detect 
any more unsatisfiable relations; like the first two 
dependences from the Forward Solve CSR example in
Section <ref>. After detecting 12 out of 83
dependences as unsatisfiable we are left with 71 dependences to use in our
evaluation. We call these 71 dependences that are satsfiable just by looking 
at their affine constraints our baseline. 


Figure <ref> categorizes the complexity of an
inspector for each dependence into 7 different classes in total. 
In this figure, nnz is number of non-zeros, and for simplicity
n denotes the number of rows or columns of the input matrix.
The black bar, “baseline”, in each class
shows the baseline number of relations with that complexity in our suite.
The  bars show how many dependences would remain
after we instantiate certain index array properties. 
The last bar in each class, the red bar, 
shows the effect of adding all the information in combination.


The main observations from analyzing Figure <ref> are as
follows:
(1) Combining the array properties and non-domain information has 
the biggest impact and helps detect significantly more unsatisfiable 
dependences than any single property.
Combining all the index array properties helped us detect 45 out of 71 relations 
as unsatisfiable, with 26 remaining as maybe satisfiable.
(2)  Monotonicity has the highest impact on detecting unsatisfiable relations
when array properties are applied independently. 
(3) The Triangular Matrix property helped detect 3 relations 
when applied independently and 11 more in combination 
with Monotonicity (not obvious in the figure). 
This property helped us detect unsatisfiability in cases 
where Monotonicity was completely handicapped; see the 
nnz and nnz*n classes in Figure <ref>.





 §.§ Simplifying Inspector Complexity Utilizing Equalities


As stated in the previous section, 
instantiating index array properties results in
45 out of 71 dependence relations being detected as unsatisfiable.
At this point, without any further simplification,
to perform a partial parallelism transformation, 
inspectors are needed for the remaining 26 dependences. 
One question we can ask about those 26 inspector is whether
their complexity is even reasonable.  We consider a runtime dependence
analysis  complexity reasonable, if it is bound by the complexity of the 
original computation.
The computations would certainly
do much more operations compared to the analysis as 
numerical algorithms
usually call these computations several time for the same sparse
matrix nonzero structure.
Thus runtime data dependence analysis is reasonable 
if it is the same complexity as the original computation. 
Nonetheless, for numerical algorithms, 
it is common to aim for a runtime data dependence analysis
that is of O(nnz), 
where nnz is the number of nonezeros in the input. 

By instantiating index array properties with expressions
from the data dependences, it is also possible to derive
equalities between some of the iterators in the dependence.
These new useful equalities can be used to eliminate extra loops in 
the runtime inspector. 
Table <ref> shows that the additional equalities increases 
the number of dependence relations with reasonable complexities (≤kernel).  For instance, the Left Cholesky code have 4 high
complexity dependence relation left. As illustrated in
Section <ref>, the additional equalities can be used to reduce
the complexity of all those relations.  Finding equalities also help reduce the
complexity of 4 dependences for Incomplete Cholesky0 and 2 dependences of
Incomplete LU0 to become reasonable.

We should also mention that in addition to these 10 complexity reductions,
the complexity of another 4 dependence relations were reduced. However,
the complexity after simplification is still higher than the kernel, and hence
these simplifications are not visible in Table <ref>.









 §.§ Impact of Utilizing Superset Relationship


The superset relations we identify uncovers dependence relations that are
redundant. We can discard the dependence relations that are found to be
subsets of another and only generate runtime inspectors for remaining
relations. As shown in Table <ref>, this results in 
fewer dependence relations to be checked at runtime. Most notably, the 
number of runtime checks were reduced from 4 to 2 for Left Cholesky, 
and both of those dependences are less comlex than the original algorithm.

As discussed in Section <ref>, the superset relation
may reveal that a relation is redundant being subset of another relation with 
lower complexity. The Incomplete Cholesky kernel were left with 4 
expensive relations even after adding equalities. 
As you can see in Table <ref>, these relations are removed 
from runtime checks by identifying the superset relations. 
For Incomplete Cholesky kernel, we have found 2 relations with less than 
orginal algorithm complexity to be superset of all the dependences that
we need to have a runtime check for. The composed effect of our proposed
technique reduces the inspector cost to 2 or fewer inexpensive tests for all of
our kernels, except for the Incomplete LU.



 §.§ Putting It All Together

We have presented a series of techniques to simplify dependence relations with
the main motivation being automatic generation of efficient inspector code.
Our approach aims to simplify the dependence relations starting from array
properties that can be succinctly specified by the experts. We show that the
array properties can be used to automatically disprove a large number of
potential dependences, as well as reduce the complexity of remaining
dependences. Combined with a method for detecting redundancies in dependence
tests, we are able to generate efficient inspectors.

Table <ref> summarizes the impact of our proposed 
approach on inspector complexity. It is interesting to note that 
Incomplete LU0 is the only kernel left with expensive inspector 
(more complex than kernel). This case is discussed further in Section <ref>.







 §.§ Discussion: Limitations


Table <ref> demonstrates that
our method significantly reduces both the number of runtime checks and their complexity. 
Nonetheless, our approach is not free of limitations, which are discussed in this section.

Two of the original kernels include dynamic index arrays and temporary arrays that
require privatization. As discussed in Section <ref>,
these kernels can be preprocessed such that it can be accepted by our compiler. This preprocessing
is currently done manually.

Using the associativity of reductions is important for Forward Solve CSC
and Incomplete Cholesky0.  We do not automate the reduction detection in this
paper, as it is a complex task on its own.  It is common for compilers and
programming models, such as openMP, to provide pragma interfaces for
programmers to signal which update should be considered as a reduction. We have
followed the same approach.

Incomplete LU0 has two dependence relations that has higher complexity
than the kernel, even with domain information.  Related work by <cit.> presents
approximation techniques that reduce the inspector complexity for these high
complexity relation to nnz × (nnz/n).  Such approximation can
potentially result in loss of some parallelism.  Nevertheless,
 <cit.> show that the approximation of dependences does not
significantly affect the performance of the partial parallelism for this code.
We have not used approximations in our work, but it would be interesting to see
how the two approaches can be combined.





§ RELATED WORK


Array data dependence analysis has been used for a variety of applications,
including automatic parallelization <cit.>, locality optimization <cit.>,
communication generation, program slicing <cit.>, detecting race
conditions <cit.>,
and high-level synthesis <cit.>.
For sparse matrix codes, this analysis is made more difficult
due to indirection through index arrays, such that the source and sink
of dependences cannot be resolved until their values are available at runtime.
For these and other situations where dependences arise that cannot
be resolved until runtime, a number of techniques for compile time and
runtime dependence analysis have been developed.




 §.§ User-Provided Assertions


<cit.> exploit  user assertions about index arrays 
to increase the precision
of dependence testing. The assertions certify common
properties of index arrays, e.g., an index array can be a permutation
array, monotonically increasing, and monotonically
decreasing.
<cit.> present a compile time analysis
for determining index array properties, such as monotonicity.
They use the analysis results for  parallelization of sparse matrix
computations.

Our approach also uses these assertions, but in
addition we use more domain-specific assertions and provide
a way to automate the general use of such assertions.
In this paper, the idea of applying constraint instantiation of universally
quantified constraints as is done in SMT solvers to find unsatisfactory
dependences is novel and the assertions about index arrays we use
are more general.




 §.§ Proving Index Arrays Satisfy the Assertions

In this work, we assume that the assertions provided by the programmer is
correct. It is useful to verify the user-provided assertions by analyzing the
code that constructs the sparse matrix data structures. There is a
large body of work in abstract interpretation that address this problem.

The major challenge in verifying the assertion about programs that manipulate
arrays is the trade-off between scalability and precision. When there is a
large number of updates to an array, keeping track of individual elements do
not scale, but approximating the whole array as a single summary significantly
degrades the precision. Many techniques to verify/infer important properties
about array contents from programs have been developed, e.g.,
<cit.>.

In the work by <cit.>, the authors present an approach
for inferring shape invariants for matrices.

While their work does not deal with sparse matrices and index arrays,
it may help generate domain-specific assertions that we could
employ to show that the data dependences are unsatisfiable.


The main subject of our work - dependence tests - does not involve array
updates, since all the index arrays, which alter the control-flow and indexing
of the data arrays, are not updated. This makes the verification of the
assertions a closely related but orthogonal problem, which we do not address in
this paper.




 §.§ More General Quantifier Elimination Techniques


The area of SMT-solving is advancing at a significant pace; the webpage for
SMT-COMP[<http://smtcomp.sourceforge.net/2017/>] provides a list
of virtually all actively developed solvers, and how they fared in each theory
category.
As these solvers are moving into a variety of domains,
quantifier instantiation and elimination has become a topic of 
central interest.
Some of the recent work in this area are: E-matching <cit.>,
Model-Based <cit.>, Conflict-Based <cit.>, and
Counter-Example Guided <cit.>.


These efforts make it clear that
quantifier instantiation is challenging, and is an area of active development.
SMT solvers often rely on 
heuristic-driven instantiations to show unsat for difficult problems.
In this context,
our work can be viewed as heuristic instantiation where 
the heuristic is inspired by decidable fragments of the array theory.

Dependence constraints with universally quantified assertions
are related to the first order theory fragments
described by <cit.> as undecidable extensions
to their array theory fragment.
However,  <cit.> claim that the proofs for undecidability
of extension theories by <cit.> are incorrect,
and declare their decidability status as an open problem.
Regardless of whether the theory fragment that encompasses our
dependence constraints is decidable or not following is true: if we
soundly prove that a relation is unsatisfiable just with compile time
information, the unsatisfiability applies in general, and having
runtime information would not change anything. However, if a dependence
detected to be satisfiable just with compile time information, we need
to have runtime tests to see if it is actually satisfiable given runtime
information, and even if it is, run time tests would determine
for what values the dependence holds.





 §.§ Dependence Analysis for Full Parallelization


Some compilation techniques have been developed to extend the dependence analysis to
sparse, or non-affine programs <cit.>. These
techniques extend to non-affine programs of various forms: while loops,
polynomial expressions, function calls, data-dependent conditions, and
indirection arrays. The outcome of such analysis is often an approximation,
which is quite pessimistic for sparse computations involving indirection
arrays. The focus of our work is not to identify (approximated) dependences,
but to reduce the cost of runtime dependence checks by disproving potential
dependences as much as possible at compile-time.

The work by <cit.> also formulate the
problem in the theory of Presburger sets with uninterpreted functions.
However, they only allow affine expressions of unquantified variables as
indexing expressions to the function symbols, excluding some of the examples in
this paper.  They propose an analysis to identify conditions for a dependence
to exist through the use of gist operator that simplifies the constraint system
given its context. The result of this analysis may involve uninterpreted
functions, and can be used to query the programmer for their domain knowledge.
This is an interesting direction of interaction that complements our work.

Several runtime approaches focus on identifying loops
we denoted fully parallel whose iterations are
independent and can safely execute in 
parallel <cit.>
or speculatively execute in parallel while testing safety <cit.>.




 §.§ Dependence Analysis for Wavefront Parallelization


For sparse codes, even when loops carry dependences, the dependences 
themselves may be sparse, and it may be possible to execute some iterations of the
loop in parallel (previously denoted partially parallel.
The parallelism is captured in a task graph, and typically executed as a parallel wavefront.
A number of prior works write specialized code to derive this task graph as part 
of their application <cit.>
or with kernel-specific code generators <cit.>.
For example, Saltz and Rothbergs worked on manual parallelization of  
sparse triangular solver codes  in the 1990s <cit.>.
There is also more recent work on optimizing sparse triangular solver
NVIDIA GPUs and Intel's multi-core CPUs <cit.>.
Even though these manual optimizations have been successful
at achieving high performance in some cases, significant programmer
effort has to be invested for each of these codes and automating
these parallelization strategies can significantly reduce this effort. 


Other approaches automate the generation of inspectors that
find task-graph, wavefront or partial parallelism. <cit.>
and others <cit.> have developed efficient and parallel inspectors
that maintain lists of iterations that read and write each memory location.
By increasing the number of dependences found unsatisfiable,
the approach presented in this paper reduces
the number of memory accesses that would need to be tracked.
For satisfiable dependences, there is a tradeoff between inspecting
iteration space dependences versus maintaining data for each memory
access. That choice could be made at runtime. There are also
other approaches for automatic generation of inspectors that
have looked at simplifying the inspector by finding equalities,
using approximation, parallelizing the inspector, and applying
point-to-point synchronization to the executor <cit.>.





 §.§ Algorithm-Specific Data Dependence Analysis


An algorithm-specific approach to represent data dependences and optimize
memory usage of sparse factorization algorithms such as Cholesky
 <cit.> uses an elimination tree, but to the
best of our knowledge, this structure is not derived automatically from source code.
When factorizing a column of a sparse
matrix, in addition to nonzero elements of the input matrix new nonzero
elements, called fill-in, might be created. Since the sparse matrices
are compressed for efficiency, the additional fills during
factorization make memory allocation ahead of factorization difficult.
The elimination tree is used to predict the sparsity pattern of the L
factor ahead of factorization so the size of the factor can be computed <cit.> or predicted <cit.>,
and captures a potential parallel schedule of the tasks.
Prior work has investigated the applicability of the elimination
tree for dependence analysis for  parallel implementation
 <cit.>.
Some techniques such 
as <cit.> 
use the elimination tree for static scheduling while
others use it for runtime scheduling.





§ CONCLUSION


In this paper, we present an automated approach for 
showing sparse code data dependences are unsatisfiable or if
not reducing the complexity for later runtime analysis.
Refuting a data dependence brings benefits to many areas of sparse matrix code
analysis, including verification and loop optimizations such as parallelization, pipelining,
or tiling by
completely eliminating the high runtime costs of
deploying runtime dependence checking.
Additionally, when a dependence remains satisfiable, our 
approach of performing constraint instantiation within
the context of the Integer Set Library (ISL)
enables equalities and subset relationships to be 
derived that simplify the runtime complexity
of inspectors for a case study with wavefront parallelism.

Parallelization of these sparse numerical methods is 
an active research area today, but one where most current approaches
require manual parallelization.
It is also worth noting that without inspector complexity reduction, most
inspectors would timeout, thus underscoring the pivotal role of  the work in
this paper in enabling parallelization and optimization of sparse codes.
Our results are established over
71 dependences extracted from 8 sparse numerical methods.

































