Enver Kayaaslan

Cevdet Aykanat

AND BORA U C ¸ar

.5D PARALLEL SPARSE MATRIX-VECTOR MULTIPLY *

Keywords: sparse matrix partitioning, parallel sparse matrix-vector multiplication, directed hypergraph model, bipartite vertex cover, combinatorial scientific computing AMS subject classifications. 05C50, 05C65, 05C70, 65F10, 65F50, 65Y05

There are three common parallel sparse matrix-vector multiply algorithms: 1D row-parallel, 1D column-parallel and 2D row-column-parallel. The 1D parallel algorithms offer the advantage of having only one communication phase. On the other hand, the 2D parallel algorithm is more scalable but it suffers from two communication phases. Here, we introduce a novel concept of heterogeneous messages where a heterogeneous message may contain both input-vector entries and partially computed output-vector entries. This concept not only leads to a decreased number of messages, but also enables fusing the input-and output-communication phases into a single phase. These findings are exploited to propose a 1.5D parallel sparse matrix-vector multiply algorithm which is called local row-column-parallel. This proposed algorithm requires a constrained fine-grain partitioning in which each fine-grain task is assigned to the processor that contains either its inputvector entry, or its output-vector entry, or both. We propose two methods to carry out the constrained fine-grain partitioning. We conduct our experiments on a large set of test matrices to evaluate the partitioning qualities and partitioning times of these proposed 1.5D methods.

1. Introduction. Sparse matrix-vector multiply (SpMV) of the form y ← Ax is a fundamental operation in many iterative solvers for linear systems, eigensystems and least squares problems. This renders the parallelization of SpMV operation an important problem. In the literature, there are three SpMV algorithms: row-parallel, column-parallel, and row-column-parallel. Row-parallel and column-parallel (called 1D) algorithms have a single communication phase, in which either the x-vector or partial results on the y-vector entries are communicated. Row-column-parallel (2D) algorithms have two communication phases; first the x-vector entries are communicated, then the partial results on the y-vector entries are communicated. We propose another parallel SpMV algorithm in which both the x-vector and the partial results on the y-vector entries are communicated as in the 2D algorithms, yet the communication is handled in a single phase as in the 1D algorithms. That is why, the new parallel SpMV algorithm is dubbed 1.5D.

Partitioning methods based on graphs and hypergraphs are widely established to achieve 1D and 2D parallel algorithms. For 1D parallel SpMV, row-wise or columnwise partitioning methods are available. The scalability of 1D parallelism is limited especially when a row or a column has too many nonzeros in the row-and columnparallel algorithms, respectively. In such cases, the communication volume is high and the load balance is hard to achieve, severely reducing the solution space. The associated partitioning methods are usually the fastest alternatives. For 2D parallel SpMV, there are different partitioning methods. Among them, those that partition matrix entries individually, based on the fine-grain model [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF], have the highest flexibility. That is why they usually obtain the lowest communication volume and achieve near perfect balance among nonzeros per processor [START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF]. However, the fine-grain partitioning approach usually results in higher number of messages; not surprisingly higher number of messages hampers the parallel SpMV performance [START_REF] Kaya | Analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF].

The parallel SpMV operation is composed of fine-grain tasks of multiply-and-add operations of the form y i ← y i + a ij x j . Here, each fine-grain task is identified with a unique nonzero and assumed to be performed by the processor that holds the associated nonzero by the owner-computes rule. The proposed 1.5D parallel SpMV imposes a special condition on the operands of the fine-grain task y i ← y i + a ij x j : the processor that holds a ij should also hold x j or should be responsible for y i (or both). The standard rowwise and columnwise partitioning algorithms for 1D parallel algorithms satisfy the condition, but they are too restrictive. The standard fine-grain partitioning approach does not necessarily satisfy the condition. Here we propose two methods for partitioning for 1.5D parallel SpMV. With the proposed partitioning methods, the overall 1.5D parallel SpMV algorithm inherits the important characteristics of 1D and 2D parallel SpMV and the associated partitioning methods. In particular, it has

• a single communication phase as in 1D parallel SpMV,

• the partitioning flexibility close to that of 2D fine-grain partitioning,

• much reduced number of messages compared to the 2D fine-grain partitioning,

• a partitioning time close to that of 1D partitioning.

We propose two methods (Section 4) to obtain a 1.5D local fine-grain partition each with a different setting and approach where some preliminary studies on these methods are given in our recent work [START_REF] Kayaaslan | Semi-two-dimensional partitioning for parallel sparse matrix-vector multiplication[END_REF]. The first method is developed by proposing a directed hypergraph model. Since current partitioning tools cannot meet 1.5D partitioning requirements, we adopt and adapt an approach similar to that of a recent work by Pelt and Bisseling. [START_REF] Pelt | A medium-grain method for fast 2d bipartitioning of sparse matrices[END_REF]. The second method has two parts. The first part applies a conventional 1D partitioning method but decodes this only as a partition of the vectors x and y. The second part decides nonzero/task distribution under the fixed partition of the input and output vectors.

The remainder of this paper is as follows. In Section 2, we give a background on parallel SpMV. Section 3 presents the proposed 1.5D local row-column-parallel algorithm and 1.5D local fine-grain partitioning. The two methods proposed to obtain a local fine-grain partition are presented and discussed in Section 4. Section 5 gives a brief review of recent related work. We display our experimental results in Section 6 and conclude the paper in Section 7.

2. Background on parallel sparse matrix-vector multiply.

2.1. The anatomy of parallel sparse matrix-vector multiply. Recall that y ← Ax can be cast as a collection of fine-grain tasks of multiply-and-add operations [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF] y i ← y i + a ij × x j .

These tasks can share input and output-vector entries. When a task a ij and the input-vector entry x j are assigned to different processors, say P and P r , respectively, P r sends x j to P , which is responsible to carry out the task a ij . An input-vector entry x j is not communicated multiple times between processor pairs. When a task a ij and the output-vector entry y i are assigned to different processors, say P r and P k , respectively, then P r performs ŷi ← ŷi + a ij × x j as well as all other multiply-and-add operations that contribute to the partial result ŷi and then sends ŷi to P k . The partial results received by P k from different processors are then summed to compute y i .

Task-and-data distributions.

Let A be an m×n sparse matrix and a ij represent both a nonzero of A and the associated fine-grain task of multiply-andadd operation [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF]. Let x and y be the input-and output-vectors of size n and This manuscript is for review purposes only. Fig. 1: A task-and-data distribution Π(y ← Ax) of matrix-vector multiply with a 2×3 sparse matrix A. m, respectively, and K be the number of processors. We define a K-way task-anddata distribution Π(y ← Ax) of the associated SpMV as a 3-tuple Π(y ← Ax) = (Π(A), Π(x), Π(y)), where Π(A) = {A (1) , . . . , A (K) }, Π(x) = {x (1) , . . . , x (K) }, and Π(y) = {y (1) , . . . , y (K) }. We can also represent Π(A) as a nonzero-disjoint summation (2)

A = A (1) + A (2) + • • • + A (K) .
In Π(x) and Π(y), each x (k) and y (k) is a disjoint subvector of x and y, respectively.

Figure 1 illustrates a sample 3-way task-and-data distribution of matrix-vector multiply on a 2×3 sparse matrix.

For given input-and output-vector distributions Π(x) and Π(y), the columns and rows of A and those of A (k) can be permuted, conformably with Π(x) and Π(y), to form K×K block structures:

(3) A =      A11 A12 • • • A1K A21 A22 • • • A2K . . . . . . . . . . . . AK1 AK2 • • • AKK      , (4) A (k) =       A (k) 11 A (k) 12 • • • A (k) 1K A (k) 21 A (k) 22 • • • A (k) 2K . . . . . . . . . . . . A (k) K1 A (k) K2 • • • A (k) KK       .
Note that the row and column orderings (4) of the individual A (k) matrices are in compliance with the row and column orderings (3) of A. Hence, each block A k of the block structure (3) of A can be written as a nonzero-disjoint summation (5)

A k = A (1) k + A (2) k + • • • + A (K) k .
Let Π(y ← Ax) be any K-way task-and-data distribution. According to this distribution, each processor P k holds the submatrix A (k) , holds the input-subvector

x (k) and is responsible for storing/computing the output subvector y (k) . The finegrain tasks (1) associated with the nonzeros of A (k) are to be carried out by P k .

An input-vector entry x j ∈ x (k) is sent from P k to P , which is called an input communication, if there is a task a ij ∈ A ( ) associated with a nonzero at column j.

On the other hand, P k receives a partial result ŷi on an output-vector entry y i ∈ y (k) from P , which is referred to as an output communication, if there is a task a ij ∈ A ( ) associated with a nonzero at row i. Therefore, the fine-grain tasks associated with the nonzeros of the column stripe A * k = [A T 1k , . . . , A T Kk ] T are the only ones that require an input-vector entry of x (k) and the fine-grain tasks associated with the nonzeros of the row stripe A k * = [A k1 , . . . , A kK ] are the only ones that contribute to the computation of an output-vector entry of y (k) .

1D parallel sparse matrix-vector multiply.

There are two main alternatives for 1D parallel SpMV, row-parallel and column-parallel.

This manuscript is for review purposes only.

E. KAYAASLAN, B. UC ¸AR AND C. AYKANAT

In the row-parallel SpMV, the basic computational units are the rows. For an output-vector entry y i assigned to processor P k , the fine-grain tasks associated with the nonzeros of A i * = {a ij ∈ A : 1 ≤ j ≤ n} are combined into a composite task of inner product y i ← A i * x which is to be carried out by P k . Therefore, for the rowparallel algorithm, a task-and-data distribution Π(y ← Ax) of matrix-vector multiply on A should satisfy the following condition:

(6) a ij ∈ A (k) whenever y i ∈ y (k) .
Then, Π(A) coincides with the output-vector distribution Π(y)-each submatrix is a row stripe of the block structure (3) of A. In the row-parallel parallel SpMV, all messages are communicated in an input-communication phase called expand where each message contains only input-vector entries.

In the column-parallel SpMV, the basic computational units are the columns. For an input-vector entry x j assigned to processor P k , the fine-grain tasks associated with the nonzeros of A * j = {a ij ∈ A : 1 ≤ i ≤ m} are combined into a composite task of "daxpy" operation ŷk ← ŷk + A * j x j which is to be carried out by P k where ŷk is the partially computed output-vector of P k . As a result, a task-and-data distribution Π(y ← Ax) of matrix-vector multiply on A for the column-parallel algorithm should satisfy the following condition:

(7) a ij ∈ A (k) whenever x j ∈ x (k) .
Here, Π(A) coincides with the input-vector distribution Π(x)-each submatrix

A (k)
is a column stripe of the block structure (3) of A. In the column-parallel SpMV, all messages are communicated in an output-communication phase called fold where each message contains only partially computed output-vector entries.

The column-net and row-net hypergraph models [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF] can be respectively used to obtain the required task-and-data partitioning for the row-parallel and column-parallel SpMV.

2.4. 2D parallel sparse matrix-vector multiply. In the 2D parallel SpMV, also referred to as the row-column-parallel, the basic computational units are nonzeros [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF][START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF]. The row-column-parallel algorithm requires fine-grain partitioning which imposes no restriction on distributing tasks and data. The row-column-parallel algorithm contains two communication and two computational phases in an interleaved manner as shown in Algorithm 1. The algorithm starts with the expand phase where the required input-subvector entries are communicated. The second step computes only those partial results that are to be communicated in the following fold phase.

In the final step, each processor computes its own output-subvector. If we have a rowwise partitioning, the steps 2, 3 and 4c are not needed and hence the algorithm reduces to the row-parallel algorithm. Similarly, the algorithm without steps 1, 2b and 4b, can be used when we have a columnwise partitioning. The row-column-net hypergraph model [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF][START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF] can be used to obtain the required task-and-data partitioning for row-column-parallel SpMV. This manuscript is for review purposes only.

Algorithm 1 The row-column-parallel sparse matrix-vector multiply

For each processor P k :

1. (expand) for each nonzero column stripe A

k ← A (k) k x (k) and (b) y ( ) k ← y ( ) k + r =k A (k) r x(
y (k) ← A (k) kk x (k) , (b) y (k) ← y (k) + A (k) k x( ) k and (c) y (k) ← y (k) + =k ŷ(k) .
In the well-known parallel SpMV, the messages are homogenous in the sense that they pertain to either x-or y-vector entries. In the proposed row-column-parallel SpMV algorithm, the number of messages are reduced with respect to the row-columnparallel algorithm by making the messages heterogenous (pertaining to both x-and y-vector entries), and by communicating them in a single expand-fold phase. If a processor P sends a message to processor P k in both of the expand and fold phases, then the number of messages required from P to P k reduces from two to one. However, if a message from P to P k is sent only in the expand phase or only in the fold phase, then there is no reduction in the number of such messages.

A Task categorization.

We introduce a two-way categorization of inputand output-vector entries and a four-way categorization of fine-grain tasks (1) according to a task-and-data distribution Π(y ← Ax) of matrix-vector multiply on A. For a task a ij , the input-vector entry x j is said to be local if both a ij and x j are assigned to the same processor; the output-vector entry y i is said to be local if both a ij and y i are assigned to the same processor. With this definition, the tasks can be classified into four groups. The task k) and y i ∈ y (k) .

y i ← y i + a ij × x j on P k is          input-output-local if x j ∈ x (k) and y i ∈ y (k) , input-local if x j ∈ x (k) and y i ∈ y (k) , output-local if x j ∈ x (k) and y i ∈ y (k) , nonlocal if x j ∈ x (
Recall that an input-vector entry x j ∈ x ( ) is sent from P to P k if there exists a task a ij ∈ A (k) at column j, which implies that the task a ij of P k is either output-local or nonlocal since x j ∈ x (k) . Similarly, for an output-vector entry y i ∈ y ( ) , P receives

This manuscript is for review purposes only.

a partial result ŷi from P k if a task a ij ∈ A (k) , which implies that the task a ij of P k is either input-local or nonlocal since y i ∈ y (k) . We can also infer from that the input-output-local tasks neither depend on the input-communication phase nor incur a dependency on the output-communication phase. However, the nonlocal tasks are linked with both communication phases.

In the row-parallel algorithm, each of the fine-grain tasks is either input-outputlocal or output-local due to the rowwise partitioning condition [START_REF] Aykanat | Patoh (partitioning tool for hypergraphs)[END_REF]. For this reason, no partial result is computed for other processors, and thus no output communication is incurred. In the column-parallel algorithm, each of the fine-grain tasks is either input-output-local or input-local due to the columnwise partitioning condition [START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF]. In the row-column-parallel algorithm, the input and output communications have to be carried out in separate phases. The reason is that the partial results on the outputvector entries to be sent are partially derived by performing nonlocal tasks that rely on the input-vector entries received.

Local fine-grain partitioning.

In order to remove the dependency between the two communication phases in the row-column-parallel algorithm, we propose the local fine-grain partitioning where "locality" refers to the fact that each finegrain task is input-local, output-local or input-output-local. In other words, there is no nonlocal fine-grain task.

A task-and-data distribution Π(y ← Ax) of matrix-vector multiply on A is said to be a local fine-grain partition if the following condition is satisfied:

(8) a ij ∈ A (k) + A ( ) whenever y i ∈ y (k) and x j ∈ x ( ) .
Notice that this condition is equivalent to

(9) if a ij ∈ A (k) then either x j ∈ y (k) , or y i ∈ x (k) , or both.
Due to (4) and ( 9), each submatrix A (k) becomes of the following form ( 10)

A (k) =          0 . . . A (k) 1k . . . 0 . . . . . . . . . . . . . . . A (k) k1 • • • A (k) kk • • • A (k) kK . . . . . . . . . . . . . . . 0 . . . A (k) Kk . . . 0          .
In this form, the tasks associated with the nonzeros of diagonal block A

A k = A (k) k + A ( ) k ,
and for each diagonal block we have

A kk = A (k) kk .
In order to clarify Equations ( 8)-( 11), we provide the following 4-way local finegrain partition on A as permuted into a 4×4 block structure.

This manuscript is for review purposes only. Fig. 2: A sample local fine-grain partition. Here, a 12 is an output-local task, a 13 is an input-output-local task, a 21 is an output-local task, and a 22 is an input-local task.

P 1 P 2 P 3 P 1 P 2 P 3 x 1 x 1 x 2 x 2 x 3 x 3 y 1 y 1 y 2 y 2 a
A =      A 11 A (1) 12 A (1) 13 A (1) 14 A (1) 21 0 0 0 A (1) 31 0 0 0 A (1) 41 0 0 0      +      0 A (2) 12 0 0 A (2) 21 A 22 A (2) 23 A (2) 24 0 A (2) 32 0 0 0 A (2) 42 0 0      +      0 0 A (3) 13 0 0 0 A (3) 23 0 A (3) 31 A (3) 32 A 33 A (3) 34 0 0 A (3) 43 0      +      0 0 0 A (4) 14
0 0 0 A (4) 24 0 0 0 A (4) 34 A (4) 41 A (4) 42 A (4) 43 A 44      . For instance, A 42 = A (2) 42 + A (4) 42 , A 23 = A (2) 23 + A (3) 23 , A 31 = A (1) 31 + A (3)
31 , . . . , etc.

Figure 2 displays a sample 3-way local fine-grain partition on the same sparse matrix used in Figure 1. In this figure, a 13 ∈ A (1) where y 1 ∈ y (2) and x 3 ∈ x (1) and thus a 13 is an input-local task of P 1 . Also, a 21 ∈ A (3) where y 2 ∈ y (3) and x 1 ∈ x (1) and thus a 21 is an output-local task of P 3 . The local row-column-parallel algorithm is composed of three steps as shown in Algorithm 2. In the first step, processors concurrently perform their input-local tasks which contribute to partially computed output-vector entries for other processors. In the expand-fold phase, for each nonzero off-diagonal block

Local row-column

A k = A (k) k + A ( ) k , P k prepares a message [x (k) , ŷ( ) k ] for P .
Here, x(k) contains the input-vector entries of x (k) that are required by the output-local tasks of P , whereas ŷ( ) k contains the partial results on the output-vector entries of y ( ) , where the partial results are derived by performing the input-local tasks of P k . In the last step, each processor P k computes output-subvector y (k) by summing the partial results computed locally by its own input-output-local tasks (step 3a) and output-local tasks (step 3b) as well as the partial results received from other processors due to their input-local tasks (step 3c).

For a message [x (k) , ŷ( ) k ] from processor P k to P , the input-vector entries of x(k) correspond to the nonzero columns of A ( ) k , whereas the partially computed output-This manuscript is for review purposes only.

Algorithm 2 The local row-column-parallel sparse matrix-vector multiply

For each processor P k :

1. for each nonzero block

A (k) k , where = k; compute y ( ) k ← A (k) k x (k) , input-local tasks of P k 2. (expand-fold) for each nonzero block A k = A (k) k + A ( )
k , where = k; (a) form vector x(k) , which contains only those entries of x (k) corresponding to nonzero columns in A

k corresponding to nonzero rows in A (k) k , (c) send vector [x (k) , ŷ( ) k ] to processor P . 3. compute output-subvector (a) y (k) ← A (k) kk x (k) , input-output-local tasks of P k (b) y (k) ← y (k) + A (k) k x( ) k and output-local tasks of P k (c) y (k) ← y (k) + =k ŷ(k) .
input-local tasks of other processors

a ij P `Pk x j y i x j ŷi ŷi ŷi + a ij ⇥ x j y i y i + ŷi P `PPk P k a ij a ij x j x j y i y i ŷi x j ŷi ŷi + a ij ⇥ x j y i y i + ŷi y i y i + a ij ⇥ x j ŷ2 a 22 ⇥ x 2 y 1 a 12 ⇥ x 2 + a 13 ⇥ x 3 y 2 a 21 ⇥ x 1 + ŷ2 P 1 P 2 P 3 x 1 x 2 x 3 y 1 y 2 a 12 a 13 a 22 a 21 [x 2 ] [x 1 , ŷ2 ] P r 3 2 2 1 1 P 1 P 1 P 2 P 2 P 3 P 1 P 2 P 2 P 3 P 3 P 2 P 1
Fig. 3: An illustration of Algorithm 2 for the local fine-grain partition in Figure 2.

vector entries of ŷ( ) k correspond to the nonzero rows of

A (k) k . That is, x(k) = [x j : a ij ∈ A ( ) k ] and ŷ( ) k = [ŷ i : a ij ∈ A (k) k ]. This message is heterogeneous if A (k) k and A ( )
k are both nonzero and homogeneous otherwise. We also note that the number of messages is equal to the number of nonzero off-diagonal blocks of the block structure

(3) of A induced by the vector distribution (Π(x), Π(y)). Figure 3 illustrates the steps of Algorithm 2 on the sample local fine-grain partition given in Figure 2. As seen in the figure, there are only two messages to be communicated. One message is homogeneous, which is from P 1 to P 2 and contains only an input-vector entry x 2 , whereas the other message is heterogeneous, which is from P 1 to P 3 and contains an input-vector entry x 1 and a partially computed output-vector entry ŷ2 .

Two proposed methods for local row-column-parallel partitioning.

We propose two methods to find a local row-column-parallel partition that is required for 1.5D local row-column-parallel SpMV. One method finds vector and nonzero distributions simultaneously, whereas the other one has two parts in which vector and nonzero distributions are found separately.

4.1.

A directed hypergraph model for simultaneous vector and nonzero distribution. In this method, we adopt the elementary hypergraph model for the fine-grain partitioning [START_REF] Uc | Revisiting hypergraph models for sparse matrix partitioning[END_REF] and introduce an additional locality constraint on partitioning in order to obtain a local fine-grain partition. In this hypergraph model

This manuscript is for review purposes only.

H 2D = (V, N ), there is an input-data vertex for each input-vector entry, an outputdata vertex for each output-vector entry and a task vertex for each fine-grain task (or per matrix nonzero) for a given matrix A. That is,

V = {v x (j) : x j ∈ x} ∪ {v y (i) : y i ∈ y} ∪ {v z (ij) : a ij ∈ A} .
The input-and output-data vertices have zero weights, whereas the task vertices have unit weights. In H 2D , there is an input-data net for each input-vector entry, and an output-data net for each output-vector entry. An input-data net n x (j), corresponding to the input-vector entry x j , connects all task vertices associated with the nonzeros at column j as well as the input-data vertex v x (j). Similarly, an output-data net n y (i), corresponding to the output-vector entry y i , connects all task vertices associated with the nonzeros at row i as well as the output-data vertex v y (i). That is

N = {n x (j) : x j ∈ x} ∪ {n y (i) : y i ∈ y} , n x (j) = {v x (j)} ∪ {v z (ij) : a ij ∈ A, 1 ≤ i ≤ m} , and
n y (i) = {v y (i)} ∪ {v z (ij) : a ij ∈ A, 1 ≤ j ≤ n}.
Note that each input-data net connects a separate input-data vertex, whereas each output-data net connects a separate output-data vertex. We associate nets with their respective data vertices.

We enhance the elementary row-column-net hypergraph model by imposing directions on the nets; this is required for modeling the dependencies and their nature.

Each input-data net n x (j) is directed from the input-data vertex v x (j) to the task vertices connected by n x (j), and each output-data net n y (i) is directed from the task vertices connected by n y (i) to the output-data vertex v y (i). Each task vertex v z (ij) is connected by a single input-data-net n x (j) and a single output-data-net n y (i).

In order to impose the locality in the partitioning, we introduce the following constraint for vertex partitioning on the directed hypergraph model H 2D : each task vertex v z (ij) should be assigned to the part that contains either input-data vertex v x (j), or output-data vertex v y (i), or both. Figure 4a displays a sample 6×7 sparse matrix. Figure 4b illustrates the associated directed hypergraph model. Figure 4c shows a 3-way vertex partition of this directed hypergraph model satisfying the locality constraint, and Fig. 4d shows the local fine-grain partition decoded by this partition.

Instead of developing a partitioner for this particular directed hypergraph model, we propose a task-vertex amalgamation procedure which will help in meeting the described locality constraint by using a standard hypergraph partitioning tool. For this, we adopt and adapt a simple-yet-effective approach of Pelt and Bisseling [START_REF] Pelt | A medium-grain method for fast 2d bipartitioning of sparse matrices[END_REF]. In our adaptation, we amalgamate each task vertex v z (ij) into either input-data vertex v x (j) or output-data vertex v y (i) according to the number of task vertices connected by n x (j) and n y (i), respectively. That is, v z (ij) is amalgamated into v x (j) if column j has a smaller number of nonzeros than row i, and otherwise it is amalgamated into v y (i), where the ties are broken arbitrarily. The result is a reduced hypergraph that contains only the input-and output-data vertices amalgamated with the task vertices where the weight of a data vertex is equal to the number of task vertices amalgamated into that data vertex. As a result, the locality constraint on vertex partitioning of the initial directed hypergraph naturally holds on any vertex partitioning on the reduced hypergraph. It so happens that after this process, the net directions become irrelevant for partitioning, and hence one can use the standard hypergraph partitioning tools.

Figure 5 illustrates how to obtain a local fine-grain partition through the described task-vertex amalgamation procedure. In Figure 5a, the up and left arrows imply that This manuscript is for review purposes only. We emphasize here that the reduced hypergraph constructed as above is equivalent to the hypergraph model of Pelt and Bisseling [START_REF] Pelt | A medium-grain method for fast 2d bipartitioning of sparse matrices[END_REF]. In that original work, the use of this model was only for two-way partitioning (of the fine grain model) which is then used for K-way fine-grain partitioning recursively. But this distorts the locality of task vertices so that a partition obtained in further recursive steps is no more a local fine-grain partition. That is why the adaptation was necessary.

Nonzero distribution to minimize the total communication vol-

ume. This method is composed of two parts. The first parts finds a vector distribution (Π(x), Π(y)). The second part finds a nonzero/task distribution Π(A) that exactly minimizes the total communication volume over all possible local fine-grain partitions which abide by (Π(x), Π(y)) of the first part. The first part can be accomplished by any conventional data partitioning method such as 1D partitioning.

Therefore, this section is devoted to the second part.

Consider the block structure (3) of A induced by (Π(x), Π(y)). Recall that in a local fine-grain partition, due [START_REF] Kaya | Analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF], the nonzero/task distribution is such that each

This manuscript is for review purposes only. k . This corresponds to assigning each nonzero of A kk to P k , for each diagonal block A kk , and assigning each nonzero of A k to either P k or P . Figure 6 illustrates a sample 10×12 sparse matrix and its block structure induced by a sample 3-way vector distribution which incurs four messages:

from P 3 to P 1 , from P 1 to P 2 , from P 3 to P 2 , and from P 2 to P 3 due to A 13 , A 21 , A 23 and A 32 , respectively.

Since diagonal blocks and zero off-diagonal blocks do not incur any communication, we focus on the nonzero off-diagonal blocks. Consider a nonzero off-diagonal block A k which incurs a message from P to P k . The volume of this message is determined by the distribution of tasks of A k between P k and P . This in turn implies that distributing the tasks of each nonzero off-diagonal block can be performed independently for minimizing the total communication volume.

In

(k) k , for a nonzero/task distribution A k = A (k) k + A ( )
k . Then, we can derive the following formula for the communication volume φ k from P to P k :

(12) φ k = n(A (k) k ) + m(A ( ) k ),
where n(•) and m(•) refer to the number of nonzero columns and nonzero rows of the input submatrix, respectively. The total communication volume φ is then computed by summing the communication volumes incurred by each nonzero off-diagonal block of the block structure. Then, the problem of our interest can be described as follows.

This manuscript is for review purposes only. (b) the induced block structure Fig. 6: A sample 10×12 sparse matrix A and its block structure induced by inputdata distribution Π(x) = {x (1) , x (2) , x (3) } and output-data distribution Π(y) = {y (1) , y (2) , y (3) }, where x (1) = {x 3 , x 4 , x 6 , x 7 , x 8 , x 9 }, x (2) = {x 2 , x 11 }, x (3) = {x 1 , x 5 , x 12 , x 10 }, y (1) = {y 4 , y 10 }, y (2) = {y 2 , y 3 , y 5 , y 6 , y 8 }, and y (3) = {y 1 , y 7 , y 9 }.

Problem 1. Given A and a vector distribution (Π(x), Π(y)), find a nonzero/task distribution Π(A) such that (i) each nonzero off-diagonal block has the form

A k = A (k) k + A ( ) k ; (ii) each diagonal block A kk = A (k)
kk in the block structure induced by (Π(x), Π(y)); and (iii) the total communication volume φ = k = φ k is minimized. A subset of vertices of a graph is called vertex cover if each of the graph edges is incident to any of the vertices in this subset. A vertex cover is minimum if its size is the least possible. In bipartite graphs, the problem of finding a minimum vertex cover is equivalent to the problem of finding a maximum matching [START_REF] Konig | Gráfok és mátrixok[END_REF]. Aschraft and Liu [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF] describe a similar application of vertex covers.

Let G k = (U k ∪ V k , E k ) be
Theorem 1. Let A k be a nonzero off-diagonal block and G k = (U k ∪ V k , E k )
be its bipartite graph representation.

For any vertex cover

S k of G k , there is a nonzero distribution A k = A (k) k + A ( ) k such that |S k | ≥ n(A (k) k ) + m(A ( ) k ),

For any nonzero distribution

A k = A (k) k + A ( ) k , there is a vertex cover S k of G k such that |S k | = n(A (k) k ) + m(A ( ) k ).
Proof. We prove the two parts of the theorem separately.

1) Take any vertex cover

S k of G k . Consider any nonzero distribution A k =
This manuscript is for review purposes only.

A (k) k + A ( ) k such that (13) a ij ∈      A (k) k if v j ∈ S k and u i ∈ S k , A ( ) k if v j ∈ S k and u i ∈ S k , A (k) k or A ( ) k if v j ∈ S and u i ∈ S k . Since v j ∈ S k for every a ij ∈ A (k) k and u i ∈ S k for every a ij ∈ A ( ) k , |S k ∩ V k | ≥ n(A (k) k ) and |S k ∩ U k | ≥ m(A ( )
k ), which in turn leads to ( 14)

|S k | ≥ n(A (k) k ) + m(A ( ) k ).
2) Take any nonzero distribution

A k = A (k) k + A ( ) k . Consider S k = {u i ∈ U k : a ij ∈ A ( ) k } ∪ {v j ∈ V k : a ij ∈ A (k) k } where |S k | = n(A (k) k ) + m(A ( ) k ). Now, consider a nonzero a ij ∈ A k and its corresponding edge {u i , v j } ∈ E k . If a ij ∈ A (k) k , then v j ∈ S k . Otherwise, u i ∈ S k since a ij ∈ A ( ) k . So, S k is a vertex cover of G k .
At this point, however, it is still not clear how the reduction from the problem of distributing the nonzeros/tasks to the problem of finding the minimum vertex cover holds. For this purpose, using Theorem 1, we show that a minimum vertex

* k | ≥ φ * k and |S * k | ≤ φ * k , respectively.
We decode an optimal nonzero/task distribution

A k = A (k) k + A ( )
k out of S * k according to [START_REF] Konig | Gráfok és mátrixok[END_REF] where one such distribution is ( 15)

A (k) k = {a ij ∈ A k : v j ∈ S * k } and A ( ) k = {a ij ∈ A k : v j ∈ S * k }.
Let φ k be the communication volume incurred by this nonzero/task distribution.

Then, |S * k | ≥ φ k due to [START_REF] Kuhlemann | Improving the communication pattern in matrix-vector operations for large scale-free graphs by disaggregation[END_REF], and Algorithm 3 gives a sketch of our method to find a nonzero/task distribution that minimizes the total communication volume based on Theorem 1. For each nonzero offdiagonal block A k , the algorithm first constructs G k , then obtains a minimum vertex cover S k , and then decodes S k as a nonzero/task distribution

φ k = φ * k since φ * k = |S * k | ≥ φ k ≥ φ * k .
A k = A (k) k + A ( ) k
according to [START_REF] Pelt | A medium-grain method for fast 2d bipartitioning of sparse matrices[END_REF]. Hence, the communication volume incurred by A k is equal to the size of the cover |S k |. In detail, each row vertex u i on the cover incurs an output

This manuscript is for review purposes only. communication of ŷi ∈ ŷ(k) , and each column vertex v j on the cover incurs an input communication of x j ∈ x( ) k . We recall that P sends ŷ(k) and x( ) k to P k in a single message in the proposed row-column-parallel sparse matrix-vector multiply algorithm.

u 2 u 3 u 5 u 8 v 3 v 4 v 6 v 7 v 8 v 9 u 6 S k`= {u 3 , u 6 , v 7 , v 8 } G kÀ k2 3
Algorithm 3 Nonzero/task distribution to minimize the total communication volume

1: procedure NonzeroTaskDistributeVolume(A, Π(x), Π (y)) 2:
for each nonzero off-diagonal block A k do See

Construct

G k = (U k ∪ V k , E k ) Bipartite graph representation 4: S k ← MinimumVertexCover(G k ) 5:
for each nonzero a ij ∈ A k do 6:

if v j ∈ S k then v j ∈ V k is a column vertex and v j ∈ S k 7: A (k) k ← A (k) k ∪ {a ij } 8: else u i ∈ U k is a row vertex and u i ∈ S k 9: A ( ) k ← A ( ) k ∪ {a ij } Figure 8
illustrates the steps of Algorithm 3 on the block structure given in Figure 6b. Figure 8a shows four bipartite graphs each corresponding to a nonzero off-diagonal block. In this figure, a minimum vertex cover for each bipartite graph is also shown. Figure 8b illustrates how to decode a local fine-grain partition from those minimum vertex covers. In this figure, the nonzeros are represented with the processor to which they are assigned. As seen in the figure, the number of entries sent from P 1 to P 2 is four, that is, φ 21 = 4, and the number of entries sent from P 3 to P 1 , from P 3 to P 2 and from P 2 to P 3 are all two, that is,

φ 13 = φ 23 = φ 32 = 2.
We note here that the objective of this method is to minimize the total communication volume under a given vector distribution. Since blocks of nonzeros are assigned, a strict load balance cannot be always maintained.

Related work.

Here we review recent related work on matrix partitioning for parallel SpMV.

Kuhlemann and Vassilevski [START_REF] Kuhlemann | Improving the communication pattern in matrix-vector operations for large scale-free graphs by disaggregation[END_REF] recognize the need to reduce the number of messages in parallel sparse matrix vector multiply operations with matrices corresponding

This manuscript is for review purposes only.

u 2 u 3 u 5 u 8 v 3 v 4 v 6 v 7 v 8 v 9 u 6 G 21 G 13 v 12 v 5 u 4 u 10 G 32 v 2 v 11 u 1 u 9 G 23 v 5 v 1 v 12 v 10 u 3 u 5 u 6 u 8 S 21 = {v 7 , v 8 , u 3 , u 6 } S 23 = {v 5 , u 5 } S 13 = {v 12 , u 10 } S 32 = {v 2 , u 9 }
(a) a minimum vertex cover for each nonzero off-diagonal block of Figure 6b. (b) a local fine-grain partition attained by optimal nonzero/task distribution.

Fig. 8: An optimal nonzero distribution minimizing the total communication volume obtained by Algorithm 3. The matrix nonzeros are represented with the processors they are assigned to. The total communication volume is 10, where P 1 sends [x 7 , x 8 , ŷ3 , ŷ6 ] to P 2 ; P 3 sends [x 12 , ŷ10 ] to P 1 ; P 3 sends [x 2 , ŷ9 ] to P 1 ; and P 3 sends [x 5 , ŷ5 ] to P 2 .

to scale-free graphs. They present methods to embed the given graph in a bigger one to reduce the number of messages. The gist of the method is to split a vertex into a number of copies (the number is determined with a simple calculation to limit the maximum number of messages per processor). In such a setting, the SpMV operations with the matrix associated with the original graph, y ← Ax, is then cast as triple sparse matrix vector products of the form y ← Q T (B(Qx)). This original work can be extended to other matrices (not necessarily symmetric, nor square) by recognizing the triplet product as a communication on x for duplication (for the columns that are split), communication of x vector entries (duplicates are associated with different destinations), multiplication, and as a communication on the output vector (for the rows that are split) to gather results. This exciting extension requires further analysis.

Boman et al. [START_REF] Boman | Scalable matrix computations on large 648 scale-free graphs using 2d graph partitioning[END_REF] propose a 2D partitioning method obtained by post-processing a 1D partition. Given a 1D partition among P processors, the method maps the P × P block structure to a virtual mesh of size P r × P c and reassigns the off-diagonal blocks so as to limit the number of messages per processor by P r + P c . The postprocessing is fast, and hence the method is as nearly efficient as a 1D partitioning method. However, the communication volume and the computational load balance obtained in the 1D partitioning phase are disturbed and the method does not have any means to control the perturbation. The proposed two-part method (Section 4.2), is similar to this work in this aspect; a strict balance cannot always be achieved; yet a finer approach is discussed in the preliminary version of the paper [START_REF] Kayaaslan | Semi-two-dimensional partitioning for parallel sparse matrix-vector multiplication[END_REF].

Pelt and Bisseling [START_REF] Pelt | A medium-grain method for fast 2d bipartitioning of sparse matrices[END_REF] propose a model to partition sparse matrices into two parts (which then can be used recursively to partition into any number of parts). The essential idea has two steps. First, the nonzeros of a given matrix A are split into two different matrices (of the same size as the original matrix), say A = A r + A c .

Second, A r and A c are partitioned together, where A r is partitioned rowwise, and A c is partitioned columnwise. As all nonzeros of A are in only one of A r or A c , the final result is a two-way partitioning of the nonzeros of A. The resulting partition on

This manuscript is for review purposes only.

in the experiments, we consider conformal partitioning of input and output vectors, by using vertex amalgamation of the input-and output-vector entries [START_REF] Uc | Revisiting hypergraph models for sparse matrix partitioning[END_REF]. We used PaToH [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF][START_REF] Aykanat | Patoh (partitioning tool for hypergraphs)[END_REF] with default parameters where the maximum allowable imbalance ratio is 3% for partitioning. We also notice that the 1.5D-V and 2D-B methods are based on 1D-H and keeps the vector distribution obtained from 1D-H intact. Hence, in the experiments, the input and output vectors for those methods are conformal as well.

Finally, since PaToH depends on randomization, we report the geometric mean of ten different runs for each partitioning instance.

In all experiments, we report the results using performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] which is very helpful in comparing multiple methods over a large collection of test cases. In a performance profile, we compare methods according to the best performing method for each test case and measure in what fraction of the test cases a method performs within a factor of the best observed performance. For example, a point (abscissa = 1.05, ordinate = 0.30) on the performance curve of a given method refers to the fact that for 30% of the test cases, the method performs within a factor of 1.05 of the best observed performance. As a result, a method that is closer to top-left corner is better.

In the load balancing performance profiles displayed in Figures 9b,9d, 10b and 10d, we compare performance results with respect to the performance of perfect balance instead best observed performance. That is, a point (abscissa = 6% and ordinate = 0.40) on the performance curve of a given method means that for 40% of the test cases, the method produces a load imbalance ratio less than or equal to 6%. As seen in Figure 9, for the set of matrices with no dense rows/columns, the relative performances of all methods are similar for K = 64 and K = 1024 in terms of both communication volume and load imbalance. As seen in Figures 9a and9c, all methods except the 1.5D-H method achieve a total communication volume at most 30% more than the best in almost 80% of the cases in this set of matrices. As seen in these two figures, the proposed 1.5D-V method performs significantly better than all other methods, whereas the 2D-H method is the second best performing method.

As also seen in the figures, 1D-H displays the third best performance, whereas 1.5D-H

shows the worst performance. As seen in Figures 9b and9d, in terms of load balance, the 2D-H method is the best performing method. As also seen in the figures, the proposed 1.5D-V method displays considerably worse performance than the others.

Specifically, all methods except 1.5D-V achieve a load imbalance below 3% in almost all test cases. In terms of the total communication volume, 2D checkerboard partitioning methods perform considerably worse than 1.5D-V, 2D-H and 1D-H methods. The first alternative 2D-B obtains better results than 2D-C. For load balance, 2D-C behaves similar to 1D-H, 2D-H and 1.5D-H methods except that 2D-C achieves a load imbalance below 5% (instead of 3%) for almost all instances. 2D-B behaves similar to 1.5D-V, and does not achieve a good load balance.

As seen in Figure 10, for the set of matrices with some dense rows/columns, all methods display a similar performance for K = 64 and K = 1024 in terms of the total communication volume. As in the previous dataset, in terms of the total communication volume, the 1.5D-V and 2D-H methods are again the best and second best methods, respectively, as seen in Figures 10a and10c. As also seen in these figures, 1.5D-H is the third best performing method in terms of the total communication vol-This manuscript is for review purposes only. row-column-parallel fine-grain [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF] row-column-parallel checker-board [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF] row-column-parallel checker-board [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF] 1D-H 2D-H 2D-B 2D-C row-column-parallel fine-grain [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF] row-column-parallel checker-board [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF] row-column-parallel checker-board [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF] 1D-H 2D-H 2D-B 2D-C In summary, the 1.5D-H method is a promising alternative for sparse matrices with dense rows/columns. It obtaines a total communication volume close to 2D-H, near-perfect balance, considerably lower message count than 2D-H, and has short partitioning time. The 1.5D-V method performs at the extremes: the best for the total communication volume, and the worst for the load balance, especially for matrices with dense rows/columns. Nevertheless, 1.5D-V could still be favorable to other methods for particular matrices due to lower communication volume. In short, if a sparse matrix contains dense rows/columns, then 1.5D-H seems to be the method of choice in general; otherwise, 1.5D-V and 1D-H are reasonable alternatives competing with each other. The 2D checkerboard based methods perform worse than the 1.5D methods, but they have good performance in terms of the message count based metrics. In particular, 2D-B is a fast method with a striking performance in reducing the latency, but load balance can be an issue. These could be deciding factors for large scale systems. On the other hand, 2D-C obtains better balance than 2D-B, but is slower.

7. Conclusion and further discussions. This paper introduced 1.5D parallelism for the sparse matrix-vector multiply (SpMV) operations. We presented the local row-column parallel SpMV that uses this novel parallelism. This multiply algorithm is the fourth one in the literature for SpMV in addition to the well-known 1D row-parallel, 1D column-parallel and 2D row-column-parallel ones. In this paper, we also proposed two methods (1.5D-H and 1.5D-V) to distribute tasks and data in accordance with the requirements of the proposed 1.5D parallel algorithm. Using a large set of matrices from the UFL sparse matrix collection, we compared the partitioning qualities of these two methods against the standard 1D and 2D methods.

The experiments suggest the use of the local row-column-parallel SpMV with a local fine-grain partition obtained by the proposed directed hypergraph model for matrices with dense rows/columns. This is because the performance of the proposed 1.5D partitioning is close to that of 2D fine-grain partitioning (2D-H) in terms of the partitioning quality, with considerably less number of messages and much faster execution.

We considered the problem mainly from a theoretical point of interest and leave the performance of 1.5D parallel SpMV algorithms in terms of the parallel multiply timings as a future work. We note that the main ideas behind the proposed 1.5D parallelism, such as heterogeneous messaging and avoiding nonlocal tasks by a locality constraint on partitioning, are of course not restricted to the parallel SpMV operation and these ideas can be extended to other parallel computations as well.

This manuscript is for review purposes only.

3. 1 .

 1 5D parallel sparse matrix-vector multiply. In this section, we propose the local row-column-parallel SpMV algorithm that exhibits 1.5D parallelism. The proposed algorithm simplifies the row-column-parallel algorithm by combining the two communication phases into a single expand-fold phase while attaining a flexibility on nonzero/task distribution close to the flexibility attained by the row-column-parallel algorithm.

(

  ) * k , where = k; (a) form vector x(k) which contains only those entries of x (k) corresponding to nonzero columns in A ( ) * k and (b) send vector x(k) to P , 2. for each nonzero row stripe A (k) * , where = k; compute (a) y ( )

  and the off-diagonal blocks of the column-stripe A (k) * k are input-output-local, output-local and input-local, respectively. Furthermore, due to (5) and (8), each off-diagonal block A k of the block structure (3) induced by the vector distribution (Π(x),Π(y)) becomes[START_REF] Kaya | Analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]

  -parallel sparse matrix-vector multiply. As there is no nonlocal tasks, the output-local tasks depend on input communication, and the output communication depends on the input-local tasks. Therefore, the tasks groups and communication phases can be arranged as: (i) input-local tasks; (ii) outputcommunication, input-communication; (iii) output-local tasks and input-output-local tasks. The input and output communication phases can be combined into the expandfold phase, and the output-local and input-output-local task groups can be combined into a single computation phase to simplify the overall execution.

  form vector ŷ( ) k , which contains only those entries of y ( )

Fig. 4 :

 4 Fig. 4: An illustration of attaining a local fine-grain partition through vertex partitioning of the directed hypergraph model that satisfies locality constraints. The input-and output-data vertices are drawn with triangles and rectangles, respectively.

Fig. 5 :

 5 Fig. 5: An illustration of local fine-grain partitioning through task-vertex amalgamations. The input-and output-data vertices are drawn with triangles and rectangles, respectively. The figure on the bottom right shows the fine-grain partition.

  the bipartite graph representation of A k , where U k and V k are the set of vertices corresponding to the rows and columns of A k , respectively, and E k is the set of edges corresponding to the nonzeros of A k . Based on this notation, the following theorem states a correspondence between the problem of distributing nonzeros/tasks of A k to minimize the communication volume φ k from P to P k and the problem of finding a minimum vertex cover of G k . Before stating the theorem we give a brief definition of vertex covers for the sake of completeness.

  cover of G k can be decoded as a nonzero distribution of A k with the minimum communication volume φ k as follows. Let S * k be a minimum vertex cover of G k and φ * k be the minimum communication volume incurred by a nonzero/task distribution of A k . Then, |S * k | = φ * k , since the first and second parts of Theorem 1 imply |S

Figure 7 k

 7 Figure 7 illustrates the reduction on a sample 5×6 nonzero off-diagonal block A k . The left side and middle of this figure respectively display A k and its bipartite graph representation G k , which contains 5 row vertices and 6 column vertices. On the middle of the figure, a minimum vertex cover S k that contains two row vertices {u 3 , u 6 } and two column vertices {v 7 , v 8 } is also shown. The right side of the figure displays how this minimum vertex cover is decoded as a nonzero/task distributionA k = A (k) k + A ( )k . As a result of this decoding, P sends [x 7 , x 8 , ŷ3 , ŷ6 ] to P k in a single message. Note that a nonzero corresponding to an edge connecting two cover vertices can be assigned to either A (k) k or A ( ) k without changing the communication volume from P to P k . The only change that may occur is in the values of partially computed output-vector entries to be communicated. For instance, in the figure, nonzero a 37 is assigned to A (k) k . Since both u 3 and v 7 are cover vertices, a 37 could be

Figures 9 and 10

 10 Figures 9 and 10 both display performance profiles of four task-and-data distribution methodsin terms of the total communication volume and the computational load imbalance. Figure 9 displays performance profiles for the set of matrices with no dense rows/columns, whereas Figure 10 displays performance profiles for the set of matrices containing dense rows/columns.

  without dense rows (K = 1024)

577

  The two 2D checkerboard methods display considerably worse performance than the 578 others (except 1D-H, which also shows a poor performance) in terms of the total com-579 munication volume. When K = 64, 2D-C shows an acceptable performance however 580 when K = 1024 its performance considerably deteriorates in terms of load balance. 581 2D-B obtains worse results. This not surprising since 2D-B is a modification of 1D-H 582 This manuscript is for review purposes only. without dense rows (K = 1024)

Fig. 10 :

 10 Fig. 10: Performance profiles comparing the total communication volume and the load balance on test matrices with dense rows/columns for K = 64 and 1024.

Figure 11d compares the

  Figure 11d compares the methods in terms of partitioning times for K = 1024.The run time of the 1.5D-V method involves the time spent for obtaining the vector distribution, which is the run time of the 1D-H method in our case. As seen in the figure, the 1D-H, 1.5D-V and 1.5D-H methods display comparable performances, whereas the 2D-H method takes significantly longer. The longer run time of 2D-H stems from the large size of the hypergraph model. 2D-B displays comparable performance (in terms of running-time) with that of 1D-H, 1.5D-V and 1.5D-H methods. Meanwhile, 2D-C is considerably slower than all others except 2D-H.

  k = {y 2 , y 3 , y 5 , y 6 , y 8 }x `= {x 3 , x 4 , x 6 , x 7 , x 8 , x 9 } Input-communication [x 7 , x 8 ]Fig.7:The minimum vertex cover model for minimizing the communication volume φ k from P to P k . According to the vertex cover S k , P sends [x 7 , x 8 , ŷ3 , ŷ6 ] to P k .

			3 4 6 7 8 9
		3
	3 4 6 7 8 9	5	A (k) k2
		6
		8
	5
	6
	8	2	3 4 6 7 8 9
		3
		5	A	(`) kÒ
		6
		8
		utput-communication [ŷ 3 , ŷ6 ]

y

1.5D PARALLEL SPARSE MATRIX-VECTOR MULTIPLY

This manuscript is for review purposes only.

A achieves load balance and reduces the total communication volume by the standard hypergraph partitioning techniques. Two-dimensional partitioning methods that bound the maximum number of messages per processor, such as the checkerboard [START_REF] Aykanat | A hypergraph-partitioning approach for coarse-grain decomposition[END_REF][START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF] and orthogonal recursive bisection [START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF] based methods, have been used in modern applications [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF][START_REF] Yoo | A scalable distributed parallel breadth-first search algorithm on bluegene/l[END_REF], sometimes without graph/hypergraph partitioning [START_REF] Yoo | A scalable eigensolver for large scale-free graphs using 2D graph partitioning[END_REF]. In almost all cases, inadequacy of 1D partitioning schemes are confirmed.

All previous work (including those that were summarized above) assumes the standard SpMV algorithm based on expanding x-vector entries, performing multiplies with matrix entries, and folding y-vector entries. Compared to all these previous work, ours has therefore a distinctive characteristic. In this work, we introduce the novel concept of heterogeneous messages where x-vector and partially computed y-vector entries are possibly communicated within the same message packet. In order to make use of this, we search for a special 2D partition on the matrix nonzeros in which a nonzero is assigned to a processor holding either the associated input-vector entry, or the associated output-vector entry, or both. The implication is that the proposed local row-column-parallel SpMV algorithm requires only a single communication phase (all the previous algorithms based on 2D partitions require two communication phases) as is the case for the parallel algorithms based on 1D partitions; yet the proposed algorithm achieves a greater flexibility to reduce the communication volume than the 1D methods.

6.

Experiments. We performed our experiments on a large selection of sparse matrices obtained from the University of Florida (UFL) sparse matrix collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF].

We used square and structurally symmetric matrices with 500-10M nonzeros. At the time of experiments, we had 904 such matrices. We discarded 14 matrices as they contain diagonal entries only, and we also excluded one matrix (kron g500-logn16) because it took extremely long to have a partition with the hypergraph partitioning tool used in the experiments. We conducted our experiments for K = 64 and K = 1024 and omit the cases when the number of rows is less than 50 × K. As a result, we had 566 and 168 matrices for the experiments with K = 64 and 1024, respectively.

We separate all our test matrices into two groups according to the maximum number of nonzeros per row/column, more precisely, according to whether the test matrix contains a dense row/column or not. We say a row/column dense if it contains at least 10 √ m nonzeros, where m denotes the number of rows/columns. Hence, for K = 64 and 1024, the first group respectively contains 477 and 142 matrices that have no dense rows/columns out of 566 and 168 test matrices. The second group contains the remaining 89 and 26 matrices, each having some dense rows/column, for K = 64 and 1024, respectively.

In the experiments, we evaluated the partitioning qualities of the local fine-grain partitioning methods proposed in Section 4 against 1D rowwise (1D-H [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF]), the 2D fine-grain (2D-H [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF]), and two checkerboard partitioning methods (2D-B [START_REF] Boman | Scalable matrix computations on large 648 scale-free graphs using 2d graph partitioning[END_REF], 2D-C [START_REF] Aykanat | A hypergraph-partitioning approach for coarse-grain decomposition[END_REF]).

For the method proposed in Section 4.1, we obtain a local fine-grain partition through the directed hypergraph model (1.5D-H) using the procedure described at the end of that subsection. For the method proposed in Section 4.2 (1.5D-V), the required vector distribution is obtained by 1D rowwise partitioning using the column-net hypergraph model. Then, we obtain a local fine-grain partition on this vector distribution with a nonzero/task distribution that minimizes the total communication volume.

The 1D-H, 2D-H, 2D-C and 1.5D-H methods are based on hypergraph models. Although all these models allow arbitrary distribution of the input-and output-vectors,

This manuscript is for review purposes only. row-column-parallel fine-grain [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel 652 sparse-matrix vector multiplication[END_REF] row-column-parallel checker-board [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF] row-column-parallel checker-board [START_REF] Ashcraft | Applications of the Dulmage-Mendelsohn decomposition and 645 network flow to graph bisection improvement[END_REF] 1D-H 2D-H 2D-B 2D-C