2.1 The MiniCP solver and GPU

Since the late Eighties several efficient constraint solvers have been developed for Constraint (Logic) Programming. For the scope of this paper we have decided to use MiniCP (Michel et al. Reference MICHEL, SCHAUS and VAN HENTENRYCK2021), a relatively new and open-source solver devised to ease development of extensions. In particular, our extensions build on MiniCPP (Gentzel et al. Reference GENTZEL, MICHEL, VAN HOEVE and Simonis2020), a C++ implementation of MiniCP that can easily include CUDA C++ code (Tardivo et al. Reference TARDIVO, DOVIER, FORMISANO and MICHEL2024a; Reference TARDIVO, DOVIER, FORMISANO, MICHEL and PONTELLI2023; Reference TARDIVO and MICHEL2024b). New global constraints can be integrated in MiniCP by creating a class under /fzn-minicpp/global_constraints/ which extends the Constraint class. Such a class, which will serve as the propagators for the new global constraints, must include, alongside a constructor, a procedure, propagate(). The propagate() method implements the propagation algorithm for the new constraint. To provide a simple mechanism to use the GPU-accelerated propagators, we rely on the MiniZinc annotation introduced by (Tardivo et al. Reference TARDIVO, DOVIER, FORMISANO, MICHEL and PONTELLI2023). Specifically, a table constraint annotated with ::gpu (i.e., a constraint declared as: constraint table(…)::gpu; in the input model) is propagated using the GPU-accelerated implementation in place of the standard one.

4.1 The propagator $\textbf{CT}_\textbf{CU}^\textbf{u}$

In this propagator the execution flow of $\mathit{enforceGAC}()$ , closely matches its serial counterpart. First, the contents of $\mathit{\_currTable}$ and $\mathit{\_s\_val}$ are copied to the device, together with the variables’ domains. The kernels responsible for updating the current table, that is $\mathit{updateTableGPU}()$ and $\mathit{reduce}()$ , are subsequently launched, see Figure 4. The outcome of this step is the same mask computed by the original $\mathit{updateTable}()$ procedure, but it is written in GPU’s global memory. It is then copied back to the host and combined bitwise with $\mathit{\_currTable}$ . Finally, the filtering procedure on the host is invoked.

Fig 1. Overview of the parallel reduction done on $\_supportsT\_dev$ , where $\_s\_val=\{x_1,x_3,\ldots \}$ . Different colors indicate which words are accessed by the threads in the block. Arrows depict the direction of the parallel reductions and $\lceil \frac {t}{32} \rceil$ is a shorthand for $\_currTableSize$ (see Section 2).

The propagator $CT_{CU}^{u}$ uses some auxiliary structures that do not have a counterpart on the serial side. They are:

• • $\mathit{\_supportsT\_dev}$ : An array of integers storing the transposed support matrix. It is used by the kernel $\mathit{updateTableGPU}()$ for updating the current table. The transpose is exploited to make threads of the same warp access adjacent memory words.

• • $\mathit{\_vars\_dev}$ : An array of $ \lceil \mathit{\_supportSize/32}\rceil$ integers. It stores the concatenation of the bitmaps representing the different variable domains. It is the counterpart of an host side variable ( $\mathit{\_vars\_host}$ ) storing variable domains. The content of $\mathit{\_vars\_host}$ is copied to the device at each propagation.

• • $\mathit{\_tmpMasks}$ : An array of integers of size $\mathit{\_currTableSize\times n}$ . It is a linearized representation of the support matrix, in which each row (of size $\mathit{\_currTableSize}$ ) is associated to a variable. Namely, for each variable $x_i \in \_s\_val$ , the row $\mathit{\_tmpMasks[i]}$ stores the $\mathit{bitwise}$ - $\mathit{OR}$ of the rows $\mathit{\_supports}[x_i,a]$ for all $a \in dom(x_i)$ . Such a structure stores the result of the kernel $\mathit{updateTableGPU}()$ and is subsequently processed by a $\mathit{reduce}()$ kernel (see below).

Fig 2. Overview of the parallel reduction done on $\mathit{\_tmpMasks}$ by $\mathit{reduce}()$ . Different colors of the blocks indicate which words are accessed by the threads in the block. Arrows show the direction of the parallel reductions while $\lceil \frac {t}{32} \rceil$ is a shorthand for $\_currTableSize$ (see Section 2).

Fig 3. Overview of the data each block considers in the kernel $\mathit{filterDomainsGPU}()$ .

As concerns the kernel $\mathit{updateTableGPU}()$ , the grid launched has $\lceil \_currTableSize/4 \rceil \times |\_s\_val|$ blocks. That is, a block of 128 threads for every four (32-bit) words used by $\mathit{\_currTable}$ and for each changed variable is launched. Note that the division of work between blocks/threads (and thus between streaming multiprocessors of the GPU) is dynamic and scales with the size of the problem (i.e. the number of rows in the table as well as the scope size).

Each row of the grid is related to a variable reduced during the last iteration, that is the $i$ -th row of the grid is related to the $i$ -th variable in $\_s\_val$ , assume it to be $x_j$ . Each block in the row is responsible for updating four (32-bit) words of the table mask accessing the transposed supports related to $x_j$ . Each block of 128 threads can be seen as a small sub-grid of $32 \times 4$ threads, where each of the 32 rows considers a (different) slice of $dom(x_j)$ . Each row of this sub-grid computes, via a $\mathit{bitwise}$ - $\mathit{OR}$ , the partial mask that needs to be added to $\mathit{\_currTable}$ . Such a mask, computed in shared memory is then copied in global memory in the $i$ -th row of $\mathit{\_tmpMasks}$ (see Figure 1).

Once $\mathit{updateTableGPU}()$ is terminated, the kernel $\mathit{reduce}()$ is launched. It considers all and only the rows of $\_tmpMasks$ related to changed variables (i.e. those in $\_s\_val$ ) and performs a $\mathit{bitwise}$ - $\mathit{AND}$ operation along the columns. To this aim, a grid of $\mathit{\_currTableSize}$ blocks of 32 threads each is launched. Each block processes a different column of $\mathit{\_tmpMasks}$ . Similarly to $updateTableGPU()$ , also for the $\mathit{reduce}()$ kernel the division of work is dynamic and scales with the number of rows in the table.

Consider the $i$ -th block of the $reduce$ grid, each thread of such a block looks at a word (same column) in different rows of $\mathit{\_tmpMasks}$ performing a $\mathit{bitwise}$ - $\mathit{AND}$ operation. Once all the rows have been considered, all the threads of the same block perform a parallel reduction (applying a $\mathit{bitwise}$ - $\mathit{AND}$ ) to produce a single 32-bit word per block (see Figure 2). After this kernel is executed, the obtained mask is copied back to the host.

4.2 The propagator $\textbf{CT}_\textbf{CU}^\textbf{f}$

In the case of the propagator $CT_{CU}^{f}$ , the update of $\mathit{\_currTable}$ is carried out by the host while the filtering procedure is carried out on the device by means of the $\mathit{filterDomainsGPU}()$ kernel. Such a kernel (see Figure 3) is responsible for providing the host, for each variable in the scope, the values to remove from their respective domains. The initialization and data structures used are the same as for $CT^u_{CU}$ , with the exception of a new array $\mathit{\_vars\_to\_remove\_host}$ on the host side. Such an array, of size $\lceil \mathit{\_supportSize}/32\rceil$ , stores the result of the $\mathit{filterDomainsGPU}()$ kernel.

As concerns the $\mathit{filterDomainsGPU}()$ kernel, the number of blocks launched is $\lceil \mathit{\_supportSize}/32\rceil$ , each consisting of 32 threads. The division of work is thus dynamic and depends on the domain size of the variables in the scope. Each block is responsible for processing a single 32-bit word of the variable domains stored in $\mathit{\_vars\_dev}$ . All threads in the same block check whether the same value is still present in the domain of the variable $x_j$ considered. If such is the case, a parallel reduction on the 32 threads is launched on the respective support row, computing a $\mathit{bitwise}$ - $\mathit{AND}$ with the current table. After the parallel reduction, the $i$ -th bit of a shared word in the block is set to 1 if the value considered must be removed from the variable, or to 0 otherwise. The kernel thus computes an array of size $\lceil \mathit{\_supportSize}/32\rceil$ , in which the bit in position $\mathit{\_supportJmp[j]}+i-\mathit{\_variablesOffsets[j]}$ is set to $1$ if and only if the $i$ -th value of the variable $x_j$ must be removed from its domain. Each block is responsible for writing one word in $\_vars\_dev$ . Such an array is then copied back to the host into $\mathit{\_vars\_to\_remove\_host}$ and the values are removed from the domains of the respective variables.

Fig 4. Overview of the execution flow of $CT_{CU}^{u}$ . Kernels are depicted by green activities.

As an example consider a table with 1000 rows and suppose the scope has two variables, $x_1$ and $x_2$ , respectively with $dom(x_1)=\{100,\ldots ,196\}$ and $dom(x_2)=\{1,\ldots ,100\}$ . A total of 7 blocks, of 32 threads each, are launched. Consider now the second block (index 1), all 32 threads of such a block check the 32 bits responsible for keeping track of whether, for the variable $x_1$ , the values in the set $\{132,\ldots , 163\}= \{\mathit{\_variableOffsets[0]}+32,\ldots ,\mathit{\_variableOffsets[0]}+63\}$ are still present in $dom(x_1)$ .

4.3 The propagator $CT_{CU}^{uf}$

This last implementation combines the kernels designed for $CT_{CU}^{u}$ and $CT_{CU}^{f}$ .

The kernel $\mathit{filterDomainsGPU}()$ is called after $\mathit{updateTableGPU}()$ and $\mathit{reduce}()$ , so that both the table update and the domain filtering are carried out by the device. The data structures used are those introduced in both $CT^u_{CU}$ and $CT^f_{CU}$ .

5.1 The benchmark suite

We report here on the experimentation we conducted using three collections of instances. The instances of the first set involve a single table constraint and a single linear equation on some of the variables in the scope of the table constraint, as in the following example:

Note that we admit variables in the linear constraint not occurring in the table constraint. We generated two batches of instances, called $LIN\_B$ and $LIN\_EB$ resp., differing in the number of variables and tuples (see Table 2). The (synthetic) instances of this benchmark are based on a simple encoding of the bounded knapsack problem. More specifically, in the bounded Knapsack problem there are $k$ items and a number representing the capacity of the knapsack. Each item is associated with a weight, a profit, and a maximum number of times the object can be included in the knapsack. The goal is to select the number of instances of each item maximizing the profit which can fit the knapsack (Assi and Haraty Reference ASSI and HARATY2018). This problem can be modeled using two constant arrays ( $\textit {profit}$ and $\textit {weight}$ ) and a variable array representing, for each object, the number of times it is put in the knapsack. A table describing all possible knapsack configurations (i.e., the variables in the scope are the elements of the variable array) is defined. Additionally, a constraint requiring the weights of the objects considered to fit the knapsack is added (this is the linear equation mentioned earlier).

Table 2. Overview of table constraint features in the first test set

The tables in the $LIN\_B$ and $LIN\_EB$ tests may not involve all variables declared in the instance and do not fix a specific solving strategy, so such tests may trigger different numbers of propagations when solved by different solvers (e,g., Gecode and $CT^{uf}_{CU}$ ). In order to have a fair comparison between our implementations and Gecode, we introduced a second set of instances where both the linear equation and the table constraint involve all variables occurring in the instance. Moreover, in each solver, the solving strategy has been fixed as follows:

Finally, we experimented with a set of instances of the $\textit {Orienteering problem (OP)}$ or Single-Agent Path Planning with Reward Maximization (Bock and Sanità Reference BOCK and SANITÀ2015). The OP seeks a path for an agent starting and ending at specified locations, imposing a bound $D$ on the length of the path and maximizing the reward accumulated by visiting the various sites along the path, up to a threshold $C$ . The problem can be formalized as a problem on an undirected graph $G = \langle V, E\rangle$ , with nodes $V=\{v_1,\ldots ,v_k\}$ and edges $E$ , such that, for each edge $(v_i,v_j)$ , $w_{i,j}$ denotes its length and for each node $v_i$ , $\pi _i$ is a profit value associated to $v_i$ . Given such $G = \langle V, E\rangle$ , a starting and an ending nodes $v_s, v_e \in V$ , a length bound $D \in \mathbb{R}$ , a set $\Pi =\{\pi _1,\ldots ,\pi _k\}$ of profit values for nodes, and a profit bound $C \in \mathbb{R}$ , an instance of OP asks to determine a path $P=\langle (v_{p_1},v_{p_2}),(v_{p_2},v_{p_3}),\ldots ,(v_{p_k-1},v_{p_k})\rangle$ (with $v_{p_1}=v_s$ and $v_{p_k}=v_e)$ such that $\sum _{(v_{p_i},v_{p_j}) \in P}w_{{p_i},{p_j}}\le D$ . Such a path must also maximize the profit collected by visiting the nodes, subject to the bound $C$ , namely, $(\sum _{(v_{p_i},v_{p_j})\in P} \pi _{p_i}) + \pi _{p_k} \leq C$ .

We used a collection of instances of OP in which $G$ is a 4-connected grid (a graph in which each node is connected to at most four other nodes) with $w_{ij} = 1$ for each $(v_i,v_j) \in E$ . Hence, a path can exit a node in one among (at most) four possible directions. We did two additional changes, w.r.t. the original definition of OP:

• • Each node can either contain a non-negative profit or be an obstacle, in which case it is not traversable by any path.

• • The profit corresponding to any node in the path of the agent is collected once.

To model this problem, we use $\lceil \frac {D}{2} \rceil$ table constraints, each describing agent position and accumulated reward across two consecutive time-steps. If $G$ is a grid of size $r\times s$ , then each table contains $2\times r \times s$ columns. In each row, only two cells are assigned values between $1$ and $D$ , representing agent position and collected reward. Instances with grid size from $7\times 7$ to $10\times 10$ have been generated (see Figure 5).

5.2 Results

We now report the results considering the solve time statistic provided by MiniZinc on the instances described in Section 5.1. Such tests have been executed on both CT $_{CU}^{f}$ and CT $_{CU}^{uf}$ implementations on a system equipped with an Intel Core i7-13700KF at base clock $3.4$ GHz (up to $5.4$ GHz in boost), 128 GB of DDR4 RAM and a NVIDIA GeForce RTX 4090 with 16K CUDA cores, grouped into 128 streaming multiprocessors, running at $2.23$ GHz (up to $2.52$ GHz). On the software side, the system runs Ubuntu 24.04 and CUDA 12.7. All tests have been executed with 15 minutes timeout. Barplots summarizing the serial ( $CT$ ) and CUDA solve times are now presented. The same results in tabular form are available in the appendix. Such tables show also the number of propagation performed for each test, allowing to calculate additional metrics, such as the performance gain in terms of propagations per second. Propagation times of some tests are analyzed more in detail in Section 5.2.1.

From Figures 6 and 7 it can be observed that the kernels responsible for the update step scale well, taking advantage of the computational power of the GPU. Indeed, for all sufficiently large instances, $CT^{uf}_{CU}$ outperformed $CT^{f}_{CU}$ , achieving an average speedup of $2.88$ compared to $2.11$ on $LIN\_B$ instances, and $4.35$ compared to $2.52$ on $LIB\_EB$ . Figure 5 shows an overview of the solving times and speedups for the Orienteering test set. Notice that while such models involve several other constraints in addition to multiple table instances, the speedups remain consistent with those of LIN_B and LIN_EB instances.

Fig 5. Barplot comparing serial and CUDA solve times for the OR instances. The grids considered are squared, the instance named $r\_C$ is an $OR$ instance with a $r\times r$ grid and profit bound $C$ .

Fig 6. Barplot comparing serial and CUDA solve times for the first batch of instances. Test instances where all implementations timed out are omitted.

Fig 7. Barplot comparing serial and CUDA solve times for the second batch of instances. Test instances where all implementations timed out are omitted.

Kernel execution profiling showed that during the tests performed, the GPU capabilities were not fully exploited. In all the tests, the average GPU utilization was roughly $45\%$ . This is because the amount of work offloaded to the GPU is often not enough to saturate all streaming multiprocessors. (To saturate the capabilities of the RTX 4090, much larger instances would be needed, involving much larger tables.) This supports the fact that the good performance achieved using the RTX 4090 is also expected when using GPUs equipped with less streaming multiprocessors.

Fig 8. Series plots of the first 20 $CT$ and $CT^{uf}_{CU}$ propagation times for $\mathit{TEST\_B\_10}$ and $\mathit{TEST\_B\_9}$ . The times are grouped for the update and filter procedures.

Note that forcing a finer division of work for the sole purpose of using all available streaming multiprocessors often results in performance degradation. In such cases, in fact, required work synchronization and additional data exchange generate overhead that outweighs the benefits of parallel computation.

As concerns memory transfers, the time taken by memory copies from the host to the device and vise versa is significant when compared to the pure processing time. In part, this is due to the fact that the $CT$ algorithm is very efficient. On average, on $LIN\_B$ and $LIN\_EB$ instances, the copy of data to the device can take up to $50\%$ of the kernel time, while the copy of the results back to the host can take up to $80\%$ of the kernel execution time.

5.2.1 $CT^{uf}_{CU}$ propagations in detail

In this section we briefly analyze the behavior of $CT^{uf}_{CU}$ executed on the first set of instances. We do so by analyzing the first 20 propagations of $LIN\_B\_9$ and $LIN\_B\_10$ , which represent the two most representative behaviors observed in experimenting with the $LIN\_B$ and the $LIN\_EB$ instances. In the first scenario, $s^{sup}$ has a stable size, making the serial filtering times quite constant. In contrast, the second case is characterized by high variability in the size of $s^{sup}$ , leading to spikes in the serial filtering times, see Figure 8. In both scenarios however, the times recorded for the filtering kernel remain almost constant, this is due to the fact that the size of the grid launched is independent from $|s^{sup}|$ . As it can be seen from Figure 8, concerning the first kind of instances, the filtering process is the procedure leading to the most sensible improvements when parallelized.

For what concerns the procedure responsible for updating $currTable$ , Figure 8 compares the sum of the times of kernels $\mathit{udpateTableGPU}()$ and $\mathit{reduce}()$ (in $CUDA\: update$ ) with those of $\mathit{updateTable}()$ . Such times, when paired with a small number of propagations, do not lead to a significative difference. Moreover, the sum of $\mathit{udpateTableGPU}()$ and $\mathit{reduce}()$ times, due to the grids not presenting a constant size, present the same time trends as the serial counterpart.

To have a direct comparison between the filtering and update steps performed on GPU and their serial counterparts, the times considered in Figure 8 do not include the copy and retrieval of the data nor the time taken to dump the domains in $\mathit{\_vars\_host}$ .

Fig 9. Series plot of the first 200 $CT$ and $CT^{uf}_{CU}$ propagation times for $\mathit{TEST\_EB\_19}$ .

Finally, we discuss the times spent by the overall propagations in $CT$ and $CT^{uf}_{CU}$ . For what concerns the serial propagations the time only includes the update of $\_currTable$ and filtering. On the other hand, for $CT^{uf}_{CU}$ the time accounts for the kernels, the dump of the domains into $\_vars\_host$ , the copy and retrieval of the data. Figure 9 shows the timings of the first 200 propagations of $LIN\_EB\_19$ . Such an instance has been chosen since already in early propagations it presents a significant backtracking phase. As can be seen, once the valid tuples and domains start to reduce, the serial propagation times become very close to the $CT^{uf}_{CU}$ ones. In some instances requiring many propagations toward the bottom part of the search tree, it can also happen that the serial propagations run faster than the parallel ones. In Figure 9 it can be noticed that around propagation no. 125 a substantial backtracking phase takes place. Such an event leads the domains and valid tuples in the table to increase again, advantaging $CT^{uf}_{CU}$ over the serial propagator.

5.3 $CT^{uf}_{CU}$ and gecode

We now present some of the results obtained by comparing $CT^{uf}_{CU}$ against Gecode (version 6.3.0). Both solvers rely on a single thread (host side) with no parallelization features enabled. The comparison involves the second set of instances of Section 5.1. The instances have a number of variables between 100 and 150, a maximum domain size of 2000, and between 10000 and 15000 tuples in the table. The number of propagations for the tests in Table 3 ranges between 35000 and 85000. On some instances we noticed a slightly different number of propagations in favor of MiniCP. This difference is less than $2\%$ and is likely due to $CT^{uf}_{CU}$ relying on the built-in method $\mathit{changed}()$ and not using lastSizes to update $s^{val}$ . However, the number of explored nodes and failures is always the same for both solvers. Note that while the instance features (such as number of tuples or domain sizes) are comparable with those of $LIN\_EB$ instances, the number of propagations (see also the tabular data in the Appendix) result very different. This is due to the fact that all variables are now involved in the table constraint. Table 3 reports an excerpt of the instances we used, along with the solve times obtained by Gecode and $CT^{uf}_{CU}$ .

Table 3. The solving times (in seconds) for gecode and $CT^{uf}_{CU}$