Olivier Benoist
email: olivier.benoist@ens.fr

Olivier Debarre
email: olivier.debarre@imj-prg.fr

À Claire

Smooth subvarieties of Jacobians

Keywords: Algebraic cycles, Jacobians, complex cobordism 2020 Mathematics Subject Classification. 14K12, 14C25, 57R77

We give new examples of algebraic integral cohomology classes on smooth projective complex varieties that are not integral linear combinations of classes of smooth subvarieties. Some of our examples have dimension 6, the lowest possible. The classes that we consider are minimal cohomology classes on Jacobians of very general curves. Our main tool is complex cobordism.

classes, which is generated by the cycle classes of codimension c algebraic subvarieties Y ⊂ X. Since the subvarieties Y may be singular, the following question going back to Borel and Haefliger [BH61, Section 5.17] naturally arises.

Question 1.1. Is the group H 2c (X, Z) alg generated by classes of smooth subvarieties of X?

The answer to Question 1.1 is obviously positive if d = 0 or if c ≤ 1. Further positive answers were obtained by Hironaka [Hir68, Theorem, Section 5, p. 50] when d ≤ min(3, c -1) and by Kleiman [Kle69,Theorem 5.8] when c = 2 and d ∈ {2, 3}; the answer is therefore positive for all n ≤ 5.

A first counterexample was then constructed by Hartshorne, Rees, and Thomas [HRT74, Theorem 1] when c = 2 and d ≥ 7 (and their method should yield counterexamples for all c ≥ 2 and d ≫ c). Other counterexamples were given by Debarre [Deb95, Théorème 6] when c = 2 and d ≥ 5, and by Benoist [Ben20, Theorem 0.3] when d ≥ c and α(c + 1) ≥ 3 (where α(m) denotes the number of ones in the binary expansion of m).

In this article, we build on the counterexample given in [START_REF] Debarre | Sous-variétés de codimension 2 d'une variété abélienne[END_REF]. There, Debarre considers the Jacobian X of a smooth projective complex curve C of genus n, polarized by its theta divisor class θ ∈ H 2 (X, Z). The minimal cohomology class θ c c! ∈ H 2c (X, Z) is the cycle class of the image W n-c (C) ⊂ X of the symmetric power C (n-c) by the Abel-Jacobi map, hence is algebraic. However, when n ≥ 2c + 2, the variety W n-c (C) is singular and there is no general reason why its class θ c c! should be a Z-linear combination of classes of smooth subvarieties of X. Debarre shows that this is indeed not the case if n ≥ 7, c = 2, and X is very general.

Our main theorem extends this result to many other values of (c, n). Recall that α(m) is the number of ones in the binary expansion of m. Theorem 1.2 (Theorem 3.7). Let (X, θ) be a very general complex Jacobian of dimension n. Let c be a nonnegative integer such that α(c + α(c)) > α(c) and n ≥ 4c -2. Then the integral classes λ θ c c! with λ odd are algebraic but are not Z-linear combinations of cycle classes of smooth subvarieties of X.

The weird condition α(c + α(c)) > α(c) springs naturally from our proof. It holds for integers c in the set {2, 4, 5, 8, 9, 12, 16, 17, . . . }. We have nothing to say when c = 3: we do not know if there exist Jacobians (X, θ) such that θ 3 3! is not a Z-linear combination of cycle classes of smooth subvarieties of X. Applied with c = 2 and n = 6, Theorem 1.2 implies the following result.

Corollary 1.3. There exists a smooth projective complex variety X of dimension 6 such that H 4 (X, Z) alg is not generated by classes of smooth subvarieties of X.

As noted above, the groups of algebraic cohomology classes of smooth projective varieties of dimension at most 5 are generated by classes of smooth subvarieties. The dimension of the variety X in Corollary 1.3 is thus the lowest possible. We do not know, however, if Question 1.1 has a positive answer if c = d = 3.

The proof of Theorem 1.2 when c = 2 and n ≥ 7 given in [START_REF] Debarre | Sous-variétés de codimension 2 d'une variété abélienne[END_REF] relies on a Barth-Lefschetz-type theorem for abelian varieties, on the Hirzebruch-Riemann-Roch theorem, and on the Serre construction. We can still rely on the same Barth-Lefschetz-type theorem under our more general hypotheses. However, computations based on the Hirzebruch-Riemann-Roch theorem become untractable in high codimension (and do not yield a proof when c = 2 and n = 6), and we cannot use the Serre construction in general because it is specific to codimension 2 subvarieties.

Both difficulties are overcome by resorting to complex cobordism, whose use in the theory of algebraic cycles was pioneered by Totaro [START_REF] Totaro | Torsion algebraic cycles and complex cobordism[END_REF]. We replace the integrality results derived from the Hirzebruch-Riemann-Roch theorem with divisibility properties of Chern numbers (Proposition 2.5). Those are obtained in Section 2 as a consequence of a detailed understanding of the structure of the complex cobordism ring. In [START_REF] Debarre | Sous-variétés de codimension 2 d'une variété abélienne[END_REF], the Serre construction was used to infer restrictions on the cohomology class of a smooth subvariety Y ⊂ X. Instead, we consider the class of Y in the complex cobordism of X and use the description of the complex cobordism of an abelian variety (in Sections 3.1-3.2). These tools are combined in Section 3.3 to prove Theorem 1.2.

In low codimension, a different method based on complex topological K-theory and on the Grothendieck-Riemann-Roch theorem, closer to the one used in [START_REF] Debarre | Sous-variétés de codimension 2 d'une variété abélienne[END_REF], gives small improvements on Theorem 1.2. We use this method in Section 4 to prove the following.

Theorem 1.4 (Theorem 4.3). Let (X, θ) be a very general complex Jacobian of dimension n. Then λ θ 4 4! is algebraic but not a Z-linear combination of classes of smooth subvarieties of X (a) if n ≥ 12 and λ is odd; (b) if n ≥ 14 and λ is not divisible by 4.

Conventions

A complex variety is a separated scheme of finite type over C. If M is a compact oriented manifold of dimension d, we let [M] ∈ H d (M, Z) be the fundamental class of M and denote by deg M :

H d (M, Z) → Z the morphism ω → deg(ω ⌢ [M]
). We let α(m) denote the number of ones in the binary expansion of m.

Congruences of Chern numbers

The Hurewicz morphism of MU

In this paragraph, we recall the structure of the Hurewicz morphism (2.1)

H : π * (MU) -→ H * (MU, Z)
of the spectrum MU representing complex cobordism. The computation of the cohomology of complex Grassmannians and the Thom isomorphism combine to show that H * (MU, Z) is a polynomial ring with integral coefficients with one generator in degree 2i for each i ≥ 1 (see [Ada74, Section I.3]). A deep theorem of Milnor shows that π * (MU) is also a polynomial ring with integral coefficients with one generator in degree 2i for each i ≥ 1 and that H is injective (see [Ada74, Theorem II.8.1 and Corollary II.8.11]). Quillen's theorem identifying π * (MU) with the coefficient ring of the universal formal group law [Ada74, Theorem II.8.2] and Hazewinkel's explicit construction of a universal formal group law [START_REF] Hazewinkel | Formal groups and applications[END_REF] allow us to be more precise.

Proposition 2.1. For each i ≥ 1, there exist u i ∈ π 2i (MU) and v i ∈ H 2i (MU, Z) such that

(a) π * (MU) = Z[u i ] i≥1 ; (b) H * (MU, Z) = Z[v i ] i≥1 ; (c) H(u i ) = λ i v i
, where λ i = p if i = p t -1 for some t ≥ 1 and some prime number p, and λ i = 1 otherwise.

Proof. Let u i ∈ π 2i (MU) be the polynomial generators of π * (MU) specified in [START_REF] Hazewinkel | Formal groups and applications[END_REF]Section 34.4.1]. Induction on i using the formula [START_REF] Hazewinkel | Formal groups and applications[END_REF](34.4.3)] shows that there exist v i ∈ H 2i (MU, Z) such that H(u i ) = λ i v i . It then follows from [Ada74, Lemmas II.7.9(iii) and II.8.10] that v i generates H 2i (MU, Z) modulo its decomposable elements. Consequently,

H * (MU, Z) = Z[v i ] i≥1 . □
For each e ≥ 1, define the ideal I e ⊂ π * (MU) to be the kernel of the composition

π * (MU) H ---→ H * (MU, Z) --→ H * (MU, Z)/2 e
of H and the reduction modulo 2 e . Working in the monomial bases associated with the u i and the v i given by Proposition 2.1 shows at once the following result.

Lemma 2.2. One has I 1 = ⟨2, u 2 t -1 ⟩ t≥1 and I e = (I 1 ) e .

Chern numbers

Let MU * be the complex cobordism ring, whose degree d elements are complex cobordism classes of d-dimensional compact stably almost complex manifolds. The Thom-Pontrjagin construction gives an identification

(2.2) ξ : MU * ∼ --→ π * (MU)
(apply [Koc96, Theorem 1.5.10] with B = BU).

Consider the polynomial ring Z[c j ] j≥1 in the Chern classes c j , graded so that c j has degree j. A degree i element P ∈ Z[c j ] j≥1 may be evaluated on a 2i-dimensional compact stably almost complex manifold M by setting P (M) deg M (P (c j (M))). This integer only depends on the complex cobordism class of M (see [Koc96, Section 4.3, p. 135], with B = BU and E the ordinary integral cohomology spectrum), so we get a morphism P : MU 2i → Z, called the Chern number associated with P . The following lemma is an immediate consequence of [Koc96, Proposition 4.3.8].

Lemma 2.3. Let i ≥ 0. A morphism MU 2i → Z is a Chern number if and only if it may be written as ψ • H • ξ for some group morphism ψ : H 2i (MU, Z) → Z, where H and ξ are as in (2.1) and (2.2).

We define the Segre classes s i ∈ Z[c j ] j≥1 to be the unique elements such that s i has degree i and

(2.3) j c j i s i = 1.
The Chern numbers associated with s i have the following properties.

Lemma 2.4. □

(a) If x ∈ MU 2i and x ′ ∈ MU 2i ′ , then s i+i ′ (xx ′ ) = s i (x)s i ′ (x ′ ). (b) For i ≥ 0 and h ≥ 1, the function s i : MU 2i → Z is divisible by 2 h if and only if α(i + h -1) > 2h -2. (c) For i ≥ 1, the function s i : MU 2i → Z

only takes even values and only takes values divisible by 4 on decomposable elements. (d) For

i = 2 t -1 ≥ 1 and u i ∈ MU 2i as in Proposition 2.1, s i (u i ) ≡ 2 (mod 4).

A congruence result for the top Segre class

The next proposition is the main goal of this section. It simultaneously generalizes the theorem of Rees and Thomas recalled in Lemma 2.4(b) (when e = 0) and [Ben20, Theorem 3.6] (when e = 1). Proposition 2.5. Fix e ≥ 0, h ≥ 1, and i ≥ 1. The following assertions are equivalent:

(i) There exists a degree i element Q ∈ Z[c j ] j≥1 such that the Chern number

s i + 2 h Q : MU 2i -→ Z only takes values divisible by 2 e+h . (ii) One has α(i + e + h -1) > e + 2h -2.
Proof. Assertion (i) implies that s i : MU 2i → Z is divisible by 2 h . So does assertion (ii) by Lemma 2.4(b). We may thus assume that the function

s i : MU 2i → Z is indeed divisible by 2 h .
Identify MU 2i and π 2i (MU) using (2.2). Consider the following statements:

(a) The function

) 2i . (f) If i = ( e ′ k=1 2 t k -1) + j for some e ′ ≤ e, then s j : MU 2j → Z is divisible by 2 h . The equivalence (i)⇔(a)

Smooth cycles on abelian varieties

Complex cobordism

We denote by MU * (X) the complex cobordism homology theory represented by the spectrum MU, evaluated on a topological space X. The group MU d (X) has a geometric interpretation as the group of complex cobordism classes of continuous maps f : M → X, where M is a d-dimensional compact stably almost complex manifold (see for instance [Swi75, Proposition 12.35]). We denote by [f ] ∈ MU d (X) the class of f . When X is a point, one recovers the ring MU * = π * (MU) which we studied in Sections 2.1 and 2.2. In general, the group MU * (X) is naturally an MU * -module.

If f : M → X is as above and if P ∈ Z[c j ] j≥1 of degree p and ω ∈ H d-2p (X, Z) are given, the integer

deg M P (c j (M)) • f * ω only depends on [f ] ∈ MU d (X) (the argument found in [Con79, Section 17, p. 54]
for unoriented cobordism also works for complex cobordism). It is called the characteristic number of f associated with P and ω. Fix a point o ∈ S 1 . The compatibility of the homology theory MU * with suspension shows that MU * (S 1 ) is free of rank 2 over MU * , generated by the classes of the inclusion {o} → S 1 and of the identity S 1 → S 1 (where S 1 is endowed with the stably almost complex structure induced by a trivialization of its tangent bundle). For N ≥ 0 and E ⊂ {1, . . . , N }, we consider the inclusion f E : T E → (S 1 ) N of the subtorus T E e∈E S 1 × e∈{1,...,N }\E {o}. Lemma 3.1. Fix N ≥ 1. When E describes the set of all subsets of {1, . . . , N },

(a) the f E, * [T E ] ∈ H |E| ((S 1 ) N , Z) form an additive basis of H * ((S 1 ) N , Z); (b) the [f E ] ∈ MU |E| ((S 1 ) N ) form an MU * -basis of MU * ((S 1 ) N ).
Proof. The first assertion follows from the Künneth formula in ordinary homology and the second assertion from the Künneth formula in complex cobordism (apply [Swi75, Theorem 13.75 i)] with E = MU). □ Lemma 3.2. Let x ∈ MU 2i , and let E ⊂ {1, . . . , N } be such that |E| = d -2i. Fix P ∈ Z[c j ] j≥1 of degree l and ω ∈ H d-2l ((S 1 ) N , Z). The characteristic number of x • [f E ] ∈ MU d ((S 1 ) N ) associated with P and ω is equal to 0 if l i and to

P (x) deg (S 1 ) N f E, * (1) • ω if l = i.
Proof. Let M be a 2i-dimensional compact stably almost complex manifold representing x ∈ MU 2i . Let g E : M ×T E → X be the composition of the second projection and f E : T E → X. We will also let h E : M ×T E → M denote the first projection. Since the stably almost complex structure on the tangent bundle on the torus T E is stably trivial, one has c j (T E ) = 0 for j > 0. It follows from the Whitney sum formula that

P (c j (M × T E )) = h * E P (c j (M)) ∈ H 2l (M × T E , Z). As a consequence, P (c j (M × T E )) • g * E ω ∈ H d (M × T E , Z) vanishes unless l = i. When l = i, we use the projection formula to compute deg M×T E P (c j (M × T E )) • g * E ω = deg M×T E h * E P (c j (M)) • g * E ω = deg M P (c j (M)) deg T E f * E ω = P (x) deg (S 1 ) N f E, * (1) • ω . □

Complex cobordism of abelian varieties

Now let X be a complex abelian variety of dimension n with a principal polarization θ ∈ H 2 (X, Z). We identify X and (S 1 ) N for N = 2n by means of a Lie group isomorphism X ≃ (S 1 ) N . By Lemma 3.1(a), there exists for each k ≥ 0 a unique Z-linear combination (3.1)

τ k = E⊂{1,...,N } |E|=2k µ E [f E ] ∈ MU 2k (X) such that |E|=2k µ E f E, * [T E ] ∈ H 2k (X, Z) is Poincaré-dual to the integral class θ n-k (n-k)! or, in other words, such that (3.2) |E|=2k µ E f E, * (1) = θ n-k (n -k)! .
Proposition 3.3. Let (X, θ) be a principally polarized complex abelian variety of dimension n. Assume that the group Hdg 2k (X, Z) of Hodge classes is generated by θ k k! for each k ≥ 0, and let τ k ∈ MU 2k (X) be as in (3.1). Let f : Y → X be a morphism of smooth projective varieties with Y of pure dimension d. Then there exists, for each i ∈ {0, . . . , d}, an element x i ∈ MU 2i such that

(3.3) [f ] = d i=0 x i • τ d-i ∈ MU 2d (X).
Proof. Let R i be the rank of the free Z-module MU 2i , and let (y i,r ) 1≤r≤R i be a Z-basis of it. Since the MU * -module MU * (X) is free with basis ([f E ]) E⊂{1,...,N } by Lemma 3.1(b), there exist unique integers ν i,r,E such that

(3.4) [f ] = d i=0 R i r=1         y i,r • |E|=2d-2i ν i,r,E [f E ]         ∈ MU 2d (X). Fix 0 ≤ i ≤ d and 1 ≤ r ≤ R i . As MU 2i ∼ --→ π 2i (MU) H --→ H 2i (MU, Z
) is an inclusion of free Z-modules of the same rank R i by Proposition 2.1, it follows from Lemma 2.3 that there exists a degree i element P ∈ Z[c j ] j≥1 such that P (y i,s ) is nonzero if and only if s = r. In view of Lemma 3.2 and the projection formula, the characteristic number of (3.4) associated with P and ω ∈ H 2d-2i (X, Z) reads

deg X f * (P (c j (Y ))) • ω = P (y i,r ) deg X         |E|=2d-2i ν i,r,E • f E, * (1) • ω         .
As ω is arbitrary, Poincaré duality on X implies

(3.5) f * (P (c j (Y ))) = P (y i,r ) |E|=2d-2i ν i,r,E f E, * (1) ∈ H 2n-2d+2i (X, Z).
The left side of (3.5) is algebraic, hence is a Hodge class. As a consequence, so is the right side. Since P (y i,r ) 0, we deduce from our hypothesis that the class |E|=2d-2i ν i,r,E f E, * (1) is an integral multiple of

θ n-d+i (n-d+i)! . Equation (3.2) now implies |E|=2d-2i ν i,r,E [f E ] = ξ i,
deg Y P (c j (Y )) • f * θ d-l (d -l)! = n d -l P (x l ).
Proof. Let us compute the characteristic number of f associated with P and θ d-l (d-l)! . Combining (3.3), (3.1), and Lemma 3.2 shows that it is

deg Y P (c j (Y )) • f * θ d-l (d -l)! = P (x l ) |E|=2d-2l µ E deg X f E, * (1) • θ d-l (d -l)! = P (x l ) deg X θ n-d+l (n -d + l)! • θ d-l (d -l)! .
Since deg X (θ n ) = n!, the proposition is proven. □

Smooth subvarieties of abelian varieties

The next proposition is an application of a Barth-Lefschetz-type theorem proved by Sommese [START_REF] Sommese | Complex subspaces of homogeneous complex manifolds. II. Homotopy results[END_REF].

Proposition 3.5. Let (X, θ) be a principally polarized complex abelian variety of dimension n such that Hdg 2k (X, Z) is generated by θ k k! for all k ≥ 0. Let f : Y → X be the inclusion of a smooth projective subvariety of pure codimension c. Assume that n ≥ 4c -2l for some l ≥ 1. Then there exist a 0 , . . . , a c-l , a c ∈ Z such that

(a) s i (Y ) = a i f * ( θ i i! ) for i ∈ {0, . . . , c -l, c}; (b) f * [Y ] is Poincaré-dual to a c θ c c! ∈ H 2c (X, Z). Proof. Since the subvariety Y ⊂ X is algebraic, the homology class f * [Y ]
bundle N Y /X is a c f * ( θ c c! ). Since the tangent bundle T X is trivial, one has c(N Y /X ) = c(Y ) -1 = s(Y ). This shows (a) for i = c.
If the abelian variety X were not simple, pulling back an ample divisor from a nontrivial quotient would produce a nonample divisor on X, contradicting the fact that Hdg 2 (X, Z) is generated by θ. We deduce that X is simple. One may thus apply [Som82, Corollary 3.5 and (3.6.1)] with B = Y to obtain π j (X, Y , y) = 0 for j ≤ n -2c + 1 and all y ∈ Y . It follows from the version [Hat02, Theorem 4.37] of the Hurewicz theorem, from the universal coefficient theorem, and from the long exact sequence of relative cohomology of the pair (X, Y ) that the restriction map H j (X, Z) → H j (Y , Z) is an isomorphism for j ≤ n -2c. For 0 ≤ i ≤ cl, one may apply this fact with j = 2i because n ≥ 4c -2l. This shows that the class s i (Y ) ∈ H 2i (Y , Z), which is Hodge because it is algebraic, is the restriction to Y of a class in Hdg 2i (X, Z). The latter is necessarily of the form a i

deg Y P (c j (Y )) • f * θ n-2c (n -2c)! = n n -2c P (x c ).
Using that Chern classes may be expressed as polynomials with integral coefficients in Segre classes by (2.3), it follows from Proposition 3.5(a) that

Q(c j (Y )) ∈ H 2c (Y , Z) is an integral multiple of f * ( θ c c! ), say Q(c j (Y )) = bf * ( θ c
c! ) for some b ∈ Z. Applying Proposition 3.5(a) again, we get P (c j (Y )) = (a c + 2b)f * ( θ c c! ). Rewriting the left side of (3.6) using the projection formula and Proposition 3.5(b), we obtain

deg X a c θ c c! • (a c + 2b) θ c c! • θ n-2c (n -2c)! = n n -2c P (x c ).
Using deg X (θ n ) = n!, we finally get (3.7) P (x c ) = 2c c a c (a c + 2b).

Our choice of Q implies that the left side of (3.7) is divisible by 2 α(c)+1 . The formula for the 2-adic valuation of the factorial given in [Rob00, Lemma, Section 5.3.1, p. 241] implies that the 2-adic valuation of 2c c is equal to α(c). We deduce that a c (a c + 2b) is even, hence that a c is even. □ Theorem 3.7. Let (X, θ) be a very general complex Jacobian of dimension n. Let c ≥ 0 be such that α(c + α(c)) > α(c) and n ≥ 4c -2. Then the classes λ θ c c! with λ odd are algebraic but are not Z-linear combinations of cycle classes of smooth subvarieties of X. (i) In Theorem 3.7, the hypothesis that (X, θ) is very general is only used to ensure that Hdg 2k (X, Z) is generated by θ k k! for k ≥ 0. As there exist Jacobians over Q whose Mumford-Tate group is the full symplectic group (see [And96, Théorème 5.2 3) and Remarque (vii) below it] which apply as all Hodge classes on abelian varieties are absolute Hodge by [Del82, Main Theorem 2.11]), one may find such an (X, θ) that is defined over Q.

Proof. The integral class

(ii) The proof of Theorem 3.7 actually shows that the class of any smooth subvariety of codimension c of X is divisible by 2 in the group H 2c (X, Z) alg .

Codimension 4 cycles

Theorem 3.7 is not optimal in several respects. When c is fixed, it is sometimes possible to give stronger restrictions on the cycle classes of smooth subvarieties of codimension c of X or results for lower values of the dimension n of X.

To obtain such improvements, one may work with complex topological K-theory instead of complex cobordism and replace the divisibility result for Chern numbers given in Proposition 2.5 with an application of the Grothendieck-Riemann-Roch theorem and an integrality property of the Chern character (cf. Lemma 4.1). This works well when c is low but becomes intractable for high values of c. We illustrate the method when c = 4.

We start with a lemma. Let X be a topological space, and let K * (X) = K 0 (X) ⊕ K 1 (X) be its Z/2graded complex topological K-theory defined in [AH61, Section 1.9]. We consider the Chern character ch :

K * (X) → H * (X, Q) as in [AH61, Section 1.10]. Lemma 4.1. For N ≥ 1, the Chern character ch : K * ((S 1 ) N ) → H * ((S 1 ) N , Q) is an isomorphism onto H * ((S 1 ) N , Z).
Proof. First, suppose that N = 1. The morphism ch : K 0 (S 1 ) → H 0 (S 1 , Q) = Q has image Z because it associates with each vector bundle its rank, and it is injective because K 0 (S 1 ) ≃ Z. That ch : K 1 (S 1 ) → H 1 (S 1 , Q) is an isomorphism onto H 1 (S 1 , Z) follows from the definition of this morphism using suspension and Bott periodicity (see [AH61, Section 1.10]) and from the fact the Chern character sends the Bott element which is a generator of K 0 (S 2 ) to a generator of H 2 (S 2 , Z) (see [AH61, Section 1.10, p. 206]).

It now follows from the Künneth formula in cohomology and in complex topological K-theory (for which see [Ati62, Lemma 1]) and from the multiplicativity of the Chern character that ch : K * ((S 1 ) N ) → H * ((S 1 ) N , Q) is an isomorphism onto H * ((S 1 ) N , Z). □ Proposition 4.2. Let (X, θ) be a principally polarized complex abelian variety of dimension n such that Hdg 2k (X, Z) is generated by θ k k! for all k ≥ 0. Let Y ⊂ X be a smooth projective subvariety of pure codimension 4. Proof. Let f : Y → X be the inclusion map. Proposition 3.5 shows the existence of integers a i such that s i (Y ) = a i f * ( θ i i! ) for i ∈ {0, 1, 2, 4} (if n ≥ 12), or i ∈ {0, 1, 2, 3, 4} (if n ≥ 14), and such that, in addition, the cohomology class of Y in X is equal to a 4

  Proof. For assertion (a), see[START_REF] Benoist | On the subvarieties with nonsingular real loci of a real algebraic variety[END_REF] Lemma 3.3]. Assertions (b) and (c) follow from a theorem of Rees and Thomas ([RT77, Theorem 3]; see [Ben20, Theorem 3.4 and Corollary 3.5] for these exact statements). The same result of Rees and Thomas ([RT77, Theorem 3] applied with r = 0 and n = 2 t -1) shows that not all the values of s i : MU 2i → Z are divisible by 4 if i = 2 t -1. This fact combined with (c) implies (d).

  is clear, and (a)⇔(b) follows from Lemma 2.3. The equivalence (b)⇔(c) holds because MU 2i is Z-free. The implication (c)⇒(d) is a consequence of the definition of I e , and the converse holds because Z/2 e is an injective Z/2 e -module [Wei94, Exercise 2.3.1]. The equivalence (d)⇔(e) is elementary. As for (e)⇔(f), it is a consequence of the description of I e given in Lemma 2.2 and of the properties of the Segre classes given in Lemma 2.4(a) and (d). Finally, applying Lemma 2.4(b) shows the equivalence (f)⇔(ii).□

  is Poincaré-dual to a Hodge class, which is necessarily of the form a c θ c c! by hypothesis. This proves (b). By the self-intersection formula [Ful98, Corollary 6.3], the top Chern class of the normal

  hypothesis. The proof is now complete. □ Proposition 3.6. Keep the hypotheses and notation of Proposition 3.5, assume that l = 1, and suppose in addition that α(c + α(c)) > α(c). Then a c is even. Proof. Let Q ∈ Z[c j ] j≥1 be the degree c homogeneous polynomial obtained by applying Proposition 2.5 with i = c, e = α(c), and h = 1. Applying Proposition 3.4 with l = c and P = s c + 2Q yields the identity (3.6)

  θ c c! is algebraic by [BL04, Poincaré's formula 11.2.1]. The hypothesis that Hdg 2k (X, Z) is generated by θ k k! for all k ≥ 0 is satisfied by [BL04, Theorem 17.5.1] and because the integral class θ k k! is primitive. One may thus combine Propositions 3.5(b) and 3.6 to show that the cycle classes of smooth codimension c subvarieties Y ⊂ X are even multiples of θ c c! . This concludes the proof. □ Remarks 3.8.

  (a) If n ≥ 12, the cohomology class of Y is an integral multiple of 2 θ 4 4! . (b) If n ≥ 14, the cohomology class of Y is an integral multiple of 4 θ 4 4! .

  As the class f * s 3 (Y ) is Hodge, we may also write f * s 3 (Y ) = b θ 7 7! for some b ∈ Z.As noted during the proof of Proposition 3.5, one has c(-N Y /X ) = s(Y ) -1 , hence c(-N Y /X ) = 1a 1 f * θ + (a 2 1a 2 /2)f * θ 2 + ((a 1 a 2a 3 1 )f * θ 3s 3 (Y )) + ((a 4 1 -3a 2 1 a 2 /2 + a 2 2 /4a 4 /24)f * θ 4 + 2a 1 s 3 (Y )f * θ) + • • • .One can then compute the Todd class (see [Ful98, Example 3.2.4]):(4.1) td(-N Y /X ) = 1a 1 f * θ/2 + (4a 2 1a 2 )f * θ 2 /24 + (a 1 a 2 -2a 3 1 )f * θ 3 /48 + (144a 4 1 -108a 2 1 a 2 + 12a 2 2 + a 4 )f * θ 4 /17 280a 1 s 3 (Y )f * θ/720 + • • • . Let χ ∈ K 0 (X)be the image in complex topological K-theory of the algebraic K-theory class f * [O Y ]. The Grothendieck-Riemann-Roch theorem [Ful98, Theorem 15.2] applied to the inclusion f and to the algebraic K-theory class [O Y ] shows that ch(χ) = f * td(-N Y /X ). Since f * f * is the cup-product with the cohomology

  s i 2 h : MU 2i → Z coincides modulo 2 e with a Chern number. (b) The function s i 2 h : MU 2i → Z/2 e factors through H 2i (MU, Z) via H. (c) The function s i 2 h : MU 2i → Z/2 e factors through H 2i (MU, Z)/2 e via H. (d) The function s i 2 h : MU 2i → Z/2 e vanishes on (I e ) 2i . (e) The function s i : MU 2i → Z only takes values divisible by 2 e+h on (I e

  r τ d-i for some ξ i,r ∈ Z.

	It finally follows from (3.4) that (3.3) holds with x i =	R i r=1 ξ i,r y i,r .

□ Proposition 3.4. Keep the hypotheses and notation of Proposition 3.3. Let P ∈ Z[c j ] j≥1 be of degree l for some 0 ≤ l ≤ d. Then,

Smooth subvarieties of Jacobians

The second-named author has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Project HyperK-grant agreement 854361).

class of Y which equals a 4 θ 4 4! , this identity may be combined with (4.1) to give (4.2) ch(χ) = a 4 θ 4 /24a 1 a 4 θ 5 /48 + (4a 2 1a 2 )a 4 θ 6 /576 + (a 1 a 2 -2a 3 1 )a 4 θ 7 /1 152 + (144a 4 1 -108a 2 1 a 2 + 12a 2 2 + a 4 )a 4 θ 8 /414 720ba 1 θ 8 /3 628 800 + • • • .

By Lemma 4.1, which applies because X is diffeomorphic to (S 1 ) 2n , the left side ch(χ) of (4.2) is an integral cohomology class, hence so is its right side.

Assume by way of contradiction that a 4 is odd. The integrality of the class ch 5 (χ) = a 1 a 4 θ 5 /48 implies that a 1 is even, and the integrality of the class ch 6 (χ) = (4a 2 1a 2 )a 4 θ 6 /576 shows that a 2 is divisible by 4. All terms appearing in the expression for ch 8 (χ) given in (4.2) are multiples of θ 8 8! with coefficients a rational number of nonnegative 2-adic valuation, with the exception of a 2 4 θ 8 /414 720 = 7a 2 4 72 θ 8

8! since a 4 is assumed to be odd. It follows that ch 8 (χ) is not an integral cohomology class, which gives a contradiction. This proves (a). Now suppose n ≥ 14. Then bθ 7 /7! = f * s 3 (Y ) = f * f * (a 3 θ 3 /3!) = a 3 a 4 θ 7 /144, so that b = 35a 3 a 4 . Assume by way of contradiction that a 4 is not divisible by 4. By the integrality of ch 5 (χ) = a 1 a 4 θ 5 /48, we see that a 1 a 4 is even. The integrality of ch 6 (χ) = (4a 2 1a 2 )a 4 θ 6 /576 shows that a 2 is even. These pieces of information imply that in the right side of the equality


Stephane Vialle

Sylvain Contassot-Vivier

Patrick Mercier

Generic algorithmic scheme for 2D stencil applications on hybrid machines

Hardware accelerators are classic scientific coprocessors in HPC machines. However, the number of CPU cores on the mother board is increasing and constitutes a non negligible part of the total computing power of the machine. So, running an application both on an accelerator (like a GPU or a Xeon-Phi device) and on the CPU cores can provide the highest performance. Moreover, it is now possible to include different accelerators in a machine, in order to support and to speedup a larger set of applications. Then, running an application part on the most suitable device allows to reach high performance, but using all unused devices in the machine should permit to improve even more the performance of that part. However, the overlapping of computations with inter-device data transfers is mandatory to limit the overhead of this approach, leading to complex asynchronous algorithms and multi-paradigm optimized codes. This article introduces our research and experiments on cooperation between several CPU and both a GPU and a Xeon-Phi accelerators, all included in a same machine.

Introduction and objectives

Hardware accelerators have become classical scientific coprocessors but they still need a CPU to control the entire machine. Also, current standard CPU have a significant number of cores and computing power. Using both CPU and accelerator cores seems an interesting way to achieve the maximal computing speed on a computing node. But using both the CPU and accelerators leads to frequently transfer intermediate results between all those devices. So, asynchronous data transfers, overlapped with computations, are required to obtain significant speedups compared to mono-device executions.

In this study, we consider hybrid computing nodes: each node hosting one or more CPU, one GPU and one Xeon-Phi, allowing one to choose the most adapted architecture for each application, or application module. However, once the most adapted device has been identified for a computation in an algorithm, it may be interesting to cooperatively use any other kind of computing device available in the machine to concurrently process that computation, and then to increase the overall performance. In order to ease the use of our multi-devices and multi-architectures machines, we developed a generic parallel algorithmic scheme achieving an overlapping of data transfers with computations between two or three devices (simultaneously), for 2D stencil applications. This generic scheme has been implemented and optimized for GPU and Xeon-Phi devices, considering the specific features of their programming models and API.

Our approach has been validated in [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF] with a Jacobi relaxation application. In this paper, we extend it to a more complex stencil application, more representative of real scientific applications: a shortest path computation. As in [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF], that application has been tested on two different hybrid machines. We were able to experiment our different kernels (on CPU, GPU and Xeon-Phi), using one, two or three devices simultaneously, and to identify the most suited combination for each application.

Application examples

Our parallel scheme (see Sec.4) aims at running 2D stencil applications, on hybrid machines. Two different applications have been used to test and validate our approach: a very classical 2D Jacobi relaxation algorithm, and a shortest path calculation on 2D+ ground (taking the elevations into account). They both work on 2D regular grids of points. The reader is invited to see [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF] for a detailed description of our Jacobi algorithm. Our shortest path algorithm computes minimal paths costs from a given position in the grid to all other positions. It is slightly different from the Jacobi algorithm as it is an adaptation of the Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]. The path cost of each point P i,j is updated according to the previous costs of its eight neighbors in the 2D grid (denoted N 8 (P i,j )). Indeed, for each neighbor, the 3D euclidian distance to the current point is computed and added to the current minimal cost of the neighbor. Then, the new cost of the current point is chosen as the minimal cost among the eight computations:

C n+1 (i, j) = min P a,b ∈N8(Pi,j ) (C n (a, b) + d(P i,j , P a,b ))
Considering h, the side length of a grid cell (between two successive points in the grid), the distance between P i,j and P a,b is: In fact, N 8 (P i,j ), the neighboring of point P i,j in the 2D grid (see Fig. 1), contains less than eight neighbors when P i,j is located on a grid edge. This algorithm requires a 2D array of costs and a 2D array of elevations (z). Initially, a target point T α,β is chosen and its cost is set to 0 while all others costs are set to +∞. Then the algorithm starts to iterate until no cost changes during an entire iteration (convergence is reached). Finally, we get a 2D array of minimal costs (distances) and paths from every point P i,j to the target point T α,β .

d(P i,j , P a,b ) = ((i -a).h) 2 + ((j -b).h) 2 + (z(i, j) -z(a, b))

Related work

Today, scientific computing on GPU accelerators is common, while using Xeon-Phi accelerators starts to be deeply investigated and some comparisons have been achieved. In [START_REF] Fang | Parallel computation of non-bonded interactions in drug discovery: Nvidia GPUs vs. Intel Xeon Phi[END_REF], authors point out the need to optimize data storage and data accesses in different ways on GPU and Xeon Phi. In [START_REF] Calandra | Evaluation of Successive CPUs/APUs/GPUs Based on an OpenCL Finite Difference Stencil[END_REF], authors optimize data storage for stencil applications on different CPU, GPU and APU (processors with both CPU and GPU cores).

Several studies like [START_REF] Gysi | MODESTO: Data-centric analytic optimization of complex stencil programs on heterogeneous architectures[END_REF][START_REF] Szustak | Adaptation of MPDATA heterogeneous stencil computation to Intel Xeon Phi coprocessor[END_REF] address the deep optimization of parallel stencil applications or other lattice based applications [START_REF] Wende | Swendsen-Wang Multi-cluster Algorithm for the 2D/3D Ising Model on Xeon Phi and GPU[END_REF], for different types of devices. However, their programs use only one device at a time. Some authors use a generic programming model and tool to develop on different architectures, like OpenCL [START_REF] Gaster | Heterogeneous Computing with OpenCL[END_REF][START_REF] Su | On the GPU-CPU Performance Portability of OpenCL for 3D Stencil Computations[END_REF]. Nevetheless, such tool does not hide to the programmer the devices hardware specificities to be taken into account to attain optimal performance. Thus, an important algorithmic effort is still required to design codes efficiently running on different architectures. Another approach consists in using the CPU memory to store data too large for the accelerator memory, and to efficiently send sub-parts of the problem to the accelerator [START_REF] Jin | Efficient utilization of memory hierarchy to enable the computation on bigger domains for stencil computation in CPU-GPU based systems[END_REF]. Here again, there is no simultaneous use of the CPU and the accelerators.

Finally, some works propose general frameworks to use different devices by performing dynamic scheduling of tasks graphs or data flows [START_REF] Courtès | C language extensions for hybrid CPU/GPU programming with StarPU[END_REF]. However, the generality of such frameworks induces additional costs compared to application specific schemes.

In [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF], we introduced our first algorithmic scheme of heterogeneous computing, successfully running a Jacobi relaxation on CPU, GPU and Xeon-Phi. In the current study, we generalize our approach to a more realistic scientific problem involving more complex and irregular computation.

Algorithmic scheme

In our context, the GPU and the Xeon-Phi are used as scientific co-processors in offload mode. The Xeon-Phi could have been used via MPI but we chose the offload mode to get similar programming paradigms between the GPU and the Phi. So, our hybrid solution uses the CPU to launch and control the computation steps on the GPU, the Xeon-Phi and its own cores, and to manage data transfers between CPU and devices.

Global strategy of the parallel scheme

Our global data distribution strategy is to partition the 2D grid in three horizontal strips (potentially of different heights), so that the first one is processed on the GPU, the second one, on the CPU, and the third one, on the Phi. Placing the CPU between the GPU and Phi in the logical organization is a strategic choice as direct data transfers between GPU and Phi are not currently possible. We name CPU boundaries the first and last lines computed by the CPU, GPU boundary the last line computed by the GPU, and MIC boundary the first line computed by the Xeon-Phi. We name corpus the other lines computed by a computing device. According to data dependencies in the application, each computing device must store, in addition to its strip of the grid, the adjacent boundary(ies) of its neighboring device(s). So, the top part of the grid is transferred to the GPU and the bottom part to the Xeon-Phi (see Fig. 2 left). Although the CPU memory may host only its associated part of the grid (and the required boundaries), for the sake of simplicity in the overall management of the system, we have chosen to store the entire current (crt) and previous (prev) relaxation grids. However, this could be easily optimized if the CPU memory were too small to store the entire grids.

At each iteration during the computation loop, the GPU boundary may be transferred to the CPU while the adjacent CPU boundary may be transferred to the GPU. Symmetrically, the CPU/Xeon-Phi boundaries may be transferred too. Indeed, as a frontier may remain unchanged during one iteration, our algorithm transfers a frontier to/from an accelerator only if it has undergone at least one modification. Also, in order to optimize the transfers, our algorithm is designed to allow direct transfers of the frontiers in their place in the local arrays on the destination device. So, no intermediate array is required to store the received frontiers coming from another device, and the CPU algorithm uses symmetric data structures and interactions for both accelerators.

Multi-device algorithm

Figure 2 (right) details the three main parts of our generic multi-device algorithm (see [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF] for a simplified version, with straightforward array allocation on accelerators and systematic frontier transfers):

Initialization step: Memory allocation on accelerators and initial data transfers from CPU to accelerators can be long, so the initialization step is parallel. Firstly, input data are prepared on the CPU. Then, the other devices are initialized and they simultaneously receive their respective part of input data. Complementary allocations and initialization are performed concurrently between all devices. At last, synchronization barriers are used to ensure that each device is ready to enter the computation loop.

Computation loop: We focused on maximizing the overlap of boundary(ies) computations, boundary(ies) transfers, and corpus computations on the different devices. Boundary transfers are performed simultaneously between devices that compute adjacent strips of the grid when boundaries are modified. On each device, boundaries computation and potential transfer are performed concurrently with the corpus computation in order to maximize the overlap of boundaries transfers with corpus computations. All the devices work concurrently thanks to asynchronous operations. Data transfers from/to the accelerators are fully efficient when two PCI express buses are present. Each iteration is ended by synchronization barriers to ensure that every device has all its newly updated data.

Final step: Results retrieving from devices to the CPU can be long and have been also overlapped. The results from one accelerator are transferred asynchronously during the synchronous transfer from the other one. Thus, transfers overlap and only one synchronization barrier is required to ensure the complete reception on the CPU. Finally, the CPU cleans up the devices and its own memory and releases them. As this last operation never appeared to be time consuming, it is done sequentially.

The design of optimized multi-core or many-core kernels is out of scope of this article as we focus on the efficient interactions between cooperating kernels. However, as we developed the computing kernels for the three types of devices by adapting the Jacobi kernels to the shortests path problem, the reader should see [START_REF] Contassot-Vivier | Algorithmic scheme for hybrid computing with CPU, Xeon-Phi/MIC and GPU devices on a single machine[END_REF] for further details. On GPU, both for Jacobi and Shortest Path applications, we designed CUDA kernels using the fast (but small) shared memory of each stream-processor, in order to load and access sub-parts of the data arrays. This is a kind of cache memory explicitly managed by the developer, which is a classical optimization on NVIDIA GPU. For the Shortest Path program, we designed a first solution that pre-compute distances between all neighbors in the grid, in order to avoid recomputing d(P i,j , P a,b ) (see Sec. 2). However, this version showed to be penalized by the induced larger data transfers between GPU memories. The adopted solution, less memory consuming, consists in taking advantage of the regular structure of the grid cells by pre-computing only X and Y axis contributions to the distances between neighbors, which represents at most 8 values instead of four times the number of grid points. This implies an important reduction of data copies into the shared memory at each iteration. In compensation, distances have to be recomputed at each iteration, but only in a partial form (thanks to X and Y precomputed contributions) that does not increase the computational load so much. This version has been used in all the experiments presented in Sec. 8.

Finally, we deployed 2D grids of 2D blocks of CUDA threads, and experimented different sizes of 2D blocks. The optimal solution depended on the GPU used: blocks of 16 × 16 (Y × X) threads run faster on the GeForce GTX Titan Black, while blocks of 16 × 32 threads run faster on the Tesla K40m (both GPU have a Kepler architecture).

On CPU and Xeon-Phi, we designed basic and cache-blocking algorithms. When blocked in cache, the grid update was processed per small sub-grids filling only the amount of L2 cache available per thread. Moreover, we attempted to guide the compiler to vectorize computational loops by using the AVX units (but we did not implement intrinsics routine calls). On the basic calculations of the Jacobi algorithm, our cache blocking mechanism has been efficient. On the opposite, on the Shortest Path algorithm our cache blocking implementation had no effect but vectorization has been efficient. Using #pragma ivdep and #pragma simd directives and AVX compiler options, and achieving minor code modifications we succeeded to significantly improve performance on Xeon CPU and on Xeon-Phi accelerator (3120 and 5100 series).

Optimized solution for Xeon-Phi and GPU accelerators

We present below the optimization of the initialization step and the computation loop.

Initial datastructure allocations and settings

Figure 3 introduces the actually implemented data-structures, not so different from the generic ones (see Fig. 2 right). Listing 1.1 illustrates how we implemented datastructures allocations and initialization by mixing GPU CUDA, CPU OpenMP and Xeon-Phi Offload semantics and syntax. This code implements and respects the initial step of our generic algorithm (see Sec. 4).

Main changes in the data-structures concern the Xeon-Phi (see Fig. 2 (right) and 3). In fact, we could observe during our developments that allocations of large arrays on the Xeon-Phi are very long when achieved from the CPU through a #pragma offload in(...) directive. However, this allocation mechanism is required to be able to transfer data between the CPU and the Xeon-Phi. So, we allocate the micPrev array and transfer CPU data with a #pragma offload directive (lines 4-5), but the _micCrt array is allocated as a pure Xeon-Phi variable (line 12), which is a fast mechanism. Then, it is initialized with a Xeon-Phi internal memcpy call (line 13). Obviously, the final results will have to be transferred from the Xeon-Phi to the CPU only through the micPrev array. This strategy significantly reduces the initialization step, but leads to allocate two additional small buffers to transfer the CPU and Xeon-Phi frontiers during each computation step, as shown in Fig. 3 (fm1 andfm2). This allocation of a transferable 2-lines array is done at line 6. All this sequence of operations is run asynchronously due to the signal clause in line 9. After launching the allocation and initialization sequence on the Xeon-Phi, the CPU immediately enters the GPU sequence at line 19. It allocates several streams to manage concurrent operations on the GPU, and several arrays on the GPU (lines 21 and 22). According to CUDA paradigm, the CPU allocates memory on the GPU calling a cudaMalloc function. At lines 24-26 the CPU runs another CUDA library function to lock in its memory the arrays that will be transferred to/from the GPU, in order to speedup the transfers and to do them asynchronously. All these memory allocations and locking operations are achieved synchronously but are fast. Then, the CPU launches long asynchronous data transfers to fill the allocated GPU arrays (lines 28-30), and an asynchronous GPU-to-GPU memory copy at lines 32-33. All these long operations are run asynchronously on stream 0 so that the CPU can concurrently perform its own initialization at lines 35-38 (running several threads to speedup a large memory copy).

Finally, we implemented the two synchronization barriers of our algorithm (Fig. 2 right). At line 40, the CPU waits for the end of all asynchronous GPU operations on stream 0, and at line 41, it waits for the end of all asynchronous Xeon-Phi operations. This implementation is very close to the initial step of our generic algorithm, save for some fast GPU memory allocation and CPU memory locking that remain synchronous. // GPU 41 #pragma offload_wait target(mic:0) wait(&micPrev) // Xeon-Phi

Asynchronous and overlapped computation loop

Code snippet in Listing 1.2 summarizes our implementation of the computation loop of the generic algorithm (Fig. 2 right). Line 1 shows that the global loop is stopped as soon as no modification occurs during one iteration. Line 4 is an asynchronous CUDA kernel launch, running frontier computation on the GPU. Indeed, a CUDA kernel launch has a specific CUDA syntax, and requires to be implemented in a .cu source file, and to be compiled with the Nvidia CUDA compiler (nvcc), that cannot compile Intel offload or OpenMP directives. So, CUDA kernels launching operations have to be encapsulated into functions located in a .cu file. Line 5 transfers a modification flag of the GPU frontier onto the CPU. This asynchronous transfer uses the same CUDA top frontier stream (TF) as the previous kernel execution, so these two asynchronous operations are sequenced: frontier computation is achieved and the modification flag set before it is transferred to the CPU. Lines 7 and 8 reproduce this kernel launch and modification flag transfer, focusing on the rest of the problem processed on the GPU, using the CUDA corpus stream.

Lines 11-16 perform a similar frontier and corpus computations on the Xeon-Phi (like on the GPU). Here, the modification flags are scalar boolean values automatically transferred between the CPU and the Xeon-Phi at the beginning and at the end of each offload directive. Moreover, when declared (at top of the source file) the heterMIC routine has been labeled to be a Xeon-Phi routine, in order the Intel compiler applies all the suited optimizations. In particular, it can take into account the vectorization directives embedded in the source code of this routine.

Lines 18-25 are the synchronous CPU computation of the two CPU frontiers (CPU/ GPU and CPU/Xeon-Phi), using OpenMP multithreading and Intel compiler vectorization on SSE or AVX units. The CellUpdate routine (lines 22-23) is compiled only for CPU, and Intel compiler can achieve vectorization optimization adapted to CPU cores. Lines 27-35 perform the asynchronous transfers of the CPU frontiers to the GPU, using a new CUDA bottom frontier stream (BF), and to the Xeon-Phi, using a new signal clause (dummyTransfer). According to our generic algorithm (see Sec. 4) these transfers are done only if the frontiers have been modified on the CPU.

At lines 37-40, the CPU waits for the end of the GPU frontier computation and modification flag transfer (TF CUDA stream) before to asynchronously transfer the GPU frontier to the CPU if it has been modified on the GPU. Lines 42-46 execute the same operations for the Xeon-Phi. Lines 48-61 correspond to the computation of the CPU part of the problem excepted its frontiers (the CPU corpus). Again, it is a multithreaded and vectorized computation, and a modification flag is raised when a modification occurs (line 60).

Lines 63-65 contain synchronization barriers for the end of the CPU/GPU frontier transfers and the end of the GPU corpus computation. Lines 67-71 perform the same operations on the Xeon-Phi. Lines 73-77 permute current and previous array pointers on the Xeon-Phi, the GPU, and the CPU. These permutations are not complex but their implementations differ: pointers on GPU memory are CPU variables managed by the CPU, while pointers on Xeon-Phi memory are stored on the Xeon-Phi and must be permuted by the Xeon-Phi. Finally, pointers on the Xeon-Phi memory have mirror pointers on the CPU (due to the offload semantics) that must be permuted too (line 77). Lines 78-81 compute global modification flags corresponding to each computing device. ----------------------------------- Finaly, mixing the Nvidia CUDA and Intel offload programming paradigm leads to a complex syntax, but we succeeded to implement our generic asynchronous parallel algorithm, including many computations / data transfers overlapping. Vectorization can be used on the CPU and Xeon-Phi kernels, but must be implemented in separated routines so that the compiler can easily generate efficient code for each architecture. Some different CUDA streams must be associated with the GPU kernel launches and data transfers, in order to concurrently run sequences of operations.

The final step of our generic algorithm mainly consists in retrieving the results from the accelerators. The mechanisms used to implement this are similar to the initialization step and consists in overlapping the results transfers from GPU and Xeon-Phi by using an asynchronous transfer from the GPU. So, we do not give code sample for that part. Instead, we mention that on the Xeon-Phi, it is necessary to copy the results into the transferable array (the one initialized with the in offload directive, line 5 in Listing 1.1) when the number of iterations is odd.

Testbeds and benchmark data for shortest path computation

The machine used at CentraleSupelec is a Dell R720 server containing two 6-cores Intel(R) Xeon(R) CPU E5-2620 at 2.10GHz, with two accelerators on separate PCIe buses. One accelerator is an Intel MIC Xeon-Phi 3120 with 57 physical cores at 1.10 GHz, supporting 4 threads each. The second accelerator is a Nvidia GPU GeForce GTX Titan Black (Kepler architecture) with 2880 CUDA cores. The machine used at Loria is a Dell R720 server containing two 8-cores Intel(R) Xeon(R) CPU E5-2640 at 2.00GHz, with two accelerators on separate PCIe buses. One accelerator is an Intel MIC Xeon-Phi 5100 with 60 physical cores at 1.05 GHz supporting 4 threads each. The second accelerator is a Nvidia GPU Tesla K40m (Kepler architecture) with 2880 CUDA cores.

In order to get representative results with the shortest path problem, we used elevation data of a large area in the French and Italian Alps from the NASA SRTM project [START_REF]Shuttle Radar Topography Mission[END_REF]. The 2D grid resolution is 30m, with 10803 columns and 18005 lines (≈324Km×540Km), and the elevations range from -87m to 4797m. The overall calculation takes 9005 iterations with a target point located at the center of the area. In this study, we did not use larger datasets to be able to compare executions on single devices.

Benchmark results with different hybrid machines

In the following paragraphs, we present the performance results of our algorithm in several steps. Firstly, we present the performance of each device alone in order to give an idea of the relative computing powers of the devices inside each machine. These relative powers are useful to deduce the theoretical optimal partitioning of the 2D grid by using a simple static load balancing scheme. Then, we consider double device mixing, where two devices cooperate (GPU-CPU, CPU-PHI, GPU-PHI). Finally, the results with the three cooperating devices are presented. All the following experiments are averages of several executions although we observed a very limited standard deviation. Standalone performance In Tab. 1 are presented the performance results for each device in both machines. Speeds are expressed in points per second instead of GFlops as the computation amount varies among the points. Relative computing powers are estimated according to the total power of each machine, deduced from the aggregation of the devices speeds. As expected, the CPU is the slowest device inside each machine and the GPU is the fastest one. However, strong differences between the two machines are observed, which is interesting as it provides significantly different hardware configurations. It allows us to see how our algorithm behaves in either cases. Absolute speeds are useful here to deduce the relative computing powers of the devices but they are omitted in the following tables as they are redundant information with gains.

Double device mixing

According to the single devices results, the relative computing powers of the devices in each machine are used to perform static load balancing. The optimal grid partitioning into two horizontal strips is computed for every couple of devices. For example, GPU is 4.7 times faster than CPU in the Loria machine, then the strip associated to the GPU will be 4.7 times larger than the CPU one, leading to a cutting line (opt. cut in Tab. 2) of 14900. The results are presented in Tab. 2 where the gains are computed relatively to the speed of the fastest device involved in the mixing.

First of all, we can see that some gains are negative and others are positive. Negative gains indicate that the mixing is less efficient than the fastest device alone. With both machines, the GPU-CPU mixing does not provide any gain. Indeed, even extremely imbalanced partitions with most of the grid on the GPU (17500 lines) provided a loss of 19.42% on the Loria machine and 8.20% on the CentralSupelec one. Those results show that the mixing overheads are never compensated by the computation gains with that grid size. Concerning the CPU-PHI mixing, we obtain significant gains on both machines, showing that our scheme is efficient for mixing those two devices. Finally, we obtain diverging results between the two machines for the GPU-PHI mixing. This mainly comes from the difference of powers between the two GPU devices. As the powers of the PHI and GPU are similar on the CS machine, their mixing is actually efficient whereas on the Loria machine the two devices have too different powers to compensate the mixing overheads. Also, it is surprising to get so different gains between GPU-CPU and CPU-PHI on the CS machine and this will require a deeper analyze. Triple device mixing For the sake of concision, we present in Tab. 3 only the percentage gains for both machine around their respective optimal cutting lines between GPU-CPU and CPU-PHI. For the CS machine, we added an extra column of GPU-CPU cutting line as the maximal performance was obtained at the limit of the initial set. It can be observed that our heterogeneous algorithm hardly obtains positive gains on the Loria machine whereas it obtains up to 37% gain over the GPU on the CS machine (better than double mixing). The bad performance with the Loria machine is quite unexpected. It probably comes from the different hardware configuration as well as the older Intel compiler used. However, results with the CS machine show that our algorithmic scheme can provide significant gains. All those results show that the efficiency of our approach is sensitive to the hardware configuration as well as to the problem size. This will deserve a more complete investigation over the behavior of our scheme according to larger grid sizes.

CPU

Conclusion

An algorithmic scheme has been presented that enables the cooperation of several computing devices of different types (CPU, GPU, PHI) to process a stencil application on a single machine. Our scheme makes an intensive use of asynchronous computations and data transfers in order to obtain an efficient overlapping of both operations.

The scheme has been evaluated with different combinations of devices cooperation on a representative elevation map of the Alps region. Results show that our algorithm can provide significant gains either with two or three devices. It has been observed that better gains are obtained when the relative computing powers of the devices are closer.

Several extensions should be interesting such as studying the behavior of our scheme with larger grid sizes, including a dynamic load balancing by re-partitioning the grid during the algorithm execution, as well as extending our scheme to cluster systems by adding explicit inter-machine communications.

Fig. 2 .

 2 Fig. 2. Data structures (left) and asynchronous and overlapping algorithm (right) to obtain efficient simultaneous computations on the three devices (CPU, GPU and Xeon-phi)

Fig. 3 .

 3 Fig. 3. Optimized data structures implemented on the three devices (CPU, Xphi and GPU)

Listing 1 . 1 . 6 in 11 /

 11611 Initial step implementation 1 // CPU allocation of data array 2 posix_memalign(&crt, ALIGN, TotCol * TotLin * sizeof(double)); 3 // ASYNCHRONOUS XEON-PHI ALLOCATIONS AND INITIALIZATION 4 #pragma offload target(mic:0) \ 5 in(micPrev:length(micTotLin * TotCol) free_if(0) align(ALIGN)) \ / Alloc and init of _micCrt array on Xeon-Phi (matching micPrev array) 12 posix_memalign(&_micCrt, ALIGN, micTotLin * TotCol * sizeof(double)); 13 memcpy(_micCrt, micPrev, micTotLin * TotCol * sizeof(double)); 14 // Init pointers on data arrays on Xeon-Phi 15 micRes = micPrev + TotCol; 16 ..

Listing 1 . 2 .

 12 Computation loop implementation 1 for(iter=0; iter<nbIters && (modifCPU || modifGPU || modifMIC); ++iter) { 2 ...

3 / 11 #

 311 / ASYNCHRONOUS GPU computations and frontier transfer to the CPU ------------4 gpuOneIterationLast(gpuPrec+..., gpuCrt+..., 1 / * 1 line * /, nbCol, ..., TF); 5 cudaMemcpyFromSymbolAsync(&modifGPULast, pt_gpu_modif_hf, ..., 6 cudaMemcpyDeviceToHost, streamTab[TF]); 7 gpuOneIteration(gpuPrev, gpuCrt, ..., nbLinsGPU-1, nbCol, ..., CORPUS); 8 cudaMemcpyFromSymbolAsync(&modifGPUCorpus, pt_gpu_modif_corpus, ..., 9 cudaMemcpyDeviceToHost, streamTab[CORPUS]); 10 // ASYNCHRONOUS Xeon-Phi computations and frontier transfer to the CPU ------pragma offload target(mic:0) nocopy(botFrontMIC, micPrev, _micCrt, ...) \ 12 signal(&botFrontMIC) 13 { heterMIC(micPrev, _micCrt, ..., botFrontMIC, &modifMICFirst); } 14 #pragma offload target(mic:0) \ 15 nocopy(micPrev, _micCrt, micTer, ..., modifMICCorpus) signal(&_micCrt) 16 { heterMIC(micPrev, _micCrt, ..., &modifMICCorpus); } 17 // Synchronous CPU frontier computations -

- 18 #pragma omp parallel num_threads(nbTHet) 19 { #pragma omp for 20 #pragma simd 21 for

 18192021 (col=1; col <= nbCol; ++col) { 22 CellUpdate(prev, crt, ..., firstLine + col, &modifCPUFirst); 23 CellUpdate(prev, crt, ..., lastLine + col, &modifCPULast); CPU frontiers transfers to accelerators ----------------------27 if (modifCPUFirst) 28 cudaMemcpyAsync(gpuCrt+..., cpuCrt+..., nbCol * sizeof(double), 29 cudaMemcpyHostToDevice, streamTab[BF]); 30 if (modifCPULast) { #pragma offload target(mic:0) \ 32 in(topFrontMIC:length(nbCol+1) alloc_if(0) free_if(0) align(ALIGN)) \ 33 in(dummyTransfert) nocopy(_micCrt) signal(&dummyTransfer) 34 { memcpy(_micCrt+1, topFrontMIC+1, nbCol * sizeof(double)); } 35 } 36 // ASYNCHRONOUS GPU frontier transfer to the CPU -----------------------------37 cudaStreamSynchronize(streamTab[TF]); 38 if (modifGPULast) 39 cudaMemcpyAsync(topFrontCPU, gpuCrt+..., nbCol * sizeof(double), 40 cudaMemcpyDeviceToHost, streamTab[TF]);

  .

	17 }
	18 // Synchronous GPU Allocations
	19 for (i=0; i<NBSTREAMS; ++i)	// 3 streams
	20	cudaStreamCreate(&(streamTab[i]));
	21 cudaMalloc(&gpuPrev, gpuTotLin * TotCol * sizeof(double));
	22 ...
	23 // Synchronous CPU memory lock (to speedup CPU-GPU data transfers)
	24 cudaHostRegister(prev, gpuTotLin * TotCol * sizeof(double), // transferred
	25	cudaHostRegisterPortable);	// to gpuPrev
	26 ...
	27 // ASYNCHRONOUS GPU INIT: data transfer from CPU to GPU
	28 cudaMemcpyAsync(GpuPrev, prev, gpuTotLin * TotCol * sizeof(double),
	29	cudaMemcpyHostToDevice, streamTab[0]);
	30 ...
	31 // ASYNCHRONOUS GPU INIT: data copy from GPU to GPU
	32 cudaMemcpyAsync(GpuCrt, GpuPrev, gpuTotLin * TotCol * sizeof(double),
	33	cudaMemcpyDeviceToDevice, streamTab[0]);
	34 // Complementary CPU init. (parallelized with OpenMP)
	35 #pragma omp parallel num_threads(nbTHet) {

36 size_t idxFirstValTh = ...; size_t nbValTh = ...; 37 memcpy(crt+idxFirstValTh, prev+idxFirstValTh, nbValTh * sizeof(double)); 38 } 39 // SYNCHRONIZATION BARRIERS on the end of the GPU and Xeon-Phi async. init. 40 cudaStreamSynchronize(streamTab[0]);

Table 1 .

 1 Performance results of each device alone and relative computing powers

		Loria machine		CentralSupelec machine
	CPU	PHI	GPU	CPU	PHI	GPU
	Speed (pts/s) 6.35E+008 1.56E+009 3.00E+009 4.46E+008 1.74E+009 1.74E+009
	Relative powers 12.20%	29.97%	57.83%	11.36%	44.32%	44.32%

Table 2 .

 2 Performance results and gains around theoretical optimal cutting lines for both machines.

  68% -23.83% -20.22% -12.83% -7.13% 11000 -27.32% -20.27% -11.68% -1.08% 19.08% 32.01%

			Loria machine			CPU	CentraleSupelec machine
	PHI		GPU-CPU cut line		PHI		GPU-CPU cut line
	cut line 9400	9900	10400	10900	11400 cut line 7000	7500	8000	8500	9000	9500
	11600	1.28% -3.38% -7.07% -11.66% -13.60%	9000 24.87% 28.94% 27.08% 25.63% 23.47%
	12100	-2.13%	1.55% -2.31% -5.70% -9.57%	9500 3.05% 23.26% 32.60% 30.65% 28.41% 22.53%
	12600 -12.73% -5.19%	3.00% -1.03% -5.78% 10000 -9.33%	1.61% 21.92% 34.83% 32.30% 26.24%
	13100 -21.84% -16.00% -9.32% -4.51% -2.01% 10500 -19.32% -10.69%	0.07% 20.51% 37.01% 30.83%
	13600 -30.

Table 3 .

 3 Performance gains for the triple mixing on both machines.

// ASYNCHRONOUS Xeon-Phi frontier transfer to the CPU - ----------------------