remapping). However, as the building of the mapping consumes only 0.9% of computation time, the latter portion is barely visible in Figure 6.3. Table 6.2 shows total elapsed time as well as the time spent in the different operations of the algorithm as an average over the 10 instances *Ta021-Ta030*.

<table><caption>Table 6.2: Average elapsed time (in seconds) and average repartition of execution time among bounding, IVM management and remapping operations. Average taken over instances *Ta021-Ta030*.</caption><thead><tr><th rowspan="2">Bound mapping method</th><th>Wall-clock time</th><th colspan="2">bounding operation</th><th colspan="2">IVM management operations</th><th colspan="2">Bound mapping operation</th></tr><tr><th>(sec)</th><th>(sec)</th><th>(%)</th><th>(sec)</th><th>(%)</th><th>(sec)</th><th>(%)</th></tr></thead><tbody><tr><td>Static</td><td>696.4</td><td>632.9</td><td>89.4</td><td>63.5</td><td>10.6</td><td>0.0</td><td>0.0</td></tr><tr><td>Remapping</td><td>395.7</td><td>329.1</td><td>82.0</td><td>63.4</td><td>17.1</td><td>3.4</td><td>0.9</td></tr></tbody></table>

The compacted mapping *remap* is clearly advantageous as it reduces the average time spent in the bound kernel by a factor 1.9. As the bounding operation amounts for more than 80% of the total execution time, the latter decreases by a factor 1.7. The overhead induced by compacting the mapping at each iteration is largely compensated by these performance gains. Indeed, thanks to the parallelization of this operation using the parallel prefix sum, the remapping operation amounts for less than 1% of the elapsed time. For comparison, using the CPU for the remapping, it amounts for about 7% of the algorithm's total execution time, mainly because of the transfer of the maps back to the device.

Using the more compact mapping *remap* instead of *static* improves the control flow efficiency<fnref n="2" /> (CFE) of the kernel. For *static* the average CFE is 0.43, meaning that for an executed instruction on average more than half of the execution slots are wasted. For the mapping *remap* the average CFE is 0.83 - the launched warps are used almost twice as efficiently. The number of warps launched at each kernel call is 960 for mapping *static*, which exceeds theoretical maximum of $13 \times 64 = 832$ resident warps for the K20m. The average number of warps launched with mapping *remap* is 300 (average per kernel call and per instance), the average maximum (per instance) being 825 warps and the minimum 4. These results show that it is a high priority optimization to adapt the configuration of the algorithm's most time-depending part to the varying workload.

<footnote n="2">²defined as $CFE = \frac{\text{not\_predicated\_off\_thread\_inst\_executed}}{32*\text{inst\_executed}}$</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 6.4: Execution time for instances *Ta021-Ta030* for different mapping choices in IVM management kernels.</figcaption></figure>

### 6.4.3 Evaluation of IVM mapping schemes

In this subsection the two mapping schemes presented in Subsection 6.3.2 are evaluated and compared to each other. The kernels concerned by these mapping schemes are the IVM management kernels (i.e work stealing, branching, selection and pruning). Figure 6.4 shows the time spent for completing the exploration with both mapping schemes. Both, version *one-thread-per-IVM* (*M1*) and version *one-warp-per-IVM* (*M2*) use the same bounding kernel (with remapping). Although the time spent managing the IVM structures is moderate compared to the bounding operation, the mapping *M2* allows a reduction of the total execution time by a factor 1.1 compared to the mapping *M1*. With respect to *M1*, mapping *M2* decreases the share of IVM management operations from 18% to 7.5%. Table 6.3 shows the average duration per call of the kernels bound (in msec), branching and selection (in µsec) and their respective share of the elapsed time (in %). The kernels pruning and work stealing amount for at less than 2% of total execution time, so they are not evaluated.

The mapping *M2* allows to use the supplementary lanes for efficient loading of the IVM structures into shared memory. In order to dissociate the impact of shared memory usage from the impact of remapping, the profiling of mapping *M2* is performed with and without shared memory usage.
### 6.4.4 Evaluation of work stealing schemes

In this subsection the behavior of the algorithm according to the instance sizes and its scalability with the number of used IVM structures ($T$) is examined. The algorithm's performance for problem instances of different sizes is compared in terms of node processing speed (#branched nodes/second), which is computed from wall-clock time. Obviously, using more explorers can only be beneficial if they can be supplied with enough work. Therefore, the relationship between instance size, the node processing speed and $T$ needs to be studied for both proposed work stealing strategies, as they strongly impact this relationship. For the experimental study of scalability only the best version of the previous subsection is considered, i.e. the one using parallel remapping for the bounding kernel and the spaced mapping $M2$ for the management kernels. Another factor that has a significant impact on the algorithm's performance is the irregularity of the explored PB&B tree, which is very hard to quantify. In order to obtain a clearer dissociation between tree size and tree irregularity, the average node processing speed for instances of similar size is considered. The instances have been grouped as shown in Table 6.5.

<table><caption>Table 6.5: Groups of similar sized flowshop instances (20 jobs × 20 machines) and the corresponding average number of nodes branched when initialized with the optimal complete permutation.</caption><thead><tr><th>Instance group</th><th>#branched nodes</th><th>Flowshop Instances</th><th>Average #branched nodes</th></tr></thead><tbody><tr><td>small</td><td>&lt; 10M</td><td>Ta028,Ta029,Ta030</td><td>5.5M</td></tr><tr><td>medium</td><td>∈ [10M, 50M]</td><td>Ta021,Ta022,Ta024,Ta025</td><td>36.3M</td></tr><tr><td>large</td><td>∈ [50M, 100M]</td><td>Ta026,Ta027</td><td>64.3M</td></tr><tr><td>huge</td><td>&gt; 100M</td><td>Ta023</td><td>140.8M</td></tr></tbody></table>

Figure 6.5 shows the average node processing speed for these four groups of instances, according to different values of $T$ and work stealing strategies 1D-Ring and 2D-Ring. $T$ is chosen as a multiple of 64, as the GK110 architecture allows up to 64 warps per SM and the management kernels reserve one warp per IVM. The number of rings $R$ in the 2D-Ring work stealing strategy is chosen such that $R$ divides $T$ while being as close to $\sqrt{T}$ as possible – the goal being to approach the ideal configuration where $R$ equals the ring-size $C$.
<figure><img src="image_21.png" /><figcaption>Figure 6.5: Average node processing speed for groups of instances *small*, *medium*, *large*, *huge* (Table 6.5), and $T = k \cdot 256$ ($k = 1, 2, ..., 6$).</figcaption></figure>

Left: 1D-Ring work stealing ; Right: 2D-Ring work stealing with dimensions $R \times C$
$= 16 \times 16, 16 \times 32, 24 \times 32, 32 \times 32, 32 \times 40$ and $32 \times 48$

From Figure 6.5 one can see that the node processing speed for *small* instances is lower than for intermediate sized instances, regardless of parameter *T* or the used work stealing strategy. This is partially due to warm-up and shut-down phases of the parallel exploration, which, for *small* instances last relatively long with respect to the total exploration time. In these phases a low overall workload limits the degree of parallelism and the work stealing mechanism must handle sharp variations of the workload. Using the 2D-Ring topology significantly improves the nodes-per-second rate for *small* instances, as the reduced diameter of the 2D-Ring accelerates the distribution of the workload. Another reason for the poorer performance for *small* instances is that the finer granularity of allocated tasks per IVM does not allow one to hide the initialization costs as efficiently as for larger instances. Similarly, the nodes-per-second performance for *medium* instances is below the performance for *large* and *huge* instances. As the groups *large* and *huge* do only count one, respectively two members, the results for these groups must be interpreted more carefully – however, they suggest that the node processing speed for large or very large instances is less influenced by the tree’s size. The size of the instance also impacts the scalability of the algorithm with *T*, especially when the 2D-Ring work stealing strategy is used. Using the 2D-Ring topology, for *small* sized problems the best value for *T* improves performance to ≈ 120k nodes/sec from ≈ 85k nodes/sec for the worst *T*, while for *large* sized problems the node processing speed doubles from ≈ 105 to ≈ 210k nodes/sec, doing the same comparison.

Comparing both work stealing strategies, one can see from Figure 6.5 that the 2D-Ring topology improves the scalability of the algorithm. Indeed, for the 1D-Ring almost no performance is gained above $T = 768$ because additional workers are left idle or inefficiently initializing. In contrast, the 2D-Ring strategy allows one to use up to $T = 1280$ IVMs efficiently. For $T = 1536$, and only when using the 2D-Ring work stealing strategy,
performance drops significantly for all instance sizes. The most likely explanation for this
performance drop is that the computation of the bounds is partially serialized due to
hardware limitations. On average, each node branching leads to the evaluation of approx-
amately 20 bounds, for all considered instance sizes. So, if all T = 1280 IVMs are busy,
25600 bounds are evaluated per average iteration, which is close to the hardware limit of
26624 concurrent threads (13 SM ×64 warps/SM ×32 threads/warp) for the GK110 ar-
chitecture. Therefore, supposing a well balanced workload, no performance improvement
can be expected from increasing T beyond 1280 on this device and at constant tree size.

<figure><img src="image_22.png" /><figcaption>Figure 6.6: *IVM-efficiency= branched nodes / iterations×#IVM, measures the overall share of IVM struc- tures in the *exploring* state. An IVM efficiency of 1 indicates that each IVM branches a node at each iteration.</figcaption></figure>

However, Figure 6.6 shows that there is some margin left to increase the percentage of busy explorers. In order to measure the overall efficiency of the work stealing strategies, the number of branched nodes is divided by (#iterations × #IVM), which gives a measure for the portion of IVMs doing meaningful work. If this *IVM efficiency* equals 1, then each IVM branches a node at each iteration – if it equals 0, then all IVMs remain idle for an infinity of iterations. An *IVM efficiency* of 0.5 indicates that an average IVM spends half of its iterations idle or initializing. The results shown in Figure 6.6 confirm that the 2D-Ring work stealing allows to keep more explorers busy. Moreover, the fact that at *T* = 1536 no sharp degradation of the IVM efficiency occurs, confirms that the performance drop at *T* = 1536 is due to hardware limitations. One can also see in Figure 6.6 that for *small* and *medium* instances the IVM efficiency is lower than for larger instances. This explains the relatively poorer performances for these instances. The smaller the instance being solved, the likelier it becomes that a given IVM frequently steals intervals that are explored within a few or zero iterations. In that case, the cost of initialization is not covered by meaningful work. Without any *a priori* knowledge concerning the amount of work contained in a given interval, it seems difficult to resolve this problem.

Concerning the algorithm's sensitivity towards the irregularity of the tree structure,
the performances obtained when solving same sized instances are compared in Figure 6.7.
Figure 6.7 zooms on the scalability with *T* of the three instances with ≈ 40M nodes. For all three instances the performance drop increasing the number of IVMs from *T* = 1280 to *T* = 1536 is clearly visible.

<figure><img src="image_23.png" /><figcaption>Figure 6.7: Scaling with T of the medium sized instances</figcaption></figure>

Figure 6.7 also shows that the instance irregularity impacts performance significantly. Exploring a roughly equal number of nodes, the peak node processing speed attained for instance *Ta024* is 202k nodes/sec while it remains below 175k nodes/sec for instances *Ta021* and *Ta025*. However, in sequential execution *Ta024* has the best nodes-per-second rate of these three instances, which means that the performance gap is not due to a higher number of nodes in the upper part of the tree (which are more costly to evaluate). The performance variety under same sized instances rather seems to be due to sharper workload variations for instances *Ta021* and *Ta025*. Although for the three instances the proposed 2D-Ring strategy scales nicely up to 1280 IVM structures, it apparently has a moment of inertia which makes the algorithm sensitive to rapid workload variations. Maybe this gap can be closed with a further improved work stealing strategy, in the sense that the work stealing strategy should deliver a faster response to those variations. This is a difficult task, in particular because IVM needs to go through initialization after the reception of a new work unit.

### 6.4.5 Comparison of PB&B@GPU and GPU LL PB&B

In this subsection the best version (according to the previous experimentations) of the PB&B@GPU is compared to the GPU-accelerated PB&B presented in [CM13] which uses a conventional LL for the storage and management of the pool of permutations. Like the PB&B@GPU, GPU LL performs a depth best-first search, retaining after each branching the better of two generated pools. Moreover, both algorithms use the same device function to compute the lower bounds for a given permutation. Table 6.6 shows the elapsed time for solving each of the ten 20-job instances, as well as the number of branched nodes per second, in order to take into account the size of the instances. The average time spent by the PB&B@GPU for exploring the ten instances is 229.7 seconds, while the GPU
LL PB&B algorithm requires on average 865.6 seconds for the same tasks. In terms of node processing speed, the PB&B@GPU algorithm branches on average 3.3 times more nodes per second than its LL counterpart. The PB&B@GPU outperforms GPU LL PB&B by at least a factor 2 for all ten instances. The highest nodes-per-second performance is attained for instance *Ta027*, branching 213,000 nodes/sec, which is almost 4 times more than the highest rate attained by the LL-based algorithm. As examined in Subsection 6.4.4, for the PB&B@GPU algorithm the node processing speed varies from one instance to another, depending on the instance's size and irregularity. In contrast, as the GPU LL PB&B algorithm regularizes the workload by dynamically adapting the size of the offloaded pools it provides an almost constant node processing speed for all instances. This performance variation is discussed in Subsection 6.4.4.

<table><caption>Table 6.6: Elapsed execution time and number of branched nodes per second (in 1000 nodes/sec) for the PB&B@GPU and the GPU LL PB&B [CM13]. The PB&B@GPU algorithm uses the 2D-Ring work stealing strategy and $T = 1280$ IVMs for all instances. Results for instances *Ta021-Ta030*.</caption><thead><tr><th rowspan="2">Flowshop instance</th><th colspan="2">GPU LL PB&B</th><th colspan="2">PB&B@GPU</th><th rowspan="2">Ratio</th></tr><tr><th>elapsed</th><th>(k nodes/sec)</th><th>elapsed</th><th>(k nodes/sec)</th></tr></thead><tbody><tr><td>Ta021</td><td>833</td><td>49.7</td><td>242</td><td>171.1</td><td>3.4</td></tr><tr><td>Ta022</td><td>415</td><td>53.3</td><td>134</td><td>163.6</td><td>3.1</td></tr><tr><td>Ta023</td><td>3089</td><td>45.6</td><td>740</td><td>189.6</td><td>4.2</td></tr><tr><td>Ta024</td><td>738</td><td>54.3</td><td>200</td><td>202.2</td><td>3.7</td></tr><tr><td>Ta025</td><td>865</td><td>47.9</td><td>239</td><td>173.1</td><td>3.6</td></tr><tr><td>Ta026</td><td>1292</td><td>55.3</td><td>348</td><td>205.6</td><td>3.7</td></tr><tr><td>Ta027</td><td>1094</td><td>52.2</td><td>268</td><td>213.3</td><td>4.1</td></tr><tr><td>Ta028</td><td>171</td><td>47.4</td><td>53</td><td>152.0</td><td>3.2</td></tr><tr><td>Ta029</td><td>125</td><td>54.4</td><td>56</td><td>121.4</td><td>2.2</td></tr><tr><td>Ta030</td><td>34</td><td>47.1</td><td>17</td><td>94.6</td><td>2.0</td></tr></tbody><tfoot><tr><td>Average</td><td>865.6</td><td>50.7</td><td>229.7</td><td>168.6</td><td>3.3</td></tr></tfoot></table>

## 6.5 Conclusions

This chapter proposes a GPU-centric PB&B algorithm which performs all PB&B operations on the GPU. During the exploration of the PB&B tree, the CPU core is only used for launching the CUDA kernels in a loop until a boolean variable, which the CPU receives at each iteration from the GPU, indicates the end of the algorithm. To the best of our knowledge, our PB&B@GPU algorithm is the first one that does not rely on the transfer of pools of permutations between host and device.

The proposed approach is based on the IVM structure, better adapted to the GPU than
LL-based data structures, which are conventionally used for the storage and management of the pool of permutations. The algorithm has two levels of parallelism. On a lower level it efficiently uses up to 1280 IVM structures to perform the branching, selection and pruning operations in parallel, exploring different parts of the PB&B tree simultaneously. For each exploring IVM the bounding operation is in turn parallelized, leading to an increased overall degree of parallelism in the bounding operation. At the junction of the two levels a remapping operation is introduced in order to adapt the configuration and the mapping of the bounding kernel to the varying workload. While the mapping for the bounding is compacted, the mapping in the management phase – characterized by a very high number of data dependent conditional instructions – is spaced, adding idle threads to the kernel. As a result, the bounding operation is accelerated on the one hand, as the control flow efficiency is improved – on the other hand, in the IVM management phase the opposed strategy is better adapted, as it reduces thread divergence and improves memory accesses.

We have proposed two work stealing strategies for work load balancing and analyzed the scalability of our algorithm with respect to both strategies. The reported experimental results show that the performance of the proposed algorithm depends crucially on the choice of the mapping, as well as on the used work stealing strategy. The proposed PB&B@GPU algorithm explores the PB&B trees of 10 Taillard's flowshop instances *Ta021-Ta030* on average with a 3.3 times higher node processing speed than a GPU LL PB&B algorithm. For all instances *Ta021-Ta030* the PB&B@GPU algorithm outperforms the GPU LL PB&B algorithm at least by a factor 2.0 and by up to a factor 4.0 on large instances, where our PB&B@GPU algorithm reaches its best performances.
CHAPTER 7

# Permutation B&B for GPU-powered clusters

## Contents

<table><tr><td>7.1 Introduction</td><td>95</td></tr><tr><td>7.2 Interval-list operators</td><td>96</td></tr><tr><td>7.2.1 Interval operators</td><td>96</td></tr><tr><td>7.2.2 Interval-list operators</td><td>97</td></tr><tr><td>7.3 Data structures and operations</td><td>99</td></tr><tr><td>7.3.1 Data structures</td><td>99</td></tr><tr><td>7.3.2 Work processing operations</td><td>101</td></tr><tr><td>7.4 Experiments</td><td>103</td></tr><tr><td>7.4.1 Experiments on a cluster of GPUs</td><td>103</td></tr><tr><td>7.4.2 Experiments on a cluster of GPUs and CPUs</td><td>108</td></tr><tr><td>7.4.3 Resolution of a big instance</td><td>109</td></tr><tr><td>7.5 Conclusions</td><td>113</td></tr></table>

## 7.1 Introduction

This chapter presents a parallel distributed hybrid version of the PB&B algorithm (PB&B@CLUSTER) for large-scale heterogeneous clusters, integrating multi-core CPU and graphics processing units (GPUs). PB&B@GPU is based on the B&B@Grid framework [MMT07c]. B&B@Grid allows to efficiently partition the PB&B tree search among distributed computing nodes, which host one or several workers. Each worker explores a portion of the search space (an interval) using a sequential PB&B. Thus, while compute nodes in B&B@Grid may be composed of multi-core processors, the latter are seen as a collection of single-core processors. Therefore, B&B@Grid is revisited with the goal of adapting it to heterogeneous computing platforms.
This chapter is divided into three main sections. Section 7.2 defines a set of operators used to handle interval-lists. Section 7.3 presents the operations of the coordinator. Section 7.4 describes the experiments carried out to validate the PB&B@CLUSTER approach.

## 7.2 Interval-list operators

In our approach, there are two types of operators, namely the interval-list and interval operators. The interval-list operators are based on interval operators. Before presenting the interval-list operators, this section first describes the interval operators.

### 7.2.1 Interval operators

Like the PB&B@CPU and PB&B@GPU approaches, a work unit, in the PB&B@CLUSTER approach, is also an interval, and each interval corresponds to an IVM. An interval can be seen as a sorted set of positive factorial integers. In order to define our approach, we have introduced some operators on these intervals. As an interval is a set of integers, these operators are inspired by operators known in set theory. The interval operators are therefore the intersection, the subtraction, the right and left divisions, the cardinality and the norm. Let's assume:

* $\mathbb{I}$: interval space;

* $\mathbb{L}$: interval-list space;

* $\mathbb{N}$ and $\mathbb{N}^+$: respectively the space of integers and strictly positive integers;

* $[a, b[$: an interval;

* $[a_1, a_2[$ and $[b_1, b_2[$, two other intervals;

* $n$: a positive integer;

The interval operators and the empty interval are defined as follows:

* Intersection operator ($\cap$):

$$ \mathbb{I} \times \mathbb{I} \rightarrow \mathbb{I} \\ [a_1, a_2[ \cap [b_1, b_2[ \mapsto [max(a_1, b_1), min(a_2, b_2)] $$

* Subtraction operator ($\setminus$):

$$ \mathbb{I} \times \mathbb{I} \rightarrow \mathbb{L} \\ [a_1, a_2[ \setminus [b_1, b_2[ \mapsto {[a_1, min(a_2, b_1)], [max(a_1, b_2), a_2]} $$
* Right division operator (÷):

$$ \mathbb{I} \times \mathbb{N}^{+} \rightarrow \mathbb{I} \\ [a, b[ \div n \mapsto \begin{cases} [a, b[, & \text{if } (b-a) \le \epsilon \\ [b - \lfloor(b-a)\frac{1}{n}\rfloor, b[, & \text{if } (b-a) > \epsilon \\ 0, & \text{otherwise} \end{cases} $$

* Left division operator (÷):

$$ \mathbb{N}^{+} \times \mathbb{I} \rightarrow \mathbb{I} \\ n \div [a, b[ \mapsto \begin{cases} [a, b[, & \text{if } (b-a) \le \epsilon \\ [a, a + \lfloor \frac{n-1}{n} (b-a) \rfloor[, & \text{if } (b-a) > \epsilon \\ 0, & \text{otherwise} \end{cases} $$

* Cardinality operator (|.|):

$$ \mathbb{I} \rightarrow \mathbb{N} \\ |[a_1, a_2[| \mapsto \begin{cases} 1, & \text{if } a_1 < a_2 \\ 0, & \text{otherwise} \end{cases} $$

* Norm operator (||.||):

$$ \mathbb{I} \rightarrow \mathbb{N} \\ \| [a_1, a_2[ \| \mapsto \begin{cases} a_2 - a_1, & \text{if } a_1 < a_2 \\ 0, & \text{otherwise} \end{cases} $$

* Emptiness value (Ø):

$$ ([a, b[=Ø) \iff (\| [a, b[ \| = 0) $$

### 7.2.2 Interval-list operators

A series of intervals form an interval-list. In our approach, each worker handles a list of IVMs and, therefore, an interval-list. In the same way, the previous subsection defines some interval operators, this section extends the definition of these operators on interval-lists, and also defines the concept of empty interval-list. Let's assume:
* $A = (A_1, ..., A_I)$: an interval-list of $I$ intervals;

* $B = (B_1, ..., B_J)$: another interval-list of $J$ intervals;

* $n$: an integer;

Interval-list operators and empty interval-list are defined as follows:

* Intersection operator ($\cap$):

$$ \mathbb{L} \times \mathbb{L} \to \mathbb{L} \\ A \cap B \mapsto \cup_{i=1, I}^{j=1, J} (A_i \cap B_j) $$

* Subtraction operator ($\backslash$):

$$ \mathbb{L} \times \mathbb{L} \to \mathbb{L} \\ A \backslash B \mapsto \cup_{i=1, I}^{j=1, J} (A_i \backslash B_j) $$

* Right division operator ($\div$):

$$ \mathbb{L} \times \mathbb{N}^{+} \to \mathbb{L} \\ (A_1, ..., A_I) \div n \mapsto (A_1 \div n, ..., A_I \div n) $$

* Left division operator ($\div$):

$$ \mathbb{N}^{+} \times \mathbb{L} \to \mathbb{L} \\ n \div (A_1, ..., A_I) \mapsto (n \div A_1, ..., n \div A_I) $$

* Cardinality operator ($|\cdot|$):

$$ \mathbb{L} \to \mathbb{N} \\ |A| \mapsto \sum_{i=1, I} |A_i| $$

* Norm operator ($\|\cdot\|$):

$$ \mathbb{L} \to \mathbb{N} \\ \|A\| \mapsto \sum_{i=1, I} \|A_i\| $$
* Emptiness value (Ø):

$$ (A = Ø) \iff (\|A\| = 0) $$

Compared to the previous subsection, this subsection also introduces a new operator called *biggest*. This operator receives, as input, an interval-list *A* and a natural integer *n*, and returns, as output, another interval-list *biggest*(*n*, *A*) made up of the *n* largest intervals of *A*. The intervals are compared with each other using the `||.` operator.

* Biggest (*biggest*):

$$ \mathbb{N} \times \mathbb{I} \to \mathbb{L} \\ \text{biggest}(n, A) \mapsto \begin{cases} (A_i)/\forall j, \|A_j\| \le \|A_i\|, & \text{if } n=1 \\ (\text{biggest}(1, A \setminus \text{biggest}(1, A)), \text{biggest}(n-1, A)), & \text{otherwise} \end{cases} $$

# 7.3 Data structures and operations

This section describes the PB&B@CLUSTER approach. Subsection 7.1 and Subsection 7.3.2 present respectively the data structures of this approach and its operations.

## 7.3.1 Data structures

<figure><img src="image_24.png" /><figcaption>Figure 7.1: The data structure of the coordinator in PB&B@CLUSTER approach.</figcaption></figure>
As shown in Figure 7.1, our approach PB&B@CLUSTER is based on a coordinator-workers
architecture. In this architecture, there are a single coordinator process and a large num-
ber of workers. The coordinator process has two roles: On the one hand, it allows PB&B
workers to exchange the best complete permutation found so far and, on the other hand,
this coordinator is used by workers for implementing a work stealing strategy. Compared
with the approaches described in previous chapters, the PB&B@CLUSTER approach intro-
duces a new data structure. This structure allows to save the best complete permutation
found so far, as well as a copy of the interval-lists being processed on each worker. This
data structure is defined as follows:

* **Best found permutation:** This variable saves the best complete permutation found so far and the cost of this permutation.

* **Work of an IVM:** The interval $W_i^j$ is the $j^{th}$ interval of the $i^{th}$ worker. In other words, $W_i^j$ is associated with the $j^{th}$ IVM of the $i^{th}$ worker. This interval $W_i^j$ is equal to $[P_i^j, E_i^j]$ [such that:

- $P_i^j$: the position vector of IVM $j$ in the worker $i$;

- $E_i^j$: the end vector of IVM $j$ in the worker $i$;

A copy of the interval $W_i^j$, denoted $W_i'^j$, is always stored at the coordinator. The copy $W_i'^j$ is equal to $W_i^j$ only when the worker sends a request to the coordinator. Between two requests, the copy $W_i'^j$ is probably not equal to the interval $W_i^j$.

* **Work of worker:** The work of a worker $W_i$ is an interval-list $(W_i^1, ..., W_i^{S_i})$, such that:

- $S_i$: #IVM in the worker $i$

- $W_i^j$ work of IVM $j$ in the worker $i$

For example, the worker 2 in Figure 7.1 has 6 intervals. A copy of the interval-list of
worker *i*, denoted $W_i'$, is stored at the coordinator. The copy $W_i'$ of the coordinator
is equal to $W_i$ only the first time the worker contacts the coordinator to get work.

* **Work of all workers:** $W'$ is the list of all copies $(W'_1, ..., W'_N)$, such that $W'_i$ is the copy of the interval-list $W_i$ of the worker $i$. For example, the coordinator of Figure 7.1 manages 3 workers. These workers, numbered 1, 2 and 3, respectively handle 7, 6 and 8 IVMs. Therefore, the copy $W'$ has three interval-lists $(W'_1, W'_2, W'_3)$, containing each a list of intervals.
7.3.2 Work processing operations

In addition to the operators explained in the previous chapters, the PB&B@CLUSTER approach essentially introduces three new operations, namely a communication operation of the worker, a checkpointing operation of the coordinator and a work stealing operation of the coordinator. Figure 7.2 shows the relation between these operations on an example of a PB&B@CLUSTER approach using three workers.

<figure><img src="image_25.png" /><figcaption>Figure 7.2: Example of a PB&B coordinator with three PB&B workers.</figcaption></figure>

*   **Worker communication operation:** The worker $i$ contacts the coordinator when one of the following three conditions occurs:

- **The rate of empty IVMs is high:** When an IVM of a worker $i$ becomes empty, it is important that the worker, responsible for managing this IVM, operates a work steal to fill the empty IVM. However, if this work stealing operation is done immediately at each empty IVM, then there is a big risk of wasting a lot of time in communication. As a result, the worker $i$ only operates a work stealing when the rate of empty of IVMs $|W_i|/S_i$ is less than a certain threshold $\alpha$, where $|W_i|$ is the number of empty IVMs, and $S_i$ is the total number of IVMs of the worker $i$.

- **New best complete permutation found:** When a new best complete permutation is found by the worker $i$, this worker communicates to the coordinator this
new permutation. The interest of this communication is to allow the workers
to be informed as soon as possible of the new permutation found. This allows
these workers to prune more nodes from the PB&B tree. By communicating this
new permutation, the worker *i* communicates its work $W_i$ to the coordinator
to update its copy $W'_i$.

- **Maximum time limit passed:** After a request to the coordinator, the work $W_i$ of the worker $i$ is equal to the copy $W'_i$ of this work at the coordinator. Between two requests to the coordinator, the more time passes, the larger the gap between $W_i$ and $W'_i$. Therefore, it is important to configure the worker by a maximum time limit. If one of the two preceding conditions did not occur within this maximum time limit, then the worker $i$ must contact the coordinator to update its copy $W'_i$. This third condition allows the coordinator to be regularly informed of the progress of the work of the worker $i$.

* **Work checkpointing operation:** As shown in Algorithm 13, the operation of updating the coordinator's work $W'$ is quite simple. This operation receives as input the work $W'$ of the coordinator, the ID $i$ of the worker and its work $W_i$. The checkpointing operation performs an intersection of the two lists $W_i$ and $W'_i$. In the algorithm, the result of this intersection, denoted $W_{tmp}$, is assigned to the two lists $W_i$ and $W'_i$.

<caption>**Algorithm 13 COORDINATOR-WORK-CHECKPOINT($W'$, $i$, $W_i$)**</caption>

1: $W_{tmp} = W_i \cap W'_i$

2: $W'_i = W_{tmp}$

3: $W_i = W_{tmp}$

* **Work stealing operation:** Algorithm 14 explains the work stealing operation of the coordinator. As this algorithm shows, this operation is executed only if the interval $W_i$ is empty. In this case, the first instruction is to seek the greatest interval-list, denoted $W'_{j}$, of all the works in $W'$. Once this interval-list is determined, the second instruction calculates the minimum size, denoted $S$, between the cardinalities of $W'_{j}$ and $W'_{i}$. The cardinality of an interval-list is equal to the number of intervals in this list. The role of the fourth instruction is to take the larger $S$ intervals, denoted $W'_{tmp}$, from the interval-list $W'_{j}$. Then, the intervals of the list $W'_{tmp}$ are divided into 2 using a right division operation. The interval-list, obtained by this division operation, is assigned to the variable $W'_{i}$, and the list of intervals $W'_{j} \setminus W'_{i}$ is assigned to the variable $W'_{j}$. This work stealing operation returns the interval-list $W_{i}$ to worker $i$.
<caption>**Algorithm 14 COORDINATOR-WORK-STEALING(W', i, Wᵢ)**</caption>

1: **if** Wᵢ = ∅ **then**
2:     **j** = arg maxᵢ(||Wᵢ||)
3:     **S** = min(|Wᵢ|, |Wⱼ|)
4:     **W'**tmp = biggest(S, Wⱼ)
5:     **W'**ᵢ = **W'**tmp ÷ 2
6:     **W'**ⱼ = **W'**ⱼ \ **W'**ᵢ
7:     **W**ᵢ = **W'**ᵢ
8: **end if**

7.4 Experiments

This section describes the two experiments performed to validate the PB&B@CLUSTER
approach, on a cluster of GPUs and a hybrid cluster of GPUs and CPUs In addition, the
section presents the obtained results when solving a large instance.

7.4.1 Experiments on a cluster of GPUs

Experimental protocol

* **Hardware testbed:** The Ouessant cluster is used for this first experiment. This cluster is located at l'Institut du Développement et des Ressources en Informatique Scientifique (IDRIS<fnref n="1" />). The Ouessant Tier1 (national) cluster is composed of 12 computing nodes. Each node consists of 2 CPUs POWER8+ (10 core CPUs, 8 threads per core, so 160 threads per node) and 4 Nvidia GPUs (Generation Pascal P100, and 16 GB of memory). In total, this cluster is composed of about 170,000 cores CPUs and 240 cores CPUs. For reasons of availability, only 9 nodes are used in our experiments.

* **Problem instance:** The ten instances of size 20 × 20 (20 jobs and 20 machines), noted *Ta021*, *Ta022*, ... and *Ta030*, are used in the previous experiments. Among all instances of this size, the *Ta023* instance is the most difficult to solve, and requires approximately 22 hours with the PB&B@CORE approach on a single CPU core. Therefore, these instances are not difficult enough for a resolution using a PB&B@CLUSTER on the Ouessant cluster. On the other hand, the ten instances of size 50 × 20, denoted *Ta051*, *Ta052*, ... and *Ta060*, are difficult for a resolution on the Ouessant cluster. So we used the ten 50 × 20 instances to generate 50 instances

<footnote n="1">¹http://www.idris.fr/</footnote>
suitable for a resolution on Ouessant. These 50 instances are denoted $Ta0i-h$, with $i \in \{51, 52, ..., 60\}$ and $h \in \{1, 2, ..., 5\}$.

To generate a $Ta0i-h$ instance, the PB&B@CORE approach is deployed with a single worker. This worker solves the $Ta0i$ instance by using a single GPU of Ouessant for $h$ hour (s). In addition, this resolution is initialized with the best known permutation of the instance $Ta0i$. At the beginning of the resolution, the work to be explored is the interval $[0, N![$, with $N = 50$ the size of the instance. At the end of these $h$ hours, it remains a list of intervals $L_i^h$ to continue exploring. Therefore, the interval-list that is already explored is $([0, N!]) \setminus L_i^h$. So the instance $Ta0i-h$ is defined as (1) the resolution of the $Ta0i$ instance (2) by initializing the PB&B with the list of intervals $([0, N!]) \setminus L_i^h$ and (3) the best known permutation of $Ta0i$.

## Obtained results

This section presents the results obtained when solving five instances $Ta0i-h$, with $i \in \{54, 55, 57, 58, 59\}$ and $h \in \{5\}$. Each of these instances is solved using six deployment configurations, namely 1, 4, 8, 16, 24, and 36 GPUs workers. As a result, 30 tests (i.e. 5 \* 6) are performed in total. The metrics measured in this experiment are, for example, percentage of redundantly explored nodes, GPU speedup, number of checkpointing operations, number of work-stealing operations, elapsed-time, etc. For each of these metrics, an average is calculated for each of the six deployment configurations.

* Redundant nodes:

<figure><img src="image_26.png" /><figcaption>Figure 7.3: Rate of redundant branched nodes.</figcaption></figure>

One of the few disadvantages of the PB&B@CORE approach is that some nodes of the PB&B tree can be explored multiple times. Therefore, it is important to estimate
the percentage of nodes explored more than one time. Figure 7.3 shows that with a single GPU worker, no node is explored twice. As the number of GPU workers increases, the percentage of redundant nodes increases. However, the percentage of redundant nodes is on average equal to 0.8% when 36 GPU workers are used. These experiments show therefore that the number of redundant nodes is not significant.

* Elapsed time and speedup:

<figure><img src="image_27.png" /><figcaption>Figure 7.4: Elapsed time and speedup over 1 GPU.</figcaption></figure>

Figure 7.4 shows the elapsed time (in minutes) in blue on the right y-Axis and the speedup with respect to a single GPU in black on the left y-Axis. The red dashed line corresponds to linear speedup with the number of GPUs. The average execution time on a single GPU is about 220 minutes. Using 36 GPUs the average execution time decreases to approximately 7.5 minutes, which corresponds to a relative speedup of 30× over a single GPU. These experiments show that the speedup remains good despite the high used computing power.

* Work checkpointing and stealing operations:

<figure><img src="image_28.png" /><figcaption>Figure 7.5: Number of checkpointing operations.</figcaption></figure>
As explained previously, the coordinator essentially performs two operations, namely work checkpointing and stealing operations. The objective is to study these two operations. In Figure 7.5 the number of checkpointing operations and work stealing operations is shown in black on the left y-Axis. The number of these two operations per second is shown in blue on the right y-Axis. Naturally, as the number of GPUs increases the number of these two operations also increases. For the considered scale, the rate of increase appears to be constant. A checkpointing operation is counted for both metrics. Therefore the number of checkpointing operations is always greater than the number of work stealing operations. One can see that the number of checkpoints increases faster than the number of work stealing operations. However, using more GPUs does not necessarily increase the number of checkpointing operations. Indeed, if all workers contact the coordinator in fixed and regular intervals, and if the elapsed time decreases linearly with the number of GPUs, then the total number of checkpointing operations would remain constant. In contrast, the results show that workers contact the coordinator more frequently. At near-linear acceleration factors, the rate at which the coordinator performs checkpointing operations increases quadratically. Many checkpointing operations are performed by replacing the coordinator's copy by the current work unit (if the copy wasn't modified remotely).

* Coordinator operations:

<figure><img src="image_29.png" /><figcaption>Figure 7.6: Exploitation of coordinator process. Lines show the time (in seconds) and bars the percentage of total elapsed time.</figcaption></figure>

Figure 7.6 focuses on the activity of the coordinator, which is split into four parts measured separately. The absolute time (in seconds) spent in these parts is represented by solid lines. Stacked bars represent this time as a percentage of the total
elapsed time. Messages are sent and received in the form of a sequence of characters (*stringstream*), which must be written/read to/from valid work units. This manipulation of the send-buffer consumes 65-90% of the coordinator's processing time (decreasing with the number of GPUs). This overhead is significant and should be reduced. Rewriting communication routines with MPI, instead of socket programming, in order to take advantage of optimized derived data-types [$S^+03$] may be a necessary modification of PB&B@CLUSTER. As the number of GPUs increases, the coordinator spends a greater portion of time actually treating the requests. For 36 GPUs, work unit intersection, division and other checkpoint operations (e. g. periodically saving all work units to the disk) consume each 10-12% of the coordinator's time. In total, the coordinator is exploited 40% of the time when 36 GPUs are used. At this rate it is likely that incoming requests from workers are queuing up, causing the coordinator to become a bottleneck.

* Scalability:

<figure><img src="image_30.png" /><figcaption>Figure 7.7: Average node processing rate and GPU efficiency</figcaption></figure>

In order to evaluate the impact of the instance size on the efficiency of PB&B@CLUSTER the scaling experiment is repeated with instances *Ta053-h*, with $h \in \{1, 2, 3, 4, 5\}$. The obtained node processing rates is shown (in Mn/s) in Figure 7.7b and GPU efficiencies in Figure 7.7a. GPU efficiency is defined analogous to the conventional parallel efficiency definition, replacing processors with GPUs. The smallest of these instances, *Ta053-1* is solved in 25 minutes on a single GPU and the largest, *Ta053-5*, lasts for 220 minutes on a single device.

Unsurprisingly, efficiency and node processing rates increase as the size of the explored tree increases. In order to exploit all 36 available GPUs efficiently (> 70%), the flowshop instance to be solved should at least require $10^9$ node branchings (*Ta053-3*,
requiring 2 hours of processing on a single GPU). For smaller instances, the repartition of the search space, represented by the interval [0, N![, among 36 × 16 384 ≈ 600 000 IVMs incurs too much overhead. Even for instance Ta053-3, solved in 4.5 minutes on 36 GPUs, each IVM branches on average only ≈ 1 500 permutations.

### 7.4.2 Experiments on a cluster of GPUs and CPUs

#### Experimental protocol

* **Hardware testbed:** In order to evaluate the scalability of PB&B@CLUSTER with GPUs and CPUs further experiments are performed on the Kepler cluster located at the Mathematics and Operations Research Department at UMONS University. Kepler is a cluster of 20 low-power system-on-a-chip (SoC) devices. Each of the 20 nodes is a Nvidia Tegra K1 SoC featuring a 32-bit quad-core ARM Cortex-A15 CPU and a Kepler GK20A GPU containing 192 CUDA cores. Tegra K1 is primarily designed for graphics intensive mobile applications, like gaming, and is used in different tablet computers. Therefore, battery lifetime is a main design objective and Tegra K1 has a TDP of less than 10 W.

* **Problem instance:** As the Kepler cluster is not very powerful, one of the 20 jobs and 20 machines instances is chosen to validate the PB&B@CLUSTER approach. The selected instance is Ta022 which can be solved in 196 minutes on a CPU core.

* **GPU-CPU efficiency definition:** We evaluate the scalability of PB&B@CLUSTER as follows. For a given problem instance, the node processing rate (in nodes/sec) on a single CPU (resp. GPU) is measured. Based on these rates the expected node processing rate on a system using n GPUs and m CPUs is deduced. The efficiency on a (n GPU+m CPU) system is expressed as a percentage of this achieved rate.

Formally, let $\alpha$ and $\beta$ be the node processing rates achieved by a single CPU (resp. GPU), and let $\tau_{m,n}$ be the processing speed measured for a system composed of m CPUs and n GPUs. The efficiency $\eta_{m,n}$ for this configuration is computed as

$$ \eta_{m,n} = \frac{\tau_{m,n}}{m\alpha + n\beta} \times 100\% $$

For instance, solving Ta022 on a single CPU (resp. GPU) a node processing rate of $\alpha = 4.90$ kn/s (resp. 14.12 kn/s) is achieved. Using 8 GPUs and 4 CPUs, the same instance is solved with an average node processing rate of 131.6 kn/s. Supposing linear scalability with CPUs and GPUs one can expect to achieve a node processing rate of $8 \times 14.12 + 4 \times 4.90 = 132.6$ kn/s. In this case, PB&B@CLUSTER reaches an efficiency of $\eta_{4,8} = 99.2\%$.
Obtained results

Table 7.1 reports the efficiency $\eta_{m,n}$ achieved for $m, n=4, 8, 12, 16, 20$ solving flowshop instance *Ta022*. The results shown in this table are averages over 5 independent runs and the relative standard deviation (RSD) is shown on the right-hand side. The initial runs on one CPU (resp. GPU) were also performed 5 times. The two tables on top (Tables 7.1a and 7.1b) show results using 128 IVMs per GPU, the two bottom tables (Tables (7.1c and 7.1d) show results for 1024 IVMs per GPU. Each CPU PB&B process is a 4-threaded (4-IVM) exploration process. The coordinator process runs on a reserved node (without concurrent exploration processes) except for the runs where \#CPUs=20 or \#GPUs=20. The node processing rates for individual workers are the following:

• (1) 4.90 kn/s for one CPU worker,

• (2) 14.12 kn/s for one GPU worker with $T = 128$

• and (3) 22.10 kn/s for one GPU worker with $T = 1024$.

For reference, the sequential processing rate on a Intel E5-2630v3 CPU is about 1.93 kn/s.
For $T = 128$, resp. $T = 1024$ IVMs, the maximal rate of the hybrid system is therefore
380 kn/s, resp. 540 kn/s. Using all 20 nodes, the achieved processing rates is 293 kn/s
($\eta_{20,20} = 77\%$), resp. 373 kn/s ($\eta_{20,20} = 69\%$), meaning that *Ta022* is solved in 75, resp.
59 seconds.

For all configurations less than 2% of nodes is explored redundantly. However, one
can observe that execution time variability is significant, especially when the number of
CPUs is high and the number of GPU is low. A more detailed analysis of exceptionally slow
explorations reveals that a high number of work allocations occur in the shutdown phase.
While this indicates a better robustness for larger instances, experimental confirmation
is needed. Although a threshold below which work units are not divided, an efficient
handling of the shutdown phase revealed challenging.

7.4.3 Resolution of a big instance

Experimental protocol

* **Hardware testbed:** The Ouessant cluster, described in Subsection 7.4.1, is used for the experiments of this subsection.

* **Problem instance:** The objective is to solve the instance *Ta056*. To the best of our knowledge, of the 10 Taillard instances defined by 50 jobs and 20 machines (*Ta051-Ta060*), *Ta056* is currently the only one for which the best known complete permutation is proven to be optimal.
<table><caption>(a) $\eta_{m,n}$ - Ta022 - 128 IVM/GPU.</caption><thead><tr><th rowspan="2"></th><th colspan="6">GPU x</th></tr><tr><th>Eff</th><th>0</th><th>4</th><th>8</th><th>12</th><th>16</th><th>20</th></tr></thead><tbody><tr><th rowspan="5">GPU x</th><th>0</th><td></td><td>95</td><td>91</td><td>86</td><td>82</td><td>79</td></tr><tr><th>4</th><td>100</td><td>103</td><td>93</td><td>87</td><td>82</td><td>78</td></tr><tr><th>8</th><td>98</td><td>99</td><td>95</td><td>87</td><td>81</td><td>77</td></tr><tr><th>12</th><td>96</td><td>91</td><td>89</td><td>89</td><td>84</td><td>79</td></tr><tr><th>16</th><td>94</td><td>89</td><td>89</td><td>84</td><td>84</td><td>78</td></tr><tr><th>20</th><td>92</td><td>86</td><td>85</td><td>81</td><td>79</td><td>77</td></tr></tbody></table>

<table><caption>(c) $\eta_{m,n}$ - Ta022 - 1024 IVM/GPU.</caption><thead><tr><th rowspan="2"></th><th colspan="6">GPU x</th></tr><tr><th>Eff</th><th>0</th><th>4</th><th>8</th><th>12</th><th>16</th><th>20</th></tr></thead><tbody><tr><th rowspan="5">GPU x</th><th>0</th><td></td><td>96</td><td>90</td><td>86</td><td>81</td><td>81</td></tr><tr><th>4</th><td>100</td><td>88</td><td>87</td><td>83</td><td>78</td><td>75</td></tr><tr><th>8</th><td>102</td><td>84</td><td>83</td><td>82</td><td>78</td><td>73</td></tr><tr><th>12</th><td>100</td><td>79</td><td>88</td><td>76</td><td>75</td><td>71</td></tr><tr><th>16</th><td>101</td><td>77</td><td>87</td><td>77</td><td>78</td><td>72</td></tr><tr><th>20</th><td>99</td><td>83</td><td>77</td><td>79</td><td>73</td><td>69</td></tr></tbody></table>

<table><caption>(b) RSD - Ta022 - 128 IVM/GPU.</caption><thead><tr><th rowspan="2"></th><th colspan="6">GPU x</th></tr><tr><th>Eff</th><th>0</th><th>4</th><th>8</th><th>12</th><th>16</th><th>20</th></tr></thead><tbody><tr><th rowspan="5">GPU x</th><th>0</th><td></td><td>0.1</td><td>0.8</td><td>0.6</td><td>0.7</td><td>2.3</td></tr><tr><th>4</th><td>0.1</td><td>3.4</td><td>0.2</td><td>2.0</td><td>1.8</td><td>1.2</td></tr><tr><th>8</th><td>0.2</td><td>4.3</td><td>5.0</td><td>2.5</td><td>4.2</td><td>1.8</td></tr><tr><th>12</th><td>0.4</td><td>8.9</td><td>6.9</td><td>3.4</td><td>6.9</td><td>5.2</td></tr><tr><th>16</th><td>0.2</td><td>11.4</td><td>7.7</td><td>8.9</td><td>7.5</td><td>4.5</td></tr><tr><th>20</th><td>0.7</td><td>6.5</td><td>11.2</td><td>3.6</td><td>9.2</td><td>5.9</td></tr></tbody></table>

<table><caption>(d) RSD - Ta022 - 1024 IVM/GPU.</caption><thead><tr><th rowspan="2"></th><th colspan="6">GPU x</th></tr><tr><th>Eff</th><th>0</th><th>4</th><th>8</th><th>12</th><th>16</th><th>20</th></tr></thead><tbody><tr><th rowspan="5">GPU x</th><th>0</th><td></td><td>0.8</td><td>1.3</td><td>0.9</td><td>1.1</td><td>0.9</td></tr><tr><th>4</th><td>0.1</td><td>2.1</td><td>2.9</td><td>2.4</td><td>2.3</td><td>1.5</td></tr><tr><th>8</th><td>0.2</td><td>2.5</td><td>7.7</td><td>2.0</td><td>4.7</td><td>4.2</td></tr><tr><th>12</th><td>0.4</td><td>1.0</td><td>1.5</td><td>4.5</td><td>4.6</td><td>4.7</td></tr><tr><th>16</th><td>0.2</td><td>9.6</td><td>3.7</td><td>4.7</td><td>3.8</td><td>2.3</td></tr><tr><th>20</th><td>0.7</td><td>13.2</td><td>11.7</td><td>3.1</td><td>8.0</td><td>3.0</td></tr></tbody></table>

Table 7.1: Mixed GPU-CPU efficiency and Relative Standard Deviation (RSD) for resolution of flowshop instance Ta022 (22.1×10⁶ nodes) using 20 Tegra K1. The upper (resp. lower) row shows results for *T* = 128 (resp. *T* = 1024) IVMs/GPU. Average efficiency and RSD over 5 runs.

## Obtained results

### Elapsed time:

<table><caption>Table 7.2: Exploration statistics for resolution of flowshop instance Ta056.</caption><thead><tr><th></th><th>Serial PB&B</th><th>B&B@Grid</th><th>PB&B@CLUSTER</th></tr></thead><tbody><tr><td>#GPUs</td><td>0</td><td>0</td><td>9 × 4 = 36</td></tr><tr><td>#CPU cores</td><td>1</td><td>Aver. 328</td><td>9 × 2 = 18</td></tr><tr><td>Elapsed time</td><td>22 years</td><td>25 days</td><td>9 hours</td></tr></tbody></table>

The optimal complete permutation of Ta056 was found and proven in 2006 using B&B@Grid [MMT07c]. The resolution required 25 days of processing, exploiting on average 328 CPU cores distributed on 9 clusters of the French experimental testbed Grid'5000 <fnref n="2" />. This result is used as a reference for the three resolutions of Ta056 which are performed under identical initial conditions. As in the B&B@Grid experiment reported in [MMT07c] the initial upper bound is set to the optimal cost plus one

<footnote n="2">²https://www.grid5000.fr/</footnote>
unit, i.e. 3680, which allows one to verify the correctness of the algorithm. Our new resolution of *Ta056* using PB&B@CLUSTER found the same optimal cost (3679), and the same optimal complete permutation.

On Ouessant and using PB&B@CLUSTER, *Ta056* is solved in 9 hours. Compared with the B&B@Grid resolution, the execution time is reduced by a factor of about 65×. Compared to the estimated sequential execution time of 22 years, the execution time is reduced by a factor of about 20,000×.

• Energy consumption:

<table><caption>Table 7.3: Exploration statistics for resolution of flowshop instance *Ta056*.</caption><thead><tr><th></th><th>B&B@Grid</th><th>PB&B@CLUSTER</th></tr></thead><tbody><tr><td>#GPUs</td><td>0</td><td>9 × 4 = 36</td></tr><tr><td>TDP of a GPU</td><td>-</td><td>300 W</td></tr><tr><td>#CPUs</td><td>Aver. 328</td><td>9 × 2 = 18</td></tr><tr><td>TDP of a CPU</td><td>30 W +</td><td>225 W</td></tr><tr><td>Elapsed time</td><td>25 days</td><td>9 hours</td></tr><tr><td>Energy consumption</td><td>≈6000 kWh</td><td>≈130 kWh</td></tr></tbody></table>

Table 7.3 indicates an approximate value for the energy consumption of each resolution. These values are based on the Thermal Design Power (TDP) of CPUs and GPUs, as listed by the respective vendors. For example, the GTX 980 GPU is listed with a TDP of 165 W and the host Xeon CPUs with 85 W, so an indicative value for the energy consumption is given by $(36 \times 300 + 18 \times 225)W \times 9h \approx 130 \text{ kWh}$.

For the 2006 resolution using B&B@Grid the energy consumption can only be roughly estimated. About $\frac{2}{3}$ of CPU cores in the computational pool exploited by B&B@Grid in 2006 are AMD Opteron dual-core CPUs, 90 nm feature size, with clock rates between 2.0 and 2.2 GHz. The remaining $\frac{1}{3}$ are Intel Pentium 4 and Celeron single-core processors with similar clock rates. The most energy efficient models of this type of CPUs are listed with TDP values above 30 W. Taking into account that most CPUs are dual-core, an optimistic estimation for the energy consumption is $328 \text{ CPU cores} \times 30W \times 25d \times 24h/d \approx 6000 \text{ kWh}$.

• Load balancing:

Figure 7.8 illustrates the workload repartition among GPU and CPU based workers, in terms of branched nodes. In addition to four GPU workers, on each Ouessant node one multi-threaded CPU based worker with 160 IVM is used $(2 \times 10 \text{ cores} \times 8 \text{ threads})$.
In this configuration, each multi-core PB&B process branches on average about $\frac{1}{10}$ the amount of nodes branched by an average GPU PB&B. Using 4 times as many GPU based workers as CPU workers, in total less than 3% of node branchings are performed on a CPU. One can see in Figure 7.8 that workers of the same type perform a roughly equal amount of work. However, a node branching represents a variable amount of work. Therefore, the number of branched nodes is only an approximative indicator for load balancing.

<table><caption>Table 7.4: Exploration statistics for resolution of flowshop instance *Ta056*.</caption><thead><tr><th></th><th>PB&B@CLUSTER</th></tr></thead><tbody><tr><td>Branched nodes</td><td>175.8 × 10<sup>9</sup></td></tr><tr><td>T<sub>coordinator</sub></td><td>38.2 min</td></tr><tr><td>Coord. exploitation (%)</td><td>7.1%</td></tr><tr><td>#checkpoints</td><td>3 568 368</td></tr><tr><td>#work allocations</td><td>33 387</td></tr></tbody></table>

<figure><img src="image_31.png" /><figcaption>Figure 7.8: Number of nodes branched per worker. Resolution of *Ta056* on 9 Minsky nodes (Ouessant), using 9 × (4 GPU based workers with 16 384 IVM + 1 CPU based worker with 160 IVM).</figcaption></figure>

Another indicator for the load imbalance throughout the execution is the number of work allocations (work stealing operations). A worker only requests new work from the coordinator when no more work is locally available. The results indicate that workers run out of work very rarely.

These metrics change significantly for the hybrid resolution on Ouessant. In order to balance the workload between 36 GPUs and 9 multi-core CPUs, more than 33 000 work unit allocations and 3 500 000 checkpointing operations are performed. In that case the farmer processor is exploited 7% of the time, i. e. 38 minutes out of 9 hours.
7.5 Conclusions

In this chapter we revisited the design of B&B@Grid to enable the integration of GPU-
and CPU-based workers. When B&B@Grid was designed in mid-2000 most components of
computational grids were single-core and dual-core CPUs. Today, more than 15 years later,
large-scale HPC platforms are becoming increasingly heterogeneous, integrating GPUs and
CPUs with larger core-counts.

The extension of B&B@Grid to PB&B@CLUSTER, the proposed PB&B for hybrid dis-
tributed HPC clusters, includes the redefinition of work units and a modification of the
communication scheme to allow asynchronous checkpointing operations that overlap with
worker computations.

A very large flowshop instance defined by 50 jobs and 20 machines, *Ta056*, was successfully solved on GPU powered clusters with a total of up to 130 000 GPU cores. A first resolution of this instance was performed in 2006, using B&B@Grid to exploit on average 328 CPU cores in a computational grid during 25 days. Using PB&B@CLUSTER the resolution of *Ta056* was performed in 9 hours on a cluster composed of 36 GPUs.

Low exploitation rates of the coordinator process and experiments performed with smaller instances are indicators for the good scalability of PB&B@CLUSTER. For a set of 50-job flowshop instances requiring 3.7 hours of computation on a single GPU, a relative speedup of 30× is achieved on 36 GPUs, solving this instance in 7.5 minutes on average.

PB&B@CLUSTER was also experimented on a mini-cluster of 20 systems-on-a-chip designed for mobile devices. Experimental results show that the availability of efficient data types for inter-node communication is a key component for the performance of the central coordinator process.
# Conclusions and perspectives

The research works, presented in this manuscript, are focused on the exact resolution of a certain type of economic problems, where the objective is to find the best permutation. A large number of methods, often serial, have been developed to solve these problems. A wide range of exact methods can be classified as PB&B algorithms. The solution space is explored by dynamically building a tree, whose size in terms of number of nodes is huge. Therefore, solving these problems often requires the use of a high performance computing system.

This last chapter recalls the main contributions, presented in this manuscript, and the conclusions to draw from these contributions. In addition, the chapter gives an overview of our other contributions in PB&B algorithms which are not described in this manuscript. And finally, the chapter indicates our main future works in the field PB&B algorithms.

## 8.1 Main contributions

Our works led to the development of four approaches, renamed in this manuscript PB&B@CORE, PB&B@CPU, PB&B@GPU and PB&B@CLUSTER. The first two methods are developed during the thesis of Rudi Leroy, and the two others in the thesis of Jan Gmys.

* The PB&B@CORE approach rethinks the conventional data structure used in PB&B. When exploring the PB&B tree, the nodes, generated but not yet processed, are often stored in a linked-list. For a single CPU core, we have developed a new data structure, called IVM, to replace the linked-list. Compared to the linked-list, the experiments, performed for the resolution of the 50 × 10 (i.e. 50 jobs and 10 machines) flowshop problem instances [Tai93], show that the management of an IVM pool requires on average ≈ 3 times less memory and ≈ 35 times less computing power.

* The PB&B@CPU approach, obtained using the IVM structure of PB&B@CORE, allows the parallelization of the PB&B on multi-core CPU processors. For such architectures, and unlike conventional parallelization, where the work unit is often a
set of nodes, the work unit of our approach is an interval of factoradics. A factoradic is an integer encoded in the factorial enumeration system. The experiments performed on a multi-core CPU accelerator, for the resolution of the ten 20 × 20 flowshop instances, show that the management of an IVM pool is on average ≈ 10 times faster than the management of linked-list pool.

* The PB&B@GPU approach, based on the factoradic intervals of PB&B@CPU, allows the parallelization of PB&B on GPU, which is a processor with a SIMD architecture. Previously, there have been many attempts in the literature to run the most time-consuming part (hot spot) of the PB&B on GPU and the rest of the algorithm on CPU. The advantage of our approach, based on a revisited work stealing method, is to offload the whole execution of PB&B, whose tree has an irregular structure, to the GPU accelerator, which has a regular SIMD architecture. When solving the ten 20 × 20 flowshop instances, experiments show that on average the resolution time is reduced on average by ≈ 3.3 times compared to a GPU linked-list PB&B.

* The PB&B@CLUSTER approach, which is based on the three previous approaches, allows the deployment of a PB&B on a cluster of CPU processors and GPU accelerators. For a distributed memory architecture, such as a cluster, our approach uses and redefines a certain number of basic operators known in set theory. This approach is validated on Ta056, which is a difficult 50 × 20 flowshop instance. The resolution of this instance required ≈ 22 days of computing using on average 328 CPU cores and a total of ≈ 6000 kWh[MMT07a]. Our new approach solved this instance in ≈ 9 hours with an energy consumption of ≈ 130 kWh.

The PB&B@CLUSTER approach is not yet published while the others are published mainly in four journals and one conference, namely two CCPE [GMMT17, GLM+16], one Parallel computing [GMMT16] and one FGCS [MGMT18a], and the 28th IEEE IPDPS conference [MLMT14a]. The PB&B@GPU approach is also awarded with the Best Paper Award at the 11th PPAM conference [GMM+15]. As future research directions for this work, we have identified some challenging perspectives summarized in the following:

## 8.2 Other contributions

In addition to these four contributions, we have several other contributions in PB&B algorithms. Some of these contributions are important, but are not presented in this manuscript. In particular, we would like to describe two of them.
* **Efficient branching and bounding operations:** In this contribution, we present a new node decomposition scheme that combines dynamic branching and lower bound refinement strategies in a computationally efficient way. To alleviate the computational burden of the two-machine bound used in the refinement stage, we propose an online learning-inspired mechanism to predict promising couples of bottleneck machines. The algorithm offers multiple choices for branching and bounding operations and can explore the search tree either sequentially or in parallel on multi-core CPUs.

In order to empirically determine the most efficient combination of these components, a series of computational experiments with 600 flowshop benchmark instances is performed. A main insight is that the problem size, as well as interactions between branching and bounding operations substantially modify the trade-off between the computational requirements of a lower bound and the achieved tree size reduction. Moreover, we demonstrate that parallel tree search is a key ingredient for the resolution of large problem instances, as strong super-linear speedups can be observed. An overall evaluation using two well-known benchmarks indicates that the proposed approach is superior to previously published PB&B algorithms.

For the first benchmark we report the exact resolution - within less than 20 minutes - of two instances defined by 500 jobs and 20 machines that remained open since more than 25 years, and for the second a total of 88 improved best known upper bounds, including proofs of optimality for 71 of them. This contribution is under a minor revision in EJOR journal.

* **PB&B on many-core processors:** On the road to exascale, coprocessors are increasingly becoming key building blocks of high performance computing platforms. In addition to their energy efficiency, these many-core devices boost the performance of multi-core processors. In this contribution, we revisit the design and implementation of PB&B algorithms for multi-core processors and Intel Xeon Phi coprocessors considering the offload mode as well as the native one. In addition, two major parallel models are considered: the master-worker and the work pool models. We address several parallel computing issues including processor-coprocessor data transfer optimization and vectorization.

The proposed approaches have been experimented using the flowshop and two hardware configurations equivalent in terms of energy consumption: Intel Xeon E5-2670 processor and Intel Xeon Phi 5110P coprocessor. The reported results show that: (1) the proposed vectorization mechanism reduces the execution time by (resp.) in the many-core (resp. multi-core) approach; (2) the offload mode allows a faster
execution on MIC than the native mode for most flowshop problem instances; (3) the many-core approach (offload or native) is in average twice faster than the multi-core approach; (4) the work pool parallel model is more suited for many/multi-core PB&B applied to flowshop than the master-worker model because of its irregular nature. This contribution is published in the FGCS journal [MGMT18b].

## 8.3 Perspectives

The contributions presented in this document open up new perspectives. These perspectives are mainly related to the validation of our approaches using other permutation problems, and a cluster powered with a larger number GPU nodes. Our future works can be summarized with the following points:

* As a short-term future work, we will validate our approaches on other single permutation permutation problems such as TSP, Job-Shop and QAP, and permutation problems with more than one permutation like Q3AP [Meh11].

* The experimental results, obtained with PB&B@CLUSTER, indicate that the algorithm is scalable on larger GPU-enhanced clusters. We plan to validate this by attempting the resolution of previously unsolved flowshop instances on a large GPU-powered supercomputer. To resolve these instances, we obtained 50,000 hours of computing on the 1,000 GPU nodes of Jean Zay supercomputer. Jean Zay is the converged platform acquired by the French Ministry of Higher Education, Research and Innovation through the intermediary of the French civil company, GENCI (Grand Equipement National De Calcul Intensif). The Jean Zay computer was installed at IDRIS, national computing centre for the French National Centre for Scientific Research (CNRS), in 2019.

* In order to further improve scalability of the approach, specially in an exascale environment [CM19], the coordinator process, usually running on a multi-core CPU, should be parallelized. Also, the checkpointing mechanism should be revisited. PB&B@CLUSTER uses the checkpointing mechanism inherited from B&B@Grid, making the approach tolerant against node failures. However, as a large portion is shifted to lower levels it becomes important to make the approach fault-tolerant against failures at the GPU and multi-core CPU level.

* The IVM data structure revealed itself particularly well suited for fine-grained permutation problems. For example, sampling methods based on the optimization of
latin hypercubes can be modeled as (multi-)permutation problems. As a future re-
search direction we plan to revisit the IVM-based algorithm to enable the resolution
of multi-permutation problems.

* Experimental results have shown strong performance variations according to the used node evaluation function. Having different bounds for the same problem matching implementations with underlying hardware. A challenging improvement of the PB&B algorithm consists in implementing a library of lower bounds for the same problem, in order to enable the different workers to use the node evaluation function which is the best fit for the underlying hardware.
Bibliography

[ABEB<sup>+</sup>16] Dabah Adel, Ahcène Bendjoudi, Didier El-Baz, Ait Zai Abdelhakim, et al. Gpu-based two level parallel b&b for the blocking job shop scheduling problem. 2016.

[ABGL02] Kurt Anstreicher, Nathan Brixius, Jean-Pierre Goux, and Jeff Linderoth. Solving large quadratic assignment problems on computational grids. *Mathe-matical Programming*, 91(3):563–588, 2002.

[ACR] ACRO. A common repository for optimizers, sandia national laboratories. https://software.sandia.gov/trac/acro/. Accessed: 2017-09-23.

[AGA99] A. Allahverdi, J.N.D Gupta, and T. Aldowaisan. A review of scheduling research involving setup considerations. *Omega*, 27(2):219–239, 1999.

[AMTM17] E. Alekseeva, M. Mezmaz, D. Tuyttens, and N. Melab. Parallel multi-core hyper-heuristic GRASP to solve permutation flow-shop problem. *Concur-rency and Computation: Practice and Experience*, 29(9), 2017.

[ASW<sup>+</sup>14] Tae-Hyuk Ahn, Adrian Sandu, Layne T Watson, Clifford A Shaffer, Yang Cao, and William T Baumann. A framework to analyze the performance of load balancing schemes for ensembles of stochastic simulations. *International Journal of Parallel Programming*, 43(4):597–630, 2014.

[Bas05] Matthieu Basseur. Conception de métaheuristiques coopératives pour l'optimisation multi-objective: Application au flowshop. *PhD Thesis from Université Lille 1*, 2005.

[BG76] M.C. Bonney and S.W. Gundry. Solutions to the constrained flowshop sequencing problem. *Operational Research Quarterly*, page 869, 1976.

[BHP05] David A Bader, William E Hart, and Cynthia A Phillips. Parallel algorithm design for branch and bound. *Tutorials on Emerging Methodologies and Applications in Operations Research: Presented at INFORMS 2004, Denver, CO*, 76, 2005.

[BL99] Robert D Blumofe and Charles E Leiserson. Scheduling multithreaded computations by work stealing. *Journal of the ACM (JACM)*, 46(5):720–748, 1999.
[BMMT18] Guillaume Briffoteaux, Nouredine Melab, Mohand Mezmaz, and Daniel Tuyttens. An adaptive evolution control based on confident regions for surrogate-assisted optimization. In *HPCS 2018 - International Conference on High Performance Computing & Simulation*, Orléans, France, July 2018.

[BMT12a] A. Bendjoudi, N. Melab, and E. G. Talbi. Hierarchical branch and bound algorithm for computational grids. *Future Gener. Comput. Syst.*, 28(8):1168–1176, October 2012.

[BMT12b] Ahcène Bendjoudi, Nordine Melab, and El-Ghazali Talbi. An adaptive hierarchical master-worker (ahmw) framework for grids—application to b&b algorithms. *Journal of Parallel and Distributed Computing*, 72(2):120–131, 2012.

[BMT14] Ahcène Bendjoudi, Nouredine Melab, and El-Ghazali Talbi. Fth-b&b: A fault-tolerant hierarchical branch and bound for large scale unreliable environments. *IEEE Transactions on Computers*, 63(9):2302–2315, 2014.

[Bob] Bobpp. Bobpp framework, université de versailles. http://www.prism.uvsq.fr/~blec/bobpp/main.html. Accessed: 2017-09-23.

[Cha13] I. Chakroun. *Parallel heterogeneous Branch and Bound algorithms for multi-core and multi-GPU environments*. PhD thesis, Université Lille 1, 2013.

[CM13] I. Chakroun and N. Melab. Operator-level gpu-accelerated branch and bound algorithms. In *Proceedings of the International Conference on Computational Science, ICCS 2013, Barcelona, Spain, 5-7 June, 2013*, pages 280–289, 2013.

[CM19] Tiago Carneiro and Nouredine Melab. Productivity-aware design and implementation of distributed tree-based search algorithms. In João M. F. Rodrigues, Pedro J. S. Cardoso, Jânio Monteiro, Roberto Lam, Valeria V. Krzhizhanovskaya, Michael H. Lees, Jack J. Dongarra, and Peter M.A. Sloot, editors, *Computational Science – ICCS 2019*, pages 253–266, Cham, 2019. Springer International Publishing.

[CMGH08] Leocadio G Casado, JA Martinez, Inmaculada García, and Eligius MT Hendrix. Branch-and-bound interval global optimization on shared memory multiprocessors. *Optimization Methods & Software*, 23(5):689–701, 2008.
[CMMB13] I. Chakroun, M. Mezmaz, N. Melab, and A. Bendjoudi. Reducing thread divergence in a GPU-accelerated branch-and-bound algorithm. *Concurrency and Computation: Practice and Experience*, 25(8):1121–1136, 2013.

[CMMT13] Imen Chakroun, Nordine Melab, Mohand Mezmaz, and Daniel Tuyttens. Combining multi-core and gpu computing for solving combinatorial optimization problems. *Journal of Parallel and Distributed Computing*, 73(12):1563–1577, 2013.

[CMNLdC11] T. Carneiro, A.E. Muritiba, M. Negreiros, and G.A. Lima de Campos. A New Parallel Schema for Branch-and-Bound Algorithms Using GPGPU. In *23rd International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD)*, pages 41–47, 2011.

[CNNdC12] Tiago Carneiro, Ricardo Holanda Nobre, Marcos Negreiros, and Gustavo Augusto Lima de Campos. Depth-first search versus Jurema search on GPU branch-and-bound algorithms: A case study. *NVIDIA's GCDF - GPU Computing Developer Forum on XXXII Congresso da Sociedade Brasileira de Computação (CSBC)*, 2012.

[CO] COIN-OR. Computational infrastructure for operations research. https://www.coin-or.org/documentation.html. Accessed: 2017-09-23.

[Cra06] Teodor Gabriel Crainic. Parallel branch-and-bound algorithms. *Parallel combinatorial optimization*, 1:1–28, 2006.

[CZ06] S. Climer and W. Zhang. Cut-and-solve: an iterative search strategy for combinatorial optimization problems, artificial intelligence. 170:714–738, 2006.

[EHP15] Jonathan Eckstein, William E Hart, and Cynthia A Phillips. Pebbl: an object-oriented framework for scalable parallel branch and bound. *Mathematical Programming Computation*, 7(4):429–469, 2015.

[EPH00] J. Eckstein, C. A. Phillips, and W. E. Hart. PICO: An object-oriented framework for parallel branch-and-bound. Research Report 40–2000, RUTCOR, 2000.

[EPH01] Jonathan Eckstein, Cynthia A Phillips, and William E Hart. Pico: An object-oriented framework for parallel branch and bound. *Studies in Computational Mathematics*, 8:219–265, 2001.
[EPS09a] Y. Evtushenko, M. Posypkin, and I. Sigal. A framework for parallel large-scale global optimization. *Computer Science-Research and Development*, 23(3):211–215, 2009.

[EPS09b] Yuri Evtushenko, Mikhail Posypkin, and Israel Sigal. A framework for parallel large-scale global optimization. *Computer Science-Research and Development*, 23(3-4):211–215, 2009.

[FRvLP10] Frank Feinbube, Bernhard Rabe, M von Löis, and Andreas Polze. Nqueens on cuda: Optimization issues. In *Parallel and Distributed Computing (ISPDC), 2010 Ninth International Symposium on*, pages 63–70. IEEE, 2010.

[GC94] B. Gendron and T.G. Crainic. Parallel Branch and Bound Algorithms: Survey and Synthesis. *Operations Research*, 42:1042–1066, 1994.

[GGS04] B. Goldengorin, D. Ghosh, and G. Sierksma. Branch and bound algorithms for the simple plant location problem. *Computers & Operations Research*, 31:241–255, 2004.

[GJS76] M.R. Garey, D.S. Johnson, and R. Sethi. The complexity of flow-shop and job-shop scheduling. *Mathematics of Operations Research*, 1:117–129, 1976.

[GLM+16] J. Gmys, R. Leroy, M. Mezmaz, N. Melab, and D. Tuyttens. Work stealing with private integer-vector-matrix data structure for multi-core branch-and-bound algorithms. *Concurrency and Computation: Practice and Experience*, 28(18):4463–4484, 2016.

[GMM+15] J. Gmys, M. Mezmaz, N. Melab, R. Leroy, and D. Tuyttens. IVM-based Work Stealing for Parallel Branch-and-Bound on GPU. In *Proc. of Int. Conference on Parallel Processing and Applied Mathematics*, Krakow, Poland, September 2015.

[GMMT16] J. Gmys, M. Mezmaz, N. Melab, and D. Tuyttens. A gpu-based branch-and-bound algorithm using integer-vector-matrix data structure. *Parallel Computing*, 59:119–139, 2016.

[GMMT17] J. Gmys, M. Mezmaz, N. Melab, and D. Tuyttens. Ivm-based parallel branch-and-bound using hierarchical work stealing on multi-gpu systems. *Concurrency and Computation: Practice and Experience*, 29(9), 2017.

[Gmy17] J. Gmys. *Heterogeneous Cluster Computing for Many-Task Exact Optimization: Application to Permutation Problems*. PhD thesis, Université de Mons et Université de Lille, 2017.
[Gol66] S. Golomb. Run-length encodings (corresp.). *Information Theory, IEEE Transactions on*, 12(3):399-401, Jul 1966.

[HS05] S. Reza Hejazi and S. Saghafian. Flowshop-scheduling problems with makespan criterion: a review. *International Journal of Production Research*, 43(14):2895-2929, 2005.

[HSH$^{+}$17] Juan F. R. Herrera, José M. G. Salmerón, Eligius M. T. Hendrix, Rafael Asenjo, and Leocadio G. Casado. On parallel branch and bound frameworks for global optimization. *Journal of Global Optimization*, Mar 2017.

[HSO07] Mark Harris, Shubhabrata Sengupta, and John D. Owens. Parallel prefix sum (scan) with cuda. In Hubert Nguyen, editor, *GPU Gems 3*. Addison Wesley, August 2007.

[IF00] A. Iamnitchi and I. Foster. A Problem-Specific Fault-Tolerance Mechanism for Asynchronous, Distributed Systems. 29th International Conference on Parallel Processing (ICPP), Toronto, Canada, August, pages 21-24, 2000.

[JAO$^{+}$11] John Jenkins, Isha Arkatkar, John D. Owens, Alok Choudhary, and Nagiza F. Samatova. Lessons learned from exploring the backtracking paradigm on the gpu. In *Proceedings of the 17th International Conference on Parallel Processing - Volume Part II*, Euro-Par'11, pages 425-437, Berlin, Heidelberg, 2011. Springer-Verlag.

[Joh54] S.M. Johnson. Optimal two and three-stage production schedules with setup times included. *Naval Research Logistis Quarterly*, 1:61-68, 1954.

[Knu97] D.E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. *Reading, Ma*, page 192, 1997. ISBN=9780201896848.

[Knu98] Donald E Knuth. The art of computer programming, volume 2: Seminumerical algorithms, addison-wesley. *Reading, Massachusetts*, 1998.

[KRR88] Vipin Kumar, V. N Rao, and K. Ramesh. Parallel depth first search on the ring architecture. Technical report, Austin, TX, USA, 1988.

[KS80] J. R. King and A. S. Spachis. Heuristics for flow-shop scheduling. *International Journal of Production Research*, 18(3):345-357, 1980.

[Lai88a] C-A. Laisant. Sur la numération factorielle, application aux permutations. *Bulletin de la Société Mathématique de France*, 16:176-183, 1888.
[Lai88b] C-A. Laisant. Sur la numération factorielle, application aux permutations. *Bulletin de la Société Mathématique de France*, 16:176–183, 1888.

[Lai88c] C-A Laisant. Sur la numération factorielle, application aux permutations. *Bulletin de la Société Mathématique de France*, 16:176–183, 1888.

[LEB12] M.E. Lalami and D. El-Baz. GPU Implementation of the Branch and Bound Method for Knapsack Problems. In *IEEE 26th Intl. Parallel and Distributed Processing Symp. Workshops PhD Forum (IPDPSW)*, pages 1769–1777, Shanghai, CHN, May 2012.

[Ler15a] R. Leroy. *Parallel Tree-based Exact Algorithms using Heterogeneous Many and Multi-core Computing for Solving Challenging Problems in Combinatorial Optimization*. PhD thesis, Université de Lille, 2015.

[Ler15b] Rudi Leroy. *Parallel Branch-and-Bound revisited for solving permutation combinatorial optimization problems on multi-core processors and coprocessors*. PhD thesis, Université Lille 1, 2015.

[LLK78] J.K. Lenstra, B.J. Lageweg, and A.H.G. Rinnooy Kan. A General bounding scheme for the permutation flow-shop problem. *Operations Research*, 26(1):53–67, 1978.

[LLW$^+$15] Liang Li, Hong Liu, Hao Wang, Taoying Liu, and Wei Li. A parallel algorithm for game tree search using GPGPU. *Parallel and Distributed Systems, IEEE Transactions on*, 26(8):2114–2127, 2015.

[MCA13] X. Meyer, B. Chopard, and P. Albuquerque. A branch-and-bound algorithm using multiple gpu-based lp solvers. In *20th Annual International Conference on High Performance Computing*, pages 129–138, Dec 2013.

[MCB14] Nouredine Melab, Imen Chakroun, and Ahcène Bendjoudi. Graphics processing unit-accelerated bounding for branch-and-bound applied to a permutation problem using data access optimization. *Concurrency and Computation: Practice and Experience*, 26(16):2667–2683, 2014.

[McC03] James McCaffrey. Using permutations in .net for improved systems security. *Microsoft Developer Network*, 2003.

[Meh11] Malika Mehdi. Parallel hybrid optimization methods for permutation based problems. *PhD Thesis from Université Lille 1 and Université du Luxembourg*, 2011.
[Men17] Tarek Menouer. Solving combinatorial problems using a parallel framework. *Journal of Parallel and Distributed Computing*, 2017.

[MGMT18a] Nouredine Melab, Jan Gmys, Mohand Mezmaz, and Daniel Tuyttens. Multi-core versus many-core computing for many-task branch-and-bound applied to big optimization problems. *Future Generation Comp. Syst.*, 82:472–481, 2018.

[MGMT18b] Nouredine Melab, Jan Gmys, Mohand Mezmaz, and Daniel Tuyttens. Multi-core versus many-core computing for many-task branch-and-bound applied to big optimization problems. *Future Generation Computer Systems*, 82:472–481, 2018.

[MLMT14a] M. Mezmaz, R. Leroy, N. Melab, and D. Tuyttens. A Multi-Core Parallel Branch-and-Bound Algorithm Using Factorial Number System. In *Proc. of 28th IEEE Intl. Parallel and Distributed Processing Symp.*, Phoenix, Arizona, USA, May 2014.

[MLMT14b] Mohand Mezmaz, Rudi Leroy, Nouredine Melab, and Daniel Tuyttens. A Multi-Core Parallel Branch-and-Bound Algorithm Using Factorial Number System. In *28th IEEE Intl. Parallel & Distributed Processing Symp. (IPDPS)*, pages 1203–1212, Phoenix, AZ, May 2014.

[MMK$^{+}$11] M. Mezmaz, N. Melab, Y. Kessaci, Y.C. Lee, E-G. Talbi, A.Y. Zomaya, and D. Tuyttens. A parallel bi-objective hybrid metaheuristic for energy-aware scheduling for cloud computing systems. *Journal of Parallel and Distributed Computing*, 71(11):1497 – 1508, 2011. ISSN 0743-7315, DOI: 10.1016/j.jpdc.2011.04.007.

[MMT07a] M. Mezmaz, N. Melab, and E-G. Talbi. A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems. In *21th IEEE Intl. Parallel and Distributed Processing Symp. (IPDPS)*, pages 1–9, Long Beach, CA, March 2007.

[MMT07b] M. Mezmaz, N. Melab, and E-G. Talbi. A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems. In *In Proc. of 21th IEEE Intl. Parallel and Distributed Processing Symp. (IPDPS)*. Long Beach, California, March 2007.

[MMT07c] M. Mezmaz, N. Melab, and E. G. Talbi. A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems.
In 2007 IEEE International Parallel and Distributed Processing Symposium, pages 1–9, March 2007.

[MMT13] M. Mezmaz, N. Melab, and D. Tuyttens. A multithreaded branch-and-bound algorithm for solving the flow-shop problem on a multicore environment, chapter 3, pages 53–70. Large Scale Network-Centric Distributed Systems. John Wiley & Sons, July 2013. ISBN-10: 0470936886, ISBN-13: 978-0470936887.

[PC04] R. Pastor and A. Corominas. Branch and win: Or tree search algorithms for solving combinatorial optimisation problems. Top, 1:169–192, 2004.

[PŽ09] R. Paulavičius and J. Žilinskas. Parallel branch and bound algorithm with combination of lipschitz bounds over multidimensional simplices for multicore computers. Parallel Scientific Computing and Optimization, pages 93–102, 2009.

[RLS04] Ted K Ralphs, Laszlo Ladányi, and Matthew J Saltzman. A library hierarchy for implementing scalable parallel search algorithms. The Journal of Supercomputing, 28(2):215–234, 2004.

[RS10] Kamil Rocki and Reiji Suda. Parallel Minimax Tree Searching on GPU, pages 449–456. Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.

[S⁺03] Xian-He Sun et al. Improving the performance of mpi derived datatypes by optimizing memory-access cost. In Cluster Computing, 2003. Proceedings. 2003 IEEE International Conference on, pages 412–419. IEEE, 2003.

[SECG11] J.F. Sanjuan-Estrada, L.G. Casado, and I. García. Adaptive parallel interval branch and bound algorithms based on their performance for multicore architectures. The Journal of Supercomputing, pages 1–9, 2011.

[SKK⁺11] Vijay. A. Saraswat, Prabhanjan Kambadur, Sreedhar Kodali, David Grove, and Sriram Krishnamoorthy. Lifeline-based global load balancing. In Proc. of the 16th symposium on Principles and practice of parallel programming, PPoPP ’11, pages 201–212, New York, NY, USA, 2011. ACM.

[Tai93] E. Taillard. Benchmarks for basic scheduling problems. Journal of Operational Research, 64:278–285, 1993.
[TFMJ13] D. Tuyttens, H. Fei, M. Mezmaz, and J. Jalwan. Simulation-based Genetic Algorithm towards an Energy-efficient Railway Traffic Control. *Mathematical Problems in Engineering, Hindawi*, 2013:12 pages, 2013. DOI: 10.1155/2013/805410.

[TRDM17] James McCracken Trauer, Romain Ragonnet, Tan Nhut Doan, and Emma Sue McBryde. Modular programming for tuberculosis control, the “autumn” platform. *BMC Infectious Diseases*, 17(1):546, Aug 2017.

[VDM13] Trong-Tuan Vu, Bilel Derbel, and Nouredine Melab. Adaptive Dynamic Load Balancing in Heterogenous Multiple GPUs-CPUs Distributed Setting: Case Study of B&B Tree Search. In *7th International Learning and Intelligent Optimization Conference (LION)*, Catania, Italy, January 2013. Lecture Notes in Computer Science.

[ZSW11] Tao Zhang, Wei Shu, and Min-You Wu. Optimization of n-queens solvers on graphics processors. In *Proceedings of the 9th International Conference on Advanced Parallel Processing Technologies*, APPT'11, pages 142–156, Berlin, Heidelberg, 2011. Springer-Verlag.