Block-space GPU Mapping for Embedded Sierpi\'nski Gasket Fractals

TL;DR

This paper introduces an efficient GPU thread mapping method for embedding Sierpiński gasket fractals in Euclidean space, significantly improving parallel space utilization and performance over traditional bounding-box approaches.

Contribution

A novel block-space mapping algorithm for GPU fractal embedding that operates in logarithmic time and reduces thread usage compared to bounding-box methods.

Findings

Achieves $ ext{O}( ext{log}_2 ext{log}_2 n)$ mapping time.

Provides up to 10x speedup for large fractals.

Reduces thread count to $ ext{O}(n^ ext{H})$ with H≈1.58.

Abstract

This work studies the problem of GPU thread mapping for a Sierpi\'nski gasket fractal embedded in a discrete Euclidean space of $n \times n$. A block-space map $\lambda: \mathbb{Z}_{\mathbb{E}}^{2} \mapsto \mathbb{Z}_{\mathbb{F}}^{2}$ is proposed, from Euclidean parallel space $\mathbb{E}$ to embedded fractal space $\mathbb{F}$, that maps in $\mathcal{O}(\log_2 \log_2(n))$ time and uses no more than $\mathcal{O}(n^\mathbb{H})$ threads with $\mathbb{H} \approx 1.58...$ being the Hausdorff dimension, making it parallel space efficient. When compared to a bounding-box map, $\lambda(\omega)$ offers a sub-exponential improvement in parallel space and a monotonically increasing speedup once $n > n_0$. Experimental performance tests show that in practice $\lambda(\omega)$ can produce performance improvement at any block-size once $n > n_0 = 2^8$, reaching approximately $10\times$ of speedup for $n=2^{16}$ under optimal block configurations.

Tables

Table 1. TABLE I: Hardware used for performance tests.

Device	Model
GPU	Titan-X Pascal, GP102, 3584 cores, 12GB
CPU	Intel i7-6950X 10-core Broadwell
RAM	128GB DDR4 2400MHz

Equations

$$\mathcal{V}({\mathcal{F}}_n^{k,s})=\sum _{i=1}^k\mathcal{V}({\mathcal{F}}_{sn}^{k,s})$$

$$\mathcal{V}({\mathcal{F}}_n^{k,s})=k^r$$

$$\mathcal{V}({\mathcal{F}}_n^{3,\frac{1}{2}})=n^{\mathcal{H}={\log }_2(3)}$$

$$\beta_{\mu}(\omega)=\Big{(}\frac{\omega_{x}((\mu+1)\mod 2)+\omega_{y}(\mu\mod 2)}{3^{\lceil\frac{\mu}{2}\rceil-1}}\Big{)}\mod 3$$

$$\displaystyle\Delta_{\mu}^{x}=\Big{\lfloor}\frac{\beta_{\mu}}{2}\Big{\rfloor},\ \ \ \Delta_{\mu}^{y}=\beta_{\mu}-\Big{\lfloor}\frac{\beta_{\mu}}{2}\Big{\rfloor}$$

$${\tau }_x^{\mu }={\Delta }_{\mu }^x{2}^{\mu -1},$$

$${\tau }_y^{\mu }={\Delta }_{\mu }^y{2}^{\mu -1}$$

$$\lambda (\omega )$$

$${\lambda }_x(\omega )$$

$${\lambda }_y(\omega )$$

$$S_{\lambda (\omega )}$$

$$\underset{n\to \infty }{\lim }{S}_{\lambda (\omega )}$$

$$=\underset{n\to \infty }{\lim }\frac{(2-\mathcal{H})n^{1-\mathcal{H}}}{\frac{1}{n{\log }_2(n))}}$$

$$=\infty$$

Full text

I Introduction

Fractals can be described as self-similar structures [14] where a similar111Depending on which fractal, the term similar can refer to exactly similar or quasi similar. geometrical pattern is found at all scales. Several natural phenomena produce fractal patterns that obey a self-similar structure [13]. Phenomena such as plant and tree growth [21, 22], terrain formation [15, 23], molecular dynamics [26], snowflake crystallization [9], blood vessels [5], morphological features of living organisms [27], among many others, display a fractal design where self-similarity is a critical feature for modeling its geometrical structure.

One well known fractal is the Sierpiński Gasket, or Sierpiński Triangle, described by Waclaw Sierpiński in 1915. Despite being over than a century old, the Sierpiński gasket remains a relevant subject as it is present in different fields such as the construction of antennas [2, 24], cellular automata [20, 28], molecular DNA self-organization [10, 26], self-assembly theory [4, 12] and phase transitions on fractal spin lattices [6, 7, 16], among others, because of its special properties. For fractal simulations such as in Cellular Automata or Monte Carlo simulation on spin models, the Sierpiński gasket is usually discretized. Is this form, the fractal is defined as a self-similar structure where level $r+1$ is composed of three repetitions of level $r$ at $1/2$ the scale, as shown in Figure 1.

Applications that involve graphical representations or data-parallel simulations with nearest neighbors interactions may benefit if the storage of the structure preserves spatial locality in memory space, i.e, the memory locations $(x\pm 1,y\pm 1)$ define a neighborhood in the actual fractal as well. One way to achieve this is to embed the fractal in a Euclidean space of $n\times n$ as shown in Figure 2.

For applications that can operate in an embedded fractal domain, computations will usually consist of accessing all the elements of the fractal and eventually perform arithmetic/logic operations that involve the data-element itself and possibly its nearest neighbors. Eventually, when the fractal is large enough to the point of containing hundreds of thousands of elements, a sequential computation can take an excessive amount of time for the practical requirements of the field. In these situations GPU computing becomes an attractive tool for accelerating the task [19].

For every GPU computation there is a stage where threads are mapped from parallel to data space. A map, defined as $f:\mathbb{Z}^{k}\rightarrow\mathbb{Z}^{m}$, transforms each $k$-dimensional point $x=(x_{1},x_{2},...,x_{k})$ in parallel space $P^{k}$ into a unique $m$-dimensional point $f(x)=(y_{1},y_{2},\cdots,y_{m})$ in data space $D^{m}$. GPU parallel spaces are defined as orthotopes $\Pi^{k}\in P^{k}$ in $k=1,2,3$ dimensions. A known way of mapping threads is to use the bounding-box (BB) approach, that builds an orthotope $\Pi^{k}$ sufficiently large to cover the whole data space and threads are mapped using the identity $f(\omega)=\omega$. Such map is highly convenient and efficient for the class of problems where data space is also defined by an orthotope; such as vectors, tables, matrices and box-shaped volumes. But for an embedded fractal such as the Sierpiński gasket, this approach is no longer efficient in terms of parallel space as many threads fall outside the domain, introducing a performance penalty to the execution time (see Figure 3).

Two research questions arise from this GPU efficiency problem on the embedded Sierpiński gasket; The first question: Is there any parallel-space efficient function, namely $\lambda(\omega)$, that can use asymptotically the same number of threads as data elements in the fractal and map blocks properly onto the embedded Sierpiński gasket? (see Figure 4).

It is important to note that Question 1 asks for a block-space map and not a thread-space one. The change from thread-space to block-space allows coalesced memory to be preserved throughout the entire domain as thread organization is not compromised inside a block.

The second question relates to performance: Will the parallel-space improvement translate into a significant GPU performance improvement?

The present work focuses on these two questions and provides positive answers for both of them. A dedicated analysis is devoted to show that an alternating unrolling strategy allows to define a parallel-space efficient $\lambda(\omega)$ that only requires $\mathcal{O}(n^{\mathbb{H}})$ threads, with $\mathbb{H}\approx 1.58...$ being the Hausdorff dimension of the Sierpiński fractal. It addition, it is shown that by taking advantage of block-parallelism, $\lambda(\omega)$ becomes computable in $\mathcal{O}(\log_{2}\log_{2}(n))$ time which is fast enough to produce monotonically increasing speedup once $n>n_{0}$ with $n_{0}$ being a constant threshold value.

The rest of the manuscript presents related work (Section II), a formal definition and analysis of $\lambda(\omega)$ (Section III) and Section (IV) presents performance results. A discussion and comments on future work is found in Section V.

II Related Work

The following related works can be classified into two categories; (1) studies on efficient GPU mapping for triangular domains and (2) general studies on the structure of the discrete Sierpiński fractal.

One of the first works that explored the possibilities of improving the GPU mapping and stage was the research of Jung et. al. [11] whom proposed packed data structures for representing triangular and symmetric matrices with applications to LU and Cholesky decomposition [8]. The strategy is based on building a rectangular box strategy for accessing and storing a triangular matrix (upper or lower). Data structures become practically half the size with respect to classical methods based on the full matrix. The strategy was originally intended to modify the data space (i.e., the matrix), however one can apply the concept analogously to the parallel space.

Ries et. al. contributed with a parallel GPU method for the triangular matrix inversion [25]. The authors identified that the parallel space indeed can be improved by using a recursive partition of the grid, based on a divide and conquer strategy. The approach takes $O(\log_{2}(n))$ time by doing a balanced partition of the structure, from the orthogonal point to the diagonal.

Q. Avril et. al. proposed a GPU mapping function for collision detection based on the properties of the upper-triangular map [1]. The map is a thread-space function $u(x)\rightarrow(a,b)$, where $x$ is the linear index of a thread $t_{x}$ and the pair $(a,b)$ is a unique two-dimensional coordinate in the upper triangular matrix. Since the map works in thread space, the map is accurate only in the range $n\in[0,3000]$ of linear problem size.

Navarro, Hitschfeld and Bustos have proposed a block-space map function for $2$-simplices and $3$-simplices [18, 17], based on the solution of an $m$ order equation that is formulated from the linear enumeration of the discrete elements. The authors report performance improvement for $2$-simplices, and for the $3$-simplex case, the mapping technique is extended to the discrete orthogonal tetrahedron, where the parallel space usage can be $6\times$ more efficient. However the authors clarify that it is difficult to translate such space improvement into performance improvement, as the map requires the computation of several square and cubic roots that introduce a significant amount of overhead to the process. From the point of view of data-reorganization, a succinct blocked approach can be combined along with the block-space thread map, producing additional performance benefits with a sacrifice of $o(n^{3})$ extra memory.

Exploring the benefits of efficient GPU mapping onto the embedded Sierpiński gasket becomes an interesting topic of research since its geometry is no longer Euclidean as in the related works, but instead it is embedded in an Euclidean one. Finding a proper efficient $\lambda(\omega)$ would produce an asymptotic improvement in parallel space and a potential performance improvement that could eventually be exploited.

III Analysis and Formulation of $\lambda(\omega)$

This Section first analyzes two important properties of the discrete embedded Sierpiński gasket, which are helpful in the formulation of the map $\lambda(\omega)$.

III-A Analysis of the Space and Packing of $\mathcal{F}_{n}^{k,s}$

The discrete Sierpiński gasket belongs to a category of discrete fractals where its structure can be built by replicating $k$ instances of itself at a scale $s$ and placed with an arbitrary spatial organization. The notation $\mathcal{F}_{n}^{k,s}$ is introduced to denote such fractal, where $n\in\mathbb{N}$ is the linear size, $k\in\mathbb{N}$ the replication factor and $0<s<1\in\mathbb{R}$ the scaling factor in terms of reduction. Since fractals have a recursive self-similar structure, their volume $\mathcal{V}(\mathcal{F}_{n}^{k,s})$ may be expressed as

[TABLE]

with $\mathcal{V}(\mathcal{F}_{1}^{k,s})=1$ being the limit condition of the recursion. Given that the replication factor $k$ is fixed, and $n$ scales by factors of $s$, the volume may be simplified into

[TABLE]

where $r=\log_{1/s}(n)$ is defined as the scale level.

The Sierpiński gasket is a particular case where $k=3$, $s=1/2$ and $r=\log_{2}(n)$. Successive steps of the fractal produce the pattern early depicted in Figure 2. The following two lemmas are introduced to support the formulation of a map $\lambda(\omega)$.

Lemma 1

The space occupied by a discrete embedded Sierpiński gasket is in correspondence with the Hausdorff dimension of the infinite Sierpiński gasket.

Proof:

The space occupied by a Sierpiński gasket of linear size $n$ is $\mathcal{V}(\mathcal{F}_{n}^{3,\frac{1}{2}})=3^{r}$. Given that $r=\log_{2}(n)$ and $3^{log_{2}(x)}=2^{\log_{2}(3)\log_{2}(x)}$, the space expression can be rearranged into

[TABLE]

where the exponent $\mathcal{H}=\log_{2}(3)\approx 1.5849...$ is the Hausdorff dimension of the original infinite Sierpiński gasket. ∎

Lemma 2

A discrete Sierpiński gasket of level $r$ packs into a $2$-orthotope $\Pi^{2}$ of dimensions $3^{\lceil\frac{r}{2}\rceil}\times 3^{\lfloor\frac{r}{2}\rfloor}$.

Proof:

Proof by induction on $r$:

• •

Base case: At scale $r=0$ the fractal has a space of $\mathcal{V}(\mathcal{F}_{1}^{3,\frac{1}{2}})=1$ element that packs into a regular $2$-orthotope of $1\times 1=3^{\lceil\frac{0}{2}\rceil}\times 3^{\lfloor\frac{0}{2}\rfloor}$ satisfying $3^{\lceil\frac{r}{2}\rceil}\times 3^{\lfloor\frac{r}{2}\rfloor}$.

• •

Induction step: It is assumed that the orthotope for $r=k$ is quasi-regular. If $k$ is even, the packing for $k+1$ will triple the horizontal dimension of the $2$-orthotope. If $k$ is odd, the packing for $k+1$ will triple the vertical dimension of the $2$-orthotope. Since even and odd must alternate, the dimensions of the packed $2$-orthotope for $k+1$ can only be $3\cdot 3^{\lceil{\frac{k}{2}}\rceil}\times 3^{\lfloor{\frac{k}{2}}\rfloor}$ or $3\cdot 3^{\lceil{\frac{k}{2}}\rceil}\times 3\cdot 3^{\lfloor{\frac{k}{2}}\rfloor}$, which is regular or quasi-regular, respectively.

∎

The packing process of Lemma (2) is illustrated in Figure 5, where the packing steps are in correspondence with the scale levels shown in Figure 5.

III-B Changing from Thread-space to Block-space

An important aspect to consider is at which level the parallel space will be mapped. Two approaches are possible; (1) thread-space mapping and (2) block-space mapping. For first one, $\lambda(\omega)$ defines $\omega$ as a unique thread location in parallel space. For the second approach, $\lambda(\omega)$ defines $\omega$ as a block coordinate in which several threads are contained. The block-space approach has three important advantages over thread-space mapping. First, in block-space, the fractal becomes a simplified version of the original, requiring less elements to be mapped. Second, since the fractal is a simplified version of itself, it is possible to work on higher sizes of $n$ before the CUDA grid maximum dimensions are reached. Third, the block-space approach allows the possibility for threads inside a block to preserve locality, which is essential for coalesced memory accesses on data.

The block-space map $\lambda(\omega)$ is now introduced, where $\omega=(\omega_{x},\omega_{y})$ denotes a two-dimensional coordinate of a block of constant size $|B|=\rho_{x}\times\rho_{y}$ threads. The change from thread-space to block-space means that blocks are mapped to a simplified version of the fractal of linear size $n_{b}=n/b$, as shown in Figure 6.

It is important to clarify that the extra green regions visible in Figure 6 do not necessarily mean unused threads. The stage of local thread mapping is detailed later in this Section, where three approaches can be used.

The formulation of $\lambda(\omega)$ continues in block-space with $\rho_{x}=\rho_{y}$ to simplify the analysis, with $n_{b}$ the new simplified linear size of the fractal with the origin $(0,0)$ located at the top-left corner for both the parallel and embedded fractal spaces, and with $y$ increasing downwards.

III-C Formulation of $\lambda(\omega)$

The function $\lambda:\mathbb{Z}_{\mathbb{E}}^{2}\mapsto\mathbb{Z}_{\mathbb{F}}^{2}$ is introduced as a mapping of block coordinates $\omega$ in parallel-space onto block coordinates in the embedded space. The intuition behind is an unrolling process applied in parallel to each $\omega\in\Pi^{2}$ through all the scale levels. At each level, different $x,y$ offsets are accumulated to form the final $(\lambda_{x}(\omega),\lambda_{y}(\omega))$ coordinate in the fractal.

Theorem 1

There exists a block-space parallel-space efficient $\lambda(\omega)$ that can map blocks in $\mathcal{O}(\log_{2}\log_{2}(n_{b}))$ time using $|B|=\mathcal{\theta}(\frac{\log_{2}(n_{b})}{\log_{2}\log_{2}(n_{b})})$ threads per block.

Proof:

By construction: let $r_{b}=2^{n_{b}}$ be the block-space scale level of the fractal, $\Pi^{2}$ the $2$-orthotope of $3^{\lceil\frac{r_{b}}{2}\rceil}\times 3^{\lfloor\frac{r_{b}}{2}\rfloor}$ blocks that maps onto the discrete Sierpiński gasket $\mathcal{F}_{n_{b}}^{3,\frac{1}{2}}$, with each block having $\rho_{x}\times\rho_{y}$ threads. By Lemma (1), $\Pi^{2}$ is parallel-space efficient in block-space.

A helper index function $\beta_{\mu}(\omega)$ is defined as

[TABLE]

to produce an index in the range $\beta_{\mu}(\omega)=0,1,2$ that identifies, within scale level $\mu\in[0..r_{b}]$, which of the three regions of the fractal does block $\omega$ belongs to. For even $\mu$, $\beta_{\mu}(\omega)$ acts on $\omega_{x}$. For odd $\mu$, it acts on $\omega_{y}$. Regions are sorted as [math] for top, $1$ for middle and $2$ for right (see Figure 1 for visual reference).

Having the $\beta_{\mu}(\omega)$ index, the weight functions

[TABLE]

compute the offset weights, [math] or $1$, for each of the $x$ and $y$ directions at scale level $\mu$. The value of the offset corresponds to $2^{\mu-1}$, which is the linear size of each region at scale level $\mu$. For a given $\mu$, the combination of the weight functions with the offset produce partial coordinates of the form

[TABLE]

that contribute to the final mapped coordinate. The summation of all partial coordinates produce the map

[TABLE]

which can be computed in $\mathcal{O}(\log_{2}\log_{2}(n))$ time (i.e., $n_{b}\in\theta(n)$) using a parallel reduction with the threads contained in the $\omega$ block. Finally, by Brent’s Theorem [3], $|B|=\mathcal{\theta}(\frac{\log_{2}(n)}{\log_{2}\log_{2}(n)})$ threads are sufficient for a block of threads to reduce efficiently in parallel. ∎

Theorem 2

$\lambda(\omega)$* requires asymptotically less work than the bounding-box approach.*

Proof:

The asymptotic work improvement factor of $\lambda(\omega)$ with respect to the bounding-box approach is the quotient of the costs of mapping all blocks using their corresponding $\Pi^{2}$ structures with the consideration $n_{b}\in\theta(n)$

[TABLE]

where $\Pi^{2}_{BB}$ and $\Pi_{\lambda(\omega)}^{2}$ are the parallel-spaces for the bounding-box and $\lambda(\omega)$ approaches, respectively. The parallel-space of $\Pi_{BB}^{2}$ corresponds to the Euclidean box of $n_{b}\times n_{b}$ blocks, and the parallel-space of $\Pi_{\lambda(\omega)}^{2}$ is $\mathcal{O}(n^{\mathcal{H}})$ by Lemma (1). Applying the limit $n\to\infty$

[TABLE]

∎

The importance of Theorem (2) is that it guarantees the existence of a $n>n_{0}$ where $\lambda(\omega)$ will start becoming each time faster than the bounding box approach. A theoretical curve is presented in Figure 7.

From the plot, one can note that space improvement is clear in the log-log scale. For the time improvement, the improvement decreases until $n_{0}\sim 7$, which is the point where it becomes a monotonically increasing function. In practice, constants could have an effect that would push $n_{0}$ further to the right.

III-D Intra-Block Mapping

Once $\lambda(\omega)$ maps a block $\omega$, all the $\rho_{x}\times\rho_{y}$ threads contained have a shared reference coordinate that is available to use for computing their individual location in the fractal. This phase of organizing the threads within a block is referred here as Intra-Block Mapping, and this subsection describes three possible approaches to accomplish this.

III-D1 Further Unrolling

In this approach threads inside their mapped block use the same $\lambda(\omega)$ map but applied to each thread. By Theorem (1), the Intra-block map is parallel-space efficient and the mapping time becomes $\mathcal{O}(\log_{2}\log_{2}(|B|))\in\mathcal{O}(1)$ as $\rho_{x},\rho_{y}$ are constant and do not grow with $n$.

III-D2 Shared Lookup Table

The second approach is to use a shared lookup table of $\rho_{x}\times\rho_{y}=\mathcal{O}(1)$ offset coordinates, available to any thread in any block. Mapping each thread would cost $\mathcal{O}(1)$ memory accesses and the extra memory introduced by the shared table is $\mathcal{O}(\rho_{x}\times\rho_{y})\in\mathcal{O}(1)$.

III-D3 Bounding Sub-boxes

The third approach consists of using small bounding-boxes in each block. This approach introduces a constant number of extra threads in each block, but allows each thread to be mapped just with $f(x)=x$ which costs $\mathcal{O}(1)$. In order to know if it is in the fractal or not, each thread evaluates if $t_{x}\&(n-1-t_{y})==0$ is true or not, respectively, with $\&$ being the bitwise AND operator.

Regardless of which Intra-block mapping approach is chosen, the final mapping time will not surpass the $\mathcal{O}(\log_{2}\log_{2}(n))$ time, as the blocks have a constant size of threads. Nevertheless, it is worth considering that Further Unrolling introduces a constant cost in mapping time. The Shared Lookup Table approach introduces a constant cost in memory and finally the Bounding Sub-boxes introduce a constant in the number of extra threads. Choosing one or another can depend on the specific application, i.e., to avoid competing with the application in the use of memory bandwidth or arithmetic operations.

IV Implementation and Performance Results

A CUDA implementation of $\lambda(\omega)$ is put under test to obtain the actual speedup for different values of $n$. The implementation uses the Bounding Sub-boxes approach to arrange threads inside each block, as it is the simplest in terms of implementation time. The parallel reduction per-block is performed using the shuffle instruction of CUDA, which allows efficient register-level communication among threads within the same warp.

The performance test consists of measuring the average time taken to write a constant value on all the elements of a Sierpiński gasket of scale level $r$, which is embedded in a $M_{n\times n}$ matrix filled with zeros. The configuration space is explored in the ranges $r=0..16$ and $\rho=1,2,4,8,16,32$ for the scale level and block-size respectively, in order to find the optimal setting that provides the best performance for both the bounding-box and $\lambda(\omega)$ approaches. The average performance measures are taken by averaging 100 sub-averages, each one being an average time of 10 consecutive synchronized kernel calls. The standard error for each mean was below $1\%$. The hardware for performance test is listed in Table I.

Figure 8 presents the speedup of $\lambda(\omega)$ over the bounding-box approach, as well as the running times for the two mapping techniques.

For values of $n<2^{8}$, one can note that only some curves offer speedup. Once $n>2^{8}$, the speedup begins to increase for all block-size configurations, reaching the higher values at $n=2^{16}=65536$, which was the highest problem size that fit in the GPU memory. An important aspect to note from the speedup curve is that for the largest possible block size, $|B|=\rho\times\rho=32\times 32$, the $\lambda(\omega)$ map runs the test approximately $6\times$ faster than the bounding-box approach. Furthermore, as blocks become smaller in $\rho$, the improvement increases dramatically, reaching up to $55\times$ of speedup.

The plot of the running times provides further insights on what configuration is the best suited for each mapping technique. By looking at the running times of the small block configurations, one can note that regardless of their high speedup, their running times are the slowest ones. For the bounding-box approach (BB) the best performance is obtained when the block-size is maximum. For $\lambda(\omega)$ the best performance is found when using a block of $|B|=16\times 16$ threads. If the curves of the best configuration for each implementation are considered, i.e., the ones with marked points on Figure 8, right, then the speedup provided by $\lambda(\omega)$ increasing further more reaching almost an order of magnitude. The other running time curves are still useful to visualize that as blocks become smaller, the value of $n_{0}$ where $\lambda(\omega)$ is more convenient moves closer to the origin, and vice versa, however the GPU with its current organization and architecture is not fully utilized by such small block configurations, thus leading to an inferior performance. Therefore, in practice large blocks should be utilized and by Theorem (2), beyond $n=2^{16}$ the speedup would keep increasing in favor of $\lambda(\omega)$.

V Discussion

The results obtained in this work have shown that a parallel-space efficient mapping function can lead to significant performance gains in applications that require processing an embedded Sierpiński gasket. The analysis and formulation of $\lambda(\omega)$ has provided three important results; (1) There exists a correspondence between a quasi-regular $2$-orthotope of discrete elements and the elements of the embedded Sierpiński gasket fractal. (2) Such correspondence from parallel to data space can be computed in $\mathcal{O}(\log_{2}\log_{2}(n))$ time using a block-space map that is based on efficient parallel reductions. (3) The total work of mapping the $2$-orthotope with $\lambda(\omega)$ is asymptotically smaller than the work generated by the bounding box approach, leading to a monotonically increasing speedup that is guaranteed to occur once a $n_{0}$ value is reached.

The experimental performance results confirm the theory as once the fractal reaches the linear size $n_{0}=2^{8}$, the speedup begins to increase monotonically for all block-sizes. Using the maximum block size of $32\times 32$ threads, $\lambda(\omega)$ reached up to $6\times$ of speedup. The maximum speedup obtained was approximately $55\times$ with blocks of $1$ thread, however such block configuration is not practical for the GPU architecture. Still, having measured the running times with different block configurations helped in understanding that reducing the block size only brings the $n_{0}$ value closer to the origin and increasing the block-size pushes it forward. Eventually, all block configurations can reach up to $55\times$ of speedup and beyond if the fractal is large enough.

The GPU map found for the embedded Sierpiński gasket may be adapted to work for other types of embedded fractals that follow a similar building scheme applying modifications to the helper, weight and offset functions. In order to obtain speedup, it is crucial to check if the overall work will be asymptotically smaller than the bounding-box approach.

Two important questions may be formulated from the results obtained in this work. The first one is: Can there exist a general $\lambda(\omega)$ that maps a family of embedded fractals who share the same building principle?, and the second question: can there exist a $\lambda(\omega)$ that maps in $\mathcal{O}(1)$ time using no more than $\mathcal{O}(1)$ extra memory?. Future research in these directions can provide important insights on the limits of efficient GPU computing for embedded fractal domains.

Acknowledgment

This project was supported by the research project FONDECYT No 3160182 from CONICYT, as well as by the Nvidia CUDA Research Center at the Department of Computer Science (DCC) from University of Chile.