Karl Schlechta
email: schcsg@gmail.com-

A Comment on Argumentation

outline of a) new theory of argumentation.

Introduction 1.Abstract description

Argumentation is about putting certain objects together.

There are three things to consider:

(1) the objects themselves, and their inner structure (if they have any) -this inner structure may be revealed successively, or be immediately present,

(2) rules about how to put them together,

(3) avoid certain results (contradictions) in the resulting pattern.

To help intuition, we picture as result of an argumentation, an inheritance network the agents can agree on.

This network may consist of strict and defeasible rules only, with no elements or sets it is applied to. Think of the argumentation going on when writing a book about medical diagnosis. This will not be about particular cases, but about strict and default rules. "Sympton x is usually a sign of illness y, but there are the following exceptions: . . . " In addition, the network might contain cycles. There is nothing wrong with cycles. Mathematics is full of cycles, equivalences and their proofs. But consider also the following: We work in the set of adult land mammals. "Most elefants weigh more than 1 ton." "Most elements (i.e. adult land mammals) which weigh more than 1 ton are elefants." There is nothing in principle wrong with this either -except, in reality, we forgot perhaps about hippopotamus etc.

Arguments need not be contradictions to what exists already. They can be confirmations, elaborations, etc. For instance, we might have the default rule that birds fly, and clarify that penguins don't fly. This is not a contradiction, but an elaboration.

1.2

The structure of the objects Facts are either so simple that a dispute seems unreasonable. Or, they are a combination of basic facts and (default) rules, like, what I see through my microscope is really there, and not an artifact of some speck of dust on the lenses. For simplicity, facts will be basic, undisputable facts.

Expert opinion may be considered a default rule, where details stay unexplained, perhaps even unexplainable by the expert himself.

Rules (classical or defaults) have three aspects:

(1) the rule itself,

(2) the application of the rule,

(3) the result of the application of the rule.

Classical rules cannot be contested. We can contest their application, i.e. one of their prerequisites, or their result, and, consequently, their application. We can confirm their result by different means, likewise, their application.

Default rules are much more complicated, but not fundamentally different. Again, we can attack their application, by showing that one of the prerequisites does not hold, or, that we are in an (known) exceptional case. We can attack the conclusion, and, consequently, the rule, or its application. In particular, we may attack the conclusion, without attacking the application or the rule itself, by arguing that we are in a surprising exceptional case -and perhaps try to find a new set of exceptions. We can attack the default rule itself, as in the case of "normally, tigers are vegans". We can confirm a rule by confirming its conclusion, or adding a new rule, which gives the same result. We can elaborate a default rule, by adding an exception set, stating that all exceptions are known, and give the list of exceptions, etc. We can stop homogenousness (downward inheritance) e.g. for Quakers which are Republicans, we stop inheriting pacifism (or its opposite). This is not a contradiction to the default itself, but to the downward inheritance of the default (or to homogenousness) by meta-default, to be precise.

Obviously, the more we add (possible) properties to the objects (here default rules), the more we can attack, elaborate, confirm.

In the following section, we describe our general picture:

(1) there is an arbiter which checks for consistency, and directs the discussion,

(2) how to handle classical arguments and resulting contradictions,

(3) how to handle default arguments.

The classical part

We suppose there is an arbiter, whose role is to check consistency, and to authorise participants to speak.

If the arbiter detects an inconsistency, then he points out the "culprits", i.e. minimal inconsistent sets. As he detects inconsistencies immediately, the last argument will be in all those sets. The last argument need not be the problem, it might be one of the earlier arguments.

He asks all participants if they wish to retract one of the arguments involved in at least one minimal inconsistent set. (They have to agree unanimously on such retraction.) If there is no minimal inconsistent set left, the argumentation proceeds with the "cleaned" set of arguments, as if the inconsistency did not arise. Of course, arguments which were based on some of the retracted arguments are now left "hanging in the air", and may be open to new attacks.

If not, i.e. at least one minimal inconsistent set is left, the participants can defend (and attack) the arguments involved in those sets. The arbiter will chose the argument to be attacked/defended. See Example 2.1 (page 3) below. Suppose α is the argument chosen, then a defense will try to prove or argue for α, an attack will try to prove or argue for Q¬α, i.e. it is possible or consistent that ¬α. In particular, an attacker might try to prove ⊥, or some other unlikely consequence of α (and some incontested β's), and he need not begin with some α → γ, it might be a more roundabout attack.

If at least one minimally inconsistent set is left with all elements defended, then there is a deadlock, and the arbiter declares failure.

Consider

Example 2.1

We argue semantically. Let A := {x, a}, B := {x, b}, C := {x, c}, Y := {a, b, c}. Let Y be the last set added. For A, B, C, the situation is symmetrical.

Let Z = Z ′ be A, B, or C, then Y ∩ Z = ∅, but Y ∩ Z ∩ Z ′ = ∅, Z ∩ Z ′ = {x}, etc. Moreover, A ∩ B ⊆ C, etc.
Thus, A and B together are an argument for C, etc., so they argue for each other, and there is no natural way to chose any of A, B, C to be attacked. Thus, it is at the discretion of the parties involved (or the arbiter) to chose the aim of any attack -apart from Y, which is not supported by any of A, B, C. Still, Y might in the end be the strongest argument.

We may add D, E, with D := {x, d}, Y := {a, b, c, d} etc., the example may be extended to arbitrarily many sets.

At any moment, any argument can be attacked, not only if an inconsistency arises. We may continue an argumentation, even if not all minimally inconsistent subsets are treated as yet, but the arbiter has to keep track of them, and of the use of their elements. They and their consequences may still be questioned.

3 Defaults

The classical part of defaults

We see defaults primarily not as rules, but as relatively complicated classical constructions, which we may see as objects for the moment. The default character is in applying those objects, not in the objects themselves.

We follow here the theory described in Chapter 11 of [START_REF] Gabbay | A New Perspective on Nonmonotonic Logics[END_REF], and for convenience of the reader, we repeat in Section 5 (page 7) (essentially) Section 11.4 of [START_REF] Gabbay | A New Perspective on Nonmonotonic Logics[END_REF].

In our view, a (semantical) default (X : Y ) says:

(1) "most" elements of X are in Y,

(2) there may be exception sets X 1 , X 2 , etc. of X, where the elements are "mostly" not in Y (but X 1 ∪X 2 ∪. . . . has to be a "small" subset of X),

(3) in addition, there may be a "very small" subset X ′ ⊆ X, which contains "surprise elements" (i.e. not previously known exceptions), which are not in Y, (4) in addition, we may require that subsets of X "normally" behave in a homogenous way.

The notions of "most", "small" etc. are left open, a numerical interpretation suffices for the intuition. These notions are discussed in depth e.g. in [START_REF] Gabbay | Logical tools for handling change in agent-based systems[END_REF] and [START_REF] Gabbay | Conditionals and modularity in general logics[END_REF].

Introducing a default has to result in a (classically) consistent theory. E.g., it must not be the case that ∀x ∈ X.x ∈ Y, this contradicts the first requirement about defaults (and any reasonable interpretation of "most").

The default part of defaults

This leads to a hierarchy as defined in Section 5.1 (page 10) below. We use the hierarchy to define the use of the defaults.

To use the standard example with birds, penguins, fly, we proceed as follows. Suppose we introduce a bird x into the discussion. We try to put x as low as possible in the hierarchy, i.e. into the set of birds, but not into any known exception set, and much less into any "surprise" set. Only (classical) inconsistency, as checked by the arbiter, may force us to climb higher. Thus, unless there is a contradiction, we let x fly.

Attacks against defaults and their conclusions

Classical rules are supposed to be always true. Thus, classical rules themselves cannot be attacked, and an attack against a classical conclusion has to be an attack against one of its prerequisites.

Attacks against defaults can be attacks against

(1) the rule itself,

(2) one of the prerequisites,

(3) membership in or not in one of the exception sets, (4) membership in or not in the surprise set,

(5) perhaps even the notions of size involved, (6) etc.

Each component of a default rule may be attacked.

Comments

We assume that there is no fundamental difference between facts and conclusions: Usually, we were told facts, remember facts, have read facts, observed facts (perhaps with the help of a telescope etc.). These things can go wrong. Situations where things are obvious, and no error seems humanly possible, will not be contradicted.

Auxiliary elements

We now introduce some auxiliary elements which may help in the argumentation.

(1) "I agree."

This makes an error in this aspect less likely, as both parties agree -but still possible!

(2) "I confirm."

I am very certain about this aspect.

(3) Expert knowledge:

Expert knowledge and its conclusions act as "black box defaults", which the expert himself may be unable to analyse. Other experts (in the same field) will share the conclusion. (This is simplified, of course.) (One way to contest an expert's conclusion is to point out that he neglected an aspect of the situation, which is outside his expertise. His "language of reasoning" is too poor for the situation.)

(4) The arbiter may ask questions.

Examples of attacks

(1) Defaults:

Normally, there is a bus line number 1 running every 10 minutes between 10 and 11 in the morning.

Attack: No, the conclusion is wrong.

Question: Why?

Elaboration:

(1.1) No, the default is wrong (e.g.: it is line number 2 running every 10 minutes).

(1.2) Yes, but this is not homogenous, i.e. does not break down to subsets, and we know more. (For instance, we know that today is Tuesday or Wednesday, and it runs that often only Monday, Thursday, Friday, Saturday, Sunday -but we do not know this, only that is does not apply to all days of the week.)

(1.3) Yes, but today is an exception, and we know this. (e.g., we know that today is Tuesday, and we know that Tuesday is an exception.) (In addition, there might be exceptional Tuesdays, Christmas market day, etc. . . . )

(1.4) Yes, but I do not know why this is an exception. (This is a surprise case, I know about different days, but today should not be an exception, still I was just informed that it does not hold today.)

We do not attack the default, nor the applicability -but agree that it fails here.

(2) Classical conclusions:

From A and B, C follows classically.

Attack: C does not hold.

Question: Why? Elaboration:

(2.1) A does not hold or B does not hold, but I do not know which.

(2.2) A does not hold.

(2.3) B does not hold.

(2.4) A does not hold, and B does not hold.

(3) Fact: A holds.

Attacks: No, A does not hold.

Question: Why? Elaboration:

(3.1) You remember incorrectly.

(3.2) You did not observe well.

(3.3) Your observation tools do not work.

(3.4) You were told something wrong.

(3.5) etc.

(4) Expert knowledge, expert concludes that A.

Attack: A does not hold.

Question: Why?

Elaboration:

The situation involves aspects where you are not an expert. It is beyond your language. (Of course, the expert can ask for elaboration . . . .)

5 Appendix -Section 11.4 in [START_REF] Gabbay | A New Perspective on Nonmonotonic Logics[END_REF] Definition 5.1

If α ∼ | φ or α ∼ | φ, we say that the default φ is attached to α.

Given any fixed default theory, let A be the set of α (or M (α)), to which some default is attached

We work here in propositional logic, and on the semantic level.

We assume a classical background theory B, and a set of classical formulas α 1 , . . . , α n to which defaults are attached,

α i ∼ | φ i,1 , . . . , α i ∼ | φ i,ji
, where some or all of the ∼ | may also be ∼ | (without being necessarily ∼ | ¬).

Condition 5.1

We assume the following consistency conditions:

(1) B is classically consistent

(2) For each α i the defaults attached to α i together with B are jointly consistent.

In particular, the theory of defaults attached to one of the α i must be consistent. For instance, α i ⊢ α j and α i ∼ | ¬α j together are inconsistent. We thus rule out default theories like {: φ/φ, : ¬φ/¬φ}.

A negated default like α ∼ | φ needs a model of α ∧ ¬φ to be consistent, so we replace α ∼ | φ by α ∧ ¬φ for the consistency check.

The ⊢-relation between the α i 's gives a specificity relation by strict inclusion. We use this for the inheritance relation, and to solve (some) conflicts.

Definition 5.2

This definition describes how to obtain a consistent default theory at every point in the universe, using a theory revision approach, with specificity solving some conflicts. Of course, modifications are possible, but the general idea seems sound.

We define the set of valid defaults at some point. This influences the relation `as defined in Definition 5.4 (page 9), but not the relation ¡ as defined in Definition 5.5 (page 10), as the latter depends only on the sets to which defaults are attached, and not which defaults are attached.

We consider now some classical formula β -it may be one of the α i 's to which defaults are attached, or not.

(1) Visible defaults at β

The defaults visible at β : All defaults attached to some α i such that β ⊢ α i are considered visible at β.

(2) Valid defaults at β (2.1) The visible defaults are ordered by the ⊢ relation between the α's to which they are attached. The more specific ones are considered stronger defaults (for this β).

Of course, we can plug in here any partial relation, if it seems more suitable, e.g. ∼ | itself, as is often done in defeasible inheritance networks.

(2.2) Consider now the set of visible defaults, together with the classical information available at β. (If β ′ ⊢ β, but not conversely, we need not consider β ′ , etc.) If this set is inconsistent:

(2.2.1) Consider first in parallel all minimal inconsistent subsets involving classical information. They must contain at least one default, as the classical information was supposed to be consistent. Eliminate simultaneously from each such set the weakest (by the ⊢ relation) defaults. (If there are several weakest ones, eliminate all the weakest ones.) As contradictions involving classical, i.e., strongest information, seem to be more serious, we do these sets first.

(2.2.2) Consider now all remaining minimal inconsistent subsets of default information only. (Note that some might already have been eliminated in the previous step, but for other reasons.) Proceed as in the previous step, i.e. eliminate the weakest information. (2.2.3) We call the remaining defaults, visible at β, the defaults valid at β.

We will now work with the valid defaults only.

We now define "relevant" sets, sets where some default may change. We work in some fixed default theory, recall from Definition 5.1 (page 7) that A is the set of α or M (α) to which at least one default is attached.

Definition 5.3
Consider the model variant of A.

Let S(A) := { A ′ : A ′ ⊆ A} and U (A) := { A ′ : A ′ ⊆ A}, Let R := {X -Y : X -Y = ∅, X ∈ S(A), Y ∈ U (A)},
where X = U (the universe), and Y = ∅ are possible.

R ∈ R is called a relevant set (or formula).
Let R f be the ⊆-minimal elements of R -"f" for finest.

Obviously, the elements of R f are pairwise disjoint. The motivation for R f is the following. Let X -Y ∈ R f , then all elements of X -Y should satisfy all defaults valid at X, but no other, as they are in no other sets to which defaults are attached.

Example 5.1

Let A, A ′ ⊆ U, A ′′ ⊆ A ′ , let defaults be attached to A, A ′ , A ′′ , let all resulting intersections and set-differences be non-empty, if possible, e.g. A ′′ -A = ∅.

See Diagram 5.1 (page 9) for illustration.

So A = {A, A ′ , A ′′ }, R = {U, U -A, U -A ′ , U -A ′′ , U -A -A ′ , A, A -A ′ , A -A ′′ , A ′ , A ′ -A, A ′ -A ′′ , A ′ -A -A ′′ , A ′′ , A ′′ -A, A ∩ A ′ , A ∩ A ′ -A ′′ , A ∩ A ′′ } and R f = {U -A -A ′ , A -A ′ , A ′ -A -A ′′ , A ′′ -A, A ∩ A ′ -A ′′ , A ∩ A ′′ }
It is useful to code the sets A, A ′ , A ′′ by 3 bits, e.g. the left bit codes A, the middle one A ′ , the right one A ′′ .

We then have

R f = {U -A -A ′ = 000, A -A ′ = 100, A ∩ A ′ -A ′′ = 110, A ∩ A ′′ = 111, A ′′ -A = 011, A ′ -A -A ′′ = 010},
where the codes 001, 101 do not exist, as A ′′ ⊆ A ′ .

We will order the models by groups, and within groups by their "quality" satisfying defaults. We first define the latter relation, to be denoted `. It will be used to order the models within the sets ω(A), see Construction 5.1 (page 11).

Diagram 5.1

U A A ′ A ′′ µ ω
The dotted line separates ω(X) from µ(X)

Example, the sets

The construction

Note that, again, the following definition is independent of the actual defaults, and uses only the fact that defaults are attached to the A ∈ A. The order to be constructed, ¡, will be used as scaffolding in the construction of ≺ in Construction 5.1 (page 11), which orders packets of models.

(See Example 5.2 (page 12) and Diagram 5.1 (page 9) and Diagram 5.2 (page 12).)

Recall Definition 5.3 (page 8).

Definition 5.5

We order R f as follows:

X ¡ Y iff {A i : X ⊆ A i } ⊂ {A i : Y ⊆ A i } for X, Y ∈ R f ,
i.e. by the subset relation, of the A i they are in.

This order expresses exceptionality of the X ∈ R f . If X ¡ Y, then Y is a subset of more (by the subset relation) A i than X is, defaults were attached to the A i , so the A i are sets of exceptions. For instance, A 0 might be the set of birds, A 1 the set of penguins, being a penguin (an element of A 0 and A 1 ) is more exceptional than being a bird (an element of A 0 only).

Recall the coding of the X ∈ R f in Example 5.1 (page 8). This example is continued in Example 5.2 (page 12), and ¡ is then (the transitive closure of)

U -A -A ′ -A ′′ (000) ¡ A -A ′ (100) ¡ A ∩ A ′ -A ′′ (110) ¡ A ∩ A ′′ (111), U -A -A ′ -A ′′ (000) ¡ A ′ -A -A ′′ (010) ¡ A ′′ -A (011) ¡ A ∩ A ′′ (111), A ′ -A -A ′′ (010) ¡ A ∩ A ′ -A ′′ (110).
See Diagram 5.2 (page 12).

We are now ready to construct the preferential relation ≺ between model sets, it is implicitely extended to their elements.

Take e.g. in Example 5.2 (page 12) the set A ∩ A ′ -A ′′ . It will "see" the defaults valid for A ∩ A ′ , but not necessarily those for A ′′ . Every element of A ∩ A ′ -A ′′ can satisfy all those defaults, or only a part of them (or none). The best elements are those which satisfy all, the worst those which satisfy none. The precise relation is described by `, see Definition 5.4 (page 9). Let us call µ(X) the set of those elements which satisfy all defaults, the set of the others ω(X), for given X. Elements which do not satisfy all defaults are doubly exceptional, as they are not in any subclass where this failure is "explained". They are the "unexcused" exceptions, they have no reason to be exceptional, they are surprises, and doubly exceptional. So they should sit higher up in the hierarchy.

We do not think that there is a unique reasonable solution. Two ideas come to mind:

(1) Put all sets ω(X) above all other elements, in their own hierarchy, defined by ¡, and ordered internally by `. This is the radical approach.

(2) A less radical idea is to put them above the immediate ¡-successors of X. E.g. ω(A ′ -A -A ′′ ) will sit above A ∩ A ′ -A ′′ and A ′ ∩ A ′′ -A -or, rather, above µ(A ∩ A ′ -A ′′ ) and µ(A ′ ∩ A ′′ -A). This is what we will do here.

Construction 5.1

(See Example 5.2 (page 12) and Diagram 5.3 (page 13).)

We define a "packetwise" order ≺ between model sets, it will then be extended to the elements.

The relation ≺ has two parts:

(

1) µ(X) ≺ µ(Y ) iff X ¡ Y for X, Y ∈ R f (2) µ(X) ≺ ω(Y ) iff X ¡ Y or X = Y or X is a direct ¡-successor of Y for X, Y ∈ R f .
Note:

(1) We do not continue the order above the m ∈ ω(A), they are maximal elements -except for the interior order `among themselves.

(2) The µ(X) are "flat", all are best possible, There is no interior order in this set µ(X).

(3) The inside structure of the ω(X) is given by `according to Definition 5.4 (page 9) for this X.

(4) If any of the sets is empty, we just omit it and close under transitivity.

More formally, for X, X ′ ∈ R f , and element wise:

m ≺ m ′ :=               
there is X and m, m ′ ∈ ω(X) and m `m′ or there are X, X ′ and m ∈ µ(X) and m ′ ∈ µ(X ′ ) and X ¡ X ′ or there are X, X ′ and m ∈ µ(X) and m ′ ∈ ω(X ′ ) and ( X ¡ X ′ or X = X ′ or X is a direct ¡ -successor of X ′ )

The final order between models is thus basically lexicographic with two parts, first by ¡, then by `.

Example 5.2

We continue Example 5.1 (page 8), see also Diagram 5.1 (page 9).

Recall

  that A = {A, A ′ , A ′′ }, and R f = {U -A -A ′ -A ′′ = 000, A -A ′ = 100, A ∩ A ′ -A ′′ = 110, A ∩ A ′′ = 111, A ′′ -A = 011, A ′ -A -A ′′ = 010}.(1) We define ¡ by Definition 5.5 (page 10):U -A -A ′ -A ′′ = 000 ¡ A -A ′ = 100 ¡ A ∩ A ′ -A ′′ = 110 ¡ A ∩ A ′′ = 111, U -A -A ′ -A ′′ = 000 ¡ A ′ -A -A ′′ = 010 ¡ A ′′ -A = 011 ¡ A ∩ A ′′ = 111, A ′ -A -A ′′ = 010 ¡ A ∩ A ′ -A ′′ =110, closed by transitivity. See Diagram 5.2 (page 12). (2) We construct the relation between models according to Construction 5.1 (page 11). as follows (considering only the model sets, and neglecting the internal structure of the ω(X)). Constructing the full relation (and closing under transitivity) is then trivial. µ(000) ≺ µ(100) ≺ µ(110) ≺ µ(111), µ(000) ≺ µ(010) ≺ µ(011) ≺ µ(111), µ(010) ≺ µ(110), µ(100) ≺ ω(000), µ(010) ≺ ω(000), µ(110) ≺ ω(100), µ(110) ≺ ω(010), µ(011) ≺ ω(010), µ(111) ≺ ω(110), µ(111) ≺ ω(011).

FinallyU

  -A -A ′ -A ′′ = 000 A -A ′ -A ′′ = 100 A ′ -A -A ′′ = 010 A ∩ A ′ -A ′′ = 110 A ′ ∩ A ′′ -A = 011 A ∩ A ′ ∩ A ′′ = 111

Definition 5.4

Suppose we are at some relevant β, with the valid defaults φ 1 , . . . , φ n . We order the models of β according to satisfaction of the φ i . There are different possibilities to define `:

(1) By subsets: if the set of φ i satisfied by m is a subset of those satisfied by m ′ , then m ′ is better than m, m ′ `m.

(2) By cardinality: if the set of φ i satisfied by m is a smaller than the set of those satisfied by m ′ , then m ′ is better than m, m ′ `m.

(3) Some more complicated order, preferring certain defaults over others. (This might be interesting for contrary-to-duty obligations. Note that we solve here the case of additional information in real-life situations, too.)

(4) In particular, we may order the valid defaults by specificity of the α to which they are attached, and satisfy the most specific ones first, then the next, etc., resulting in a lexicographic order.

Note that the construction is more adapted to sets with rare exceptions, than to a classification, like vertebrates into fish, amphibians, reptiles, birds, and mammals, where all are exceptions -and those outside these subsets are the real exceptions.


E Agullo

L Giraud

A Guermouche

S Nakov

J Roman

Task-based Conjugate Gradient: from multi-GPU towards heterogeneous architectures

Keywords: High Performance Computing (HPC), multi-GPUs, heterogeneous architectures, task-based model, runtime system, sparse linear systems, Conjugate Gradient Calcul Haute performance (HPC), multi-GPUs, architectures hétérogènes, modéle en tâches, support d'exécution, systèmes linéaires creux, Gradient conjugué

Whereas most parallel High Performance Computing (HPC) numerical libaries have been written as highly tuned and mostly monolithic codes, the increased complexity of modern architectures led the computational science and engineering community to consider more modular programming paradigms such as task-based paradigms to design new generation of parallel simulation code; this enables to delegate part of the work to a third party software such as a runtime system. That latter approach has been shown to be very productive and efficient with compute-intensive algorithms, such as dense linear algebra and sparse direct solvers. In this study, we consider a much more irregular, and synchronizing algorithm, namely the Conjugate Gradient (CG) algorithm. We propose a task-based formulation of the algorithm together with a very fine instrumentation of the runtime system. We show that almost optimum speed up may be reached on a multi-GPU platform (relatively to the mono-GPU case) and, as a very preliminary but promising result, that the approach can be effectively used to handle heterogenous architectures composed of a multicore chip and multiple GPUs. We expect that these results will pave the way for investigating the design of new advanced, irregular numerical algorithms on top of runtime systems.

Introduction

In the last decade, the architectural complexity of High Performance Computing (HPC) platforms has strongly increased. To cope with this complexity, programming paradigms are being revisited. Among others, one major trend consists of writing the algorithms in terms of task graphs and delegating to a runtime system both the management of the data consistency and the orchestration of the actual execution. This paradigm has been first intensively studied in the context of dense linear algebra [START_REF] Agullo | LU factorization for accelerator-based systems[END_REF][START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF][START_REF] Agullo | Faster, Cheaper, Better -a Hybridization Methodology to Develop Linear Algebra Software for GPUs[END_REF][START_REF] Bosilca | Flexible Development of Dense Linear Algebra Algorithms on Massively Parallel Architectures with DPLASMA[END_REF][START_REF] Buttari | Parallel tiled QR factorization for multicore architectures[END_REF][START_REF] Kurzak | Scheduling dense linear algebra operations on multicore processors[END_REF][START_REF] Quintana-Ortí | Scheduling of QR factorization algorithms on SMP and multi-core architectures[END_REF][START_REF] Quintana-Ortí | Solving dense linear systems on platforms with multiple hardware accelerators[END_REF] and is now a common utility for related state-of-the-art libraries such as PLASMA, MAGMA, FLAME, DPLASMA and Chameleon. Dense linear algebra algorithms were indeed excellent candidates for pioneering in this direction. First, their regular computational pattern allows one to design very wide task graphs so that many computational units can execute tasks concurrently. Second, the building block operations they rely on, essentially level-three Basic Linear Algebra Subroutines (BLAS), are compute intensive, which makes it possible to split the work in relatively fine grain tasks while fully benefiting from GPU acceleration. As a result, these algorithms are particularly easy to schedule in the sense that state-of-the-art greedy scheduling algorithms may lead to a performance close to the optimum, including on platforms accelerated with multiple Graphics Processing Units (GPUs). Because sparse direct methods rely on dense linear algebra kernels, a large effort has been made to turn them into task-based algorithms [START_REF] Lacoste | Taking advantage of hybrid systems for sparse direct solvers via task-based runtimes[END_REF][START_REF] Agullo | Implementing multifrontal sparse solvers for multicore architectures with Sequential Task Flow runtime systems[END_REF].

In this paper, we tackle another class of algorithms, the Krylov subspace methods, which aim at solving large sparse linear systems of equations of the form Ax = b, where A is a sparse matrix. Those methods are based on the calculation of approximated solutions in a sequence of embedded spaces, that is intrinsically a sequential numerical scheme. Second, their unpreconditioned versions are exclusively based on non compute intensive kernels with irregular memory access pattern, Sparse Matrix Vector products (SpMV) and level-one BLAS, which need very large grain tasks to benefit from GPU acceleration. For these reasons, designing and scheduling Krylov subspace methods on a multi-GPUs platform is extremely challenging, especially when relying on a task-based abstraction which requires to delegate part of the control to a runtime system. We discuss this methodological approach in the context of the Conjugate Gradient (CG) algorithm on a shared-memory machine accelerated with multiple GPUs using the StarPU runtime system [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] to process the designed task graph. The CG solver is a widely used Krylov subspace method, which is the numerical algorithm of choice for the solution of large linear systems with symmetric positive definite matrices [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF].

The objective of this study is not to optimize the performance of CG on an individual GPU, which essentially consists of optimizing the matrix layout in order to speed up SpMV. We do not either consider the opportunity of reordering the matrix in order to improve the SpMV. Finally, we do not consider numerical variants of CG which may exhibit different parallel patterns. These three techniques are extremely important but complementary and orthogonal to our work. Instead, we rely on routines from vendor libraries (NVIDIA cuSPARSE and cuBLAS) to implement individual GPU tasks, we assume that the ordering is prescribed (we do not apply permutation) and we consider the standard formulation of the CG algorithm [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. On the contrary, the objective is to study the opportunity to accelerate CG on multiple GPUs by designing an appropriate task flow where each individual task is processed on one GPU and all available GPUs are exploited to execute these tasks concurrently. We first propose a natural taskbased expression of CG. We show that such an expression fails to efficiently accelerate CG. We then propose successive improvements on the task flow design to alleviate the synchronizations, exhibit more parallelism (wider graph) and reduce the volume of exchanges between GPUs.

The rest of the paper is organized as follows. We first propose a natural task-based expression of CG in Section 2. We then present the experimental set up in Section 3. We then show how the

1: r ← b 2: r ← r -Ax 3: p ← r 4: δ new ← dot(r, r) 5: δ old ← δ new 6: for j = 0, 1, ..., until b-Ax b ≤ eps do 7: q ← Ap /* SpMV */ 8:
α ← dot(p, q) /* dot */ baseline task-based expression can be enhanced for efficiently pipelining the execution of the tasks in Section 4. We present a performance analysis of a multi-GPU execution in Section 5. Section 6 presents concluding remarks together with preliminary experiments in the fully heterogeneous case.

Baseline Sequential Task Flow (STF) Conjugate Gradient algorithm

In this section, we present a first task-based expression of the CG algorithm whose pseudo-code is given in Algorithm 1a. This algorithm can be divided in two phases, the initialization phase (lines 1-5) and the main iterative loop (lines 6-16). Since the initialization phase is executed only once, we only focus on an iteration occurring in the main loop in this study. Three types of operations are used in an iteration of the algorithm: SpMV (the sparse matrixvector product, line 7), scalar operations (lines 9, 13, 14) and level-one BLAS operations (lines 8, 10, 11, 12, 15). In particular three different level-one BLAS operations are used: scalar product (dot, lines 8 and 12), linear combination of vectors (axpy, lines 10, 11 and 15) and scaling of a vector by a scalar (scal , line 15). The scal kernel at line 15 is used in combination with an axpy. Indeed, in terms of BLAS, the operation p ← r + βp consists of two successive operations: p ← βp (scal ) and then p ← r + p (axpy). In our implementation, the combination of these level-one BLAS operations represents a single task called scale-axpy. The key operation in an iteration is the SpMV (line 7) and its efficiency is thus critical for the performance of the whole algorithm.

According to our STF programming paradigm, data need to be decomposed in order to provide opportunities for executing concurrent tasks. We consider a 1D decomposition of the Inria sparse matrix, dividing the matrix in multiple block-rows. The number of non-zero values per block-rows is balanced and the rest of the vectors follows the same decomposition.

After decomposing the data, tasks that operate on those data can be defined. The tasks derived from the main loop of Algorithm 1a are shown in Figure 1b, when the matrix is divided in six block-rows. Each task is represented by a box, named after the operation executed in that task, and edges represent the dependencies between tasks.

The first instruction executed in the main loop of Algorithm 1a is the SpMV . When a 1D decomposition is applied to the matrix, dividing it in six parts implies that six tasks are submitted for this operation (the green tasks in Figure 1b): q i ← A i p, i ∈ [START_REF] Agullo | LU factorization for accelerator-based systems[END_REF][START_REF] Bosilca | Flexible Development of Dense Linear Algebra Algorithms on Massively Parallel Architectures with DPLASMA[END_REF]. For these tasks, a copy of the whole vector p is needed (vector p is unpartitioned). But in order to extract parallelism of other level-one BLAS operations where vector p is used (lines 8 and 15 in Algorithm 1a), in respect with our programming, the vector p needs to be partitioned. The partitioning operation is a blocking call; it thus represents a synchronization point in this task flow. Once vector p is partitioned, both vectors p and q are divided in six parts. Thus six dot tasks are submitted. Each dot operation accesses α in read-write mode, which induces a serialization of the operation. This sequence thus introduces new synchronizations in the task flow each time we need to perform a dot operation. The twelve axpy tasks (six at line 10 and six at line 11) can then all be executed in parallel. Another dot operation is then performed (line 12) and induces another serialization point. After the scalar operations at lines 13 and 14 in Algorithm 1a, the last scale-axpy operation of the loop is executed, which updates the vector p. At this stage, the vector is partitioned in six pieces. In order to perform the SpMV tasks for the next iteration, an unpartitioned version of the vector p is needed. This is done with the unpartition operation, similar to the partition operation, which is a blocking call.

All in all, this task flow contains four synchronization points per iteration, two for the partition/unpartition operation and two issued from the dot operations. The task flow is also very thin. Section 4.1 exhibits the induced limitation in terms of pipelining, while Sections 4.2, 4.3 and 4.4 propose successive improvements allowing us to alleviate the synchronizations and design a wider task flow, thus increasing the concurrency and the performance.

Experimental setup

All the tests presented in Section 5 have been run on a cache coherent Non Uniform Memory Access (ccNUMA) machine with two hexa-core processors Intel Westmere Xeon X5650, each one having 18GB of RAM, for a total of 36GB. It is equipped with three NVIDIA Tesla M2070 GPUs, each one equipped with 6GB of RAM memory. The task-based CG algorithm proposed in Section 2 is implemented on top of the StarPU v1.2. We use the opportunity offered by StarPU to control each GPU with a dedicated CPU core.To illustrate our discussion we consider the matrices presented in Table 1. All needed data is prefetched to the target GPU before the execution and assement of all the results presented in this paper.

Scheduling and mapping strategy. As discussed in Section 2, the task flow derived from Algorithm 1a contains four synchronization points per iteration and is very thin, ensuring only a very limited concurrency. Pipelining this task flow efficiently is thus very challenging. In particular, dynamic strategies that led to close to optimum scheduling in dense linear algebra [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] are not well suited here. We have indeed experimented such a strategy (Minimum Completion Time (MCT) policy), but all studied variants failed to achieve a very high performance. We have thus implemented a static scheduling strategy. We perform a cyclic mapping of the block-rows on the available GPUs in order to ensure load balancing. Building block operations. In order to explore the potential parallelism of the CG algorithm, we first study the performance of its building block operations, level-one BLAS and SpMV . The granularity does not penalize drastically the performance for SpMV operation. Additionally when three GPUs are used, a speed-up of 2.93 is obtained. On the other hand, in order to efficiently exploit multiple GPUs, vector with sizes of at least few millions are needed.

Achieving efficient software pipelining

In accordance with the example discussed in Section 2, the matrix is split in six block-rows and three GPUs are used. We pursue our illustration with matrix 11pts-128.

Assessment of the proposed task-based CG algorithm

Figure 3 shows the execution trace of one iteration of the task flow (Figure 1b) derived from Algorithm 1a with respect to the mapping proposed in Section 3. An iteration starts with the execution of a SpMV operation (line 7 in Algorithm 1a) which corresponds time interval [t 0 ,t 1 ] in Figure 3. Following the cyclic mapping strategy presented in Section 3, each GPU is thus in charge of two SpMV tasks. At time t 1 , the vector p is unpartitioned. The vector p is partitioned into six p i pieces, i ∈ [START_REF] Agullo | LU factorization for accelerator-based systems[END_REF][START_REF] Bosilca | Flexible Development of Dense Linear Algebra Algorithms on Massively Parallel Architectures with DPLASMA[END_REF], with respect to the blockrow decomposition of the matrix. However, this data partitioning operation is a blocking call (see Section 2) which means that no other task can be submitted until it is completed at time t 1 (the red vertical bar after the SpMV tasks in Figure 1b). Once vector p is partitioned, tasks for all remaining operations (lines 8-15) are submitted. The dot tasks are executed sequentially with respect to the cyclic mapping strategy. The reason for this, as explained in Section 2, is that the scalar α is accessed in read-write mode. In addition, α needs to be moved to GPUs between each execution of a dot task (time interval [t 1 , t 2 ] in Figure 3). Once the scalar product at line 8 is computed, the scalar division follows (line 9) executed on GPU 1 (respecting the task flow in Figure 1b). The execution of the next two instructions follows (lines 10 and 11). But before the beginning of the execution of the axpy tasks on GPU 2 and GPU 3, the new value of α is sent (the purple period at t 2 in Figure 3). The axpy tasks (yellow tasks in Figure 1b) are then executed during the period [t 2 , t 3 ] in parallel. The scalar product at line 12 is then executed during the time interval [t 3 , t 4 ], following the same sequence as explained above for line 8. Next, β and δ old are computed on GPU 1 at time t 4 in Figure 3, representing the scalar operations from lines 13 and 14 of Algorithm 1a. Tasks related to the last operation of the iteration (scale-axpy tasks in Figure 1b) are then processed during the time interval [t 4 ,t 5 ]. When all the new vector blocks p i are calculated, the vector p is unpartitioned (red vertical bar after the scale-axpy tasks in Figure 1b). As explained in Section 2, this data unpartition is another synchronization point and may only be executed in the RAM. All blocks p i of vector p are thus moved by the runtime system from the GPUs to the RAM during the time interval [t 5 ,t 6 ] for building the unpartitioned vector p. This vector is then used to perform the q i ← A i × p tasks related to the first instruction of the next iteration (SpMV at line 7). We now understand why the iteration starts with an active waiting of the GPUs (purple parts before time t 0 ): vector p is only valid in the RAM and thus needs to be copied on the GPUs.

During the execution of the task flow derived from Algorithm 1a (Figure 1b), the GPUs are idle during a large portion of the time (red and purple parts in Figure 3). In order to achieve more efficient pipelining of the algorithm, we present successive improvements on the design of the task flow: relieving synchronization points (Section 4.2), reducing volume of communication that is achieved using a packing data mechanism (Section 4.3) and relying on a 2D decomposition (Section 4.4).

Relieving synchronization points

Alternatively to the sequential execution of the scalar product, each GPU j can compute locally a partial sum (α j ) and perform a StarPU reduction to compute the final value of the scalar (α = n_gpus j=1 α j ). Figure 4a illustrates the benefit of this strategy. The calculation of the scalar product, during the time interval [t 0 , t 1 ] is now performed in parallel. Every GPU is working on its own local copy of α and once they have finished, the reduction is performed on GPU 1 just after t 1 .

t 0 t 1 t 2 t 3 t 4 RAM{ GPU 1{ GPU 2 GPU 3 { { (a) t 0 t 1 t 2 t 3 RAM { GPU 1{ GPU 2 GPU 3 { { (b)
Figure 4: Execution trace of one iteration when the dot is performed in reduction mode (left) and after furthermore avoiding data partitioning and unpartitioning (right).

The partition (after instruction 7 of Algorithm 1a) and unpartition (after instruction 15) of vector p, that are special features of StarPU, represent two of the four synchronization points within each iteration. They furthermore induce extra management and data movement costs. Indeed, after instruction 15, each GPU owns a valid part of vector p. For instance, once GPU 1 has computed p 1 , StarPU moves p 1 to the RAM and then receives it back. Second, vector p has to be fully assembled in the main memory (during the unpartition operation) before prefetching a copy of the fully assembled vector p back to the GPUs (after time t 3 in Figure 4a). We have designed another scheme where vector p is kept by StarPU in a partitioned form all along the execution (it is thus no longer needed to perform partitioning and unpartitioning operations at each iteration). Instead of building and broadcasting the whole unpartitioned vector p, each GPU gathers only the missing pieces. This enables us to "remove" the two synchronization points related to the partition and unpartition operations, since they are not called anymore, and decrease the overall traffic. Figure 4b illustrates the benefits of this policy. Avoiding the unpartitioning operation allows us to decrease the time required between two successive iterations from 8.8 ms to 6.6 ms. Furthermore, since the partitioning operation is no longer needed, the corresponding synchronization in the task flow control is removed. The corresponding idle time (red part at time t 0 in Figure 4a) is removed and instructions 7 and 8 are now pipelined (period [t 0 ,t 1 ] in Figure 4b).

Coming back to Figure 4a, one may notice that GPUs are idle for a while just before time t 1 and again just before time t 2 . This is due to the reduction that finalizes each dot operation (dot(p, q) at instruction 8 and dot(r, r) at instruction 12, respectively). In Algorithm 1a, vector x is only used at lines 10 (in read-write mode) and 6 (in read-only mode). The execution of instruction 10 can thus be moved anywhere within the iteration as long as the other input data of instruction 9, i.e. p and α have been updated to the correct values. In particular, instruction 10 can be moved after instruction 12. This delay enables StarPU to overlap the final reduction of the dot occurring at instruction 12 with the computation of vector x. The red part before t 2

Inria

in Figure 4a becomes (partially) green in Figure 4b. The considered CG formulation does not provide a similar opportunity to overlap reduction finalizing the dot operation at instruction 8.

Reducing communication volume by packing data

By avoiding data partition and data unpartition operations, the broadcast of vector p has been improved (from period [t 2 , t 4 ] in Figure 4a to period [t 3 , t 4 ] in Figure 4b), but still the communication time remains the large performance bottleneck (time interval [t 3 , t 4 ] in Figure 4b). This volume of communication can be decreased. Indeed, if a column within the block-row A i is zero, then the corresponding entry of p is not involved in the computation of the task q i ← A i p. Therefore, p can be pruned. We now explain how we can achieve a similar behavior with a task flow model. Instead of letting StarPU broadcast the whole vector p on every GPU, we can define tasks that only transfer the required subset. Before executing the CG iterations, this subset is identified with a symbolic preprocessing step. Based on the structure of the block A i,j , we determine which part of p j is needed to build q i . If p j is not fully required, we do not transfer it directly. Instead, it can be packed into an intermediate data, p i,j . StarPU provides an elegant support for implemented all these advanced techniques through the definition of new data types. We rely on that mechanism for implementing this packing scheme. Furthermore, the packing operation may have a non negligible cost whereas sometimes the values of p i,j that needs to be sent are almost contiguous. In those cases, it may thus be worth sending an extra amount of data in order to directly send the contiguous superset of p i,j ranging from the first to the last index that needs to be transferred. We have implemented such a scheme. A preliminary tuning is performed for each matrix and for each p i,j block to choose whether p i,j is packed or transferred in a contiguous way. Although StarPU can perform automatic prefetching, the prefetching operation is performed once all the dependencies are satisfied. In the present context, with the static mapping, this may be too late and further anticipation may be worthy. Therefore, we help the runtime system in performing data prefetching as soon as possible performing explicit prefetching within the callback of the scale-axpy task. We also do so after the computation of the α and β scalar values (lines 9 and 13 in Algorithm 1a) for broadcasting them on all GPUs.

RAM { GPU 1 { GPU 2 { GPU 3 { t 0 t 1 t 2 RAM { GPU 1{ GPU 2 GPU 3 { {
Figure 5 depicts the execution trace. The time interval [t 3 , t 4 ] in Figure 4b needed for the broadcasting of the vector p has been reduced to the interval [t 0 , t 1 ] in Figure 5. In the rest of the paper we refer to as the f ull algorithm when all the blocks are transferred, or to as the packed algorithm if this packing mechanism is used.

2D decomposition

The 1D decomposition scheme requires that for each SpMV task, all blocks of vector p (packed or not packed) are in place before starting the execution of the task. In order to be able to overlap the time needed for broadcasting the vector p (time interval [t 0 , t 1 ] in Figure 5), a 2D decomposition must be applied to the matrix. The matrix is first divided in block-rows, and then the same decomposition is applied to the other dimension of the matrix. Similarly as for a 1D decomposition, all the tasks SpMV associated with the entire block-row will be mapped on the same GPU. Contrary to the 1D decomposition, where we had to wait for the transfer of all missing blocks of the vector p, with the 2D decomposition, time needed for the transfer of the vector p can be overlapped with the execution of the SpMV tasks for which the blocks of the vector p are already available on that GPU. On the other hand, the 2D SpMV tasks yield lower performance then 1D (see Figure 6b and 2b), since they are executed on lower granularity. The result of the impact of a 2D decomposition is shown in Figure 6a. During the time interval [t 1 , t 2 ] in Figure 5 there is no communication, while in Figure 6a communications are overlapped with the execution of the SpMV tasks. In the rest of the paper we refer to either 1D or 2D depending on the data decomposition used. The trade-off between large task granularity (1D) and increased pipeline (2D) will be discussed in Section 5.

Performance analysis

We now perform a detailed performance analysis of the task-based CG algorithm designed above. We propose to analyze the behavior of our algorithms in terms of speed-up (S) and parallel efficiency (e) with respect to the execution occurring on one single GPU. In order to understand in more details the behavior of the proposed algorithms, we decompose the parallel efficiency into three effects, following the methodology proposed in [START_REF] Nakov | On the design of sparse hybrid linear solvers for modern parallel architectures[END_REF]: the effect of operating at a lower granularity (e granularity ), the effect of concurrency on the performance of individual tasks (e tasks ) and the effect of achieving a high task pipeline (e pipeline ). As shown in [START_REF] Nakov | On the design of sparse hybrid linear solvers for modern parallel architectures[END_REF], the following equality holds: e = e granularity × e tasks × e pipeline .

Inria

We observed that the efficiency of the task is maximum (e tasks = 1). Indeed, in a multi-GPU context, the different workers do not interfere with each other (they do not share memory or caches) and hence do not deteriorate the performance of one another. In the sequel, we thus only report on e granularity and e pipeline .

1D

2D 1D 2D 1.

Table 2 presents the performance achieved for all matrices in our multi-GPU context. The optimal performance is represented for each scheme in bold value. The first thing to be observed is that for all the matrices, the pack version of our algorithm where only just the needed part of the vector p is broadcasted, yields the optimal performance. Broadcasting entire sub-blocks is too expensive and thus considerably slows down the execution of the CG algorithm. For the matrices that have a regular distribution of the non zeros, i.e. the 11pts-256, 11pts-128 and the af_0_k101 matrices, the 1D algorithm outperforms the 2D algorithm. On the other hand, in the case of the Audi_kw matrix that has an unstructured pattern, the 2D algorithm which exhibits more parallelism, yields the best performance. Table 3 allows for analyzing how the overall efficiency is decomposed according to the metrics proposed above. Dividing the 11pts-256 matrix in several block-rows does not induce a penalty on the task granularity (e granularity = 0.99 ≈ 1). Furthermore, thanks to all the improvements of the task flow proposed in Section 4, a very high pipeline efficiency is achieved (e pipeline = 0.97), leading to an overall efficiency of the same (very high) magnitude. For the 11pts-128 matrix, the matrix decomposition induces a similar granularity penalty e granularity = 0.98. The slightly lower granularity efficiency is a direct consequence of the matrix order. For smaller matrices, the tasks are performed on smaller sizes, thus the execution time per task is decreased. This makes our algorithm more sensitive to the overhead created by the communications induced by the dot-products and the broadcasting of the vector p, ending up with a less optimal (but still high) pipeline efficiency (e pipeline = 0.86). The overall efficiency for this matrix is e = 0.84. This phenomenon is amplified when the matrix order is getting lower, such as in the case of the af_0_k101 matrix, resulting in a global efficiency of e = 0.65. The Audi_kw matrix yields optimal performance with the 2D algorithm (see Section 4.4). Although the 2D algorithm requires to split the matrix in many more blocks inducing a higher penalty on granularity (e granularity = 0.87), it allows for a better overlap of communication with computation ensuring that a higher pipeline (e pipeline = 0.91) is achieved. With this trade-off, an overall efficiency equal to e = 0.79 is obtained.

6 Towards a fully heterogeneous CG solver One advantage of relying on task-based programming is that the architecture is fully abstracted. We prove here that we can benefit from this design to run on an heterogeneous node composed of all available computational resources. Because the considered platform has 12 CPU cores and three GPUs, but that each GPU has a CPU core dedicated to handle it, we can only rely on 9 CPU workers and 3 GPU workers in total.

Figure 7 presents execution traces of preliminary experiments that rely on two different strategies for balancing the load between CPU cores and GPU. These traces show the ability of task-based programming in exploiting heterogeneous platforms. However, they also show that more advanced load balancing strategies need to be designed in order to achieve a better occupancy. This question has not been fully investigated yet and will be further investigated in a future work.

Inria

Figure 1 :

 1 Figure 1: Conjugate Gradient (CG) linear solver.

  Performance of the axpy kernel. The rest of the BLAS-1 kernels follow the same behavior.

  Performance of the SpMV kernel on the Audi_kw matrix.

Figure 2 :

 2 Figure 2: Performance of the building block operations used in the CG algorithm. All data is prefetched before execution and performance assessment.

Figure 3

 3 Figure3shows the execution trace of one iteration of the task flow (Figure1b) derived from Algorithm 1a with respect to the mapping proposed in Section 3. Figure3can be interpreted as follows. The top black bar represents the state of the CPU RAM memory during the execution. Each GPU is represented by a pair of bars, one for the state of the GPU and the black bar which

Figure 3 :

 3 Figure 3: Execution trace of an iteration with the CG task flow of Figure 1b using three GPUs.

Figure 5 :

 5 Figure 5: Execution trace when furthermore the vector p is packed.

Figure 6 :

 6 Figure 6: Execution trace when relying of a 2D decomposition of the matrix (left) and the performance of SpMV kernel when 2D decompostion is applied to the matrix (right).

  (a) nnz-based load balancing. (b) Performance model-based load balancing.

Figure 7 :

 7 Figure 7: Traces of an execution of one iteration of the CG algorithm in the heterogeneous case (9 CPU and 3 GPU workers) with different partitioning strategies for the Audi_kw matrix. In 7a, the nnz is balanced per block-row (33µs). In 7b a feed-back from a history based performance model is used for the partitioning of the matrix (16µs).

Table 1 :

 1 Overview of sparse matrices used in this study. The 11pts-256 and 11pts-128 matrices are obtained from a 3D regular grid with 11pt discretization stencil. The Audi_kw and af_0_k101 matrices come from structural mechanics simulations on irregular finite element 3D meshes.

	Matrix name	nnz	N	nnz/ N flop / iteration
	11pts-256	183 M	17 M	11	2 G
	11pts-128	23 M	2 M	11	224 M
	Audi_kw	154 M 943 K	163	317 M
	af_0_k101	18 M 503 K	34	38 M

Table 2 :

 2 Performance (in Gflop/s) of our CG algorithm for the matrix collection presented in Table

	# GPUs	full	pack.	full	pack. full pack. full pack.
			11pts-256		11pts-128
	1		9.74			9.58
	2	12.33 19.10 16.66 17.24 11.5	17.6	14.3 16.1
	3	11.70 28.39 13.26 23.17 9.01 24.2	9.22 20.6
			Audi_kw		af_0_k101
	1		10.0		9.84
	2	15.6	15.6	16.3	16.7 12.1	16.3	13.6 15.0
	3	17.7	20.0	22.0	23.6 11.1 19.4 12.5 18.2

Table 3 :

 3 Obtained speed-up (S), overall efficiency (e), effects of granularity on efficiency e granularity and effects of pipeline on efficiency e pipeline for matrices presented in Table1on 3 GPUs.

	Matrix	11pts-256 11pts-128 Audi_kw af_0_k101
	S	2.91	2.52	2.36	1.97
	e	0.97	0.84	0.79	0.65
	e granularity	0.99	0.98	0.87	0.96
	e pipeline	0.97	0.86	0.91	0.68

RR n°8912

http://dip.inria.fr RR n°8912

Inria

avenue de la Vielle Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

ACKNOWLEDGEMENT

The authors acknowledge the support by the INRIA-TOTAL strategic action DIP 1 and especially Henri Calandra who closely followed this work.


No El Nguyen
email: nnguyen@fapse.unige.ch

EPG Bidimensional Data Reduction

Keywords: electropalatography, data reduction, articulatory parameter

published or not. The documents may come

Introduction

Electropalatography (EPG) is a technique for gathering data about the dynamic variations of the tongue-palate contacts in the production of speech [START_REF] Hardcastle | New developments in electropalatography: a state-of-the-art report[END_REF][START_REF] Marchal | La Palatographie[END_REF]. In this domain, a great deal of attention has been recently devoted to the problem of data reduction [START_REF] Fontdevila | The contact index method of electropalatographic data reduction[END_REF][START_REF] Hardcastle | EPG data reduction methods and their implications for studies of lingual coarticulation[END_REF]. EPG provides us with a vast amount of information which most often needs to be compressed. In the Reading EPG system for example, which we have used in the present work, tongue-palate contact patterns are recorded using a 62-electrode arti cial plate, at a common sampling frequency of 200 Hz. When an experiment involves a large corpus and/or a large number of speakers, EPG recordings necessarily have to be transformed into a more manageable set of measurements. Moreover, besides its practical importance, EPG data reduction has also signi cant theoretical implications. Thus, previous studies on lingual movements in speech have led to the assumption that the tongue is characterized by a limited number of degrees of freedom. This hypothesis can be partially assessed by showing that it is possible to subdivide empirically the spatial arrangement of linguo-palatal contacts into a small number of independent components, by means of a statistical analysis [START_REF] Nguyen | Modeling tongue-palate contact patterns in the production of lingual consonants[END_REF]. In this perspective, EPG data reduction methods can allow us not only to cut down the amount of information to be processed, but also to better understand how tongue movements are controlled.

Many methods of EPG data reduction have already been proposed (see Hardcastle

Gibbon & Nicolaidis, 1991, for a review). To our knowledge however, no attempt has yet been made to apply to EPG data the techniques that are commonly used in the domain of image processing, although palatograms are images of a certain kind. In this paper, we present a new method for reducing EPG data based on such a technique, namely the bidimensional discrete cosine transform.

2 Representing an EPG pattern as a weighted sum of basis pictures

In this work, it is assumed that an EPG pattern can be regarded as a binary image, in which activated electrodes are represented as black dots and unactivated electrodes as white dots. Our method consists in extracting the main features of this image by means of a two-dimensional spectral transform. It is well-known that 1-D spectral transforms, such as the discrete Fourier transform (DFT), are commonly applied to speech sounds, to determine how the acoustic energy is distributed on a frequency scale. In such transforms, the signal is converted into a weighted sum of orthogonal basis vectors (sampled sine and cosine functions in the case of the DFT).

Analogously, an EPG image can be submitted to a 2-D transform which gives the possibility of characterizing the spatial arrangement of the tongue-palate contacts.

In two-dimensional processing, the signal is also represented as a set of orthogonal components, which are called basis pictures (hereafter BPs). Interestingly, it is possible to interpret some of these pictures in articulatory terms, as will be shown below.

In our approach, EPG images are characterized using the so-called discrete cosine transform (DCT). This transform is easier to handle than the DFT, as it only involves real numbers. When applied to a two-dimensional array of 8 8 elements, such as an EPG image 1 , the DCT is computed as follows [START_REF] Clarke | Transform Coding of Images[END_REF]:

C pq = C 0 p 2 8 7 X k=0 7 X l=0
f kl cos((2l + 1)q 16 ) cos((2k + 1)p 16 ) p = 0; : : :; 7; q = 0; : : :; 7 and more interestingly, Figure 2 shows that the three highest coe cients correspond to BPs that have been identi ed above as having articulatory signi cance, namely picture 1, picture 3, and picture 6. Conversely, note that picture 2, that we have also interpreted in articulatory terms (lateral asymmetry), was associated with a rather small DCT coe cient. This might indicate that, for our speaker at least, the linguopalatal contact patterns were relatively symmetrical on average. Finally, one other coe cient, namely coe cient 8, appears to account for a rather signi cant part of the image energy. Figure 1 indicates that this coe cient might re ect the presence of a semi-circular occlusion, such as that involved in the production of an alveolar stop, for example (although this will be true only if this coe cient has a negative value, since a positive value will correspond to a semi-circular pattern oriented in the opposite direction). In general, it seems plausible to assume that high DCT coe cients are related to articulatorily-relevant dimensions in an EPG pattern.

C 0 = 8 > > < > > : 1 p 2 if p = q = 0 1 otherwise In this formula, f
3

Conclusion

The main advantage of the method proposed in this work resides in the fact that the DCT BPs are completely orthogonal. In other words, multiplying these BPs with one another always gives an entirely white image as a result 2 . Therefore, the four indices that we have devised are related to dimensions which are completely independent from each other. In that sense, the information they provide us with about the spatial con guration of linguo-palatal contacts is maximized.

Another interesting feature is that our indices are quite robust, and are de ned in a way which does not depend on any speaker-speci c characteristic. In this respect, our method could be usefully employed in studies on defective utterances, for which more classical EPG parameters might prove to be more di cult to measure.

Footnotes 1 The Reading electropalate has 62 electrodes arranged in eight horizontal rows, with six electrodes only in the frontmost row. For the sake of simplicity however, an EPG \image" is considered here as having eight rows and eight columns (the two supernumary dots being always white).

Figure captions

1. First nine basis pictures for the discrete cosine transform (after [START_REF] Clarke | Transform Coding of Images[END_REF].

The pictures are ordered according to spatial frequency. White corresponds to numerical values comprised between 1 and 0; black corresponds to values between 0 and 1.

  Figure2gives rise to several remarks. First, the main part of the image energy resides in the lowest-order coe cients. Not surprisingly, this means that the con guration of the tongue-palate contacts most often formed a simple pattern, composed of one connected component only, with no gaps or irregular boundaries. Second,

  Figure 1, respectively, and they are computed as follows: SSU = 1 8 7 X k=0

2.

  Root mean square values of the DCT rst 48 coe cients calculated over a set of 2524 EPG frames. 3. EPG patterns for /k/ onset (upper panel) and for /l/ release (lower panel) in the production of the sequence /akla/. The normalized values of the indices are indicated above the corresponding basis pictures. Each picture has been weighted by the associated index value. 4. EPG patterns for /s/ in /asa/ (upper panel) and for / R / in /a R a/, with the corresponding index normalized values.

In mathematical terms, this means that the product of the two corresponding matrices is a matrix only composed of zeros.

Acknowledgements

This work was supported by the ESPRIT Working Group # 7098. Thanks are due to Paola Merlo for useful stylistic suggestions. We are also most grateful to Fiona Gibbon for helpful comments.