Venmugil Elango

Fabrice Rastello

Louis-Noël Pouchet

J Ramanujam

P Sadayappan

On Characterizing the Data Movement Complexity of Computational DAGs for Parallel Execution

Keywords: red-blue pebble game, I/O complexity, I/O lower bound, Computational DAG, parallel architectures, data movement jeux des pions rouges et noirs, complexité I/O, nombre minimal de communications, graphe de tâches acyclique, architectures distribuées, communications

Technology trends are making the cost of data movement increasingly dominant, both in terms of energy and time, over the cost of performing arithmetic operations in computer systems. The fundamental ratio of aggregate data movement bandwidth to the total computational power (also referred to the machine balance parameter ) in parallel computer systems is decreasing. It is therefore of considerable importance to characterize the inherent data movement requirements of parallel algorithms, so that the minimal architectural balance parameters required to support it on future systems can be well understood. In this paper, we develop an extension of the well-known red-blue pebble game to develop lower bounds on the data movement complexity for the parallel execution of computational directed acyclic graphs (CDAGs) on parallel systems. We model multi-node multi-core parallel systems, with the total physical memory distributed across the nodes (that are connected through some interconnection network) and in a multi-level shared cache hierarchy for processors within a node. We also develop new techniques for lower bound characterization of non-homogeneous CDAGs. We demonstrate the use of the methodology by analyzing the CDAGs of several numerical algorithms, to develop lower bounds on data movement for their parallel execution.

Introduction

Recent technology trends have resulted in much greater rates of improvement in computational processing rates of processors than the bandwidths for data movement across nodes or within the memory/cache hierarchies within nodes in a parallel system. This mismatch between maximum computational rate and peak memory bandwidth means that data movement and communication costs of an algorithm will be increasingly dominant determinants of performance. Although hardware techniques for data pre-fetching and overlapping of computation with communication can alleviate the impact of memory access latency on performance, the mismatch between maximum computational rate and peak memory bandwidth is much more fundamental; the only solution is to limit the total rate of data movement between components of a parallel system to rates that can be sustained by the interconnects at different components and levels of a parallel computer system.

It is therefore of considerable importance to develop techniques to characterize lower bounds on the data movement complexity of parallel algorithms. We address this problem in this paper. We formalize the problem by developing a parallel extension of the red-blue pebble game model introduced by Hong and Kung in their seminal work [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF] on characterizing the data access complexity (called I/O complexity by them) for sequential execution of computational directed acyclic graphs (CDAGs). Our extended pebble game abstracts data movement in scalable parallel computers today, which are comprised of multiple nodes interconnected by a high-bandwidth interconnection network, with each node containing a number of cores that share a hierarchy of caches and the node's physical main memory.

In contrast to some other prior efforts that have modeled lower bounds for data movement in parallel computations, we focus on relating data movement lower bounds to the critical architectural balance parameter of the ratio of peak data movement bandwidth (in GBytes/sec) to peak computational throughput (in GFLOPs) at different levels of a parallel system. We develop techniques for deriving lower bounds for data movement for CDAGs under the parallel red-blue pebble game, and use the techniques to analyze a number of numerical algorithms. Interesting insights are provided on architectural bottlenecks that limit the performance of the algorithms.

This paper makes several contributions:

• It develops an extension of the red-blue pebble game that effectively models essential characteristics of scalable parallel computers with multi-level parallelism; (i) multiple nodes with local physical memory that are interconnected via a high-speed interconnection network like Infiniband or a custom interconnect (e.g., IBM BlueGene system [START_REF] Blue | The ibm blue gene project[END_REF], or Cray XE6 [10]), and (ii) many cores at each node, that share a hierarchy of caches and the node's physical main memory. • It develops a lower bound analysis methodology that is effective for analysis of nonhomogeneous CDAGs using a decomposition approach. • It develops new parallel lower-bounds analysis for a number of numerical algorithms.

• It presents insights into implications on different architectural parameters in order to achieve scalable parallel execution of the analyzed algorithms.

2 Background

Computational Model

The model of computation we use is a computational directed acyclic graph (CDAG), where computational operations are represented as graph vertices and the flow of values between operations is captured by graph edges. Two important characteristics of this abstract form of representing a computation are that [START_REF] Aggarwal | A model for hierarchical memory[END_REF] there is no specification of a particular order of execution of the op-erations: although the program executes the operations in a specific sequential order, the CDAG abstracts the schedule of operations by only specifying partial ordering constraints as edges in the graph; [START_REF] Aggarwal | The input/output complexity of sorting and related problems[END_REF] there is no association of memory locations with the source operands or result of any operation. We use the notation of Bilardi & Peserico [START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF] to formally describe CDAGs. We begin with the model of CDAG used by Hong & Kung:

Definition 1 (CDAG-HK) A computational directed acyclic graph (CDAG) is a 4-tuple C = (I, V, E, O) of finite sets such that: (1) I ⊂ V is the input set and all its vertices have no incoming edges; (2) E ⊆ V × V is the set of edges; (3) G = (V, E) is a directed acyclic graph;(4) V -I is called the operation set and all its vertices have one or more incoming edges; (5) O ⊆ V is called the output set.

The Red-Blue Pebble Game

Hong & Kung used this computational model in their seminal work [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]. The inherent I/O complexity of a CDAG is the minimal number of I/O operations needed while optimally playing the Red-Blue pebble game. This game uses two kinds of pebbles: a fixed number of red pebbles that represent the small fast local memory (could represent cache, registers, etc.), and an arbitrarily large number of blue pebbles that represent the large slow main memory. Starting with blue pebbles on all inputs nodes in the CDAG, the game involves the generation of a sequence of steps to finally produce blue pebbles on all outputs. A game is defined as follows.

Definition 2 (Red-Blue pebble game [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]) Given a CDAG C = (I, V, E, O) such that any vertex with no incoming (resp. outgoing) edge is an element of I (resp. O), S red pebbles and an arbitrary number of blue pebbles, with a blue pebble on each input vertex. A complete game is any sequence of steps using the following rules that results in a final state with blue pebbles on all output vertices: R1 (Input) A red pebble may be placed on any vertex that has a blue pebble (load from slow to fast memory), R2 (Output) A blue pebble may be placed on any vertex that has a red pebble (store from fast to slow memory), R3 (Compute) If all immediate predecessors of a vertex of V -I have red pebbles, a red pebble may be placed on that vertex (execution or "firing" of operation), R4 (Delete) A red pebble may be removed from any vertex (reuse storage).

The number of I/O operations for any complete game is the total number of moves using rules R1 or R2, i.e., the total number of data movements between the fast and slow memories. The inherent I/O complexity of a CDAG is the smallest number of such I/O operations that can be achieved, among all possible valid red-blue pebble games on that CDAG. The optimal red-blue pebble game is a game achieving this minimal number of I/O operations.

S-partitioning for Lower Bounds on I/O Complexity

This red-blue pebble game provides an operational definition for the I/O complexity problem. However, it is not practically feasible to generate all possible valid games for large CDAGs. Hong & Kung developed a novel approach for deriving I/O lower bounds for CDAGs by relating the red-blue pebble game to a graph partitioning problem defined as follows.

Definition 3 (S -partitioning of CDAG [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]) Given a CDAG C. An S-partitioning of C is a collection of h subsets of V such that:

P1 ∀i = j, V i ∩ V j = ∅, and h i=1 V i = V P2 there is no circuit between subsets Inria P3 ∀i, ∃D ∈ Dom(V i ) such that |D| ≤ S P4 ∀i, |Min(V i )| ≤ S
where a dominator set of V i , D ∈ Dom(V i ) is a set of vertices such that any path from I to a vertex in V i contains some vertex in D; the minimum set of V i , Min(V i ) is the set of vertices in V i that have all its successors outside of V i ; and |Set| is the cardinality of the set Set. We say that there is a circuit between two sets V i and V j , if there is an edge from any vertex in V i to a vertex in V j and vice-versa.

Hong & Kung showed a construction for a 2S -partition of a CDAG, corresponding to any complete red-blue pebble game on that CDAG using S red pebbles, with a tight relationship between the number of vertex sets h in the 2S -partition and the number of I/O moves q in the pebble-game shown next.

Theorem 1 (Pebble game, I/O and 2S -partition [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]) Any complete calculation of the redblue pebble game on a CDAG using at most S red pebbles is associated with a 2S-partition of the CDAG such that S × h ≥ q ≥ S × (h -1), where q is the number of I/O moves in the game and h is the number of subsets in the 2S-partition.

The tight association from the above theorem between any pebble game and a corresponding 2S -partition provides the following key lemma that served as the basis for Hong & Kung's approach to deriving lower bounds on the I/O complexity of CDAGs.

Lemma 1 (Lower bound on I/O [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]) Let H(2S) be the minimal number of vertex sets for any valid 2S-partition of a given CDAG (such that any vertex with no incoming -resp. outgoing -edge is an element of I -resp. O). Then the minimal number Q of I/O operations for any valid execution of the CDAG is bounded by: Q

≥ S × (H(2S) -1)
This key lemma has been useful in proving I/O lower bounds for several CDAGs [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF] by reasoning about the maximal number of vertices that could belong to any vertex-set in a valid 2S -partition.

The following corollary is generally useful while deriving I/O lower bounds through 2Spartitioning.

Corollary 1 Let U(2S) be the largest vertex-set of any valid 2S-partition of a given CDAG C = (I, V, E, O). Let V ′ = V \ I. Then the minimal number Q of I/O operations for any valid execution of the CDAG is bounded by: Q

≥ S × |V ′ | |U(2S)| -1
3 Parallel Red-Blue-White Pebble Game Application codes are typically constructed from a number of sub-computations using the fundamental composition mechanisms of sequencing, iteration and recursion. For instance, the conjugate gradient method, described in Sec. 5.2, consists of sequence of sparse matrix-vector product, vector dot-product and SAXPY operations, for every iteration. Applying the I/O lower bounding techniques directly on the CDAG of such composite application codes can produce very weak lower bounds. For instance, consider the following code segment.

1 I n p u t s : p, q, r, s : V e c t o r s o f s i z e N 2 Output : sum : S c a l a r 3

A = p × q T 4 B = r × s T 5 C = AB 6 sum = N i=1 N j=1 Cij
The computational complexity of this calculation can be simply obtained by adding together the computational costs of the constituent steps, i.e., N 2 + N 2 + 2N 3 + N 2 arithmetic operations. In contrast, the data movement complexity for this computation cannot so simply be obtained by adding together the data movement lower bounds for the individual steps. Let us consider data movement costs in a two-level memory hierarchy with unbounded main memory and a limited number of words (S) in fast storage -this might represent the number of registers in the processor, or scratchpad memory or cache memory. It is known [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF][START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF][START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF] that an asymptotic lower bound on data movement between (arbitrarily large) slow memory and fast memory for matrix multiplication of N × N matrices is N 3 /2 √ 2S. An outer-product of two vectors of size N requires 2N input operations from slow memory and output of the N 2 results back to slow memory, i.e., total I/O of 2N +N 2 , independent of the fast memory capacity S. Similarly, the last step has a data movement complexity of N 2 + 1 I/O operations between slow and fast memory. But a lower bound on the data movement complexity of the total calculation cannot be obtained by simply adding together contributions for the steps. It is not even possible to assert that the maximum among them is a valid lower bound on the data movement complexity of the total calculation. The reason is that data from a previous step could possibly be passed to a later step in fast storage without having to be stored in main memory. With 4N + 4 fast memory locations, it is feasible to perform the above computation with a total of only 4N + 1 I/O operations, 4N to bring in the four input vectors into fast memory, and repeatedly recompute elements of A and B to contribute to an element of C, and when ready, accumulate it into sum. The I/O complexity of the composite multi-step computation is thus lower than that of the matrix multiply step contained in it. This motivates us to split the CDAG based on individual sub-computations, determine the lower bound for each sub-CDAG separately, and finally compose the result to obtain the I/O lower bound of the whole computation. However, using the original red/blue pebble game model of Hong & Kung, as elaborated below, it is not feasible to analyze the I/O complexity of sub-computations and simply combine them by addition.

The Hong & Kung red/blue pebble game model places blue pebbles on all CDAG vertices without predecessors, since such vertices are considered to hold inputs to the computation, and therefore assumed to start off in slow memory. Similarly, all vertices without successors are considered to be outputs of the computation, and must have blue pebbles at the end of the game. If the vertices of a CDAG corresponding to a composite application are disjointly partitioned into sub-DAGs, the analysis of each sub-DAG under the Hong & Kung red/blue pebble game model will require the initial placement of blue pebbles on all predecessor-free vertices in the sub-DAG, and final placement of blue pebbles on all successor-free vertices in the sub-DAG. The optimal pebble game for each sub-DAG will require at least one load (R1) operation for each input and a store (R2) operation for each output. But in playing the red/blue pebble game on the full composite CDAG, clearly it may be possible to pass values in a red pebble between vertices in different sub-DAGs, so that the I/O complexity is less than the sum of the I/O costs for the optimal games for each sub-DAG. In fact, it is not even possible to assert that the maximum among the I/O lower bounds for sub-DAGs of a CDAG is a valid lower bound for the composite CDAG.

In order to enable such decomposition, a modified game called the Red-Blue-White pebble game [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF] was defined, with the following changes to the Hong & Kung pebble game model (the Red-Blue-White pebble game is formally defined in Sec. vertices without predecessors must be input vertices, and all vertices without successors must be output vertices, the RBW model allows flexibility in indicating which vertices are labeled as inputs and outputs. In the modified variant of the pebble game, predecessor-free vertices that are not designated as input vertices do not have an initial blue pebble placed on them. However, such vertices are allowed to fire using rule R3 at any time, since they do not have any predecessor nodes without red pebbles. Vertices without successors that are not labeled as output vertices do not require placement of a blue pebble at the end of the game. However, all compute vertices in the CDAG are required to have fired for any complete game. 2. Prohibition of multiple evaluations of compute vertices: The RBW game disallows recomputation of values on the CDAG, i.e., each non-input vertex is only allowed to evaluate once using rule R3. Several other efforts [START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF][START_REF] Ballard | Graph expansion and communication costs of fast matrix multiplication[END_REF][START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF][START_REF] Scquizzato | Communication lower bounds for distributed-memory computations[END_REF][START_REF] Ranjan | Strong I/O lower bounds for binomial and FFT computation graphs[END_REF][START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF][START_REF] Savage | Models of computation -exploring the power of computing[END_REF][START_REF] Savage | Cache-optimal algorithms for option pricing[END_REF][START_REF] Cook | An observation on time-storage trade off[END_REF][START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF][START_REF] Ranjan | Upper and lower I/O bounds for pebbling rpyramids[END_REF][START_REF] Ranjan | Vertex isoperimetric parameter of a computation graph[END_REF] have also imposed such a restriction on the pebble game model. While such a model is indeed more restrictive than the original Hong & Kung model, the restriction in the model enables the development of techniques to form tighter lower bounds [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF].

The Red-Blue-White Pebble Game

Definition 4 (Red-Blue-White (RBW) pebble game) Given a CDAG C = (I, V, E, O), S red pebbles and an arbitrary number of blue and white pebbles, with a blue pebble on each input vertex, a complete game is any sequence of steps using the following rules that results a final state with white pebbles on all vertices and blue pebbles on all output vertices: R1 (Input) A red pebble may be placed on any vertex that has a blue pebble; a white pebble is also placed along with the red pebble, unless the vertex already has a white pebble on it. R2 (Output) A blue pebble may be placed on any vertex that has a red pebble. R3 (Compute) If a vertex v does not have a white pebble and all its immediate predecessors have red pebbles on them, a red pebble along with a white pebble may be placed on v. R4 (Delete) A red pebble may be removed from any vertex (reuse storage).

In the modified rules for the RBW game, all vertices are required to have a white pebble at the end of the game, thereby ensuring that the entire CDAG is evaluated. Non-input vertices without predecessors do not have an initial blue pebble on them, but they are allowed to fire using rule R3 at any time -since they have no predecessors, the condition in rule R3 is trivially satisfied. But if all successors of such a node cannot be fired while maintaining a red pebble, "spilling" and reloading using R2 and R1 is forced because the vertex cannot be fired again using R3.

Definition 3 is adapted to this new game so that Theorem 1 and thus Lemma 1 can hold for the RBW pebble game.

Definition 5 (S -partitioning of CDAG -RBW pebble game) Given a CDAG C. An Spartitioning of C is a collection of h subsets of V -I such that:

P1 ∀i = j, V i ∩ V j = ∅, and h i=1 V i = V -I P2 there is no circuit between subsets P3 ∀i, |In(V i )| ≤ S P4 ∀i, |Out(V i )| ≤ S where the input set of V i , In(V i ) is the set of vertices of V \ V i that have at least one successor in V i ; the output set of V i , Out(V i )
is the set of vertices of V i also part of the output set O or that have at least one successor outside of V i .

The proof of Theorem 1 under the RBW pebble game is provided in [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF].

RR n°8522

For (sub-)graphs without input/output sets, the application of S-partitioning will however lead to a trivial partition with all vertices in a single set (e.g., h = 1). A careful tagging of vertices as virtual input/output nodes will be required for better I/O complexity estimates, as described below.

Decomposition

Definition 4 allows the partitioning of a CDAG C into sub-CDAGs C 1 , C 2 , . . . , C p , to compute lower bounds on the I/O complexity of each sub-CDAG IO(C 1 ), IO(C 2 ), . . . , IO(C p ) independently and simply add them to bound the I/O complexity of C. This is stated in the following decomposition theorem, whose proof may be found in [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF].

Theorem 2 (Decomposition) Let C = (I, V, E, O) be a CDAG. Let V 1 , V 2 , . . . , V p be an arbitrary (not necessarily acyclic) disjoint partitioning of V (i = j ⇒ V i ∩ V j = ∅ and 1≤i≤p V i = V ) and C 1 , C 2 , . . . , C p be the induced partitioning of C (I i = I ∩ V i , E i = E ∩ V i × V i , O i = O ∩ V i ). Then 1≤i≤p IO(C i ) ≤ IO(C). In particular, if Q i is a lower bound on the I/O cost of C i , then 1≤i≤p Q i is a lower bound on the I/O cost of C.
We state the following corollary and theorem, which are useful in practice for deriving tighter lower bounds. The complete proofs can be found in [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF].

Corollary 2 (Input/Output Deletion) Let C and C ′ be two CDAGs: C ′ = (I ∪ dI, V ∪ dI ∪ dO, E ′ , O ∪ dO), C = (I, V, E ′ ∩ V × V, O).
Then IO(C ′ ) can be bounded by a lower bound of IO(C) as follows:

IO(C) + |dI| + |dO| ≤ IO(C ′ ) (1)
There are cases where separating input/output vertices leads to very weak lower bounds. This happens when input vertices have high fan out such as for matrix-multiplication: if we consider the CDAG for matrix-multiplication and remove all input and output vertices, we get a set of independent chains that can each be computed with no more than 2 red pebbles. To overcome this problem, the following theorem allows us to compare the I/O of two CDAGs: a CDAG C ′ = (I ′ , V, E, O ′ ) and another C = (I, V, E, O) built from C ′ by just transforming some vertices without predecessors into input vertices, and some others into output nodes so that I ′ ⊂ I and O ′ ⊂ O. In contrast to the prior development above, instead of adding/removing input/output vertices, here we do not change the vertices of a CDAG but instead only change the labeling (tag) of some vertices as inputs/outputs in the CDAG. So the CDAG remains the same, but some input/output vertices are relabeled as standard computational vertices, or vice-versa.

Theorem 3 (Input/Output (Un)Tagging -RBW) Let C and C ′ be two CDAGs of the same DAG G = (V, E): C = (I, V, E, O), C ′ = (I ∪ dI, V, E, O ∪ dO).
Then, IO(C) can be bounded by a lower bound on IO(C ′ ) as follows (tagging):

IO(C ′ ) -|dI| -|dO| ≤ IO(C) (2)
Reciprocally, IO(C ′ ) can be bounded by a lower bound on IO(C) as follows (untagging):

IO(C) ≤ IO(C ′ ) (3)
Some algorithms will benefit from decomposing their CDAGs into non-disjoint vertex sets. For instance, when we have computations that are surrounded by an outer time loop, a common technique to derive their lower bound is to decompose the CDAG, where vertices computed Inria during each outer loop iteration are placed in separate sub-CDAGs. In such cases, when the vertices, V , computed in iteration t are used as inputs for iteration t + 1, by placing V in the sub-DAGs corresponding to both iterations t and t + 1, we could obtain a lower bound that is tighter by atleast a constant factor.

Theorem 4 (Non-disjoint Decomposition) Consider an optimal game P of C with S red pebbles. We let Q L1 be the number of R1 transitions (loads) in C associated to a vertex of C -[D x + x]. We let Q S1 be the number of R2 transitions (stores) in C associated to a vertex of C -D x . We let Q 2 be the number of R1 and R2 transitions (loads/stores) in C associated to a vertex in D x . We have that

IO S (C) >= Q L1 + Q 2 + Q S1 .
Proof. The idea of the proof is to show that Q L1 + Q S1 >= IO S+1 (C) and that Q2 >= IO S (C2).

Let us start with

Q L1 + Q S1 >= IO S+1 C.
We consider the restriction of C to the vertex of C 1 . This is not a valid game for C 1 yet, as the predecessors of any vertex in In(D x ) are not the same in C 1 than in C. But we have one more red pebble that we dedicate to stay on x. As soon as it is computed: we remove any R1 (loads) and R4 (delete) transitions associated to vertex x. This gives a valid game for C1:

• as C -D x is a sub-graph of C all transitions associated to a vertex of C -[D x + x] plus the transition R3 (compute) of x are valid (this part of the game has been unchanged). • for a vertex in D x the only kept transitions are R3 / R2(compute / store) and is valid as all its predecessors are in C -D x which associated transitions are unchanged (apart from x which keeps a red pebble as soon as it is computed). The cost of this valid game (with S + 1 red pebbles) for C1 is Q L1 + Q S1 which proves the inequality.

Let us now prove that

Q L2 + Q S2 >= IO S (C2).
We consider the restriction of C to the vertex of C 2 = D x . This is a valid game for C 2 of cost Q2 which proves the second inequality.

Min-Cut for I/O Complexity Lower Bound

In [START_REF] Elango | Data access complexity: The red/blue pebble game revisited[END_REF], we developed an alternative lower bounding approach. It was motivated from the observation that the Hong & Kung 2S-partitioning approach does not account for the internal structure of a CDAG, but essentially focuses only on the boundaries of the partitions. In contrast, the min-cut based approach captures internal space requirements using the abstraction of wavefronts. This section describes the approach. Definitions: We first present needed definitions. Given a graph G = (V, E), a cut is defined as any partition of the set of vertices V into two parts S and T = V -S. An s-t cut is defined with respect to two distinguished vertices s and t and is any (S, T ) cut satisfying the requirement that s ∈ S and t ∈ T . Each cut defines a set of cut edges (the cut-set), i.e., the set of edges (u, v) where u ∈ S and v ∈ T . The capacity of a cut is defined as the sum of the weights of the cut edges. The minimum cut problem (or min-cut) is one of finding a cut that minimizes the capacity of the cut. We define vertex u as a cut vertex with respect to an (S, T ) cut, as a vertex u ∈ S that has a cut edge incident on it. A related problem of interest for this paper is the vertex min-cut problem which is one of finding a cut that minimizes the number of cut vertices.

We consider a convex cut (S x , T x ) associated to x as follows: S x includes x ∪ Anc(x); T x includes Desc(x); in addition, S x and T x must be constructed such that there is no edge from T x to S x . With this, the sets S x and T x partition the graph G into two convex partitions. We define the wavefront induced by (S x , T x ) to be the set of vertices in S x that have at least one outgoing edge to a vertex in T x .

Schedule Wavefront: Consider a pebble game instance P that corresponds to some scheduling (i.e., execution) of the vertices of the graph G = (V, E) that follows the rules R1-R4 of the Red-Blue-White pebble game (see Definition 4 in Sec. 3.1). We view this pebble game instance as a string that has recorded all the transitions (applications of pebble game rules). Given P, we define the wavefront W P (x) induced by some vertex x ∈ V at the point when x has just fired (i.e., a white pebble has just been placed on x) as the union of x and the set of vertices u ∈ V that have already fired and that have an outgoing edge to a vertex v ∈ V that have not fired yet. Viewing the instance of the pebble game P as a string, W P (x) is the set of vertices x and those white-pebbled vertices to the left of x in the string associated with P that have an outgoing edge in G to not-white-pebbled vertices that occur to the right of x in P. With respect to a pebble game instance P, the set W P (x) defines the memory requirements at the time-stamp just after x has fired. Correspondence with Graph Min-cut Note that there is a one-to-one correspondence between the wavefront W P (x) induced by some vertex x ∈ V and the (S x , T x ) partition of the graph G. For a valid convex partition (S x , T x ) of G, we can construct a pebble game instance P in which at the time-stamp when x has just fired, the subset of vertices of V that are white pebbled exactly corresponds to S x ; the set of fired (white-pebbled) nodes that have a successor that is not white-pebbled constitute a wavefront W P (x) associated with x. Similarly, given wavefront W P (x) associated with x in a pebble game instance P, we can construct a valid (S x , T x ) convex partition by placing all white pebbled vertices in S x and all the non-white-pebbled vertices in T x .

A minimum cardinality wavefront induced by x, denoted W min G (x) is a vertex min-cut that results in an (S x , T x ) partition of G defined above. We define w max G as the maximum value over the size of all possible minimum cardinality wavefronts associated with vertices, i.e., define

w max G = max x∈V W min G (x) .
Lemma 2 Let C = (∅, V, E, O) be a CDAG with no inputs. For any x ∈ V , 2 W min G (x) -S ≤ IO(C). In particular, 2 (w max G -S) ≤ IO(C).

Parallel Red-Blue-White (P-RBW) Pebble Game

In this section, we extend the RBW pebble game to the parallel environment. P-RBW assumes that multiple nodes are connected in a distributed environment, with each node containing multiple cores. Further, the memory within each machine is organized in a hierarchical way.

More formally, we have: (1) N L main memories (storage of level L) connected through ethernet;

(2) P processors, each of them having exactly S 1 registers (storage of level 1); (3) for each level (1 < l < L), N l overall caches of size S l each; (4) one given cache of level l has a unique (parent) cache of level l + 1 to which it is connected. We consider that the bandwidth between a storage of level l and its children of level l -1 is shared between all its children. In other words, the I/Os of the P l = P/N l processors associated to a given level l storage instance are to be done sequentially. Note that those P/N l processors have S l-1 × N l-1 /N l storage available at level l -1. Fig. 1 illustrates this setup.

Definition 6 (Parallel RBW (P-RBW) pebble game) Let C = (I, V, E, O) be a CDAG. Given for each level

1 ≤ l ≤ L, N l ×S l number of red pebbles of different shades R 1 l , R 2 l , • • • , R N l l
, respectively, and unlimited blue and white pebbles, with a blue pebble on each input vertex, a complete game is any sequence of steps using the following rules that results in a final state with white pebbles on all vertices and blue pebbles on all output vertices: Inria L can be placed on any vertex that has a blue pebble; a white pebble is also placed along with the shade of red pebble, unless the vertex already has a white pebble on it. R2 (Output) A blue pebble can be placed on any vertex that has a level-L pebble on it. R3 (Remote get) A level-L pebble, R i L can be placed on any vertex that has another level-L shade pebble R j L . R4 (Move up) For 1 ≤ l < L, a level-l red pebble, R i l can be placed on any vertex that has a level-l + 1 pebble R j l+1 where R i l is in a cache that is a child of the cache that holds R j l+1 . R5 (Move down) For 1 < l ≤ L, a level-l red pebble, R j l can be placed on any vertex that has a level-l -1 pebble R i l-1 where R i l-1 is in a cache that is a child of the cache that holds R j l . R6 (Compute) If a vertex v does not have a white pebble and all its immediate predecessors have level-1 red pebbles on them, then a level-1 red pebble R p 1 along with a white pebble may be placed on v; here p is the index of the that comuptes vertex v. R7 (Delete) Any shade of red pebble may be removed from any vertex (reuse storage).

I/O Lower Bound for parallel machines

In this section, we provide the necessary tools to analyze the lower bounds for the parallel case. In particular, we consider two distinct cases of data movement:

1. data movement along the memory hierarchy within a processor, which we call the vertical data movement; 2. data movement across processors, which we call the horizontal data movement.

I/O Lower Bound for Vertical Data Movement

The hierarchical memory can enforce either the inclusion or exclusion policy. In case of inclusive hierarchical memory, when a copy of a value is present at a level-l, it is also maintained at all the levels l + 1 and higher. These values may or may not be consistent with the values held at the lower levels. The exclusive cache, on the other hand, does not guarantee that a value present in the cache at level-l will be available at the higher levels ≥ l + 1. The following result is derived for the inclusive case. But, they also hold true for the exclusive case, where the difference lies only in the number of red pebbles that we consider in the corresponding two-level pebble game.

Theorem 5 (Vertical I/O Cost)

Let C = (I, V, E, O) be a CDAG. Consider any valid P-RBW game on C; for this valid game, consider the level-l storage j with the maximum number of R5 transitions placing a R j l shade red pebble. The corresponding amount of move down transitions to this shade is at least IO 1 (C, S l-1 × N l-1 )/N l , where IO 1 (C, S) is the I/O lower bound of C for a single processor with local memory of size S.

Proof. Consider a P-RBW game of C that minimizes the overall amount of I/O between levels k < l and level l storage. This amount of I/O will be bounded by IO 1 (C, S l-1 ×N l-1 ). Consider one of the N l caches with the maximum amount of I/O. It will be bounded by IO 1 (C, S l-1 × N l-1 )/N l .

Theorem 6 (Vertical I/O Cost) Let C = (I, V, E, O) be a CDAG. Consider any valid P-RBW game on C; for this valid game, consider the level-l storage j with the maximum number of R5 transitions placing a R j l shade red pebble. The corresponding amount of move down transitions to this shade is at least

[|V |/(U (C, 2S l-1 ) × N l ) -N l-1 /N l ] × S l-1 ≈ |V |×S l-1 U (C,2S l-1 )×N l where |V | is the total amount of work, U (C, 2S) is the largest 2S-partition of the CDAG C.
Proof. Consider a P-RBW game of C that minimizes the overall amount of I/O between levels k < l and level-l storage. Consider the group of P/N l processors that do the more computation in this game. They do at least |V |/N l amount of work. Let us consider the partition of those P/N l processors into N l-1 /N l sets of P/N l-1 processors that share the same level-l -1 storage unit. Each set of processor (that we denote

P i with 0 < i ≤ N l-1 /N l ) does at least α i × |V | N l
amount of work where i α i = 1. We let V i be the subset of nodes of C fired by P i .

Let us denote S l-1 by S to simplify the notations. The goal is to show that each P i performs at least |V | i /U (C, 2S) -1 × S I/O to its level-l storage where U (C, 2S) is the largest 2Spartition (RBW pebble game) of CDAG C. Consider an RBW game of C with S red pebbles. Consider the partitioning of the game into C 1 , • • • , C h used in the proof of Theorem 1 for RBW. We let V i j be the set of vertices of

V i fired in C j ( V i j = V i ; V i j ∩ V ′i j = ∅ for j = j ′ ).
With the usual reasoning we can prove that |In(V i j )| ≤ 2S and |Out(V i j )| ≤ 2S ie each V i j is a 2S-partition of C. Thus for each j, |V i j | ≤ U (C, 2S). Now from a valid P-RBW game, we can build a valid RBW game where the restriction to V i matches the P-RBW game. By construction, each V i j is associated to at least S I/O to level-l storage in the P-RBW game. Thus the total amount of I/O for P i is at least

|V i |/|V i j | -1 × S ≥ W i /U (C, 2S) -1 × S.
If we sum up the I/O of each set of processors with our level-l storage unit associated to the P/N l processors that do the more computation in the game we get

[W/N l -N l-1 /N l × U (C, 2S)]× S/U (C, 2S) = [W/(U (C, 2S).N l ) -N l-1 /N l ] × S

I/O Lower Bound for Horizontal Data Movement

The following theorem extends the S -partitioning technique to the horizontal case.

Theorem 7 (Horizontal I/O Cost)

Let C = (I, V, E, O) be a CDAG. Consider any valid P-RBW game on C; for this valid game, consider the level-L storage i whose group of processors P i perform the maximum number of R6 (compute) transitions. The corresponding amount of remote get transitions is atleast

|V | U (C,2S L ).Pi -1 × S L Proof.

Inria

We let V i be the subset of nodes of C fired by P i . Let us denote S L by S for simplicity. Consider an RBW game of C with S red pebbles. Consider the partitioning of the game into C 1 , • • • , C h used in the proof of Theorem 1 for RBW. We let V i j be the set of vertices of

V i fired in C j ( V i j = V i ; V i j ∩ V ′i j = ∅ for j = j ′ ).
With the usual reasoning we can prove that |In(V

i j )| ≤ 2S and |Out(V i j )| ≤ 2S ie each V i j is a 2S-partition of C. Thus for each j, |V i j | ≤ U (C, 2S
). Now from a valid P-RBW game, we can build a valid RBW game where the restriction to V i matches the P-RBW game. By construction, each V i j is associated to at least S I/O operations in the P-RBW game. Thus the total amount of I/O for P i is at least

|V i |/|V i j | -1 × S ≥ W i /U (C, 2S) -1 × S.
Since the group P i performs maximum number of computations, W i ≥ W/P . Hence, the total amount of remote get of processors

P i is atleast ((|V | /(U (C, 2S L ).P i )) -1) × S L

Evaluation

A processor's machine balance is the ratio of the peak memory bandwidth to the peak floatingpoint performance. Lower and upper bound analysis of the algorithms can help us identify whether an algorithm is bandwidth bound at different levels of memory hierarchy.

The lower bound results can be related to the architectural machine balance as below: Consider a multi-node/multi-core system with P processors. Let N l be the total number of memory units available at level l. Consider a memory unit at level l, M i l , that incurs the maximum communication. M i l , is shared by the processor set P i l , such that P i l = P/N l . Let B i l denote the total available memory bandwidth between M i l and all its children at level l -1. Let C = (I, V, E, O) be the CDAG of the algorithm being analyzed and C l,i ⊂ C be the sub-CDAG executed by the processors P i l . The time taken for execution of C is given by

T ≥ max(T i l , T comp )
where, T i l denotes the communication time at M i l and T comp denotes the computation time for C.

For the algorithm to be not bound by memory at level l,

T i l ≤ T comp (4)
Let IO i l denote the amount of data transferred between M i l and all its children at level l -1 for the execution of a C i l . Then,

T i l = IO i l B i l ≥ LB i l B i l ( 5
)
where, LB i l denotes the lower bound on the amount of data transfer at memory unit M i l for any valid execution of C i l . The computation time of C is given by (F below indicates FLOPs)

T comp ≥ |V | P × 1 F (6)
From Equations ( 4), ( 5) and ( 6), we have,

LB i l B i l ≤ |V | P × 1 F or LB i l |V | ≤ B i l P × 1 F

RR n°8522

As P = P i l × N i l ,

LB i l × N i l |V | ≤ B i l P i l × F (7)
The term at the right-hand side of equation 7 is the machine balance value for the machine. Hence, any algorithm that fails to satisfy the condition 7, will be bandwidth bound at level l irrespective of any optimizations we do to the code.

Through similar argument, given that U B i l is the upper bound on the minimum amount of data transfer required by the algorithm at memory unit M i l , we can show that if the algorithm is communication bound, then it definitely satisfies the condition,

U B i l × N i l |V | ≥ B i l P i l × F (8)
Hence, if an algorithm fails to satisfy condition 8, we can safely conclude that there is atleast one execution order of C that is not constrained by the memory bandwidth at level l.

In particular, we were interested in understanding the memory bandwidth requirements (1) between the main memory and L2 cache within each processor, and, (2) between different processors for various algorithms. For simplicity, we assume that the L2 cache is shared by all the cores within a node, which is common in practice. Our claim is that the vertical data movement between the main memory and L2 cache will be the major bottleneck in the future machines, compared to the inter-node data movement. We show that this is true for various algorithms by comparing the lower bound for vertical data movement and the upper bound for horizontal data movement against the machine balance values of different machines.

Considering the particular case of data movement between L2 cache and the main memory, equation 7 becomes,

LB vert × N nodes |V | ≤ B vert N cores × F (9)
where, LB vert is the vertical data movement lower bound, B vert is the total bandwidth between DRAM and L2 cache, N nodes represents the number of nodes in the system, and N cores represents the number of cores within each node. Similarly, considering the inter-node communication, equation 8 becomes,

U B horiz × N nodes |V | ≥ B horiz N cores × F (10)
where, U B horiz and B horiz represent the upper bound on the horizontal data movement cost and inter-processor communication bandwidth, respectively. Specifications for some of the powerful computing systems are shown in table 1. We plan to use this list to compare various algorithms in the following sections to determine their memory requirement constraints.

Before we present the results for various numerical algorithms, we provide a brief introduction to the type of problem solved by these numerical solvers, in the following section.

Brief introduction on discretization

Many real world problems involve solving partial differential equations (PDEs). As an example, consider the heat flow on a long thin bar of unit length, of uniform material and insulated, so that the heat can enter and exit only at the boundaries (refer Fig. 2 where, α is the thermal diffusivity of the bar. (For mathematical treatment, it is sufficient to consider α = 1). Since the problem is continuous, to numerically solve the heat equation, it needs to be discretized (through finite difference approximation) to reduce it to a finite problem. In the discretized problem, the values of u(x, t) are only computed at discrete points at regular intervals of the bar, called the computational grid or mesh. The state variables at these grid points are given by u(x(i), t(m)), where x(i) = i × h, 0 ≤ i ≤ n + 1 = 1/h and t(m) = m × k; h and k are the grid spacing and timestep, respectively. Fig. 2(b) shows an example grid obtained by discretizing the one-dimensional bar.

x=0 x=1 i=0 i=4 i=8 (a) (b) boundary boundary
The governing equation, after discretization, yields the following equation at grid point i and timestamp m + 1.

-a 2 × U (i -1, m + 1) + (1 + a) × U (i, m + 1) - a 2 × U (i + 1, m + 1) = a 2 × U (i -1, m) + (1 -a) × U (i, m) + a 2 × U (i + 1, m)
where, U (p, q) = u(x(p), t(q)) and a = k/h 2 . Hence, the solution to the problem involves solving a linear system of n -1 equations at each timestamp till convergence. Each timestamp m + 1 is dependant on values of the previous timestamp m. This linear system can be represented in tridiagonal matrix form as follows:

     1+a -a 2 -a 2 1+a -a 2 -a 2 1+a -a 2 . . . . . . . . . -a 2 1+a -a 2 -a 2 1+a      ×      U (1,m+1) U (2,m+1) U (3,m+1)
. . .

U (n-1,m+1) U (n,m+1)      =      b(1,m) b(2,m) b(3,m) . . . b(n-1,m) b(n,m)      (11)
where, b(i, m) represents the right-hand side of the i-th equation at timestamp m+1. Solving this linear system of the form Ax = b for vector x provides the solution to the original problem. In general, for a d-dimensional problem, the coefficient matrix is of size n d -by-n d , while the vectors are of size n d . In practice, the elements of the matrix are not explicitly stored. Instead, their values are directly embedded in the program as constants thus eliminating the space requirement and the associated I/O cost for the matrix.

In practice, the prohibitive problem size prevents direct solution of [START_REF] Czechowski | Balance principles for algorithm-architecture co-design[END_REF]. Hence, various iterative methods were developed to efficiently solve such large linear systems. The following section derives the vertical and horizontal data movement bounds for some of these iterative linear system solvers using the results from Sections 3 and 4 and compares it against the machine balance values.

Conjugate Gradient (CG)

The Conjugate Gradient method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] is suitable for solving symmetric positive-definite linear systems. CG maintains 3 vectors at each timestep -the approximate solution x, its residual r = Axb, and a search direction p. At each step, x is improved by searching for a better solution in the direction p.

Each iteration of CG involves one sparse matrix-vector product, three vector updates, and three vector dot-products. The complete pseudocode is shown in Fig. 3.

1 f u n c t i o n CG 2
x i s t h e i n i t i a l g u e s s Recursively applying theorem 4 on the complete CDAG, C, provides us T sub-CDAGs, C 1 , C 2 , . . . , C T , corresponding to each outer loop iteration. (Vertices of vector p due to line 11 are shared between neighboring sub-CDAGs). Further, non-disjointly sub-dividing each of these sub-CDAGs, C i , into C i |x andC i |y (vertices of vector r new from line 9 are shared between C i |x andC i |y ), to decompose the effects of wavefronts W min G (υ x ) and W min G (υ y ), we obtain a lower bound of,

Q ≥ T × (2(2n d -S)) + T × (2(n d -S)) = T × (2(3n d -2S))
which tends to 6n d T as n becomes ≫ S. Finally, application of theorem 5 provides a lower bound of 6n d T /P for the parallel case.

Horizontal data movement cost

Let us assume that the input grid is block partitioned among the processors. Hence, each processor holds the input data corresponding to its local grid points and computes the data needed by those grid points. Let B = n/N 1/d nodes be the size of the block along each dimension. Computating the sparse matrix-vector product, at line 6 in Fig. 3, majorly contributes to the communication cost. This involves getting the values of the ghost cells from the neighboring processors. This value is given by (

B + 2) d -B d . If Q is the minimum I/O cost for executing CG, then, Q ≤ ((B + 2) d -B d ) × T = (B d + d 1 B d-1 2 1 + d 2 B d-2 2 2 + • • • + d d -1 B 1 2 d-1 + d d B 0 2 d B d ) × T = O(2dB d-1 T )

Analysis

Equations 9 and 10 provided us conditions to determine the vertical and horizontal memory constraints of the algorithms. We will use them to show that the running time of CG is mainly constrained by the vertical data movement. Consider a 3D-grid (d = 3), with n = 1000. The total operation count (FLOP) is 20n 3 T .

The I/O lower bound per node is given by 6n

3 T P × N cores = 6n 3 T N nodes . Hence, LB vert × N nodes |V | = 6n 3 T /N nodes × N nodes 20n 3 T = 6 20 = 0.3
This value is higher than the machine balance value of any machine (refer table 1), leaving condition 9 unsatisfied. This shows that CG will be unavoidably bandwidth bound along the vertical direction for the problems that cannot fit into the cache. The only way to improve the performance would be to increase the main memory bandwidth.

On the other hand, let us consider the horizontal data movement cost.

U B horiz × N nodes |V | = 6B 2 T × N nodes 20n 3 T = 6 n/N (1/3) nodes 2 N nodes 20n 3 = 6N (1/3) nodes 20n
This value easily falls below the machine balance values of various machines, indicating that the inter-node communication is not a bottleneck for the execution of CG algorithm.

Generalized Minimal Residual (GMRES)

The Generalized Minimum Residual (GMRES) [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] method is designed to solve non-symmetric linear systems. The most popular form of GMRES is based on the modified Gram-Schmidt procedure. The method arrives at the solution after m iterations, where m represents the dimension of a linear subspace-called Krylov subspace-formed by an orthogonal basis. At each step, x is improved by searching for a better solution along the direction of one of these basis vectors. Each outer loop iteration i of GMRES involves one sparse matrix-vector product, (i + 1) vector dot-products and i vector updates. The complete pseudocode is shown in Fig. 4. which tends to 6n d m as n grows. Finally, application of theorem 5 provides a lower bound of 6n d m/P for the parallel case.

Horizontal data movement cost

The horizontal data movement trend for GMRES is similar to CG (Sec. 5.2.2). Hence, upon applying similar analysis on the GMRES algorithm, we obtain an upper bound of

Q = O(2dB d-1 m)
where, B is the block size along each dimension.

Analysis

Consider a 3D-grid (d = 3), with n = 1000. The total number of operations for GMRES is 20n 3 m + n 3 m 2 . The vertical data movement cost per FLOP,

LB vert × N nodes |V | = 6 m + 20
For smaller values of m, this value stays higher than the machine balance value for current systems. But as m gets higher, the computational time begins to dominate the vertical data movement cost.

The Horizontal data movement cost per FLOP is given by,

U B horiz × N nodes |V | = 6N
(1/3) nodes nm This value is orders of magnitude smaller than the machine balance values of current systems, showing that the algorithm is not inter-node bandwidth bound. On the other hand, for the vertical data movement, as the lower and upper bounds do not match, it is not possible to draw a decisive conclusion unless the rate of convergence value, m is known for the problem.

Ranjan et al. [START_REF] Ranjan | Strong I/O lower bounds for binomial and FFT computation graphs[END_REF] develop better bounds than Hong & Kung for FFT using a specialized technique adapted for FFT-style computations on memory hierarchies. Ranjan et al. [START_REF] Ranjan | Upper and lower I/O bounds for pebbling rpyramids[END_REF] derive lower bounds for pebbling r-pyramids under the assumption that there is no recomputation. Recently, Ranjan et al. [START_REF] Ranjan | Vertex isoperimetric parameter of a computation graph[END_REF] develop a technique for binomial graphs.

Very recent work from U.C. Berkeley [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs That Reference Arrays -Part 1[END_REF] has developed a very novel approach to developing parametric I/O lower bounds applicable/effective for a class of nested loop computations but is either inapplicable or produces weak lower bounds for other computations (e.g., stencil computations, FFT, etc.).

The P-RBW game developed in this paper extends the parallel model for shared-memory architectures by Savage and Zubair [START_REF] Savage | A unified model for multicore architectures[END_REF] to also include the distributed-memory parallelism present in all scalable parallel architectures. The works of Irony et al. [START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF] and Ballard et al. [START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF] model communication across nodes of a distributed-memory system. Bilardi and Preperata [START_REF] Bilardi | Processor -Time Tradeoffs under Bounded-Speed Message Propagation: Part II, Lower Bounds[END_REF] develop lower bound results for communication in a distributed-memory model specialized for multidimensional mesh topologies. Our model in this paper differs from the above efforts in defining a new integrated pebble game to model both horizontal communication across nodes in a parallel machine, as well as vertical data movement through a multi-level shared cache hierarchy within a multi-core node.

Czechowski et al. [START_REF] Czechowski | Balance principles for algorithm-architecture co-design[END_REF][START_REF] Czechowski | A theoretical framework for algorithm-architecture co-design[END_REF] consider the relationship between the ratio of an algorithm's data movement cost to arithmetic work and the machine balance ratio of memory bandwidth to peak performance. Our analysis of algorithms in this paper involves a very similar theme as theirs and is inspired by their work, but we develop new lower bounds analysis to perform the analysis. Further, we compare and contrast data movement demands for horizontal across-node communication versus vertical within-node data movement and observe that the latter is often the more constraining factor.

Conclusion

Characterizing the parallel data movement complexity of a program is a cornerstone problem, that is particularly important with current and emerging power-constrained architectures where the data transfer cost will be the dominant energy and performance bottleneck. In this paper we presented an extension to the Hong and Kung red-blue pebble game model to enable development of lower bounds on data movement for parallel execution of CDAGs. The model distinguishes horizontal data movement between nodes of a distributed-memory parallel system from vertical data movement within the multi-level memory/cache hierarchy within a multi-core node. The utility of the model and the developed lower bounding techniques was demonstrated by analysis of several numerical algorithms and the garnering of interesting insights on the relative significance of horizontal versus vertical data movement for different algorithms.

1 .

 1 3.1):Inria Flexible input/output vertex labeling: Unlike the Hong & Kung model, where all

Figure 1 :

 1 Figure 1: Distributed-memory system

Figure 2 :

 2 Figure 2: One-dimensional heat flow problem

Figure 3 :

 3 Figure 3: Conjugate Gradient method

Figure 4 :

 4 Figure 4: Basic GMRES

Table 1 :

 1 (a)). Let u(x, t) represent Specifications of various computing systems

	Inria

RR n°8522

Jacobi Method

Jacobi's method involves stencil computations, that begins with an initial guess for the unknown vector x and iteratively replaces the current approximate solution at each grid point by a weighted average of its nearest neighbors on the grid. Hence, the information at one grid point can only propogate to its adjacent grid points in one iteration. Thus, it takes atleast n steps to propogate the information throughout the grid and reach to the solution.

Vertical data movement cost

In this section, we derive the I/O lower bound for Jacobi computation on a d-dimensional grid. We provide the proof for a 2D-grid below, which can be generalized to a grid of d-dimensions as shown later.

Theorem 10 (I/O lower bound for Jacobi) For the 9-points Jacobi of size n × n with T -1 time steps, the minimum

Proof. The CDAG of Jacobi computation has the property that all inputs can reach all outputs through vertex-disjoint paths. These vertex-disjoint paths will be called

), where L is the total number of vertices on the lines. From the structure of the CDAG for 2D-Jacobi computation, it can be seen that

With the similar reasoning, the I/O lower bound can be extended to higher dimensions, leading to the I/O cost of Q ≥ n d T /4.P.(2S) 1/d , for a d-dimensional grid.

It could be seen that this lower bound is tight as the tiled stencil computation algorithm has the I/O cost that matches this bound.

Horizontal data movement cost

The horizontal data movement cost is due to the communication of the ghost cells. This amounts to the I/O cost of 4BT , where B is the block size along each dimension.

Analysis

From [START_REF] Bilardi | Processor -Time Tradeoffs under Bounded-Speed Message Propagation: Part II, Lower Bounds[END_REF] and Theorem 6, and since the lower bound derived in Theorem 10 is tight, for the computation to be not bandwidth-bound along the vertical direction, the following relation has to be satisfied:

Inria

On Characterizing the Data Movement Complexity of Computational DAGs for Parallel Execution21

From Theorem 10, for a d-dimensional Jacobi, U (C, 2S l-1 ) = 4S l-1 (2S l-1 ) 1/d . Hence,

From Table 1, for IBM BG/Q, the vertical machine balance parameter for the data movement between main memory and L2 cache is 0.052. Hence, 1/4(2S 2 ) 1/d ≤ 0.052 or, d ≤ 0.21 log(2S 2 ). Substituting the value of S 2 = 4 MWords, we get, d ≤ 4.83.

By following similar reasoning and considering the machine parameters for the data movement between L2 and L1 caches, for the computation to be not bandwidth-bound, d ≤ 96.

This shows that the data movement between main memory and L2 cache is critical for the preformance and the algorithm is bandwidth bound only for higher dimensional stencils of dimension d ≥ 5, which are not common in practice.

Related Work

Hong & Kung provided the first characterization of the I/O complexity problem using the red/blue pebble game and the equivalence to 2S-partitioning of CDAGs [START_REF] Hong | I/O complexity: The red-blue pebble game[END_REF]. Their 2S-partitioning approach uses dominators of incoming edges to partitions but does not account for the internal structure of partitions. In this paper, in addition to using the 2s-partitioning technique, we also use an alternate lower bound approach that models the internal structure of CDAGs, and uses graph mincut as the basis. In addition, Hong & Kung's original model does not lend itself easily to development of effective lower bounds for a CDAG from bounds for component sub-graphs. With a change of the pebble game model to the RBW game, we were able to use CDAG decomposition to develop tight composite lower bounds for inhomogeneous CDAGs.

Several works followed Hong & Kung's work on I/O complexity in deriving lower bounds on data accesses [START_REF] Aggarwal | The input/output complexity of sorting and related problems[END_REF][START_REF] Aggarwal | A model for hierarchical memory[END_REF][START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF][START_REF] Bilardi | On the space and access complexity of computation dags[END_REF][START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF][START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF][START_REF] Savage | Models of computation -exploring the power of computing[END_REF][START_REF] Ranjan | Strong I/O lower bounds for binomial and FFT computation graphs[END_REF][START_REF] Ranjan | Upper and lower I/O bounds for pebbling rpyramids[END_REF][START_REF] Valiant | A bridging model for multi-core computing[END_REF][START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF][START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF][START_REF] Ballard | Graph expansion and communication costs of fast matrix multiplication[END_REF][START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs That Reference Arrays -Part 1[END_REF][START_REF] Solomonik | Minimizing communication in all-pairs shortest paths[END_REF][START_REF] Savage | Cache-optimal algorithms for option pricing[END_REF]. Aggarwal et al. provided several lower bounds for sorting algorithms [START_REF] Aggarwal | The input/output complexity of sorting and related problems[END_REF]. Savage [START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF][START_REF] Savage | Models of computation -exploring the power of computing[END_REF] developed the notion of S-span to derive Hong-Kung style lower bounds and that model has been used in several works [START_REF] Ranjan | Strong I/O lower bounds for binomial and FFT computation graphs[END_REF][START_REF] Ranjan | Upper and lower I/O bounds for pebbling rpyramids[END_REF][START_REF] Savage | Cache-optimal algorithms for option pricing[END_REF]. Irony et al. [START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF] provided a new proof of the Hong-Kung result on I/O complexity of matrix multiplication and developed lower bounds on communication for sequential and parallel matrix multiplication. More recently, Demmel et al. have developed lower bounds as well as optimal algorithms for several linear algebra computations including QR and LU decomposition and all-pairs shortest paths problem [START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF][START_REF] Ballard | Graph expansion and communication costs of fast matrix multiplication[END_REF][START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF][START_REF] Solomonik | Minimizing communication in all-pairs shortest paths[END_REF]. Bilardi et al. [START_REF] Bilardi | On the space and access complexity of computation dags[END_REF][START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF] develop the notion of access complexity and relate it to space complexity. Bilardi and Preparata [START_REF] Bilardi | Processor -Time Tradeoffs under Bounded-Speed Message Propagation: Part II, Lower Bounds[END_REF] developed the notion of the closed-dichotomy size of a DAG G that is used to provide a lower bound on the data access complexity in those cases where recomputation is not allowed. Our notion of schedule wavefronts is similar to the closed-dichotomy size in their work; but, unlike the work of [START_REF] Bilardi | Processor -Time Tradeoffs under Bounded-Speed Message Propagation: Part II, Lower Bounds[END_REF], we use it do develop an effective automated heuristic to compute lower bounds for CDAGs. Extending the scope of the Hong & Kung model to more complex memory hierarchies has been the subject of some research. Savage provided an extension together with results for some classes of computations that were considered by Hong & Kung, providing optimal lower bounds for I/O with memory hierarchies [START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF]. Valiant proposed a hierarchical computational model [START_REF] Valiant | A bridging model for multi-core computing[END_REF] that offers the possibility to reason in an arbitrarily complex parameterized memory hierarchy model.

Unlike Hong & Kung's original model, several models have been proposed that do not allow recomputation of values (also referred to as "no repebbling") [START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF][START_REF] Ballard | Graph expansion and communication costs of fast matrix multiplication[END_REF][START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF][START_REF] Scquizzato | Communication lower bounds for distributed-memory computations[END_REF][START_REF] Ranjan | Strong I/O lower bounds for binomial and FFT computation graphs[END_REF][START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF][START_REF] Savage | Models of computation -exploring the power of computing[END_REF][START_REF] Savage | Cache-optimal algorithms for option pricing[END_REF][START_REF] Cook | An observation on time-storage trade off[END_REF][START_REF] Irony | Communication lower bounds for distributed-memory matrix multiplication[END_REF][START_REF] Ranjan | Upper and lower I/O bounds for pebbling rpyramids[END_REF][START_REF] Ranjan | Vertex isoperimetric parameter of a computation graph[END_REF]. Savage [START_REF] Savage | Extending the Hong-Kung model to memory hierarchies[END_REF] develops results for FFT using no repebbling. Bilardi and Peserico [START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF] explore the possibility of coding a given algorithm so that it is efficiently portable across machines with different hierarchical memory systems, without the use of recomputation. Ballard et al. [START_REF] Ballard | Minimizing communication in numerical linear algebra[END_REF][START_REF] Ballard | Graph expansion and communication costs of fast matrix multiplication[END_REF] assume no recomputation is allowed in deriving lower bounds for linear algebra computations.


Victor Y Pan
email: victor.pan@lehman.cuny.edu

Elias P Tsigaridas
email: elias.tsigaridas@inria.fr

Accelerated Approximation of the Complex Roots of a Univariate Polynomial Extended Abstract

Highly efficient and even nearly optimal algorithms have been developed for the classical problem of univariate polynomial root-finding (see, e.g., [START_REF] Pan | Optimal (up to Polylog Factors) Sequential and Parallel Algorithms for Approximating Complex Polynomial Zeros[END_REF], [START_REF] Pan | Univariate Polynomials: Nearly Optimal Algorithms for Factorization and Rootfinding[END_REF], [START_REF] Mcnamee | Numerical Methods for Roots of Polynomials[END_REF], and the bibliography therein), but this is still an area of active research. By combining some powerful techniques developed in this area we devise new nearly optimal algorithms, whose substantial merit is their simplicity, important for the implementation.

We first recall the basic concept of the isolation ratio, central also for [START_REF] Pan | Optimal (up to Polylog Factors) Sequential and Parallel Algorithms for Approximating Complex Polynomial Zeros[END_REF], [START_REF] Pan | Univariate Polynomials: Nearly Optimal Algorithms for Factorization and Rootfinding[END_REF]. Assume a real or complex polynomial p = p(x) = d i=0 pix i = pn d j=1 (x -zj), p d = 0, of degree d, an annulus A(X, R, r) = {x : r ≤ |x -X| ≤ R} on the complex plane with a center X and the radii r and R of the boundary circles. Then the internal disc D(X, r) = {x : |x -X| ≤ r} is R/r-isolated and R/r is its isolation ratio if the polynomial p has no roots in the annulus. Next we reproduce [START_REF] Renegar | On the worst-case arithmetic complexity of approximating zeros of polynomials[END_REF]Corollary 4.5]. It shows that Newton's iteration converges quadratically to a single simple root of p if is initiated at the center of a 5d 2 -isolated disc that contains just this root.

Theorem 1. Suppose both discs D(c, r) and D(c, r/s) for s ≥ 5d 2 contain a single simple root α of a polynomial p = p(x) of degree d. Then Newton's iteration

x k+1 = x k -p(x k )/p ′ (x k ), k = 0, 1, . . . (1)
converges quadratically to the root α right from the start provided x0 = c.

Now suppose that we are given a disc with a single simple zero of p having an isolation ratio 1 + η for a fixed constant η > 0. Can we increase the ratio to 5d 2 ? Yes, we just need to apply a technique already used in [START_REF] Schönhage | The Fundamental Theorem of Algebra in Terms of Computational Complexity[END_REF] for the computation of the power sums of the roots lying inside such a disc. In our case this is a single root, the power sum is the root itself, and we just need its approximation c within an error at most ∆ such that rη/∆ ≥ 5d 2 . Indeed in this case ∆ ≤ 0.2rη/d 2 , and so the disc D(c, ∆) is 5d 2 -isolated.

We can shift and scale the variable x, and so wlog we assume dealing with a (1 + t) 2 -isolated disc D(0, r) for r = 1/(1 + t) for a fixed t > 0, and with polynomial p having a single simple root z1 in this disc. Recall the Laurent expansion,

p ′ (x) p(x) = d j=1 1 x -zj = - ∞ k=1 S k x k-1 + ∞ k=0 s k x -k-1 = ∞ h=-∞ c h x h . (2)
Here

|x| = 1, s0 = 1, s k = z k 1 , S k = d i=2 z -k i , k = 1, 2, . . . Consequently s k = z k
1 , whereas S k is the kth power sum of the zeros of the reverse polynomial prev(x) that lie in the disc D(0, r). The leftmost equation of ( 2) is verified by the differentiation of p(x) = pn d j=1 (x -zj). The middle equation is implied by the decompositions 1

x-z 1 = 1 x ∞ h=0 z 1 x h and 1 x-z i = -1 z i ∞ h=0 x z i h for i > 1, provided |x| = 1 for all i.
We cover the case of any positive integer k, although we only need the case where k = 1. For a fixed positive integer q we compute the approximations s * k ≈ s k as follows,

s * k = 1 q q-1 j=0 ω j(k+1) p(ω j )/p ′ (ω j ), k = 1, 2, . . . , q -1. ( 3
)
Here ω = ωq = exp(2π √ -1/q) is a primitive qth root of unity. Then the evaluation of the polynomial p(x) at the qth roots of unity amounts to the same task for a polynomial pq(x) of degree at most q -1 with the coefficients pq,i = l j=0 pi+jq for l = ⌊d/q⌋ obtained by means of less than d additions of the coefficients of p.

Having computed the polynomial pq(x) we reduce the evaluation of all the desired approximations s * k for k = 1, . . . , q -1 essentially to performing three DFTs, each on q points, that is to a total of O(q log(q)) ops. Namely, we apply two DFTs to compute p(ω i ) and p ′ (ω i ) for i = 0, 1, . . . , q -1 and a single DFT to multiply the DFT matrix Ω = [ω hi ] q-1 h,i=0 by the vector v = [p(ω i )/p ′ (ω i ] q-1 i=0 .

Let us estimate the approximation errors. Equations ( 2) and (3) imply that

s * k = +∞ l=-∞ c -k-1+lq . Moreover, (2) for h = -k -1, k ≥ 1 implies that s k = c -k-1 , whereas (2) for h = k -1, k ≥ 1 implies that S k = -c k-1 . Consequently s * k -s k = ∞ l=1 (c lq-k-1 + c -lq-k-1 ).
We assumed in (3) that 0 < k < q -1. It follows that c -lq-k-1 = s lq+k and c lq-k-1 = -S lq-k for l = 1, 2, . . ., and we obtain

s * k -s k = ∞ l=1 (s lq+k -S lq-k ). ( 4
)
On the other hand

|s h | ≤ z h , |S h | ≤ (d -1)z h , h = 1, 2, . . . where z = max 1≤j≤d min(|zj |, 1/|zj |)
, and so z ≤ 1 1+t in our case. Substitute these bounds into (4) and obtain |s * k -s k | ≤ (z q+k + (d -1)z q-k )/(1 -z q ). Therefore it is sufficient to choose q of order log(d) to decrease the error of the approximation to the root z1 by a factor of gd h for any pair of constants g and h, and so we can ensure the desired error bound ∆. To support this computation we only need less than d additions, followed by O(log(d)) evaluations of the polynomial pq(x) of degree q -1 at the lth roots of unity for l = O(log(d)). This involves O(log(d) log(log(d))) ops overall. (Here and hereafter "ops" stand for "arithmetic operations".) Summarizing we obtain the following estimates. Corollary 4. Suppose that we are given d discs, each containing a single simple root of a polynomial p = p(x) of degree d and each being (1 + η) 2 -isolated for a fixed η > 0. Then we can approximate all d roots of this polynomial within a fixed positive error bound ǫ < 1 by using O(d log 2 (d)(1 + log(log(1/ǫ)))) ops.

Proof: Apply the same algorithm that supports Corollary 3 concurrently in all d given discs, but instead of the qth roots of unity use q equally spaced points at the boundary circle of each input disc (that is dq = O(d log d) points overall) and instead of FFT apply the Moenck-Borodin algorithm for multipoint polynomial evaluation. Also use it at the stage of performing concurrent Newton's iteration initialized at the centers of the 5d 2 -isolated subdiscs of the d input discs, each subdisc computed by the algorithm that supports Theorem 2. Here we work with the dth degree polynomial p rather than with the qth degree polynomials pq because to support transition to polynomials pq of the degree q for d discs we would need to perform d shifts and scalings of the variable x. Instead we employ the Moenck-Borodin algorithm, which still enables us to obtain a nearly optimal root-refiner. Technically, in a relatively minor change of our algorithm, we replace the matrix Ω = [ω j(k+1) ] j,k in (3) by the matrix [c + ω j(k+1) ] j,k = c[1] j,k + Ω where c is invariant in j and k. The multiplication of the new matrix by a vector v is still reduced to multiplication of the matrix Ω by a vector v with the additional 3d ops for computing the vector c[1] j,k v and adding it to the vector Ωv.

The Moenck-Borodin algorithm uses nearly linear arithmetic time, and [START_REF] Kirrinnis | Polynomial Factorization and Partial Fraction Decomposition by Simultaneous Newton's Iteration[END_REF] proved that this algorithm supports multipoint polynomial evaluation at a low Boolean cost as well (see also [START_REF] Van Der Hoeven | Fast composition of numeric power series[END_REF], [START_REF] Pan | Proc. International Symp. on Symbolic and Algebraic Computations (ISSAC 2013)[END_REF], [START_REF] Kobel | Fast Approximate Polynomial Multipoint Evaluation and Applications[END_REF], [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF], [START_REF] Pan | Transformations of Matrix Structures Work Again, accepted by Linear Algebra and Its Applications[END_REF], [START_REF] Pan | Fast Approximation Algorithms for Computations with Cauchy Matrices and Extensions[END_REF]). Consequently our algorithm supporting Corollary 4 can be extended to support a nearly optimal Boolean cost bound for refining all simple isolated roots of a polynomial.

We can immediately relax the assumption that the roots are simple because our proof of Theorem 2 applies to a multiple root as well. Furthermore deduce from the Lucas theorem that the isolation ratio of the basic discs in our algorithms does not decrease when we shift from a polynomial to its derivative and higher order derivatives. Therefore we can just apply Newton's iteration to the derivative or to a higher order derivative to approximate a double or multiple root, respectively.

Theorem 2 .

 2 Suppose the unit disc D(0, r) = {x : |x| ≤ 1} is (1 + η) 2 -isolated for (1 + η)r = 1 and a fixed η > 0 and contains a single simple root z of a polynomial p = p(x) of a degree d. Then it is sufficient to apply less than d additions and O(log(d) log(log(d))) other ops to compute a 5d 2 -isolated subdisc of D(0, r) containing this root. Combine Thm. 1 and 2 and obtain the following result. Corollary 3. Under the assumptions of Theorem 2 we can approximate the root z of the polynomial p(x) within a fixed positive error bound ǫ < 1 by using O(log(d) log(log(d)) + d log(log(1/ǫ))) ops.

Acknowledgments. VP is supported by NSF Grant CCF-1116736. ET is partially supported by GeoLMI (ANR 2011 BS03 011 06), HPAC (ANR ANR-11-BS02-013) and an FP7 Marie Curie Career Integration Grant.

Boolean cost bounds

Hereafter OB denotes the bit or Boolean complexity ignoring logarithmic factors. To estimate it we apply some results from [START_REF] Pan | Proc. International Symp. on Symbolic and Algebraic Computations (ISSAC 2013)[END_REF]- [START_REF] Pan | Nearly Optimal Computations with Structured Matrices[END_REF], which hold in the general case where the coefficients of the polynomials are known up to an arbitrary precision. In our case the input polynomial is known exactly; the parameter λ to be specified in the sequel could be considered as the working precision.

Let p be given as a λ-approximation, ie lg p -p ∞ ≤ -λ. We compute pq by using d additions. This produces a polynomial such that lg pq ∞ ≤ τ +lg d, and lg pq -

Similar bounds hold for p ′ q , ie lg p ′ q ∞ ≤ τ + 2 lg d, and [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF]Lemma 16], and similar bounds hold for pq(ω i ).

The divisions pq(ω i )/p ′ q (ω i ) output complex numbers such that |pq(ω)/p ′ q (ω)| ≤ τ +2 lg d+lg lg d+2 with the logarithm of the error ≤ -λ+τ lg(4d)+3/2 lg 2 d+9/2 lg d+2 lg lg d+11.

The final DFT produces numbers such that the logarithms of their magnitudes are not greater than τ +2 lg d+2 lg lg d+ 4 and the logarithms of their approximation errors are at most -λ + τ lg(8d) + 3/2 lg 2 d + 13/2 lg d + 4 lg lg d + 18, [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF]Lemma 16].

To achieve an error within 2 -ℓ in the final result, we perform all the computations with accuracy λ = ℓ + τ lg(8d

We perform d additions at the cost OB(dλ) and perform the rest of computations, that is the 3 DFTs, at the cost OB(lg d lg lg d µ(λ)) or OB(d(ℓ + τ )) [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF]Lemma 16].

If the root that we want to refine is not in the unit disc, then we replace τ in our bounds with dτ .

We apply a similar analysis from [10, Section 2.3] to the Newton iteration (see also [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF]Section 2.3]) and arrive at the same asymptotic bounds on the Boolean complexity. Only the overhead constants change because now we perform computations with complex numbers.

The overall complexity is OB(d 2 τ + dℓ) and the working precision is O(dτ + ℓ).

Here we assume the exact input, that is assume the coefficients of the input polynomials known up to arbitrary precision. For the refinement of the root up to precision of L bits, we arrive at an algorithm with the complexity in OB(d 2 τ + dL).

If we are interested in refining all complex roots, we cannot work anymore with the polynomial pq of degree q = O(lg d) unless we add the cost of d shifts of the initial approximations to the origin. Instead we rely on fast algorithms for multipoint evaluation. Initially we evaluate the polynomial p of degree d at O(d lg d) points, and we assume that lg p ∞ ≤ τ . These d points approximate the roots of p, and so their magnitude is at most ≤ 2 τ .

We use the following result of [START_REF] Pan | Nearly Optimal Computations with Structured Matrices[END_REF]Lemma 21]. Similar bounds appear in [START_REF] Kirrinnis | Polynomial Factorization and Partial Fraction Decomposition by Simultaneous Newton's Iteration[END_REF][START_REF] Kobel | Fast Approximate Polynomial Multipoint Evaluation and Applications[END_REF][START_REF] Van Der Hoeven | Fast composition of numeric power series[END_REF].

Lemma 5 (Modular representation). Assume m+1 polynomials, F ∈ ( C[x] of degree 2mn and Pj ∈ ( C[x] of degree n, for j = 1, . . . , m such that F ∞ ≤ 2 τ 1 and all roots of the polynomials Pj for all j have magnitude of at most 2 ρ . Furthermore assume λ-approximations of F by F and of Pj by Pj such that F -

Then we can compute an ℓ-approximations Fj of Fj = F mod Pj for j = 1, . . . , m such that Fj -Fj ∞ ≤ 2 -ℓ in OB(m n (ℓ + τ1 + m n ρ)).

Using this lemma we bound the overall complexity of multipoint evaluation by OB(d(L+dτ )). The same bounds holds at the stage where we perform Newton's iteration. We need to apply Newton's operator O(1) for each root. Each application of the operators consists of two polynomial evaluations. We perform the evaluations simultaneously and apply Lemma 5 to bound the complexity. On similar estimates for the refinement of the real roots see [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF].

Extensions

The algorithm of [START_REF] Mehlhorn | From Approximate Factorization to Root Isolation with Application to Cylindrical Algebraic Decomposition[END_REF] computes at nearly optimal cost 64disolated initial discs for all d roots of a polynomial p(x). By combining this algorithm with ours we obtain a distinct alternative algorithm, which like the one of [START_REF] Pan | Optimal (up to Polylog Factors) Sequential and Parallel Algorithms for Approximating Complex Polynomial Zeros[END_REF], [START_REF] Pan | Univariate Polynomials: Nearly Optimal Algorithms for Factorization and Rootfinding[END_REF], supports the record nearly optimal bounds on the Boolean complexity of the approximation of all complex polynomial roots, but has the advantage of allowing substantially simpler implementation.

Finally the same algorithm of [START_REF] Schönhage | The Fundamental Theorem of Algebra in Terms of Computational Complexity[END_REF] approximate the power sums of any number m of roots (forming, e.g., a single cluster or a number of clusters) in an isolated disc. The algorithm runs at about the same cost, already stated and depending just on the isolation ratio. Having the power sums available we can readily compute the coefficients of the factor f of p of degree d f , whose roots are exactly the roots of p in this disc: this is a numerically stable algorithm using O(m log(m)) ops (cf. [1, pages 34-35]).