-0.5cm

# Maximal line-free sets in $\mathbb{F}_{p}^{n}$

Christian Elsholtz  Christian Elsholtz and Jakob Führer, Institute of Analysis
and Number Theory, Graz University of Technology, Kopernikusgasse 24/II, 8010
Graz, Austria<EMAIL_ADDRESS>,  Jakob Führer
<EMAIL_ADDRESS>,  Erik Füredi  Erik Füredi, ELTE Eötvös Loránd
University Faculty of Science, 1117 Budapest, Pázmány Péter sétány 1/A,
Hungary<EMAIL_ADDRESS>,  Benedek Kovács  Benedek Kovács, ELTE Linear
Hypergraphs Research Group, Eötvös Loránd University, 1117 Budapest, Pázmány
Péter sétány 1/A, Hungary<EMAIL_ADDRESS>,  Péter Pál Pach  Péter
Pál Pach, Department of Computer Science and Information Theory,
Budapest University of Technology and Economics,
Műegyetem rkp. 3., H-1111 Budapest, Hungary;
MTA-BME Lendület Arithmetic Combinatorics Research Group, ELKH,
Műegyetem rkp. 3., H-1111 Budapest, Hungary;
Extremal Combinatorics and Probability Group (ECOPRO), Institute for Basic
Science (IBS), South Korea<EMAIL_ADDRESS>,  Dániel Gábor Simon  Dániel Gábor
Simon, Alfréd Rényi Institute of Mathematics, Reáltanoda street 13-15, H-1053
Budapest, Hungary<EMAIL_ADDRESS>and  Nóra Velich  Nóra Velich,
University of Cambridge, Lucy Cavendish College, Lady Margaret Road, Cambridge
CB3 0BU, UK<EMAIL_ADDRESS>
###### Abstract.

We study subsets of $\mathbb{F}_{p}^{n}$ that do not contain progressions of
length $k$. We denote by $r_{k}(\mathbb{F}_{p}^{n})$ the cardinality of such
subsets containing a maximal number of elements.

In this paper we focus on the case $k=p$ and therefore sets containing no full
line. A trivial lower bound $r_{p}(\mathbb{F}_{p}^{n})\geq(p-1)^{n}$ is
achieved by a hypercube of side length $p-1$ and it is known that equality
holds for $n\in\\{1,2\\}$. We will however show
$r_{p}(\mathbb{F}_{p}^{n})\geq(p-1)^{3}+p-2\sqrt{p}$, which is the first
improvement in the three dimensional case that is increasing in $p$.

We will also give the upper bound $r_{p}(\mathbb{F}_{p}^{3})\leq
p^{3}-2p^{2}-(\sqrt{2}-1)p+2$ as well as generalizations for higher
dimensions.

Finally we present some bounds for individual $p$ and $n$, in particular
$r_{5}(\mathbb{F}_{5}^{3})\geq 70$ and $r_{7}(\mathbb{F}_{7}^{3})\geq 225$
which can be used to give the asymptotic lower bound $4.121^{n}$ for
$r_{5}(\mathbb{F}_{5}^{n})$ and $6.082^{n}$ for $r_{7}(\mathbb{F}_{7}^{n})$.

###### Key words and phrases:

line-free sets, blocking sets, sets without arithmetic progressions, extremal
combinatorics, combinatorics in $\mathbb{F}_{p}^{n}$

###### 1991 Mathematics Subject Classification:

51E21, 11B25, 05D05

##  1\. Introduction

In the intersection of finite geometry and extremal combinatorics numerous
problems of finding maximal subsets of affine or projective spaces avoiding
certain configurations have been studied. One natural question asks for bounds
on the cardinality of subsets of the $n$-dimensional affine space over a
finite field $\mathbb{F}_{q}$ that do not contain a full line.

We denote by $r_{k}(\mathbb{F}_{p}^{n})$ the cardinality of a subset
$S\subseteq\mathbb{F}_{p}^{n}$, containing a maximal number of elements such
that $S$ contains no $k$ points in arithmetic progression. Note that in the
case when $k=p$ is a prime, $k$-progressions in $\mathbb{F}_{p}^{n}$
correspond to lines in the $n$-dimensional affine space and we are therefore
interested in bounds on $r_{p}(\mathbb{F}_{p}^{n})$.

When $p=3$ the problem coincides with the cap set problem, a well studied area
where one can use the fact that $x,y,z$ form a line exactly when they fulfil a
non-trivial linear equation $ax+by+cz=0$ where $a+b+c=0$. Ellenberg and
Gijswijt [9] gave the first exponential improvement to the trivial upper bound
of $3^{n}$ with $r_{3}(\mathbb{F}_{3}^{n})<2.756^{n}$, for large enough $n$
which was further improved by Jiang [16] by a factor of $\sqrt{n}$. The best
lower bound was given by Tyrrell [23] with
$r_{3}(\mathbb{F}_{3}^{n})>2.218^{n}$, for large enough $n$. The exact values
of $r_{3}(\mathbb{F}_{3}^{n})$ are known up to $n=6$, where
$r_{3}(\mathbb{F}_{3}^{6})=112$ was proven by Potechin [22].

For the general case surprisingly few results on $r_{p}({\mathbb{F}}_{p}^{n})$
are known. There is the trivial lower bound
$r_{p}(\mathbb{F}_{p}^{n})\geq(p-1)^{n}$ achieved by a hypercube of side
length $p-1$. Jamison [15] and Brouwer and Schrijver [5] independently proved
that this is sharp for $n=2$. For $n=3$ the only improvement to this
construction was by a single point described in the post of Zare in a
mathoverflow thread [1]. We will prove the following lower bounds:

######  Theorem 1.1.

Let $p\geq 5$ be a prime then

$r_{p}(\mathbb{F}_{p}^{3})\geq(p-1)^{3}+p-2\sqrt{p}=p^{3}-3p^{2}+4p-2\sqrt{p}-1.$

This can be improved in some special cases.

######  Theorem 1.2.

Let $p$ be a prime with $p\equiv 7$ (mod $24$), then

$r_{p}(\mathbb{F}_{p}^{3})\geq(p-1)^{3}+(p-1)=p^{3}-3p^{2}+4p-2.$

Moreover, $r_{7}(\mathbb{F}_{7}^{3})\geq 225$.

The simple upper bound $r_{p}(\mathbb{F}_{p}^{n})\leq
p^{n}-\frac{p^{n}-1}{p-1}$ was given by Aleksanyan and Papikian [2] and is
achieved by removing at least one point from each line going through a fixed
point. In particular $r_{p}(\mathbb{F}_{p}^{3})\leq p^{3}-p^{2}-p-1$. This was
improved by Bishnoi et al. [4] to $p^{n}-2p^{n-1}+1$ and $p^{3}-2p^{2}+1$,
respectively. We will give the following new bounds:

######  Theorem 1.3.

Let $p\geq 3$ be a prime, $k\in\\{3,4,...,p\\}$ and $n\in\mathbb{N}$ then

${r_{k}(\mathbb{F}_{p}^{n+1})}\leq\frac{2(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})+p^{n}-\sqrt{4(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})(p^{n}-r_{k}(\mathbb{F}_{p}^{n}))+p^{2n}}}{2p^{n}},$

where the three-dimensional case gives the following corollary.

######  Corollary 1.4.

Let $p\geq 3$ be a prime then

$r_{p}(\mathbb{F}_{p}^{3})\leq\frac{2p^{5}-4p^{4}+2p^{3}-p^{2}+4p-2-\sqrt{8p^{6}-20p^{5}+17p^{4}-12p^{3}+20p^{2}-16p+4}}{2p^{2}},$

in particular,

$r_{p}(\mathbb{F}_{p}^{3})\leq p^{3}-2p^{2}-(\sqrt{2}-1)p+2.$

For other dimensions, there is the lower bound $r_{p}(\mathbb{F}_{p}^{2p})\geq
p(p-1)^{2p-1}$ due to Frankl et al. [13], using large sunflower-free sets.

We found a $70$ point $5$-progression-free set in $\mathbb{F}_{5}^{3}$ via a
branch and cut approach (see Figure 4) and we will show the following upper
bounds for small primes.

######  Theorem 1.5.

$r_{5}(\mathbb{F}_{5}^{3})<74.$

######  Theorem 1.6.

$r_{7}(\mathbb{F}_{7}^{3})<243.$

One can use the tensor product $S_{1}\times S_{2}$ of two line-free sets
$S_{1}\in\mathbb{F}_{p}^{n_{1}}$, $S_{2}\in\mathbb{F}_{p}^{n_{2}}$ to get a
line-free set in the higher dimension $n_{1}+n_{2}$. This construction also
provides us the lower bound $|S_{1}|^{1/{n_{1}}}$ for
$\alpha_{p}:=\lim\limits_{n\rightarrow\infty}(r_{p}(\mathbb{F}_{p}^{n}))^{1/n}$
and therefore the asymptotic lower bound $(|S_{1}|^{1/{n_{1}}}-o(1))^{n}$ for
$r_{p}(\mathbb{F}_{p}^{n})$ (see e.g. [7], [18]). The strongest known lower
bound for general $p$ is $\alpha_{p}\geq p^{1/{2p}}(p-1)^{(2p-1)/{2p}}$ using
the results of Frankl et al. [13], however for small primes the new three-
dimensional lower bounds $r_{5}(\mathbb{F}_{5}^{3})\geq 70$ and
$r_{7}(\mathbb{F}_{7}^{3})\geq 225$ give better lower bounds, namely,
$\alpha_{5}\geq 4.121$ and $\alpha_{7}\geq 6.082$.

We will also show the following explicit lower bound for arbitrary dimension.

######  Theorem 1.7.

Let $p\geq 3$ be a prime then
$r_{p}(\mathbb{F}_{p}^{n})\geq(p-1)^{n}+\frac{n-2}{2}(p-1)(p-2)^{n-3}.$

##  2\. Related results

  * •

Davis and Maclagan [7] studied the card game SET, where the cards can be
described as points in $\mathbb{F}_{3}^{4}$ and one is interested whether the
displayed cards form a cap set. The results of Tyrell [23] are building on the
construction of Edel [8], who gave the previously best lower bound for cap
sets. Elsholtz and Lipnik [11] and Elsholtz and Pach [12] studied cap sets in
other spaces than $\mathbb{F}_{3}$.

  * •

Croot et al. [6] gave an upper bound for $3$-progression-free sets in
$\mathbb{Z}_{4}^{n}$ that is exponentially smaller than $4^{n}$. Their methods
also led to the results of Ellenberg and Gijswijt [9]. Petrov and Pohoata [20]
gave an improved upper bound for $3$-progression-free sets in
$\mathbb{Z}_{8}^{n}$, Pach and Palincza [19] gave both upper and lower bounds
for $6$-progression-free sets in $\mathbb{Z}_{6}^{n}$. Elsholtz et al. [10]
studied the general case of $k$-progression-free sets in $\mathbb{Z}_{m}^{n}$.
An overview on known bounds is given by Pach [18].

  * •

Moser [17] asked for the maximal size of a subset of $\\{1,2,...,k\\}^{n}$
without a geometric line. Similarly, Hales and Jewett asked for a subset
without a combinatorial line. The result of Furstenberg and Katznelson [14]
also known as the density Hales–Jewett theorem implies that in both cases
these sets have to be asymptotically smaller than $k^{n}$ as $n$ tends to
infinity. Polymath [21] gave some explicit bounds for special cases.

  * •

Sets that intersect every affine subspace of codimension $s$ are called
$s$-blocking sets. The complement of a line-free set in a finite
$n$-dimensional affine space is therefore also called an $(n-1)$-blocking set.
It is known that the union of any $n$ independent lines intersecting in a
single point form a $1$-blocking set in $\mathbb{F}_{p}^{n}$ which is optimal
(see e.g. [3], [5], [15]). However, for $(n-1)$-blocking sets, the union of
$n$ independent hyperplanes, which seems to be the obvious algebraic
construction, are not optimal, as will be shown in this paper. Bishnoi et al.
[4] gave several upper bounds for the size of $s$-blocking sets.

##  3\. Notation

We write $\mathbb{Z}_{n}$ for $\mathbb{Z}/n\mathbb{Z}$ and
$\mathbb{F}_{p}=\mathbb{Z}_{p}$ is the field with $p$ elements whenever $p$ is
a prime.

We write $[k,\ell]$ for the set $\\{k,k+1,...,\ell\\}$ either as a subset of
$\mathbb{Z}$ or of $\mathbb{F}_{p}$.

We use both row and column vectors for the elements of $\mathbb{F}_{p}^{n}$
and we call these elements points.

Given a subset $S\subseteq\mathbb{F}_{p}^{3}$ we call the image of
$S\cap(\\{j\\}\times\mathbb{F}_{p}^{2})$ under the projection
$\phi\colon\mathbb{F}_{p}^{3}\longrightarrow\mathbb{F}_{p}^{2}$,
$(a,b,c)\mapsto\left(b,c\right)$ the $j$-layer of $S$.

##  4\. Proofs of the upper bounds

###### Proof of Theorem 1.3.

Let $A\subseteq\mathbb{F}_{p}^{n+1}$ be $k$-progression-free. We count the
number of the point pairs on every $n$-dimensional hyperplane

$s=\left|\\{(\\{a,b\\},S)\ |\ a,b\in A,\ a\neq b,\ a,b\in S,\ \text{$S$ is an
$n$-dimensional hyperplane}\\}\right|.$

On every hyperplane, the number of points is at most
$r_{k}(\mathbb{F}_{p}^{n})$. Firstly, we assume
$r_{k}(\mathbb{F}_{p}^{n+1})\geq(p-1)r_{k}(\mathbb{F}_{p}^{n})$, then the sum
of number of point pairs for $p$ parallel hyperplanes is maximal if there are
$p-1$ hyperplanes with $r_{k}(\mathbb{F}_{p}^{n})$ points and one with
$r_{k}(\mathbb{F}_{p}^{n+1})-(p-1)r_{k}(\mathbb{F}_{p}^{n})$.

There are $\frac{p^{n+1}-1}{p-1}$ disjoint sets of parallel hyperplanes, so

$s\leq\Big{(}\frac{p^{n+1}-1}{p-1}\Big{)}\Big{(}(p-1)\binom{r_{k}(\mathbb{F}_{p}^{n})}{2}+\binom{r_{k}(\mathbb{F}_{p}^{n+1})-(p-1)r_{k}(\mathbb{F}_{p}^{n})}{2}\Big{)}.$

Note here that this inequality still holds if
$r_{k}(\mathbb{F}_{p}^{n+1})<(p-1)r_{k}(\mathbb{F}_{p}^{n})$, as in this case
the number of point pairs is clearly less than

$(p-1)\binom{r_{k}(\mathbb{F}_{p}^{n})}{2}$

and

$\binom{r_{k}(\mathbb{F}_{p}^{n+1})-(p-1)r_{k}(\mathbb{F}_{p}^{n})}{2}\geq 0.$

On the other hand, every point pair defines a line that is included in exactly
$\frac{p^{n}-1}{p-1}$ $n$-dimensional hyperplanes, so

$s=\frac{p^{n}-1}{p-1}\binom{r_{k}(\mathbb{F}_{p}^{n+1})}{2}.$

We get the quadratic inequality

$p^{n}\big{(}r_{k}(\mathbb{F}_{p}^{n+1})\big{)}^{2}-\big{(}p^{n}+2(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})\big{)}r_{k}(\mathbb{F}_{p}^{n+1})+(p^{n+2}-p)\big{(}r_{k}(\mathbb{F}_{p}^{n})\big{)}^{2}\geq
0$

with roots

$\frac{2(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})+p^{n}\pm\sqrt{4(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})(p^{n}-r_{k}(\mathbb{F}_{p}^{n}))+p^{2n}}}{2p^{n}}.$

As

${r_{k}(\mathbb{F}_{p}^{n+1})}\leq p(r_{k}(\mathbb{F}_{p}^{n}))$

but

$\begin{split}&\frac{2(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})+p^{n}+\sqrt{4(p^{n+1}-1)r_{k}(\mathbb{F}_{p}^{n})(p^{n}-r_{k}(\mathbb{F}_{p}^{n}))+p^{2n}}}{2p^{n}}\\\
>\;&p(r_{k}(\mathbb{F}_{p}^{n}))+\frac{1}{2}-\frac{r_{k}(\mathbb{F}_{p}^{n})}{p^{n}}+\frac{\sqrt{p^{2n}}}{2p^{n}}\geq
p(r_{k}(\mathbb{F}_{p}^{n}))+\frac{1}{2}-1+\frac{1}{2}=p(r_{k}(\mathbb{F}_{p}^{n})),\end{split}$

the theorem follows. ∎

###### Proof of Corollary 1.4.

The first statement follows immediately from Theorem 1.3 using
$r_{p}(\mathbb{F}_{p}^{2})=(p-1)^{2}$. For the second statement we are using
that $8p^{6}-20p^{5}+17p^{4}-12p^{3}+20p^{2}-16p+4$ can be bounded by
$(2\sqrt{2}p^{3}-5/\sqrt{2}p^{2})^{2}$ from below for $p\geq 3$ and we get

$\begin{split}r_{p}(\mathbb{F}_{p}^{3})\leq&\;p^{3}-2p^{2}+p-\frac{1}{2}+\frac{2}{p}-\frac{1}{p^{2}}-\sqrt{2}p+\frac{5}{2\sqrt{2}}\\\
\leq&\;p^{3}-2p^{2}-(\sqrt{2}-1)p-\frac{1}{2}+\frac{2}{3}+\frac{5}{2\sqrt{2}}\\\
\leq&\;p^{3}-2p^{2}-(\sqrt{2}-1)p+2.\end{split}$

∎

###### Proof of Theorem 1.5.

Assume that $S\subseteq\mathbb{F}_{5}^{3}$ is a $5$-progression-free set of
size $74$. We will compute a weighted sum over all lines containing $4$ points
to reach a contradiction.

Let us call a line containing exactly $r$ points an $r$-line. Let $l$ be a
$4$-line in $S$ and let $H_{1},H_{2},...,H_{6}$ be the planes containing $l$.
Then $\sum_{i=1}^{6}|H_{i}\cap S|=(74-4)+6\cdot 4=94$. Note that
$r_{5}(\mathbb{F}_{5}^{2})=16$ and $r_{4}(\mathbb{F}_{5}^{2})=11$, which can
be easily checked by computer search. Therefore, $|H_{i}\cap A|\geq 94-5\cdot
16=14$ for all $i$ and there is no plane in $\mathbb{F}_{5}^{3}$ containing
$12$ or $13$ points. Hence, there are five different distributions for the
number of points in five parallel planes:

  * (a)

$\\{10,16,16,16,16\\}$

  * (b)

$\\{11,15,16,16,16\\}$

  * (c)

$\\{14,14,14,16,16\\}$

  * (d)

$\\{14,14,15,15,16\\}$

  * (e)

$\\{14,15,15,15,15\\}.$

Denote by $a,b,c,d,e$ the number of classes of parallel planes having these
distributions. Note that

(4.1) $a+b+c+d+e=31.$

If we compare the number of pairs of points in each plane with the total
number of pairs we get
$(\binom{10}{2}+4\binom{16}{2})a+(\binom{11}{2}+\binom{15}{2}+3\binom{16}{2})b+(3\binom{14}{2}+2\binom{16}{2})c+(2\binom{14}{2}+2\binom{15}{2}+\binom{16}{2})d+(\binom{14}{2}+4\binom{15}{2})e=6\binom{74}{2}$

(4.2) $\Leftrightarrow 525a+520b+513c+512d+511e=16206,$

since each pair lies in exactly six planes.

Now denote by $A$, $B$, and $C$ the number of pairs $(\ell,H)$ where $H$ is a
hyperplane containing $16$, $15$ and $14$ points, respectively and
$\ell\subseteq H$ is a $4$-line. Again let $\ell$ be a $4$-line and let
$H_{1},H_{2},...,H_{6}$ be the planes containing $\ell$. Then

$\\{\left|H_{i}\cap
S\right||\;i\in[1,6]\\}\in\\{\\{14,16,16,16,16,16\\},\\{15,15,16,16,16,16\\}\\}$

as multisets and therefore

(4.3) $A-2B-5C=0$

To bound the size of $A$, $B$ and $C$ we need the following claims.

######  Claim 1.

Every plane containing $16$ points contains at least twelve $4$-lines.

###### Proof of Claim 1.

Consider a plane $H$ containing $16$ points and let $x_{i}$ be the number of
$i$-lines in $H$ for $i\in\\{1,2,3,4\\}$. By double counting the points in $H$
we get $x_{1}+2x_{2}+3x_{3}+4x_{4}=6\cdot 16=96$ and by double counting the
pairs of points in $H$ we get $x_{2}+3x_{3}+6x_{4}=\binom{16}{2}=120$. By
taking the difference of the two equations we get $-x_{1}-x_{2}+2x_{4}=24$,
implying that $2x_{4}\geq 24$. ∎

######  Claim 2.

For $m\in\\{14,15\\}$, every plane containing $m$ points contains at most $m$
$4$-lines.

###### Proof of Claim 2.

As $5\cdot 3+1=16>m$, every point in $S$ can be contained in at most four
$4$-lines and therefore the number of $4$-lines in the plane is bounded from
above by $\frac{4m}{4}=m.$ ∎

Finally combining (4.1), (4.2) and (4.3) we obtain the following system of
linear equations and inequalities.

$\begin{split}&a+b+c+d+e=31\\\ &525a+520b+513c+512d+511e=16206\\\
&A-2B-5C=0\\\ &A\geq 48a+36b+24c+12d\\\ &B\leq 15b+30d+60e\\\ &C\leq
42c+28d+14e\\\ &a,b,c,d,e,A,B,C\geq 0,\end{split}$

which does not have any integral solution, a contradiction to $|S|=74$.

∎

###### Proof of Theorem 1.6.

Assume that $S\subseteq\mathbb{F}_{7}^{3}$ is a $7$-progression-free set of
size $243$. Note that we have the following bounds.

######  Claim 3.

Every plane containing $36$ points contains at least $18$ $6$-lines and every
plane containing $35$, $34$ or $33$ contains at most $33$, $30$, $28$
$6$-lines, respectively. Moreover, $r_{7}(\mathbb{F}_{7}^{2})=36$ and
$r_{6}(\mathbb{F}_{7}^{2})=29$.

###### Proof of Claim 3.

Consider a plane $H$ containing $m$ points and let $x_{i}$ be the number of
$i$-lines in $H$ for $i\in[0,6]$. There are $56$ lines in the plane and
therefore

(4.4) $x_{0}+x_{1}+x_{2}+x_{3}+x_{4}+x_{5}+x_{6}=56$

By double counting the points in $H$ we get

(4.5) $x_{1}+2x_{2}+3x_{3}+4x_{4}+5x_{5}+6x_{6}=8m$

and by double counting the pairs of points in $H$ we get

(4.6) $x_{2}+3x_{3}+6x_{4}+10x_{5}+15x_{6}=\binom{m}{2}.$

If $m=36$ we take the difference of the (4.6) and two times (4.5) and get
$-2x_{1}-3x_{2}-3x_{3}-2x_{4}+3x_{6}=54$, implying that $x_{6}\geq 18$. If
$m\in\\{33,34,35\\}$, then by taking three times (4.4) minus two times (4.5)
plus (4.6) we get $3x_{0}+x_{1}+x_{4}+3x_{5}+6x_{6}=168-16m+\binom{m}{2}$ and
therefore $6x_{6}\leq 168-16m+\binom{m}{2}$ which gives the desired bounds.

The last two claims can be easily checked by computer search. ∎

If we now proceed analogously to the proof of Theorem 1.5 we again arrive at a
contradiction. ∎

##  5\. Proofs of the lower bounds

###### Proof of Theorem 1.7.

We consider three different types of $2$-dimensional layers:

  * •

$A:=[0,p-2]^{2}$,

  * •

$B:=[0,p-1]^{2}\setminus\\{(i,i)\ |\
i\in[0,p-1]\\}\setminus\big{(}\\{p-1\\}\times[0,\frac{p-3}{2}]\big{)}\setminus\big{(}[0,\frac{p-3}{2}]\times\\{p-1\\}\big{)}$,

  * •

$C:=\\{(i,i)\ |\ i\in[0,\frac{p-3}{2}]\\}$,

and three disjoint subsets of $\mathbb{F}_{p}^{n-2}$:

  * •

$\mathcal{A}:=[0,p-3]^{n-2}$,

  * •

$\mathcal{B}:=[0,p-2]^{n-2}\setminus[0,p-3]^{n-2}$,

  * •

$\mathcal{C}:=\bigcup_{j\in[1,n-2]}\\{x\in\mathbb{F}_{p}^{n-2}\;|\;(x_{j}=p-1)\land(x_{i}\in[0,p-3]\;\forall
i\neq j)\\}.$

We show that $S:=(\mathcal{A}\times A)\cup(\mathcal{B}\times
B)\cup(\mathcal{C}\times C)$ is $p$-progression-free.

First consider the case $n=3$. Let
$L:=\\{(a_{1},a_{2},a_{3})+(b_{1},b_{2},b_{3})i\ |\ i\in[0,p-1]\\}$ be a
$p$-progression in $\mathbb{F}_{p}^{3}$ with
$a_{1},a_{2},a_{3},b_{1},b_{2},b_{3}\in\mathbb{F}_{p}$.

  * •

Case 1: $b_{1}=0$ and $a_{1}\neq p-2:$

$[0,p-2]^{2}$ is $p$-progression-free and $|\\{(i,i)\ |\
i\in[0,\frac{p-3}{2}]\\}|<p$, therefore $L$ is not contained in $S$.

  * •

Case 2: $b_{1}=0$ and $a_{1}=p-2:$

$L^{\prime}:=\\{(a_{2},a_{3})+(b_{2},b_{3})i\ |\ i\in[0,p-1]\\}$ and
$\\{(i,i)\ |\ i\in[0,p-1]\\}$ are both lines in $\mathbb{F}_{p}^{2}$. If they
are not parallel or they are equal, they do intersect, and $L$ is not
contained in $S$. Otherwise we can rewrite $L^{\prime}=\\{(i,c+i)\ |\
i\in[0,p-1]\\}$ with $c\in[1,p-1].$ If $c\in[1,\frac{p-1}{2}]$ then
$c+(p-1)\in[0,\frac{p-3}{2}]$ (choose $i=p-1$) and $(p-2,p-1,c+(p-1))\in
L\setminus S.$

Similarly if $c\in[\frac{p+1}{2},p-1]$, then $p-1-c\in[0,\frac{p-3}{2}]$
(choose $i=p-1-c$) and $(p-2,p-1-c,p-1)\in L\setminus S.$ Therefore, $L$ is
not contained in $S$.

  * •

Case 3: $b_{1}\neq 0$:

Without the loss of generality let $b_{1}=1$ and $a_{1}=p-2$. If
$b_{2}=b_{3}=0$ then $L$ is not contained in $S$ because the $(p-2)$-layer and
$(p-1)$-layer of $S$ have no common point. Otherwise, without the loss of
generality, let $b_{2}\neq 0$ and therefore $\\{a_{2}+b_{2}i\ |\
i\in[0,p-1]\\}=[0,p-1].$ Assume that $L\subseteq S$. Then $a_{2}=p-1$ and
$a_{3}\in[\frac{p-1}{2},p-2]$ because the $(p-2)$-layer is the only layer
containing points with the coordinate $p-1$. Since the $(p-1)$-layer does not
have coordinates in $[\frac{p-1}{2},p-2]$, also $b_{3}\neq 0$ and consequently
$\\{a_{3}+b_{3}i\ |\ i\in[0,p-1]\\}=[0,p-1].$ As before it follows that
$a_{3}=p-1$ contradicting that $L\subseteq S$. Thus, $L$ is not contained in
$S$ and $S$ is $p$-progression-free.

Now consider $n>3$. We have already seen that every layer is $p$-progression-
free, so we only consider progressions $L:=\\{a+bi\ |\ i\in[0,p-1]\\}$
visiting $p$ non-empty layers. Let $m$ be the number of non-zero entries in
the first $n-2$ coordinates of $b$. Since only layers of type $C$ are placed
where one of the first $n-2$ coordinates is $p-1$ and all layers where two of
the first $n-2$ coordinates are $p-1$ are empty, $m$ is also the number of
type $C$ layers visited by $L$ and $m\leq p$.

  * •

If $m=1$, $L$ is not contained in $S$, analogously to the $3$-dimensional
case.

  * •

If $m\geq 2$ the last two coordinates of every point in $L$ are equal, since
the projection of $L$ in the last two coordinates is a line containing two
points in the main diagonal, or it is a single point in the main diagonal. Now
since only layers of type $B$ are placed where one of the first $n-2$
coordinates is $p-2$, $L$ also visits a layer of type $B$. Therefore $L$ is
not contained in $B$ because layers of type $B$ contain no points on the main
diagonal.

Finally, note that layers of type $A$ and $B$ contain $(p-1)^{2}$ points and
layers of type $C$ contain $(p-1)/2$ points and thus

$|S|=(p-1)^{2}(p-1)^{n-2}+\frac{p-1}{2}(n-2)(p-2)^{n-3}=(p-1)^{n}+\frac{n-2}{2}(p-1)(p-2)^{n-3}.$

∎

Layer ALayer BLayer CAAAAAAAAABBBBBBBCCCCCC Figure 1. A description of the
line-free set in Theorem 1.7 for $p=5$ and $n=4$. Figure 2. The last two
layers of the line-free set in Theorem 1.1 for $p=11$.

###### Proof of Theorem 1.1.

Let $k=\lfloor\sqrt{p}\rfloor$, $t=\lfloor p/k\rfloor$, $K:=[0,k-1]$ and
$T:=\\{jk-1\ |\ j\in[1,t]\\}$. Consider the set

$\begin{split}S:=&[0,p-3]\times[0,p-2]^{2}\\\
\cup&\\{p-2\\}\times([0,p-1]^{2}\setminus\\{(j,j)\ |\
j\in[0,p-1]\\}\setminus((K\cup\\{p-1\\})\times(T\cup\\{p-1\\})))\\\
\cup&\\{p-1\\}\times K\times T.\end{split}$

we will show that $S\subseteq\mathbb{F}_{p}^{3}$ is $p$-progression-free.

Let $L:=\\{(a_{1},a_{2},a_{3})+(b_{1},b_{2},b_{3})i\ |\ i\in[0,p-1]\\}$ be a
$p$-progression in $\mathbb{F}_{p}^{3}$ with
$a_{1},a_{2},a_{3},b_{1},b_{2},b_{3}\in\mathbb{F}_{p}$.

  * •

Case 1: $b_{1}=0$ and $a_{1}\neq p-2:$

$[0,p-2]^{2}$ is $p$-progression-free and $|K\times T|=kt<p$, therefore $L$ is
not contained in $S$.

  * •

Case 2: $b_{1}=0$ and $a_{1}=p-2:$

$L^{\prime}:=\\{(a_{2},a_{3})+(b_{2},b_{3})i\ |\ i\in[0,p-1]\\}$ and
$\\{(i,i)\ |\ i\in[0,p-1]\\}$ are both lines in $\mathbb{F}_{p}^{2}$. If they
are not parallel or they are equal, they do intersect, and $L$ is not
contained in $S$. Otherwise we can rewrite $L^{\prime}=\\{(i,c+i)\ |\
i\in[0,p-1]\\}$ with $c\in[1,p-1].$ $\\{(i,c+i)\ |\
i\in[0,k-1]\\}\cap(K\times(T\cup\\{p-1\\}))\neq\emptyset$ and therefore $L$ is
not contained in $S$.

  * •

Case 3: $b_{1}\neq 0$:

Without the loss of generality, let $b_{1}=1$ and $a_{1}=p-2$. If
$b_{2}=b_{3}=0$ then $L$ is not contained in $S$ because the $(p-2)$-layer and
$(p-1)$-layer of $S$ have no common point. Else, if $b_{2}\neq 0$ and
$b_{3}\neq 0$ then $\\{a_{2}+b_{2}i\ |\ i\in[0,p-1]\\}=\\{a_{3}+b_{3}j\ |\
j\in[0,p-1]\\}=[0,p-1].$ Since the $(p-2)$-layer is the only layer with $p-1$
entries but $(p-2,p-1,p-1)\not\in S$, $L$ is not contained in $S$. Finally, if
either $b_{2}=0$ or $b_{3}=0$ but not both, one of the last two coordinates is
constant, and the other one visits every possible value. Now again the
$(p-2)$-layer is the only layer with $p-1$ entries but the $(p-1)$-layer has
empty rows and columns wherever the $(p-2)$-layer has $p-1$ entries and
therefore $L$ is not contained in $S$.

Note that since $p$ is a prime and $k\geq 2$, $t\leq\frac{p-1}{k}$ and that
from the definition of $k$ it follows that

$\begin{split}&k\in[\sqrt{p}-1,\sqrt{p}+1]\\\ \Leftrightarrow\
&k^{2}-2\sqrt{p}k+p-1\leq 0\\\ \Leftrightarrow\ &k+\frac{p-1}{k}\leq
2\sqrt{p},\end{split}$

and therefore $k+t\leq 2\sqrt{p}$. Hence,

$\begin{split}|S|=&(p-2)(p-1)^{2}+(p^{2}-p-(kt-1)-k-t)+kt\\\
=&(p-2)(p-1)^{2}+p^{2}-p+1-k-t\\\ \geq&(p-2)(p-1)^{2}+p^{2}-p+1-2\sqrt{p}\\\
=&(p-1)^{3}+p-2\sqrt{p}\end{split}$

∎

Figure 3. The line-free set in Theorem 1.2 for $p=7$.

###### Proof of Theorem 1.2.

Let $p=7+24\ell$ for $\ell\in\mathbb{Z}_{\geq 0}$, let $A$ be the set of
quadratic residues, that is, $A=\\{a^{2}\ |\ a\in\mathbb{F}_{p}^{*}\\}$ and
$B:=\mathbb{F}_{p}^{*}\setminus A$. Note that $|A|=|B|=\frac{p-1}{2}$ and the
law of quadratic reciprocity yields

$\genfrac{(}{)}{}{}{-1}{p}=(-1)^{\frac{p-1}{2}}=(-1)^{3+12\ell}=-1,$

$\genfrac{(}{)}{}{}{2}{p}=(-1)^{\frac{p^{2}-1}{8}}=(-1)^{6+42\ell+72\ell^{2}}=1,$

$\genfrac{(}{)}{}{}{3}{p}=(-1)^{\frac{p-1}{2}\frac{3-1}{2}}\genfrac{(}{)}{}{}{p}{3}=(-1)^{3+12\ell}\genfrac{(}{)}{}{}{1}{3}=-1,$

and therefore $2\in A$ and $\\{-1,3\\}\subseteq B$. Note here that $A$ is a
subgroup of $\mathbb{F}_{p}^{*}$ and this means that multiplication by $2$ or
$\frac{1}{2}$ leaves elements of $A$ or $B$ in the same set, while
multiplication by $-1$, $3$ or $\frac{1}{3}$ changes the set. For instance,
$3a\in B$ for all $a\in A$ and $-\frac{3b}{2}=(-1)\cdot 3\cdot\frac{1}{2}\cdot
b\in B$ for all $b\in B$. Let

$\begin{split}S:=&[1,p-1]^{3}\cup\big{(}\\{(a,0,a)|a\in
A\\}\cup\\{(0,a,a)|a\in A\\}\big{)}\\\ \setminus&\big{(}\\{(a,a,a)|a\in
A\\}\cup\\{(a/2,a/2,a)|a\in A\\}\big{)}\\\ \cup&\big{(}\\{(3b/2,0,b)|b\in
B\\}\cup\\{(0,3b/2,b)|b\in B\\}\cup(3b,0,b)|b\in B\\}\cup\\{(0,3b,b)|b\in
B\\}\big{)}\\\ \setminus&\big{(}\\{(b,b,b)|b\in B\\}\cup\\{(3b/2,3b/2,b)|b\in
B\\}\cup\\{(b/3,b/3,b)|b\in B\\}\big{)}\\\
\setminus&\big{(}\\{(3b,-3b/2,b)|b\in B\\}\cup\\{(-3b/2,3b,b)|b\in
B\\}\big{)}\\\ \cup&\big{(}\\{(b,b,0)|b\in B\\}\cup\\{(2a,-a,0)|a\in
A\\}\cup(-a,2a,0)|a\in A\\}\big{)}.\end{split}$

We will show that $S$ is $p$-progression-free.

Note that $S$ is symmetric in the first two coordinates. We will therefore, in
this proof, skip one of two symmetric cases, whenever possible.

Let $L:=\\{(c_{1},c_{2},c_{3})+(d_{1},d_{2},d_{3})i\ |\ i\in[0,p-1]\\}$ be a
$p$-progression in $\mathbb{F}_{p}^{3}$ with
$c_{1},c_{2},c_{3},d_{1},d_{2},d_{3}\in\mathbb{F}_{p}$ and assume that
$L\subseteq S$.

First, assume that $d_{3}=0$.

  * •

Case 1: $c_{3}=0$:

Since $S$ contains no points where the third and one of the first two
coordinates is $0$, $L$ is not contained in $S$.

  * •

Case 2: $c_{3}\in A$:

Let $a:=c_{3}$. Since $(a,0,a)$ and $(0,a,a)$ are the only points where the
third coordinate is $a$ and one of the first two coordinates is $0$, we can
assume $(a,0,a)\in L$. If $d_{1}=0$ then $(a,a,a)\in L$, a contradiction. If
$d_{1}\neq 0$ then also $(0,a,a)\in L$ and consequently
$(\frac{a}{2},\frac{a}{2},a)=\frac{1}{2}(a,0,a)+\frac{1}{2}(0,a,a)\in L$,
again a contradiction.

  * •

Case 3: $c_{3}\in B$:

Let $b:=c_{3}$. First, assume $d_{1}\neq 0$ and $d_{2}\neq 0$. Since
$(\frac{3b}{2},0,b)$, $(0,\frac{3b}{2},b)$, $(3b,0,b)$ and $(0,3b,b)$ are the
only points where the third coordinate is $b$ and one of the first two
coordinates is $0$, we only have to consider the following cases:

If $(\frac{3b}{2},0,b)\in L$ and $(0,\frac{3b}{2},b)\in L$, then also
$(-\frac{3b}{2},3b,b)=(-1)(\frac{3b}{2},0,b)+2(0,\frac{3b}{2},b)\in L$, if
$(\frac{3b}{2},0,b)\in L$ and $(0,3b,b)\in L$, then also
$(b,b,b)=\frac{2}{3}(\frac{3b}{2},0,b)+\frac{1}{3}(0,3b,b)\in L$ and if
$(3b,0,b)\in L$ and $(0,3b,b)\in L$, then also
$(\frac{3b}{2},\frac{3b}{2},b)=\frac{1}{2}(3b,0,b)+\frac{1}{2}(0,3b,b)\in L$.
Consequently, we arrived at a contradiction. Now, if $d_{1}=0$ or $d_{2}=0$,
again $L$ has to contain one of the points $(\frac{3b}{2},0,b)$,
$(0,\frac{3b}{2},b)$, $(3b,0,b)$, $(0,3b,b)$ and therefore $L$ also contains
one of the points $(\frac{3b}{2},\frac{3b}{2},b)$, $(3b,\frac{-3b}{2},b)$,
$(\frac{-3b}{2},3b,b)$, again a contradiction.

Now assume that $d_{3}\neq 0$. If $d_{1}=d_{2}=0$, $L$ contains a point with a
zero last coordinate. We get that either $(b,b,0)$ and therefore also
$(b,b,b)$ is in $L$ for some $b\in B$ or $(2a,-a,0)\in L$ for some $a\in A$
and therefore also $(3b,\frac{-3b}{2},b)\in L$ for the unique $b\in B$ such
that $3b=2a$, both a contradiction.

In the remaining case $d_{3}\neq 0$ and at least one of $d_{1}$ and $d_{2}$ is
non-zero. Since there is no point in $S$ where both the third and one of first
two coordinates is zero, $L$ has to include a point with third coordinate
being zero and a different point where one of the other two coordinates is
zero. We are therefore left with checking the following cases, where $L$ is
given by a pair of two points in $S$. For some of these cases it is important
to note that $S$ contains no points where one of the first two coordinates is
$0$ and the other is in $B$.

  * •

$(b,b,0)\in L$ and $(0,a,a)\in L$:

$L=\bigg{\\{}\begin{pmatrix}b\\\ b\\\ 0\end{pmatrix}+k\begin{pmatrix}-b\\\
a-b\\\ a\end{pmatrix}\bigg{|}\ k\in\mathbb{F}_{p}\bigg{\\}}$

Since $a\neq b$, setting $k:=\frac{b}{b-a}$, we get
$(x,y,z):=(-\frac{ab}{b-a},0,\frac{ab}{b-a})\in L$. Now $x=-z$ and therefore
$z\in B$ and $x\in A$, so $-1=\frac{3}{2}$ or $-1=3$, a contradiction.

  * •

$(b,b,0)\in L$ and $(0,\frac{3b^{\prime}}{2},b^{\prime})\in L$:

$L=\bigg{\\{}\begin{pmatrix}b\\\ b\\\ 0\end{pmatrix}+k\begin{pmatrix}-b\\\
\frac{3b^{\prime}}{2}-b\\\ b^{\prime}\end{pmatrix}\bigg{|}\
k\in\mathbb{F}_{p}\bigg{\\}}$

Since $\frac{3b^{\prime}}{2}\neq b$, setting
$k:=\frac{b}{b-\frac{3b^{\prime}}{2}}$, we get
$(x,y,z):=(-\frac{3bb^{\prime}}{2(b-\frac{3b^{\prime}}{2})},0,\frac{bb^{\prime}}{b-\frac{3b^{\prime}}{2}})\in
L$. Now $x=-\frac{3}{2}z$ and therefore $x,z\in A$, so $-\frac{3}{2}=1$, a
contradiction.

  * •

$(b,b,0)\in L$ and $(0,3b^{\prime},b^{\prime})\in L$:

$L=\bigg{\\{}\begin{pmatrix}b\\\ b\\\ 0\end{pmatrix}+k\begin{pmatrix}-b\\\
3b^{\prime}-b\\\ b^{\prime}\end{pmatrix}\bigg{|}\
k\in\mathbb{F}_{p}\bigg{\\}}$

Since $3b^{\prime}\neq b$, setting $k:=\frac{b}{b-3b^{\prime}}$, we get
$(x,y,z):=(-\frac{3bb^{\prime}}{b-3b^{\prime}},0,\frac{bb^{\prime}}{b-3b^{\prime}})\in
L$. Now $x=-3z$ and therefore $x,z\in A$, so $-3=1$, a contradiction.

  * •

$(2a,-a,0)\in L$ and $(0,a^{\prime},a^{\prime})\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\ 0\end{pmatrix}+k\begin{pmatrix}-2a\\\
a^{\prime}+a\\\ a^{\prime}\end{pmatrix}\bigg{|}\ k\in\mathbb{F}_{p}\bigg{\\}}$

Since $a^{\prime}\neq-a$, setting $k:=\frac{a}{a+a^{\prime}}$, we get
$(x,y,z):=(\frac{2aa^{\prime}}{a+a^{\prime}},0,\frac{aa^{\prime}}{a+a^{\prime}})\in
L$. Now $x=2z$ and therefore $x,z\in A$, so $2=1$, a contradiction.

  * •

$(2a,-a,0)\in L$ and $(a^{\prime},0,a^{\prime})\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\
0\end{pmatrix}+k\begin{pmatrix}a^{\prime}-2a\\\ a\\\
a^{\prime}\end{pmatrix}\bigg{|}\ k\in\mathbb{F}_{p}\bigg{\\}}$

Assume $a^{\prime}\neq 2a$, then setting $k:=\frac{2a}{2a-a^{\prime}}$, we get
$(x,y,z):=(0,\frac{aa^{\prime}}{2a-a^{\prime}},\frac{2aa^{\prime}}{2a-a^{\prime}})\in
L$. Now $2y=z$ and therefore $y,z\in A$, so $1=2$, a contradiction. If
$a^{\prime}=2a$, then setting $k:=3$, we get $(x,y,z):=(2a,2a,6a)\in L$, a
contradiction since $6a\in B$.

  * •

$(2a,-a,0)\in L$ and $(\frac{3b}{2},0,b)\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\
0\end{pmatrix}+k\begin{pmatrix}\frac{3b}{2}-2a\\\ a\\\ b\end{pmatrix}\bigg{|}\
k\in\mathbb{F}_{p}\bigg{\\}}$

Assume $\frac{3b}{2}\neq 2a$, then setting $k:=\frac{2a}{2a-\frac{3b}{2}}$, we
get
$(x,y,z):=(0,\frac{3ab}{2(2a-\frac{3b}{2})},\frac{2ab}{2a-\frac{3b}{2}})\in
L$. Now $y=\frac{3z}{4}$ and therefore $z\in B$ and $x\in A$, so
$\frac{3}{4}=\frac{3}{2}$ or $\frac{3}{4}=3$, a contradiction. If
$\frac{3b}{2}=2a$, then setting $k:=3$, we get $(x,y,z):=(2a,2a,4a)\in L$, a
contradiction since $4a\in A$.

  * •

$(2a,-a,0)\in L$ and $(0,\frac{3b}{2},b)\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\ 0\end{pmatrix}+k\begin{pmatrix}-2a\\\
\frac{3b}{2}+a\\\ b\end{pmatrix}\bigg{|}\ k\in\mathbb{F}_{p}\bigg{\\}}$

Assume $\frac{3b}{2}\neq-3a$, then setting $k:=\frac{3a}{\frac{3b}{2}+3a}$, we
get
$(x,y,z):=(\frac{3ab}{\frac{3b}{2}+3a},\frac{3ab}{\frac{3b}{2}+3a},\frac{3ab}{\frac{3b}{2}+3a})\in
L$. Now $x=y=z$, a contradiction. If $\frac{3b}{2}=-3a$, then setting
$k:=-\frac{1}{2}$, we get $(x,y,z):=(3a,0,a)\in L$, so $1=3$, a contradiction.

  * •

$(2a,-a,0)\in L$ and $(3b,0,b)\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\
0\end{pmatrix}+k\begin{pmatrix}3b-2a\\\ a\\\ b\end{pmatrix}\bigg{|}\
k\in\mathbb{F}_{p}\bigg{\\}}$

Since $b\neq a$, then setting $k:=\frac{a}{a-b}$, we get
$(x,y,z):=(\frac{ab}{a-b},\frac{ab}{a-b},\frac{ab}{a-b})\in L$. Now $x=y=z$, a
contradiction.

  * •

$(2a,-a,0)\in L$ and $(0,3b,b)\in L$:

$L=\bigg{\\{}\begin{pmatrix}2a\\\ -a\\\ 0\end{pmatrix}+k\begin{pmatrix}-2a\\\
a+3b\\\ b\end{pmatrix}\bigg{|}\ k\in\mathbb{F}_{p}\bigg{\\}}$

since $3b\neq-a$, setting $k:=\frac{a}{a+3b}$, we get
$(x,y,z):=(\frac{6ab}{a+3b},0,\frac{ab}{a+3b})\in L$. Now $x=6z$ and therefore
$z\in B$ and $x\in A$, so $6=\frac{3}{2}$ or $6=3$, a contradiction.

Finally, note that the layer with third coordinate $0$ contains
$|B|+2|A|=\frac{3}{2}(p-1)$ points, layers with the third coordinate in $A$
contain $(p-1)^{2}$ points and layers with the third coordinate in $B$ contain
$(p-1)^{2}-1$ points and thus

$|S|=\frac{p-1}{2}(p-1)^{2}+\frac{p-1}{2}((p-1)^{2}-1)+\frac{3}{2}(p-1)=(p-1)^{3}+(p-1).$

In the special case of $p=7$, the layers with third coordinate in $B$ actually
contain $(p-1)^{2}=36$ points, since $\frac{3}{2}=\frac{1}{3}$, thus giving
the lower bound of $225$.

∎

### Acknowledgements

C.E. was supported by a joint FWF-ANR project ArithRand (I 4945-N and
ANR-20-CE91-0006). J.F. was supported by the Austrian Science Fund (FWF) under
the project W1230. D.G.S. was supported by the ERC Advanced Grant ”Geoscape”.
B.K. was supported by the ÚNKP-22-3, New National Excellence Program of the
Ministry for Culture and Innovation from the source of the National Research,
Development and Innovation Fund. P.P.P. was supported by the Lendület program
of the Hungarian Academy of Sciences (MTA) and by the National Research,
Development and Innovation Office NKFIH (Grant Nr. K124171).

## References

  * [1] Douglas Zare (https://mathoverflow.net/users/2954/douglas-zare)  “Blocking sets in three dimensional finite affine spaces”, MathOverflow  URL: https://mathoverflow.net/q/229442
  * [2] Ara Aleksanyan and Mihran Papikian  “On Blocking Sets of Affine Spaces”  arXiv, 1999  DOI: 10.48550/ARXIV.MATH/9910084
  * [3] Noga Alon  “Tools from higher algebra”  In _Handbook of combinatorics, Vol. 1, 2_ Elsevier Sci. B. V., Amsterdam, 1995, pp. 1749–1783
  * [4] Anurag Bishnoi, Jozefien D’haeseleer, Dion Gijswijt and Aditya Potukuchi  “Blocking sets, minimal codes and trifferent codes”  arXiv, 2023  DOI: 10.48550/ARXIV.2301.09457
  * [5] A.. Brouwer and A. Schrijver  “The blocking number of an affine space”  In _J. Combinatorial Theory Ser. A_ 24.2, 1978, pp. 251–253  DOI: 10.1016/0097-3165(78)90013-4
  * [6] Ernie Croot, Vsevolod F. Lev and Péter Pál Pach  “Progression-free sets in $\mathbb{Z}^{n}_{4}$ are exponentially small”  In _Ann. of Math. (2)_ 185.1, 2017, pp. 331–337  DOI: 10.4007/annals.2017.185.1.7
  * [7] Benjamin Lent Davis and Diane Maclagan  “The card game SET”  In _Math. Intelligencer_ 25.3, 2003, pp. 33–40  DOI: 10.1007/BF02984846
  * [8] Yves Edel  “Extensions of generalized product caps”  In _Des. Codes Cryptogr._ 31.1, 2004, pp. 5–14  DOI: 10.1023/A:1027365901231
  * [9] Jordan S. Ellenberg and Dion Gijswijt  “On large subsets of $\mathbb{F}^{n}_{q}$ with no three-term arithmetic progression”  In _Ann. of Math. (2)_ 185.1, 2017, pp. 339–343  DOI: 10.4007/annals.2017.185.1.8
  * [10] Christian Elsholtz, Benjamin Klahn and Gabriel F. Lipnik  “Large subsets of $\mathbb{Z}_{m}^{n}$ without arithmetic progressions”  In _Designs, Codes and Cryptography_ , 2022  DOI: 10.1007/s10623-022-01145-w
  * [11] Christian Elsholtz and Gabriel F. Lipnik  “Exponentially larger affine and projective cap”  In _Mathematika_ 69.1, 2023, pp. 232–249  DOI: 10.1112/mtk.12173
  * [12] Christian Elsholtz and Péter Pál Pach  “Caps and progression-free sets in $\mathbb{Z}_{m}^{n}$”  In _Des. Codes Cryptogr._ 88.10, 2020, pp. 2133–2170  DOI: 10.1007/s10623-020-00769-0
  * [13] P. Frankl, R.. Graham and V. Rödl  “On subsets of abelian groups with no $3$-term arithmetic progression”  In _J. Combin. Theory Ser. A_ 45.1, 1987, pp. 157–161  DOI: 10.1016/0097-3165(87)90053-7
  * [14] H. Furstenberg and Y. Katznelson  “A density version of the Hales-Jewett theorem”  In _J. Anal. Math._ 57, 1991, pp. 64–119  DOI: 10.1007/BF03041066
  * [15] Robert E. Jamison  “Covering finite fields with cosets of subspaces”  In _J. Combinatorial Theory Ser. A_ 22.3, 1977, pp. 253–266  DOI: 10.1016/0097-3165(77)90001-2
  * [16] Zhi Jiang  “Improved explicit upper bounds for the Cap Set Problem”  arXiv, 2023  DOI: 10.48550/ARXIV:2103.06481
  * [17] Leo Moser  “Problems for Solution, P.170”  In _Canadian Mathematical Bulletin_ 13.2 Cambridge University Press, 1970, pp. 267–272  DOI: 10.1017/S0008439500031660
  * [18] Péter Pál Pach  “Bounds on the size of progression-free sets in $\mathbb{Z}^{n}_{m}$”  In _Unif. Distrib. Theory_ 17.1, 2022, pp. 1–10  DOI: doi:10.2478/udt-2022-0005
  * [19] Péter Pál Pach and Richárd Palincza  “Sets avoiding six-term arithmetic progressions in $\mathbb{Z}_{6}^{n}$ are exponentially small”  In _SIAM J. Discrete Math._ 36.2, 2022, pp. 1135–1142  DOI: 10.1137/21M1413766
  * [20] Fedor Petrov and Cosmin Pohoata  “Improved bounds for progression-free sets in $C_{8}^{n}$”  In _Israel J. Math._ 236.1, 2020, pp. 345–363  DOI: 10.1007/s11856-020-1977-0
  * [21] D… Polymath  “Density Hales-Jewett and Moser numbers”  arXiv, 2010  DOI: 10.48550/ARXIV.1002.0374
  * [22] Aaron Potechin  “Maximal caps in ${\rm AG}(6,3)$”  In _Des. Codes Cryptogr._ 46.3, 2008, pp. 243–259  DOI: 10.1007/s10623-007-9132-z
  * [23] Fred Tyrrell  “New lower bounds for cap sets”  arXiv, 2022  DOI: 10.48550/ARXIV.2209.10045

##  6\. Appendix

Figure 4. A line-free set showing $r_{p}(\mathbb{F}_{p}^{3})\geq 70$


11institutetext:  Old Dominion University, Norfolk, VA 23529, USA
22institutetext: NVIDIA, Santa Clara, CA 95051-0952, USA 33institutetext:
Fermi National Accelerator Laboratory, Batavia, IL 60510

# Porting numerical integration codes from CUDA to oneAPI: a case
study††thanks: code available at https://github.com/marcpaterno/gpuintegration

Ioannis Sakiotis 11 Kamesh Arumugam 22 Marc Paterno 1133 Desh Ranjan 11 Balša
Terzić 11 Mohammad Zubair 11

(January 2023)

###### Abstract

We present our experience in porting optimized CUDA implementations to oneAPI.
We focus on the use case of numerical integration, particularly the CUDA
implementations of PAGANI and $m$-Cubes. We faced several challenges that
caused performance degradation in the oneAPI ports. These include differences
in utilized registers per thread, compiler optimizations, and mappings of CUDA
library calls to oneAPI equivalents. After addressing those challenges, we
tested both the PAGANI and $m$-Cubes integrators on numerous integrands of
various characteristics. To evaluate the quality of the ports, we collected
performance metrics of the CUDA and oneAPI implementations on the Nvidia V100
GPU. We found that the oneAPI ports often achieve comparable performance to
the CUDA versions, and that they are at most $10\%$ slower.

##  1 Introduction

Historically, general-purpose GPU programming has been characterized by
divergent architectures and programming models. A lack of widely adopted
common standards led to the development of different ecosystems comprised of
compilers and tools that were practically exclusive to specific GPU
architectures. Most importantly, the emergent architectures themselves were
not compatible with all ecosystems. Portability could only be achieved through
the maintenance of multiple code bases. Traditionally, the proprietary CUDA
programming model has been the most popular but is exclusively targeted to
Nvidia GPUs.

In the absence of universally adopted standards, a viable solution for
achieving general portability is to rely on platform-agnostic programming
models that target multiple architectures via a unifying interface. This
enables the execution of a single code base across various architectures.
These programming models would ideally enable the utilization of platform-
specific low-level features on their native hardware. This would allow highly-
optimized implementations in such portable programming models to remain
competitive with platform-specific alternatives. Without these capabilities,
use cases with extreme performance requirements would disqualify the use of
such portable models.

The need for performant multi-platform execution is only increasing with the
emergence of exascale supercomputers such as Frontier and Aurora that do not
carry Nvidia GPUs. Projects requiring computing cores at that scale must
develop new software solutions compatible with non-Nvidia GPUs or port
existing CUDA implementations without significant loss of performance.

Portable programming models such as RAJA, Kokkos, and oneAPI have been in
development and are already available for use. These portable alternatives
lack maturity when compared to proprietary alternatives. As such, applications
requiring portable solutions must be evaluated to quantify any necessary
concessions.

In this paper, we discuss the porting process of two numerical integration
implementations, PAGANI and m-Cubes, from CUDA to oneAPI’s Data Parallel C++
language. The oneAPI ecosystem provides a suite of compilers, libraries, and
software tools, including the Data Parallel Conversion tool (DPCT) that
automates the majority of the porting process. Reliance on the C++ and SYCL
standards as well as the capability to quickly port large CUDA
implementations, places oneAPI at the forefront of the portability initiative.

We faced challenges during the porting process due to the lack of support for
certain libraries utilized by the CUDA implementation. For example, the CUDA
implementation of PAGANI uses the Nvidia Thrust library to perform common
parallel operations on the host side, such as inner product and min-max. Even
though there is a multitude of library options in oneAPI, we encountered
difficulties with the DPCT mapping of Nvidia Thrust library calls, which were
not fully supported on all backends.

We also observed performance degradation for the ported oneAPI
implementations. We conducted numerous experiments with integrands of various
characteristics to identify the issues. Most of these issues pertained to
optimization differences between the NVCC and Clang compilers, and time
differences when executing mathematical functions. After addressing these
challenges, the oneAPI ports were at most $10\%$ slower than the optimized
CUDA versions. We observe that the cases with the highest performance
penalties for the oneAPI ports, require significantly more registers than the
CUDA originals. This decreases the occupancy in the oneAPI implementation and
causes performance degradation. When the number of registers is similar to the
CUDA version, we observe penalties lower than $5\%$.

The remainder of this paper is structured as follows. First, we provide
background information on oneAPI and other portability solutions in section
$2$. Then, we discuss the two numerical integration CUDA implementations in
section $3$. Section $4$ details the porting process and challenges we faced
using DPCT and the oneAPI platform. In section $5$, we present a performance
comparison of the CUDA and oneAPI implementations of PAGANI and m-Cubes. We
finish in section $6$ with a discussion of our conclusions regarding the
oneAPI platform’s viability and ease of use. We demonstrate that the oneAPI
implementation does not induce significant performance penalties and that it
is a viable platform for attaining performance on Nvidia GPUs.

##  2 Background

There are multiple programming models targeting different architectures. Among
the most prominent, are OpenCL [24], [16], OpenACC [5], OpenMP [1], RAJA,
Alpaka [30], and Kokkos [10]. The Khronos group was the first to address
portability by developing the OpenCL standard to target various architectures.
The same group later followed with the SYCL standard. SYCL is a higher-level
language that retained OpenCL features but significantly improved ease of use
with the utilization of C++ and the adoption of a single-source model. There
are multiple implementations of SYCL such as DPC++, ComputeCpp, HipSYCL, and
triSYCL [28]. DPC++ is conformant to the latest SYCL and C++ standards and is
integrated into the oneAPI ecosystem [2].

###  2.1 Intel oneAPI and DPC++

The oneAPI toolkit was released in 2020, providing a programming platform with
portability across multiple architectures at the core of its mission. The
toolkit included various tools along with the DPC++ language which was based
on the SYCL and C++ standards [8]. The reliance on these open standards that
are intended to evolve over time is one of the most attractive features of
DPC++. Such evolution is facilitated by DPC++ extensions with various features
that can be later introduced to the standards after periods of
experimentation. Such examples include the use of _Unified Memory_ and
_filtered Device selectors_ , which were missing from SYCL 1.2.1 but were
later included in the SYCL 2020 standard. DPC++ achieves execution platform
portability through its use of SYCL and various backends (implemented as
shared libraries) that interface with particular instruction sets such as PTX
for Nvidia GPUs and SPIR-V for Intel devices. It is worth noting that there is
no reliance on OpenCL, which is instead one of several available backends. As
such, DPC++ implementations can target various CPUs, GPUs, and FPGAs. This is
a similar approach to Kokkos, Alpaka, and RAJA.

###  2.2 CUDA-backend for SYCL

While CUDA is the native and most performant programming model for Nvidia
GPUs, Nvidia provided support to the OpenCL API [25]. As a result, non-CUDA
implementations could be executed on Nvidia GPUs. The ComputeCpp
implementation of SYCL by CodePlay, provided such functionality through
OpenCL, but its performance was not comparable to native CUDA as not all
functionality was exposed [3].

As such, CodePlay developed the CUDA backend for DPC++, which is part of the
LLVM compiler project. CUDA support is not enabled by default and is at an
experimental stage. To enable the backend, we must build the LLVM compiler
project for CUDA. This can be achieved through easy-to-follow instructions
that involve CUDA-specific flags, and the use of _clang++_ instead _dpcpp_ to
compile source code. As a result, DPC++ code can generate PTX code by using
CUDA directly instead of relying on the OpenCL backend. This approach not only
enables the use of Nvidia libraries and profiling tools with DPC++ but also
the capability to theoretically achieve the same performance as CUDA.

###  2.3 Related Work

The oneAPI programming model may not be as mature as CUDA but the literature
already includes several examples of utilizing DPC++. The authors of [12]
validated the correctness of a DPC++ tsunami simulator ported from CUDA. A
Boris Particle Pusher port from an openMP version was discussed in [27], where
a DPC++ implementation was $10\%$ slower than the optimized original. In [14],
CUDA and DPC++ implementations of a matrix multiplication kernel were compared
on different matrix sizes; the execution time on an Nvidia GPU was slower with
DPC++ code by $7\%$ on small problem sizes but as much as $37\%$ on larger
ones. On the contrary, [15] and [13] included experiments where a DPC++
biological sequence alignment code showed no significant performance penalty
compared to CUDA, and even a case of $14\%$ speedup. Spare matrix-vector
multiplication kernels and _Krylov_ solvers in [26] reached $90\%$ of a CUDA
version’s bandwidth. There were also cases with non-favorable performance for
DPC++ ports. In [17] a bioinformatics-related kernel performed twice as fast
in CUDA and HIP than in DPC++. In [11] DPC++ versions generally reported
comparable performance to CUDA but there were multiple cases where the penalty
ranged from $25-190\%$.

There seems to be a deviation in the attainable performance. This is
reasonable due to the variety of applications and the relatively early stage
of development for the oneAPI ecosystem. We also expect that the level of
optimization in CUDA implementations is an important factor. In our
experience, highly optimized codes typically yield performance penalties in
the range ($5-10\%$). There are multiple cases displaying approximately $10\%$
penalty compared to native programming models. This indicates that DPC++ can
achieve comparable performance to CUDA, though careful tuning and additional
optimizations may be needed.

##  3 Numerical Integration Use Case

Numerical integration is necessary for many applications across various fields
and especially physics. Important examples include the simulation of beam
dynamics and parameter estimation in cosmological models [21] [6] [9]. Even
ill-behaving integrands (oscillatory, sharply peaked, etc.) can be efficiently
integrated with modest computational resources, as long the integration space
is low dimensional (one or two variables). On the contrary, solving medium to
high-dimensional integrands is often infeasible on standard computing
platforms. In such cases, we must execute on highly parallel architectures to
achieve performance at scale. There are a few GPU-compatible numerical
integration algorithms [22] [23] [7] [19] [29]. Unfortunately, exploration of
execution-platform portability has been limited, with CUDA being the most
common choice. Since CUDA is a proprietary language, such optimized
implementations cannot be executed on non-Nvidia GPUs. To our knowledge, the
only mentions of potential portability in numerical integration libraries are
found in [22] where a Kokkos implementation of the PAGANI integrator is
briefly mentioned to be in development and in [23] which compares the CUDA
implementation of $m$-Cubes with an experimental Kokkos version.

###  3.1 PAGANI

PAGANI is a deterministic quadrature-based algorithm designed for massively
parallel architectures. The algorithm computes an integral by evaluating the
quadrature rules, which are a series of weighted summations of the form
$\sum_{i=1}^{f_{eval}}w_{i}\cdot f(x_{i})$. The computation involves an
integrand function $f$ which we invoke at the $d$-dimensional points $x_{i}$.
Each point $x_{i}$ has a corresponding weight $w_{i}$ and there are $f_{eval}$
such points in the summation. PAGANI computes an initial integral and error
estimate, and it progressively improves its accuracy until reaching a user-
specified threshold. The accuracy improvement is achieved by applying the
quadrature rules in smaller regions of the integration space and accumulating
those values to get the integral estimate.

PAGANI, is a memory-hungry algorithm due to its rapid sub-division of the
integration space; high-accuracy executions typically generate millions of
regions for ill-behaving integrands. This is a design choice that is greatly
suited for parallel systems such as GPUs. A significant component of the
algorithm involves the careful classification of the regions, with the goal of
eliminating certain classes of regions from the limited available memory.
While this is a critical step for the progression of the algorithm, it is not
the most computationally intensive.

Algorithm 1 PAGANI

1:procedure pagani-kernel($f$, $n$, $d$, $R[n]$, $I[n]$, $E[n]$, $K[n]$, $q$)

2: for each region $i$ to $n$ parallel do

3: local[5] $\leftarrow 0$ $\triangleright$ thread-local storage for
evaluations of the five quad. rules

4: if thread 0 then

5: shared [$f_{eval}$] $\leftarrow 0$ $\triangleright$ store function
evalutions in shared mem.

6: $S\leftarrow$ Init-shared-mem($R$) $\triangleright$ collection of info. on
region processed

7: for each $f_{id}<f_{eval}$ parallel do

8: x[$d$] $\leftarrow$ compute-eval-point($S$, $d$, $q$)

9: shared[$t_{id}$] $\leftarrow f(x)$

10: for $i\leftarrow 0$ to $4$ do $\triangleright$ per quad-rule

11: local[$i$] $\mathrel{+}=f(x)\cdot q.weights[t_{id}]$

12: for $i\leftarrow 0$ to $4$ parallel do $\triangleright$ per quad-rule

13: local[$i$] $\leftarrow$ Reduce($local[i]$)

14: if thread 0 then

15: $S.k\leftarrow$ Compute-split-axis($shared,q,d$)

16: ${S.e}\leftarrow$ find-max-err($local,q$)

17: $S.i\leftarrow$ local[$0$]

18: $I,E,K\leftarrow S.i,S.e,S.k$ $\triangleright$ final results

In fact, it is the Evaluate method (listed in Algorithm $2$ of [22]) that
consistently takes more than $90\%$ of total execution time. Its function is
to compute an integral/error estimate for each region and select one of the
dimensional axes for splitting. As such, it can be viewed as the core of
PAGANI, both from an algorithmic and performance standpoint.

We describe the Evaluate method in Listing 1, which we will refer to as
Pagani-kernel. The input consists of the integrand $f$, which will be
repeatedly evaluated within the boundaries of each region, the number of
regions $n$, and the dimensionality $d$. Then, we supply the region
coordinates $R$ in the form of two lists $R.c$ (low boundary of a region on
each axis) and $R.l$ (length of a region on each axis). Next, we use the lists
$I$ and $E$ to store the integral and error estimate of each region, while the
list $K$ contains an index representing the axis that should be split in the
future to improve the estimates. Finally, we supply a collection of read-only
lists through the parameter $q$. Each list in $q$ contains constant values
used to generate and permute over the evaluation points and weights associated
with the five quadrature rules.

Each thread group processes a single region $i$ and uses all threads in a
group to parallelize the integrand function evaluations. On lines $3$ to $6$
each thread initializes its variables, while thread zero initializes variables
shared in the group. The $local$ list at line $3$ has a single entry for each
of the five quadrature rules. All threads compute their assigned function
evaluations for each quadrature rule and accumulate the results in the $local$
list. Each of the five entries in $local$, contains the partial results
computed by the thread for the respective quadrature rule summations. $shared$
is an array in shared memory. It is used for the temporary storage of function
evaluations that need to be retained and not just accumulated.

At line $6$, thread zero of each thread group will cache the sub-region
assigned to its thread group into the $S$ shared memory location. $S$ contains
the region boundaries, global boundaries of integration space, volume, longest
dimensional axis in the boundaries of the region, integral/error-estimate, and
axis to split. All of these values are associated with a single region.

Then the for-loop at line $7$ executes each function evaluation $f_{id}$ in
parallel. At line $8$, each thread computes the integrand $f$ at a
$d$-dimensional point $x$. This involves reading $d$ values from one of the
lists in $q$ (generators list) and scaling those values to the region
boundaries of each axis (bounds stored in $S$). Then the threads evaluate the
integrand at line $9$ and store the result at the appropriate index in the
list $shared$. Once that is done, the threads use the weights associated with
each quadrature rule (stored in $q.w$) to multiply by the function evaluation
at line $11$.

Once the threads exit the for-loop of line $7$, all function evaluations have
been completed. At line $13$, the algorithm accumulates the results from all
threads in the group, which yield five integral estimates per region. At the
end, thread zero will store the reduced data as the final estimates in line
$18$. Thread zero also computes the split axis at line $15$. This requires
reading two values from a list in $q$, and the serial iteration of several
function evaluations which are stored in $shared$. Thread zero performs a
series of simple operations and then stores the index of the axis to split at
the shared location $S.k$. The error estimate is computed at line $16$ and
involves determining the maximum difference between the estimates reported by
four of the five quadrature rules. This requires reading from certain lists in
$q$ and the produced result is the error estimate for the integral estimate of
that region.

The CUDA implementation was optimized for the Nvidia _V100_ GPU. The kernel is
launched in groups of $64$ threads and the loop at line $7$ in iterated in a
strided fashion. This allows the threads to coalesce accesses to global memory
(reads in $q$) needed for the function evaluations. Additionally, the kernel
relies on the “ldg” intrinsic for reading the arrays within $q$, suggesting to
the compiler their placement in the read-only cache.

It is worth noting that this is a compute-bound kernel with hundreds and often
thousands (depending on the integrand) of floating-point operations per byte.
While the algorithm typically minimizes branch divergence, it cannot be
avoided entirely due to thread zero performing operations on shared memory and
storing the final results in global memory. Due to fast-access shared memory
utilization, the implementation uses synchronization points in the following
cases: prior to the loop in line $7$, between the main loop strides, and
before reducing the data at line $13$. The reduction utilizes warp-level
primitives to accumulate results within the warp, but needs shared memory to
reduce the values among the warps.

###  3.2 $m$-Cubes

$m$-Cubes is a probabilistic Monte Carlo algorithm based on the VEGAS
integrator [20]. It operates by randomizing the sample generation across the
integration space to solve integrands and relies on the standard deviation of
the Monte Carlo estimate to produce error estimates for the computation. Just
like VEGAS, $m$-Cubes utilizes importance and stratified sampling to
accelerate the Monte Carlo rate of convergence. Importance sampling works by
progressively approximating a step-wise function that mirrors the shape of the
integrand. The step-wise function, allows the algorithm to sample from a
probability distribution that concentrates the function evaluations where the
integrand varies the most. This approach achieves the adaptive behavior that
is necessary for numerical integrators to be efficient in high-dimensional
spaces.

In terms of stratified sampling, $m$-Cubes partitions the integration space
into $m$ sub-cubes that are sampled separately. This improves the convergence
rate and eliminates statistical errors and inefficient sampling that could
occur on sharp-peaked integrands. The sub-division resolution is dictated by a
user-specified number of samples per iteration. $m$-Cubes partitions the
integration space such that a $2$-point Monte Carlo integration can be
executed on each partition while generating approximately $n$ samples.

We describe the main kernel in Listing 2, though it is also detailed in [23].
The implementation relies on many thread groups of size $x$, requiring a total
$\frac{m}{x}$ such groups where $m$ is the number of sub-cubes resulting from
the initial sub-division. Each thread processes its own sub-cube batch of size
$s$.

Algorithm 2 $m$-Cubes

1:procedure mcubes-kernel($f$, $d$, $m$, $s$, $p$, $B[500\cdot d]$,
$C[500\cdot d]$, $r$)

2: for $m$/($b$ total threads) parallel do

3: Set-Random-Generator($seed$)

4: $I$, $E\leftarrow 0$ $\triangleright$ cumulative estimates of thread

5: for $t=0$ to $s$ do

6: $I_{t}$, $E_{t}\leftarrow 0$ $\triangleright$ estimates of sub-cube t

7: for $k\leftarrow 1$ to $p$ do

8: $x[1:d]\leftarrow$ Generate()

9: $I_{k},E_{k}\leftarrow$ Evaluate($f,x$)

10: $I_{t}\leftarrow I_{t}+I_{k}$ $\triangleright$ Accumulate sub-cube
contributions

11: $E_{t}\leftarrow E_{t}+E_{k}$

12: $b[1:d]\leftarrow$ Get-Bin-ID($x$)

13: for $j\leftarrow 1$ to $d$ do $\triangleright$ Store bin contributions

14: AtomicAdd($C[b[j]]$, $I_{k}^{2}$)

15: $E_{t}\leftarrow$ UpdateVariance($E_{t}$, $I_{t}$, $p$)

16: $I\leftarrow I+I_{t}$ $\triangleright$ update cumulative values

17: $E\leftarrow E+E_{t}$

18: $I\leftarrow$ Reduce($I$)

19: $E\leftarrow$ Reduce($E$)

20: if thread 0 within group then

21: AtomicAdd($r[0]$, $I$)

22: AtomicAdd($r[1]$, $E$)

The input of mcubes-kernel, consists of the integrand $f$ of dimensionality
$d$, the number of sub-cubes $m$, sub-cube batch size $s$, number of samples
per sub-cube $p$, bin bounds $B[500\cdot d]$, bin contributions $C[500\cdot
d]$, and result $r$. The mcubes-kernel outputs an estimate for the integral,
variance, and the accumulated bin contributions $C$ from all samples.

All $m/b$ threads are launched in parallel at line $2$. At line $3$ the random
number generator is initialized. Each thread processes its assigned sub-cubes
serially through the loop at line $5$, and takes $p$ samples per sub-cube at
line $7$. The results from each sample are stored in $I_{k}$ and $E_{k}$ (line
$9$). We accumulated those results to $I_{t}$ and $E_{t}$ (line $11$) which
correspond to the sub-cube estimates. Finally, we accumulate the results from
all sub-cubes in $I$ and $E$.

When evaluating each sample, we first generate a $d$-dimensional point $x$ in
line $8$. The kernel evaluates the integrand $f$ at that point(line $9$) .
Then, at line $12$, the kernel identifies the IDs of the $d$ bins used to to
generate the point $x$, and stores their indices in the list $C$, at the local
array $b$. Those indices are used to update the bin contribution at line $14$
through atomic addition.

The bin contribution is the square of the integral estimate $I_{k}$. Then, the
variance is updated at line $15$ and the error and integral estimate for the
$s$ sub-cubes is accumulated in lines $16$ and $17$.

The kernel then uses block-reduction to accumulate $I$ and $E$ from each
thread at lines $18$ and $19$. Finally, two atomic additions from thread $0$
of each thread-group leads to the computation of the total integral and error
estimate from all thread-groups. These are stored in memory location $r$ at
lines $21$ and $22$.

The CUDA implementation was optimized for the _V100_ GPU. The kernel consisted
of $128$ threads per block and utilized atomic addition for accumulating
results originating from different thread blocks. This is particularly
important for updating the random-access bin contributions, which refer to
$500$ bins per dimensional axis. The reduction operations at lines $18/19$
operate on local memory and utilize warp-level primitives, though shared
memory is used to accumulate the values from the different warps.

##  4 Porting Process

The maturity of the CUDA programming model along with the more widespread
utilization of highly performant Nvidia GPUs make CUDA an intuitive choice for
high-performance applications. As such, PAGANI and $m$-Cubes were designed and
optimized for CUDA on a V100 Nvidia GPU [22] [23]. This makes DPCT the most
appropriate tool to facilitate the porting process from CUDA to DPC++.

###  4.1 DPCT

DPCT is intended to automate the majority of CUDA code migration to DPC++,
instead of performing a total conversion [4]. In our experience as well as
those reported in [14], [18] and many others, DPCT functions exactly as
intended. An easy-to-complete conversion process requires few manual code
insertions. When manual editing is needed, DPCT displays helpful suggestions
in the form of comment blocks to guide the user. There were certain code
segments that were functional but needed simple fixes to improve performance.
e.g. use of 3D _nd_item_ instead of 1D equivalent. In our experience, those
cases were few and we expect that such effects will be less pronounced as
oneAPI and DPCT evolve. Expert users are anticipated to produce higher-quality
implementations than automated tools, but even then DPCT greatly facilitates
the porting process by automating the tedious and often error-prone
translation of API calls and indexing schemes.

###  4.2 Challenges

####  4.2.1 Porting Issues with Nvidia Thrust Library

PAGANI uses Thrust to perform common parallel operations on the host side,
such as reduction, dot-product, prefix sum, and finding the minimum/maximum
within a list. oneAPI provides an array of libraries to accommodate all of our
required operations. While the CUDA implementation defaults to cuRAND to
generate the random samples in $m$-Cubes, we use a simple custom-made random
generator in both the CUDA and oneAPI versions. This facilitates the
comparison between the two implementations and avoids a comparison of random-
number-generator libraries.

The usage of oneAPI libraries is more limited on the CUDA-backend. In our
experience, certain calls to library functions inserted by DPCT worked on
Intel GPUs such as the P630 but were not viable on the CUDA-backend, yielding
compile-time errors in certain cases such as when using _mkl::stats::min_max_.
Additionally, DPCT headers provided library calls with CUDA-like syntax in the
_dpct_ namespace. This worked on Intel GPUs but not on the CUDA-backend, e.g
_dpct::inner_product_ , _dpct::get_device_pointer_. Still, we found
alternatives from the oneDPL and oneMKL libraries, which allowed us to execute
successfully on both the CUDA-backend and on Intel devices.

####  4.2.2 Performance Degradation

We encountered more difficulties when attempting to achieve comparable
performance to the original CUDA implementation. We observed that the NVCC and
Clang compilers optimized differently when using the O3 optimization flag; we
needed to set inline thresholds and to manually enable/disable loop-unrolling
in certain locations. Furthermore, code elimination in some of our
benchmarking kernels that invoked mathematical functions on the device was not
equivalent in the two compilers. Additionally, we found that for our
implementation, work-group reduction caused approximately 10% slowdown in the
oneAPI version compared to manual reduction using shared memory and warp-
intrinsics.

Another challenge in our attempt to achieve comparable performance to CUDA was
deviations in the performance of SYCL and CUDA mathematical functions.
Exponential functions displayed comparable performance on benchmark kernels.
On the contrary, we observed a slowdown of various degrees in oneAPI when
using _power_ or trigonometric functions. This is most likely attributed to
the compilers utilizing different optimizations. We did not use any _fast-
math_ flags, since high accuracy is critical in numerical integration use
cases.

Finally, the use of atomic addition in $m$-Cubes caused orders of magnitude
slowdown on both the mcubes-kernel and benchmark kernels. This was attributed
to the lack of an architecture-specific flag that must be set to enable
efficient atomics when supported. After setting the Volta architecture flag,
atomic addition was as performant as in the native CUDA implementation.

####  4.2.3 Software Engineering Issues

Utilization of the Catch2 testing framework and CMake was largely successful
but more error-prone in oneAPI than in the case of CUDA. Header-inclusion
issues often caused non-intuitive compilation errors if headers for the oneDPL
library were placed after oneMKL headers. Setting up an environment to compile
for the CUDA backend through CMake was not as easy as a build for fully-
supported architectures since we require separate CMake instructions and flags
when executing on Intel or Nvidia GPUs. As the CUDA-backend for SYCL support
continues, we expect that such features will be addressed.

##  5 Experimental Results

We conducted a series of experiments to evaluate the performance and
correctness of the oneAPI ports relative to the optimized CUDA implementations
of PAGANI and $m$-Cubes. We used a single node with a V100 Nvidia GPU and a
$2.4$ GHz Intel Xeon R Gold $6130$ CPU. We also used the Devcloud environment
to verify that the DPC++ implementations were portable and could be executed
on a P630 Intel GPU. Due to the V100 GPU having significantly more computing
cores than the P630, we do not make any performance comparisons between the
two GPUs. Instead, we focus on the attainable performance of DPC++ on Nvidia
hardware.

When executing the CUDA implementations, we used gcc $8.5$ and CUDA $11.6$.
For the CUDA-backend execution, we used the same environment but compiled with
clang $15$, an inline threshold of $10000$, and the following compilation
flags: “-fsycl -fsycl-targets=nvptx64-nvidia-cuda -Xsycl-target-backend –cuda-
gpu-arch=sm_70”. We verified the correctness of our ports, by comparing the
results on both the Nvidia (V100) and Intel (P630) GPUs, to the results
generated by the CUDA originals on a V100 GPU.

In terms of evaluating performance, we chose the same benchmark integrands
originally used to evaluate PAGANI and $m$-Cubes in [22] and [23]. These
functions belong to separate integrand families with features that make
accurate estimation challenging. We list those integrands in equations $1$ to
$6$. All experiments use the same integration bounds $(0,1)$ on each
dimensional axis. Similar to [22] and [23], we perform multiple experiments
per integrand.

We deviate from [22] and [23] in that we do not execute the PAGANI and
$m$-Cubes methods in their entirety. Instead, we execute their main kernels
pagani-kernel and mcubes-kernel, which is where more than $90\%$ of execution
is spent. With this approach, we can evaluate the effectiveness of each
programming model in terms of offloading workloads to the device. It allows us
to separate kernel evaluation from memory management operations (allocations,
copies, etc.) and library usage. This comparison of custom kernel
implementations is a better indicator of performance implications when porting
CUDA codes to DPC++.

$f_{1,d}\left(x\right)=\cos\left(\sum_{i=1}^{d}i\,x_{i}\right)$ (1)

$f_{2,d}\left(x\right)=\prod_{i=1}^{d}\left(\frac{1}{50^{2}}+\left(x_{i}-1/2\right)^{2}\right)^{-1}$
(2)

$f_{3,d}\left(x\right)=\left(1+\sum_{i=1}^{d}i\,x_{i}\right)^{-d-1}$ (3)

$f_{4,d}\left(x\right)=\exp\left(-625\sum_{i=1}^{d}\left(x_{i}-1/2\right)^{2}\right)$
(4)

$f_{5,d}\left(x\right)=\exp\left(-10\sum_{i=1}^{d}|x_{i}-1/2|\right)$ (5)

$f_{6,d}\left(x\right)=\begin{cases}\exp\left(\sum_{i=1}^{d}\left(i+4\right)x_{i}\right)&\text{
if }x_{i}<\left(3+i\right)/10\\\ 0&\text{otherwise}\end{cases}$ (6)

###  5.1 Offloading Mathematical Computations to Kernels

A critical stage in pagani-kernel and m-Cubes-kernel is the invocation of the
integrand at various $d$-dimensional points. Integrands with trigonometric or
exponential functions and table look-ups will have larger execution times
compared to other simple integrands that only contain basic mathematical
operations. To attain satisfactory performance, both the invocation of the
integrand functions and the remaining operations within the kernels must
achieve comparable performance to the CUDA implementation.

We tested the efficiency of the integrand oneAPI implementations with a simple
kernel that performs a series of invocations on many $d$-dimensional points.
The points are randomly generated on the host and then copied to device
memory. Each thread invokes the integrand serially $1$ million times and
writes its accumulated results to global memory. Writing the results prevents
the NVCC and Clang compilers from disregarding the integrand computations due
to optimization.

We first tested simple integrands that contained only a particular function
such as _sin, pow, powf, sycl::exp, sycl::pow, sycl::pown_. We invoked these
mathematical functions with $d$ arguments that comprise each $d$-dimensional
point. We did not use fast-math flags as accuracy is critical in numerical
integration. We observed small but consistent penalties of at most $2\%$ when
invoking the _power_ and _exponential_ functions. On the contrary,
trigonometric functions are approximately $40\%$ slower on the CUDA backend.

We performed the same experiment on the six benchmark integrands for
dimensions $5$ to $8$. We summarize the results in Table 1. The timings in
CUDA and oneAPI columns are the means of 10 kernel executions per integrand.
The ratio of those timings shows that the oneAPI version is at most $4\%$
slower. The largest penalty is observed in the $f1$ integrand which makes use
of the $cos$ function. The remaining integrands only make use of _exponential_
and _power_ functions and yield small penalties.

These experiments on the execution time of the integrand invocations
demonstrate that the user-defined computations do not display significant
performance penalties. The one exception is the extended use of trigonometric
functions. None of the benchmark integrands make extended use of trigonometric
functions ($f1$ has one call to $cos$ per invocation). As such, we do not
expect any slowdown larger than $5\%$ in either PAGANI or $m$-Cubes to be
attributed to the integrand implementations.

Table 1: mean ($\mu$) and standard deviation ($\sigma$) of execution times for invoking $5-8D$ benchmark integrands id |  $\mu$ CUDA (ms) |  $\mu$ oneAPI (ms) |  $\sigma$ CUDA |  $\sigma$ oneAPI |  $\frac{\mu\;oneAPI}{\mu\;CUDA}$
---|---|---|---|---|---
f1 |  1866.4 |  1952.4 |  13.3 |  21.4 |  1.04
f2 |  8413.9 |  8487.3 |  5012.5 |  5042.9 |  1.009
f3 |  1812.4 |  1828.3 |  18.5 |  27.1 |  1.009
f4 |  11416.1 |  11410.1 |  2184.9 |  2148.1 |  0.99
f5 |  634.3 |  654.4 |  73.5 |  67.3 |  1.03
f6 |  300.4 |  300.8 |  32.05 |  32.6 |  1.001

###  5.2 Benchmark Integrands Performance Comparison

Another set of experiments involved the invocation of the pagani-kernel and
mcubes-kernel on the benchmark integrands. To address different degrees of
computational intensity, we vary the number of thread-blocks used to launch
the kernels. For the mcubes-kernel, we achieve this effect by varying the
required number of samples per iteration in the range ($1e8,3e9$). This leads
to different block sizes per kernel. For pagani-kernel, the number of thread
blocks corresponds to the number of regions being processed. We perform high-
resolution uniform splits to generate region lists of different sizes and
supply them to the pagani-kernel for evaluation.

We report the penalty of using oneAPI for the benchmark integrands, in the
ratio columns of Tables 2 and 3. We used four thread-block sizes for each
integrand for the kernel executions. Each kernel configuration (number of
thread groups) was repeated $100$ times to provide a statistical mean and
standard deviation for the execution times.

Across our experiments, the average execution time ratio
($\frac{oneAPI}{CUDA}$) is in the range ($0-10\%$). The $f2$ and $f4$
integrands which make repeated use of the _power_ function display the largest
performance penalties for both PAGANI and $m$-Cubes. It is worth noting that
both $f2$ and $f4$ display the largest execution times among the benchmark
integrands for both integrators.

Table 2: $m$-Cubes: mean ($\mu$) and standard deviation ($\sigma$) of execution times for 8D benchmark integrands id |  $\mu$ CUDA (ms) |  $\mu$ oneAPI (ms) |  $\sigma$ CUDA |  $\sigma$ oneAPI |  $\frac{\mu\;oneAPI}{\mu\;CUDA}$
---|---|---|---|---|---
f1 |  286.7 |  286.7 |  2.1 |  0.9 |  1.0
f2 |  402.1 |  443.1 |  2.6 |  0.9 |  1.1
f3 |  284.5 |  285.8 |  1.6 |  1.4 |  1.0
f4 |  385.7 |  423.5 |  2.4 |  0.5 |  1.1
f5 |  284.3 |  285.9 |  2.1 |  1.7 |  1.0
f6 |  283.8 |  285.4 |  1.9 |  1.6 |  1.0
Table 3: PAGANI: mean ($\mu$) and standard deviation ($\sigma$) of execution times for 8D benchmark integrands id |  $\mu$ CUDA (ms) |  $\mu$ oneAPI (ms) |  $\sigma$ CUDA |  $\sigma$ oneAPI |  $\frac{\mu\;oneAPI}{\mu\;CUDA}$
---|---|---|---|---|---
f1 |  172.3 |  177.5 |  0.9 |  1.2 |  1.02
f2 |  1500.4 |  1651.0 |  0.3 |  2.1 |  1.1
f3 |  286.4 |  290.7 |  0.8 |  0.4 |  1.01
f4 |  1434.7 |  1524.9 |  0.4 |  1.9 |  1.06
f5 |  166.5 |  170.7 |  0.6 |  0.4 |  1.03
f6 |  136.8 |  139.4 |  0.4 |  0.2 |  1.02

###  5.3 Simple Integrands Performance Comparison

In addition to the benchmark integrands, we also evaluate integrands that only
perform a summation of the arguments ($\sum_{i=1}^{d}x_{i}$) where $d$ is the
number of dimensions. This avoids any bias in the comparison by avoiding
mathematical functions that could either call different implementations, cause
differences in register usage or lead to different optimizations. The ratios
in Tables 4 and 5, display timings on addition integrands for dimensions five
to eight. Once more, we observe penalties smaller than $10\%$ and for both
integrators these penalties decrease on higher dimensionalities.

Table 4: $m$-Cubes: mean ($\mu$) and standard deviation ($\sigma$) of execution times for addition integrands ($\sum_{i=1}^{d}x_{i}$) id |  $\mu$ CUDA (ms) |  $\mu$ oneAPI (ms) |  $\sigma$ CUDA |  $\sigma$ oneAPI |  $\frac{\mu\;oneAPI}{\mu\;CUDA}$
---|---|---|---|---|---
5D |  206.1 |  214.5 |  2.1 |  1.7 |  1.04
6D |  214.1 |  217.2 |  2.2 |  1.0 |  1.01
7D |  234.1 |  235.2 |  1.8 |  0.9 |  1.005
8D |  284.7 |  285.7 |  1.9 |  1.9 |  1.005
Table 5: PAGANI: mean ($\mu$) and standard deviation ($\sigma$) of execution times for addition integrands ($\sum_{i=1}^{d}x_{i}$) id |  CUDA (ms) |  oneAPI (ms) |  Std. CUDA |  Std. oneAPI |  $\frac{oneAPI}{CUDA}$
---|---|---|---|---|---
5D |  1.5 |  1.7 |  0.05 |  0.06 |  1.1
6D |  24.8 |  26.7 |  0.3 |  1.4 |  1.1
7D |  129.8 |  131.6 |  0.7 |  0.2 |  1.01
8D |  137.4 |  137.6 |  1.3 |  1.0 |  1.001

###  5.4 Factors Limiting Performance

Both pagani-kernel and mcubes-kernel, are compute bound, performing thousands
of computations for each byte of accessed memory. For compute-bound kernels,
the number of registers per thread is a factor limiting the number of
concurrent threads that can be executed; the amount of shared memory and
registers per thread limit warp/work-group occupancy, which in turn degrades
performance.

In most cases, the oneAPI implementations assigned more registers to each
thread compared to their CUDA equivalents. We illustrate the magnitude of this
difference in registers per thread in Figures 2 and 2. We observe the largest
difference in integrands $f2$ and $f4$, which make extended use of the _power
function_. It is the same functions that display the two largest execution
time penalties for the benchmark integrands in Tables 2 and 3.

We observe a similar pattern on the simple addition integrands (Table 4 and
5). In those cases, there are no mathematical functions (_pow_ , _exp_ , etc.)
and the integrands only perform a summation. The difference in registers
decreases on higher dimensions, leading to degraded performance on low
dimensions. This is evident in tables 4 and 5 where higher-dimensional
integrands have smaller values in the $\frac{oneAPI}{CUDA}$ column. The same
pattern is observed for the benchmark integrands, where the high dimensional
versions perform better than the low dimension equivalents. It can be seen in
Figure 2, that this effect is more prominent in $m$-Cubes, since it displays a
larger deviation across all dimensions. These observations lead us to believe
that register difference and its effect on occupancy is the main reason behind
the performance degradation.

Figure 1: Register difference on simple addition integrands
($\sum_{i=1}^{d}x_{i}$). The y-axis displays the number of additional
registers per thread in the DPC++ implementation.

Figure 2: Register difference on the benchmark integrands. The y-axis displays
the number of additional registers per thread in the DPC++ implementation.

##  6 Conclusion

We presented our experience of porting two numerical integration
implementations, PAGANI and $m$-Cubes, from CUDA to DPC++ . We utilized
Intel’s DPCT to automate the conversion process from CUDA to DPC++ and
successfully attained the capability to execute the same implementation on
both Intel and Nvidia GPUs. We experimented with various workloads consisting
of different mathematical functions. We found that the assigned registers per
thread can deviate in oneAPI and CUDA codes. This affects occupancy which in
turn can negatively impact performance, particularly in compute-bound kernels.
We faced additional challenges with mapping library calls to oneAPI
equivalents, matching compiler optimizations of NVCC with Clang, and using
build and testing libraries like CMake and Catch2. We addressed those
challenges and demonstrated that the performance penalty of using oneAPI ports
instead of optimized CUDA implementations can be limited to $10\%$ on Nvidia
GPUs. Additionally, numerous cases exhibited comparable performance to the
original CUDA implementations, with execution time differences in the $1-2\%$
range. We compared oneAPI and CUDA implementations on the same Nvidia V100
GPU. We were able to execute on an Intel P630 GPU but we did not compare these
timings with those on the V100 GPU due their significant difference in
computing power. In the future, we plan to execute on the high end Intel Ponte
Vecchio GPU and compare performance metrics with Nvidia high end GPUs such as
A100.

The vast array of libraries, ease of portability, and small margin of
performance degradation, make oneAPI an appropriate software solution for the
use case of numerical integration.

##  7 Acknowledgements

The authors would like to thank Intel Corporation and Codeplay for providing
technical support in the conversion process. The authors are also grateful for
the support of the Intel oneAPI Academic Center of Excellence at Old Dominion
University.

Work supported by the Fermi National Accelerator Laboratory, managed and
operated by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359
with the U.S. Department of Energy. The U.S. Government retains and the
publisher, by accepting the article for publication, acknowledges that the
U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide
license to publish or reproduce the published form of this manuscript, or
allow others to do so, for U.S. Government purposes. FERMILAB-
CONF-23-007-LDRD-SCD.

We acknowledge the support of Jefferson Lab grant to Old Dominion University
16-347. Authored by Jefferson Science Associates, LLC under U.S. DOE Contract
No. DE-AC05-06OR23177 and DE-AC02- 06CH11357.

## References

  * [1] https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5.0.pdf
  * [2] Argonne leadership computing facility, https://www.alcf.anl.gov/support-center/aurora/sycl-and-dpc-aurora#:~:text=DPC%2B%2B%20(Data%20Parallel%20C,versions%20of%20the%20SYCL%20language
  * [3] Computecpp™ community edition, https://developer.codeplay.com/products/computecpp/ce/2.11.0/guides/#computecpp
  * [4] Migrate cuda* to dpc++ code: Intel® dpc++ compatibility tool, https://www.intel.com/content/www/us/en/developer/tools/oneapi/dpc-compatibility-tool.html#gs.lx007q
  * [5] What is openacc?, https://www.openacc.org/
  * [6] et al., G.: Dark energy survey year 3 results: Redshift calibration of the maglim lens sample from the combination of sompz and clustering and its impact on cosmology (2022)
  * [7] Arumugam, K., Godunov, A., Ranjan, D., Terzic, B., Zubair, M.: A memory efficient algorithm for adaptive multidimensional integration with multiple gpus. In: 20th Annual International Conference on High Performance Computing. pp. 169–175. IEEE (2013)
  * [8] Ashbaugh, B., Bader, A., Brodman, J., Hammond, J., Kinsner, M., Pennycook, J., Schulz, R., Sewall, J.: Data parallel c++: Enhancing sycl through extensions for productivity and performance. In: Proceedings of the International Workshop on OpenCL. IWOCL ’20, Association for Computing Machinery, New York, NY, USA (2020). https://doi.org/10.1145/3388333.3388653, https://doi.org/10.1145/3388333.3388653
  * [9] Bridle, S., Dodelson, S., Jennings, E., Kowalkowski, J., Manzotti, A., Paterno, M., Rudd, D., Sehrish, S., Zuntz, J.: Cosmosis: a system for mc parameter estimation. Journal of Physics: Conference Series 664(7), 072036 (dec 2015). https://doi.org/10.1088/1742-6596/664/7/072036, https://dx.doi.org/10.1088/1742-6596/664/7/072036
  * [10] Carter Edwards, H., Trott, C.R., Sunderland, D.: Kokkos: Enabling manycore performance portability through polymorphic memory access patterns. Journal of parallel and distributed computing 74(12), 3202–3216 (2014)
  * [11] Castaño, G., Faqir-Rhazoui, Y., García, C., Prieto-Matías, M.: Evaluation of intel’s dpc++ compatibility tool in heterogeneous computing. Journal of Parallel and Distributed Computing 165, 120–129 (2022). https://doi.org/https://doi.org/10.1016/j.jpdc.2022.03.017, https://www.sciencedirect.com/science/article/pii/S0743731522000727
  * [12] Christgau, S., Steinke, T.: Porting a legacy cuda stencil code to oneapi. In: 2020 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). pp. 359–367 (2020). https://doi.org/10.1109/IPDPSW50202.2020.00070
  * [13] Costanzo, M., Rucci, E., García-Sánchez, C., Naiouf, M., Prieto-Matías, M.: Migrating cuda to oneapi: A smith-waterman case study. In: Bioinformatics and Biomedical Engineering, pp. 103–116. Lecture Notes in Computer Science, Springer International Publishing, Cham (2022)
  * [14] Costanzo, M., Rucci, E., Sanchez, C.G., Naiouf, M.: Early experiences migrating cuda codes to oneapi (2021)
  * [15] Costanzo, M., Rucci, E., Sánchez, C.G., Naiouf, M., Prieto-Matías, M.: Assessing opportunities of sycl and intel oneapi for biological sequence alignment (2022)
  * [16] Doerfert, J., Jasper, M., Huber, J., Abdelaal, K., Georgakoudis, G., Scogland, T., Parasyris, K.: Breaking the vendor lock-performance portable programming through openmp as target independent runtime layer. Tech. rep., Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States) (2022)
  * [17] Haseeb, M., Ding, N., Deslippe, J., Awan, M.: Evaluating performance and portability of a core bioinformatics kernel on multiple vendor gpus. In: 2021 International Workshop on Performance, Portability and Productivity in HPC (P3HPC). pp. 68–78 (2021). https://doi.org/10.1109/P3HPC54578.2021.00010
  * [18] Jin, Z., Vetter, J.: Evaluating cuda portability with hipcl and dpct. In: 2021 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). pp. 371–376 (2021). https://doi.org/10.1109/IPDPSW52791.2021.00065
  * [19] Kanzaki, J.: Monte carlo integration on gpu. The European physical journal. C, Particles and fields 71(2), 1–7 (2011)
  * [20] Peter Lepage, G.: A new algorithm for adaptive multidimensional integration. Journal of Computational Physics 27(2), 192–203 (1978). https://doi.org/https://doi.org/10.1016/0021-9991(78)90004-9, https://www.sciencedirect.com/science/article/pii/0021999178900049
  * [21] Ranjan, N., Terzić, B., Krafft, G., Petrillo, V., Drebot, I., Serafini, L.: Simulation of inverse compton scattering and its implications on the scattered linewidth. Physical review. Accelerators and beams 21(3), 030701 (2018)
  * [22] Sakiotis, I., Arumugam, K., Paterno, M., Ranjan, D., Terzić, B., Zubair, M.: PAGANI: A Parallel Adaptive GPU Algorithm for Numerical Integration. Association for Computing Machinery, New York, NY, USA (2021), https://doi.org/10.1145/3458817.3476198
  * [23] Sakiotis, I., Arumugam, K., Paterno, M., Ranjan, D., Terzić, B., Zubair, M.: m-cubes: An efficient and portable implementation of multi-dimensional integration for gpus. In: High Performance Computing, pp. 192–209. Lecture Notes in Computer Science, Springer International Publishing, Cham (2022)
  * [24] Stone, J.E., Gohara, D., Shi, G.: Opencl: A parallel programming standard for heterogeneous computing systems. Computing in science & engineering 12(3), 66–73 (2010)
  * [25] Su, C.L., Chen, P.Y., Lan, C.C., Huang, L.S., Wu, K.H.: Overview and comparison of opencl and cuda technology for gpgpu. In: 2012 IEEE Asia Pacific Conference on Circuits and Systems. pp. 448–451 (2012). https://doi.org/10.1109/APCCAS.2012.6419068
  * [26] Tsai, Y.M., Cojean, T., Anzt, H.: Porting sparse linear algebra to intel gpus. In: Euro-Par 2021: Parallel Processing Workshops, pp. 57–68. Lecture Notes in Computer Science, Springer International Publishing, Cham (2022)
  * [27] Volokitin, V., Bashinov, A., Efimenko, E., Gonoskov, A., Meyerov, I.: High performance implementation of boris particle pusher on dpc++. a first look at oneapi. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), pp. 288–300. Lecture Notes in Computer Science, Springer International Publishing, Cham (2021)
  * [28] Wong, M., Liber, N., Bassini, S., Richards, A., Butler, M., McVeigh, J., Cook, B., Sugimoto, H., Cordoba, C., Fahringer, T., et al.: Sycl - c++ single-source heterogeneous programming for acceleration offload (Jan 2014), https://www.khronos.org/sycl/
  * [29] Wu, H.Z., Zhang, J.J., Pang, L.G., Wang, Q.: Zmcintegral: A package for multi-dimensional monte carlo integration on multi-gpus. Computer Physics Communications 248, 106962 (2020). https://doi.org/https://doi.org/10.1016/j.cpc.2019.106962, https://www.sciencedirect.com/science/article/pii/S0010465519303121
  * [30] Zenker, E., Worpitz, B., Widera, R., Huebl, A., Juckeland, G., Knüpfer, A., Nagel, W.E., Bussmann, M.: Alpaka - an abstraction library for parallel kernel acceleration. In: arXiv.org. Cornell University Library, arXiv.org, Ithaca (2016)