# Multi-level Parallel Branch-and-Bound Algorithms for Solving Permutation Problems on GPU-accelerated Clusters

Mohand Mezmaz

► To cite this version:

Mohand Mezmaz. Multi-level Parallel Branch-and-Bound Algorithms for Solving Permutation Problems on GPU-accelerated Clusters. Complexité [cs.CC]. Université Lille 1, 2020. tel-03100303

HAL Id: tel-03100303

https://hal.science/tel-03100303

Submitted on 6 Jan 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
HABILITÉ À DIRIGER DES RECHERCHES
DISCIPLINE : INFORMATIQUE

présentée et soutenue publiquement,
le 18 septembre 2020, par

Mohand Mezmaz

Multi-level Parallel Branch-and-Bound Algorithms
for Solving Permutation Problems on
GPU-accelerated Clusters

Algorithmes Parallèles Multi-niveaux de Séparation et Évaluation pour la Résolution de
Problèmes de Permutation sur des Clusters Accélérés par GPUs

Président : Pierre Boulet Jury
Professeur, Université Lille Nord de France

Rapporteurs : Enrique Alba Professeur, Université de Málaga
Frédéric Saubion Professeur, Université d'Angers
Pierre Sens Professeur, Sorbonne Université

Examinateurs : Pascal Bouvry Professeur, Université du Luxembourg
Amir Nakib MdC HDR, Université Paris-Est Créteil

Garant : Nouredine Melab Professeur, Université Lille Nord de France
Acknowledgments

Je tiens tout particulièrement à exprimer ma profonde gratitude à Nouredine Melab, Garant de mon HDR: tannemirt! C'est une grande chance d'avoir continué, après ma thèse, de bénéficier de la rigueur de Nouredine, de sa motivation et de ses conseils.

Mes remerciements sincères vont également à Enrique Alba, Frédéric Saubion et Pierre Sens de m'avoir fait l'honneur de rapporter cette HDR, ainsi qu'a Pierre Boulet, Pascal Bouvry et Amir Nakib d'avoir accepté de faire partie de mon jury. Leur grande implication dans cette soutenance m'a permis de présenter mon HDR, même dans des conditions particulières.

Je voudrais aussi adresser un grand merci à Daniel Tuyttens, Chef du Service MARO de l'Université de Mons. Par sa confiance et ses conseils, Daniel a beaucoup contribué à mes recherches depuis toutes ces années.

Mes remerciements sincères vont également à Rudi Leroy et Jan Gmys, que j'ai co-
encadré durant leurs thèses. Leurs contributions constituent les ingrédients secrets du
succès de cette HDR.

Merci à tous les collègues et anciens collègues du Service MARO d'avoir participé à cet esprit d'équipe formidable.

Mes derniers remerciements et non les moindres vont à ma famille, et en particulier
à mes parents. Leurs efforts, leurs patientes et leurs encouragements sont déterminants
pour cette réussite.
# Contents

<table><tr><td>1</td><td>Main research activities since my PhD</td><td>1</td></tr><tr><td>2</td><td>Introduction</td><td>9</td></tr><tr><td>3</td><td>Background and related works</td><td>15</td></tr><tr><td>3.1</td><td>Introduction</td><td>15</td></tr><tr><td>3.2</td><td>Permutation problems</td><td>16</td></tr><tr><td>3.2.1</td><td>Permutations</td><td>16</td></tr><tr><td>3.2.2</td><td>Flowshop problem</td><td>17</td></tr><tr><td>3.3</td><td>Linked-list based PB&B</td><td>18</td></tr><tr><td>3.3.1</td><td>Branching</td><td>20</td></tr><tr><td>3.3.2</td><td>Bounding</td><td>21</td></tr><tr><td>3.3.3</td><td>Pruning</td><td>22</td></tr><tr><td>3.3.4</td><td>Selection</td><td>23</td></tr><tr><td>3.4</td><td>Related works</td><td>24</td></tr><tr><td>3.4.1</td><td>Parallel CPU B&B</td><td>24</td></tr><tr><td>3.4.2</td><td>Parallel GPU B&B</td><td>26</td></tr><tr><td>3.4.3</td><td>Cluster B&B</td><td>28</td></tr><tr><td>3.5</td><td>Conclusions</td><td>29</td></tr><tr><td>4</td><td>Single-core IVM-based permutation B&B</td><td>31</td></tr><tr><td>4.1</td><td>Introduction</td><td>31</td></tr><tr><td>4.2</td><td>IVM data structure</td><td>32</td></tr><tr><td>4.3</td><td>PB&B Operations</td><td>34</td></tr><tr><td>4.3.1</td><td>Selection</td><td>35</td></tr><tr><td>4.3.2</td><td>Branching</td><td>37</td></tr><tr><td>4.3.3</td><td>Bounding</td><td>37</td></tr><tr><td>4.3.4</td><td>Pruning</td><td>39</td></tr><tr><td>4.3.5</td><td>Initialization and finalization</td><td>39</td></tr><tr><td>4.4</td><td>Experiments</td><td>40</td></tr><tr><td>4.4.1</td><td>Experimental protocol</td><td>40</td></tr><tr><td>4.4.2</td><td>Obtained results</td><td>42</td></tr><tr><td>4.5</td><td>Conclusions</td><td>44</td></tr></table>
<table><tr><td>Contents</td><td>7</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><tr><td>8 Conclusions and perspectives</td><td>115</td></tr><tr><td>Bibliography</td><td>121</td></tr></table>
CHAPTER 1

# Main research activities since my
PhD

The research works done in my PhD thesis, defended at the end of 2007, focused on the parallelization of combinatorial optimization algorithms. During this thesis, we mainly presented a new approach, called B&B@Grid [MMT07b], for the parallelization of Branch-and-Bound (B&B) algorithms.

The results of our experiments show that the coding of B&B@Grid divides the size of the communicated information by an average of 766. These experiments have also shown that the coding of the communicated information, performed using B&B@Grid, is on average 92 smaller than the coding of the token approach of PICO [EPH00], and 465 smaller than the coding of the variable approach published in [IF00]. These experiments demonstrate the advantage of using B&B@Grid in clusters to reduce the size of the communicated information, and therefore the communication delays.

On the other hand, despite the high number of CPU cores in our clusters, about 1111
on average, the experiments, carried out with the *master-worker* paradigm, showed that
the master exploits its processor at 1.29% on average, and that workers spend an average
of 99.94% of their time in computing. These two percentages are good indicators of the
quality of B&B@Grid's load balancing strategy and its ability to scale up.

The B&B@Grid approach is used to solve an instance of the flowshop problem, known
as *Ta056* (50 jobs and 20 machines), and published in [Tai93] which has never been
optimally solved before. In terms of the used computing power, the *Ta056* resolution ranks
second among the large-scale challenges tackled in combinatorial optimization field before
2007. On average, 328 processors were used for more than 25 days, and a peak of about
1200 CPU cores was recorded during this resolution.

In the recent years, GPU-powered clusters appeared to be as more and more inter-
esting alternatives for the parallelization of some algorithms. GPU accelerators are often
used in today's largest high performance computing systems for regular and data parallel
applications. However, B&B algorithms are highly irregular in terms of workload, control
flow and memory access patterns. The research works, presented in this document, are
a continuation of B&B@Grid to adapt B&B algorithms, having an irregular structure, to
GPU-powered clusters, having accelerators with a highly regular architecture.

Before presenting our contributions, in the parallelization of B&B algorithms, this chapter describes our main industrial and academic research activities, done after my PhD and which are not presented in this manuscript. This industrial and academic research is presented in Subsection 1 and Subsection 1, respectively.

## 2.1 Industrial research

Our industrial research is mainly carried out in the railway field. These research works allowed us to have several contributions, such as the development of a train simulator, the development of a train autopilot, reducing the energy consumption of a rail network, or automatic learning of the technical parameters of a train. In general, our work for the railway industry is not intended for scientific publications. One of our main industrial partners, Alstom, prefers to keep confidentiality for promising results.

### G-Drive : an automatic train operation

* **Problem description:** The cost of energy, for the traction of passenger trains, is estimated at ~100 million euros during 2006 in Belgium, and at ~850 million euros during 2008 in France. Some studies show that, in the same journey and the same driving conditions, a train can consume up to 50 % more energy than another train. Therefore, the cost of the wasted energy for train traction is relatively high.

* **Conventional approach:** To reduce this cost, railway operators use two types of approaches. The first approach is to give the driver a list of expected times at some track positions during journey. And the second is to equip trains with a driver assistance system.

* **Our approach:** In collaboration with Alstom, we developed an automatic train control system. This system automatically controls and drive train by optimizing energy consumption and respecting the constraints.

* **Obtained results:** Simulator tests show that our system saves from 15% to 25% energy compared to drivers of the National Railway Company of Belgium.

* **Publications:** We decided to do not publish the used optimization method because Alstom plans to market this system in the coming months. An international patent is under submission. In this patent, I am one of the authors of the developed system.
## Energy-efficient railway traffic control

* **Problem description:** In modern railway system, most of the energy required by trains is supplied by the electric network. In recent years, the reduction of energy consuming has become one of the main concerns of the railway managers.

* **Conventional approach:** Although railway traffic control is an important topic in the modern railway management system, few results about the studies on railway traffic control, especially real-time energy-efficient traffic control for general railway system, have been published in the literature. Most of the studies simplify the problem by using average velocities or constant speeds to estimate the energy consumption of the trains along their journeys.

* **Our approach:** Our method proposes suitable driving profiles for the trains that run in the same railway network in order to reduce the total energy consumption of the whole railway network. To minimize the total energy consumption of a targeted railway network, it is important to not only propose ecodriving driving profile for each train but also to synchronize the operations of different trains. This synchronization introduces a high reutilization rate of the energy regenerated by dynamic braking operations of the trains.

* **Obtained results:** Experimental data used in to evaluate our approach are supplied by the Belgian railway infrastructure manager.

* **Publications:** This work is published in Mathematical Problems in Engineering journal [TFMJ13]. According to the experimental results, our proposed method can generate energy-efficient traffic control solution, where the driving profile of a train is defined by a suitable mono-train driving strategy that can be different from the driving strategy used by another train part in the same railway network.

## Learning the parameters of a train

* **Problem description:** In some algorithms, it is important to modelize the train. The obtained train model is defined by a set of parameters. These parameters are related to: (1) the train traction system; (2) its braking system; (3) its aerodynamic; (4) the different masses of the train; and (5) its length; In total, any realistic train model is defined by more than 100 parameters.

* **Conventional approach:** For each train, it is important to know the values of these parameters. To find these values, the conventional approach is (1) to contact the
manufacturer of the train to get the values of some parameters, and (2) to measure
the values of the other parameters. However, this approach has two disadvantages:
(1) measuring the value of a parameter can be an expensive operation; and (2) some
parameters, such as the total mass of a train, are constantly changing;

* **Our approach:** Our learning algorithm deduces the values of the parameters by only observing the behavior of the train. The train is seen as a system which receives input commands. These commands can be traction, braking, or no control. A command changes the state of the train. This state is defined by the position of train, its velocity, and its acceleration. The state of the train can be measured by an odometer, and this state is measured with a certain noise. The role of the learning algorithm is to deduce the parameters of the train using a great number of pair of values. Each pair is defined by (1) an input command and (2) the obtained state.

* **Obtained results:** With our method, we are able to learn the exact values of all parameters when the states of train is known without noise. If the states contain noise, the algorithm is able to learn these values with an average error of less than 2%.

* **Publications:** These results are obtained in January 2017. Next weeks, we will discuss with our partners the possibility of publishing this work.

## Developing a train simulator

* **Problem description:** In the optimization algorithm of the railway sector, the cost function is sometimes based on the simulation of a train.

* **Conventional approach:** At the beginning of a project, we received a professional simulator, where the simulation is based on a time discretization. However, we observed that the obtained simulations are not sufficiently accurate, and some values do not match those of the train's mathematical model.

* **Our approach:** So we developed a train simulator, based on a time discretization, and with a better accuracy. This simulator takes into account, not only, the control system, the traction and braking systems, and the aerodynamic of the train, but also, the track topology, namely gradients and curves.

* **Obtained results:** With negligible maximum errors in the positions, velocities and accelerations, the obtained simulations correspond to those of the train's mathematical model.
* **Publications:** This simulator is more a development work than a research one. As the scientific contribution is not important, we do not plan to publish this work.

## Multi-objective optimized railway timetable

* **Problem description:** A timetable can be improved according to different criteria, such as robustness, which measures the sensitivity of a timetable to disturbances, and total energy consumption. As the Belgian network is relatively dense, it is important to improve the timetable by taking into account these two criteria.

* **Conventional approach:** To the best of our knowledge, there is no work in the literature which optimize the timetable according to the robustness and the energy consumption, using a realistic energy model. In addition, some works try to optimize the timetable by building a new timetable from scratch. However, the users do not like having lot of changes in the timetable.

* **Our approach:** The objective of our algorithm is not to generate a new timetable but to improve the current one, transparently to the user. Between two successive stations, a timetable defines the moments at which the train must pass at certain intermediate positions of the track (i.e. between two stations). Our algorithm does not modify the arrivals and departures times at the stations, but updates only the times of these intermediate positions.

* **Obtained results:** The evaluation of the energy consumed by a network requires the use of a special simulator. In this project, our role was not to develop this simulator, but only the optimization algorithm. Our optimization software was provided to our partners approximately six months before the end of the implementation of the simulator.

* **Publications:** This work is not published.

## Energy-aware railway timetable re-scheduling during disturbances

* **Problem description:** Disturbance management can be seen as a three-step process. The first step is the detection of a disturbance. Then, a decision is made to manage this event. And finally, this decision is implemented.

* **Conventional approach:** In our work, we assume that another system is responsible for detecting disturbances. Decision-making, such as the removal of a train, its deviation, or its slowdown, is a complicated process, which is hard to automate. When a
decision is taken to change the path of a train, the Belgian railway operator, Infrabel,
manually generates a local timetable for this path.

* **Our approach:** The software, which we have developed, allows to generate this local timetable automatically by optimizing at the same time the tardiness and the energy consumption.

* **Obtained results:** The evaluation of the energy consumed by a network requires the use of a special simulator. In this project, our role was not to develop this simulator, but only the optimization algorithm. Our optimization software was provided to our partners approximately six months before the end of the implementation of the simulator.

* **Publications:** This work is not published.

## 2.2 Academic research

Simultaneously with industrial research, we also had contributions in the academic re-
search. This section presents only our main academic works which are not related to B&B
algorithms. These contributions belong to three different domains, namely combinatorial
optimization, cloud computing and machine learning.

### Hyper-heuristic GRASP

* **Problem description:** GRASP is a well-known two-phase metaheuristic. First, a construction phase builds a complete solution iteratively by component by component by a greedy randomized algorithm. After that, a local search phase improves this solution. The basic GRASP configuration is defined by (1) a cost function, (2) a probabilistic parameter of greediness, and (3) a neighborhood structure.

* **Conventional approach:** Traditionally, the choice of a GRASP configuration (greediness, cost function, and neighborhood) is guided by the theoretical studies or by manually testing the different parameter values. However, hyper-heuristic framework exists to design and adapt automatically heuristics to solve hard computational search problems. The feature of this framework is a protocol (also called a high-level heuristic or mechanism) that manages with the low-level heuristics (tabu search, genetic algorithms, etc.) looking for the best configuration of their parameters.

* **Our approach:** Therefore, in our work we propose a hyper-heuristic that runs GRASP configurations in accordance with several predefined configurations. Each
configuration manages a set of one-iteration GRASPs with different parameter. The goal is to automatically test multiple GRASPs and collect performance measures in order to detect the leading configurations for given test instances. We consider 7 well-known neighborhood structures, 9 greediness values, and 5 cost functions (including a cost function based on a bounding which is integrated for the first time in GRASP).

* **Obtained results:** Our computational experiments have not revealed one leading configuration, however the winning one is a stochastic mix of several configurations.

* **Publications:** This work is published in CCPE journal [AMTM17] (Concurrency and Computation: Practice and Experience).

## Optimization of precedence-constrained parallel application problem

* **Problem description:** Parallel programs can be represented by a directed acyclic graph where each node represents a task and each edge from a node A toward a node B means that the task B starts after the end of the task A. In this graph, The weight of a node A represents the computation cost of the task A, and the weight of an edge (A,B) represents the communication cost from the task A to the task B. These parallel tasks must be deployed on several computing units. Each computing unit can work with different voltage levels. For each voltage level, a computing unit has a processing speed and an energy consumption. Precedence-constrained parallel applications are one of the most typical application models used in scientific and engineering fields.

* **Conventional approach:** In the literature, usually the goal is to deploy this application in order to minimize the total processing time only.

* **Our approach:** Our objective was to minimize not only the total processing time but also the total energy consumption.

* **Obtained results:** Experiments show that our method improves on average the results of the literature for about 10,000 standard instances. The energy consumption is reduced by about 50% and the processing time by 10%.

* **Publications:** This work is published in JPDC [MMK+11] (Journal of Parallel and Distributed Computing) which is an A*-journal. According to Google Scholar, this paper is the third most cited paper in JPDC during the last 5 years (more than 120 citations).
Solving Tuberculosis Transmission Control (TTC) problem using parallel simulation-
and surrogate-assisted optimization

* **Problem description:** Tuberculosis is the most lethal infectious disease with 10 million cases of active tuberculosis recognized in 2017. The World Health Organization (WHO) proposed the End Tuberculosis Targets, with the objective to cut new cases by 90% between 2015 and 2035.

* **Conventional approach:** The AuTuMN simulator [TRDM17] implements a tuberculosis transmission dynamic model and facilitates predictions of future tuberculosis epidemic trajectories across diverse scenarios. Besides, the coupling with an economic model allows to compare and analyze different control policies, which then helps decision makers in their effort to decrease local tuberculosis burden. However, solving TTC problem using AuTuMN simulator, like many simulation-based optimization resolutions, are computationally expensive.

* **Our approach:** In the PhD thesis of Guillaume Briffoteaux, which I co-supervise, we investigate Evolution Controls (ECs) strategies that define the alternation between the simulator and the surrogate within the optimization process. The challenge is to find the best trade-off between the quality (in terms of precision) and the efficiency (in terms of execution time) of the resolution. We consider several ECs from the literature including the mechanism built on the distance-based concept of Hyper-sphere Confident Regions (HCR) and a mechanism derived from Monte Carlo dropout (MCDropout), a technique originally dedicated to quantify uncertainty in deep learning. The investigation of this latter uncertainty-aware EC is uncommon in surrogate-assisted optimization.

* **Obtained results:** We identify some characteristics common to several ECs and then derive from them different variants that are compared considering a pioneering application to the TTC problem. The reported results show that HCR- and MCDropout-based EC coupled with massively parallel computing outperforms the literature strategies.

* **Publications:** A first part of this work is published in the HPCS'2018 conference [BMMT18], the second part is under minor revision for the Future Generation Computer Systems (FGCS) journal (Impact Factor of 5.768), and the third part is submitted to the Swarm and Evolutionary Computation (SWEVO) journal (Impact Factor of 6.330).
CHAPTER 2

Introduction

This manuscript gives an overview of my research activities focused mainly on two PhD theses that I co-supervised, namely those of Rudi Leroy [Ler15a] and Jan Gmys [Gmy17] defended respectively in 2015 and 2017. Other activities conducted within the context of academic and industrial collaborations are summarized in the CV part of this manuscript (Chapter 1).

Many industrial and economic problems, like flowshop scheduling, are permutation
combinatorial optimization problems. Solving these problems consists in finding an opti-
mal permutation of elements among a large finite set of possible permutations. A wide
range of these problems is known to be large in size and NP-hard to be solved. The Per-
mutation Branch-and-Bound (Permutation B&B or PB&B) algorithms are one of the most
used exact methods to solve these permutation optimization problems. These algorithms
are based on an implicit enumeration of all the feasible permutations of the problem to
be tackled. Building and exploring the PB&B tree are performed using four main opera-
tions: branching, bounding, selection and pruning. In PB&B algorithms, if the lower bound
for some tree node A is greater than the best permutation found so far, then A may be
discarded from the search. This key idea of the PB&B algorithms significantly reduces
the number of explored nodes. However, the execution time of a PB&B notably increases
with the size of the tackled problem instance, and often only small or moderately sized
instances can be practically solved.

For this reason, over the last decades, parallel computing has been revealed as an attrac-
tive way to deal with larger instances of combinatorial optimization problems. Recently,
multi-core processors (CPUs), Graphics processing units (GPUs) and computing clusters
have been used for the parallelization of several algorithms. To the best of our knowledge,
all parallel PB&B algorithms developed in the literature are based on using one centralized
or distributed pool which stores nodes [GC94]. In these conventional approaches, PB&B
threads cooperate by adding nodes to or removing them from this or these pool(s), usually
implemented as a linked list (LL) [Cra06]. A parallel PB&B algorithm stops when this or
these pool(s) is/are empty. In the two theses, our work has led to the development of
four approaches, renamed in this manuscript PB&B@CORE, PB&B@CPU, PB&B@GPU and
PB&B@CLUSTER.
## 2.1 PB&B on single-core processors

A single-core processor is a microprocessor able to run one single thread at any time. Processors remained single-core until it was difficult to improve their computing power by increasing clock speed or the number of transistors. Therefore even processors of smartphones are nowadays no longer single-core processors.

However before developing a parallel approach of an algorithm, it is often important to optimize the serial approach of this algorithm. A serial algorithm can be optimized in order to run it more rapidly or use less memory. If code optimization can be done using new procedures or functions, memory optimization often requires the development of a new data structure.

In our work, we propose an original and pioneering single-core PB&B@CORE algorithm, based on a new data structure, called Integer-Vector-Matrix (IVM), and the factorial number system. This special numbering system [Lai88a], also called factoradic number system, is a mixed radix numeral system adapted to numbering permutations. The objective of our new IVM-based PB&B approach is to accelerate the management of the PB&B pools and to reduce the size of the memory used to store these pools.

Our new PB&B@CORE is compared with a typical LL-based approach [MMT13] in terms of memory usage and CPU time used to manage the PB&B pool. This comparison shows that our PB&B@CORE approach outperforms this typical LL-based PB&B approach.

## 2.2 PB&B on multi-core processors

As indicated before, for decades it was possible to improve performance of a general-purpose CPU by increasing its operating frequency. However in about 2012, manufacturers encounter two main technical barriers to improve the CPU performance, namely the power barrier (i.e. the power consumption exponentially increases with each factorial increase of operating frequency), and the memory barrier (i.e. the gap between memory and processor speeds is increasing). In order to continue delivering regular performance improvements for general-purpose processors, manufacturers such as Intel and AMD have turned to multi-core designs. Multi-core processors embed two or more independent execution CPU cores into a single processor.

By providing multiple execution cores, each sequence of instructions, or thread, has a hardware execution environment entirely to itself. This enables each thread run in a truly parallel manner. Nowadays, most processors can be considered as parallel machines. When designing an algorithm for a multi-core processor, it is therefore important to take into account a certain number of issues, such as the definition of work unit and work
sharing/stealing strategies. A major advantage of multi-core processors is the possibility to parallelize using threads instead of processes. Unlike processes, which have their own virtual memory, the threads of a process share the same virtual memory. Therefore, communications between these threads are faster than between processes. Many programming methods (such as [EPS09a], [PŽ09], and [SECG11]) based on the use of threads have been developed.

Our new method PB&B@CPU is based on a new approach to manage the pools of nodes of a PB&B parallel algorithm using five different work stealing strategies. This approach aims at reducing the CPU time used to manage the thread private PB&B pools, the memory behavior of the algorithm, the number of performed work stealing operations and finally the total execution time. It is based on the use of the IVM data structure to manage the thread private PB&B pools. In our stealing strategies, work units stolen by threads are intervals of factoradics instead of sets of nodes. The IVM-based PB&B@CPU is compared with the LL-based PB&B approach by solving some hard flowshop scheduling problem instances using the five work stealing strategies.

## 2.3 PB&B on GPU accelerators

Commodity Graphics Processing Units (GPUs) have rapidly evolved to become high performance accelerators for data-parallel computing. Recent GPUs, like Tesla V100, contain more than 5000 of processing units, capable of achieving up to 14 TFLOPS for single-precision arithmetic, and 7 GFLOPS for double-precision calculations. In addition, modern high-performance computing optimized GPUs contain up to 32GB of on-board memory, and are capable of sustaining memory bandwidths up to 900GB/sec. GPU-based accelerators are also becoming popular because of their extraordinary energy efficiency, as illustrated by The Green500 List<fnref n="1" />. The maximum power consumption of Tesla V100 is 300W. However, the parallelization of an application on a GPU can not be done directly. It is important to take into account a certain number of issues, such as the optimization of data transfer between CPU and GPU, the reduction of thread divergence, an efficient mapping of data on the GPU cores, and the best location of data on the different memories of the GPU. Because of all these issues, the usage of GPU is often restricted to regular and data parallel applications.

The irregular nature, in terms of workload, control flow and memory access patterns, of applications such as PB&B may seriously degrade the performance of the GPU. The acceleration of PB&B algorithms using GPUs is therefore a challenging task which is

<footnote n="1">¹http://www.green500.org</footnote>
addressed by only a few works in the literature, such as [CM13], using flowshop as a test case, [CMNLdC11], applied to the traveling salesman problem and [LEB12], applied to the knapsack problem where the search tree is binary. All these approaches use linked lists (deques, stacks, etc.) to store and manage the pool of nodes, likewise most parallel PB&B algorithms in the literature. Such data structures are very difficult to handle on the GPU and often induce prohibitive performance penalties. For this reason all GPU accelerated PB&B algorithms at our knowledge perform the management of the pool of nodes at least partially on the CPU, requiring costly data transfers between host and device. In [MCB14] it is shown that the bounding operation for flowshop consumes 97 – 99% of the execution time of a sequential PB&B and that the GPU based parallelization of this operation can provide a substantial acceleration of the algorithm. However, as the management of a list of pending nodes is performed on the CPU, the transfer of data between CPU and GPU constitutes a bottleneck for GPU accelerated B&B algorithms.

This manuscript describes our new PB&B@GPU approach, which is based on the par-
allelization of PB&B on GPU. This parallel PB&B@GPU algorithm is, to the best of our
knowledge, the first one that implements all PB&B operations on the GPU, requiring vir-
tually no interaction with the CPU during the exploration process. It is based on the IVM
data structure, which allows the efficient storage and management of the pool of nodes
in permutation based combinatorial optimization problems. The IVM structure provides
some regularization as it allows to store and manage the pool of nodes with data structures
of constant size. However, the IVM-based parallel PB&B is still highly irregular in terms
of workload, control flow and memory access patterns. None of these three issues can be
ignored when implementing the PB&B algorithm on the GPU and all three are addressed
in PB&B@GPU. The focus is put on the reduction of thread divergence which arises in
CUDA's SIMD execution model as a consequence of control flow irregularities.

## 2.4 PB&B on GPU-powered clusters

Large-scale GPU clusters are gaining popularity in high performance computing commu-
nity. In November 2019, nearly 40% of the total compute power on the TOP500 clusters<fnref n="2" />
(i.e. 626 petaflops) comes from GPU-accelerated systems. Just over a decade ago, no
supercomputers and clusters on the list were accelerated.

By gathering GPU processing power distributed across multiple computing nodes, it
is possible to run advanced, large-scale applications efficiently, reliably, and quickly. This
acceleration delivers a dramatic boost in throughput and cost savings, paving the way

<footnote n="2">²https://www.top500.org</footnote>
to exascale science applications. The massively parallel hardware architecture and high
performance of floating point arithmetic and memory operations on GPUs make them
particularly well-suited to many of the same scientific and engineering workloads that
occupy HPC clusters, leading to their incorporation as HPC accelerators.

To distribute an algorithm on the computing nodes of a GPU-powered cluster, it is
necessary to rethink the conventional parallel mechanisms in order to adapt them to the
characteristics of this environment. Taking into account of such characteristics can be
done by the resolution of the issues related to communication delays (i.e. optimizing the
communication cost) and processor's heterogeneity (i.e. optimizing the load-balancing
strategy).

In our work, we have developed a new method, called PB&B@GPU, based on the
coordinator-worker B&B@Grid approach [MMT07c]. In order to enable the use of a hybrid
CPU-GPU cluster in PB&B@GPU, a redefinition of work units is proposed using interval
lists. A detailed description of the coordinator and worker new operations is provided,
focusing on the revisited communication scheme. Finally, PB&B@GPU is experimented on
three different GPU-enhanced clusters with up to 36 GPUs. The experimental evaluation
includes scalability and stability analysis on solving very large flowshop instances.

## 2.5 Document contents

In addition to Chapter 3, which gives a state of the art of B&B algorithms, this manuscript is organized in four main chapters. Chapter 4 proposes an original and pioneer serial PB&B@CORE algorithm, based on a new data structure. Chapter 5 explains our PB&B@CPU, designed to work on multi-core CPUs. Chapter 6 describes the new PB&B@GPU approach, which is based on the parallelization of PB&B on GPUs. And finally, Chapter 7 presents our new PB&B@GPU approach for GPU-powered clusters.
CHAPTER 3

# Background and related works

## Contents

<table><tr><td>3.1 Introduction</td><td>15</td></tr><tr><td>3.2 Permutation problems</td><td>16</td></tr><tr><td>3.2.1 Permutations</td><td>16</td></tr><tr><td>3.2.2 Flowshop problem</td><td>17</td></tr><tr><td>3.3 Linked-list based PB&B</td><td>18</td></tr><tr><td>3.3.1 Branching</td><td>20</td></tr><tr><td>3.3.2 Bounding</td><td>21</td></tr><tr><td>3.3.3 Pruning</td><td>22</td></tr><tr><td>3.3.4 Selection</td><td>23</td></tr><tr><td>3.4 Related works</td><td>24</td></tr><tr><td>3.4.1 Parallel CPU B&B</td><td>24</td></tr><tr><td>3.4.2 Parallel GPU B&B</td><td>26</td></tr><tr><td>3.4.3 Cluster B&B</td><td>28</td></tr><tr><td>3.5 Conclusions</td><td>29</td></tr></table>

## 3.1 Introduction

The exact resolution methods used in combinatorial optimization are often B&X algorithms. These methods are mainly available in three variants: Branch-and-Bound (B&B), Branch-and-Cut (B&C) and Branch-and-Price (B&P). There are other less known B&B variants like Branch-and-Peg [GGS04], Branch-and-Win [PC04] and Branch-and-Cut-and-Solve [CZ06]. This list is certainly not exhaustive. It is also possible to consider a simple tree-based algorithm like Divide-and-Conquer as a base for the B&B algorithm. It is enough to remove the pruning operation of the B&B, explained later, to obtain the Divide-and-Conquer algorithm. Some authors consider the B&C, B&P algorithms and other variants as separate B&B algorithms. In the remainder of the manuscript, the B&B algorithm designates the simple B&B itself or any other variant of this algorithm. In addition, permutation
B&B (PB&B) designates any B&B algorithm or its variants when it solves a permutation problem.

Many problems, solved by B&B algorithms, are permutation problems, where the goal is to find the optimal scheduling in a set of elements. To the best of our knowledge, these permutation problems are solved using B&B without taking into account the permutation aspect of these problems. With taking into account this aspect, it is possible to significantly improve B&B algorithms. Therefore, the remainder of this manuscript focuses primarily on PB&B. In other words, our manuscript focuses on solving permutation problems using B&B algorithms.

In addition to this introduction and a conclusion, this chapter is divided into three sections. In Section 3.2 and Section 3.3, a background is given on respectively permutation problems and the conventional PB&B algorithm based on linked-list. Finally, Section 3.4 gives a state of the art of B&B algorithms.

## 3.2 Permutation problems

This section presents a brief background on permutations focusing on the flowshop problem considered as a test case to validate our approaches.

### 3.2.1 Permutations

Let's assume a permutation problem where the objective is to find the best permutation in a set of *N* elements. These elements can be jobs, cities, locations, and so on. It is always possible to assign a number to each of these *N* elements. The first element can be designated by 1, the second element by 2, ..., and the last by *N*. Therefore, any permutation problem of *N* elements can be represented as a permutation problem of the first *N* positive natural numbers.

A permutation is obtained after the assignment of numbers to positions. Before making this assignment, all numbers are free and the positions are empty. At the end of the assignment, all numbers are assigned and positions are occupied. During the assignment, some numbers are free, others assigned, some positions are occupied and others empty. In the rest of the manuscript, a number can be said free or assigned, and a position can be said occupied or free.

After the assignment of all numbers to all positions, the permutation is called complete. For example, (2, 3, 1, 4) is a complete permutation where the numbers 2, 3, 1 and 4 are respectively assigned to positions 1,2, 3 and 4. Before starting the assignment, the permutation is said to be empty. For example, {1, 2, 3, 4} is an empty permutation where
all numbers are free and all positions are empty. During an assignment, the permutation is said to be partial. For example, $(2, 3, \{1, 4\})$ is a partial permutation, where numbers 2 and 3 are assigned, numbers 1 and 4 are free, positions 1 and 2 are occupied, and positions 3 and 4 are free. A partial permutation can also have the form $(2, \{1, 4\}, 3)$. In this example, positions 1 and 4 are occupied and the other positions are empty.

Among the empty positions, it is possible to specify **the first, the second, ... and the last empty position**. For example, in $(2, \{1, 4\}, 3)$, the first empty position is 2 and the last empty position is 3. Among the assigned numbers, it is possible to specify also **the first, the second, ... and the last number assigned**. For example, in $(2, \{3, 1, 4\})$, the last number assigned is 2. On the other hand, looking at the $(2, \{1, 4\}, 3)$, it is not possible to find the first and the last number assigned. This manuscript defines **the size of a permutation** as the number of free numbers or empty positions. For example, the permutations $(2, \{1, 4\}, 3)$, $\{1, 2, 3, 4\}$ and $(2, 3, 1, 4)$ respectively have as sizes 2, 4 and 0.

Any permutation, whether empty, partial or complete, contains a certain set of complete permutations. For example, the partial permutation $(2, \{1, 4\}, 3)$ contains the two permutations $(2, 1, 4, 3)$ and $(2, 4, 1, 3)$. The empty permutation $\{1, 2, 3, 4\}$ contains all 4! possible permutations. In the same way, the complete permutation $(2, 3, 1, 4)$ contains only one complete permutation. **The space of a permutation P**, denoted $Space(P)$, is defined as the set of the complete permutations contained in P.

### 3.2.2 Flowshop problem

The four approaches presented in this document are validated using the flowshop problem. In manufacturing environments, it is common to find permutation flowshop scheduling problems [BG76, KS80, AGA99] where *n* jobs have to be processed on *m* machines where the goal is to optimize an objective function. The objective of the flowshop is to schedule a set of *n* jobs on a set of *m* machines where each job $J_i$, $J_2$, ..., $J_n$ is processed on the machines $M_1, M_2, ..., M_m$ organized in the line. Each job $J_i$ with $i = 1, 2, ..., n$ is made of a sequence of *m* operations $O_{i1}, O_{i2}, ..., O_{im}$ where operation $O_{ik}$ is the processing of job $J_i$ on machine $M_k$ for a processing time $p_{ik}$ that can not be interrupted. The objective of the flowshop is to find a processing order on each machine $M_k$ which minimizes the time necessary to complete all jobs, also known as the *makespan*. In this manuscript, each reference to the flowshop is actually a reference to the permutation flowshop [AGA99, HS05]. Using Johnson's algorithm [Joh54], it is possible to find an optimal schedule for the flowshop in $O(n \log n)$ steps when $m = 2$. The problem is NP-hard when $m \ge 3$ [GJS76]. This is why it is often tackled using metaheuristics [Bas05] to deal with large problem instances. Figure 3.1 shows an example of a flowshop instance where $n = 3$ and
$m = 4$, it also shows the optimal complete permutation.

<table><thead><tr><th></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th><th>M<sub>4</sub></th></tr></thead><tbody><tr><th>J<sub>1</sub></th><td>5</td><td>3</td><td>4</td><td>1</td></tr><tr><th>J<sub>2</sub></th><td>2</td><td>2</td><td>1</td><td>4</td></tr><tr><th>J<sub>3</sub></th><td>1</td><td>3</td><td>5</td><td>2</td></tr></tbody></table>

Processing Times

<figure><img src="image_1.png" /></figure>

Optimal Solution

Figure 3.1: Illustration of a permutation flowshop where $n = 3$ and $m = 4$. The table shows the processing times of the jobs on each machine. The Gantt diagram shows the optimal complete permutation for this particular instance.

These are the constraints that a valid flowshop complete permutation should satisfy:

* A machine can not start processing a job before the preceding machines have finished the processing of that job. In other words, machine $M_j$ can not process operation $O_{ij}$ before it is completed on machine $M_{j-1}$.

* An operation can not be interrupted, and since a machine processes one job at a time, the machines are critical resources.

* The sequence of jobs must be the same on all machines, e.g. if job $J_3$ is processed in second position on the first machine, job $J_3$ must also be processed in second position on all the other machines.

The lower bound proposed by Lageweg et al. [LLK78] is used in our bounding operation. This bound is known for its good results and has complexity of $O(m^2nlog(n))$, where $n$ is the number of jobs and $m$ the number of machines. This lower bound is mainly based on Johnson's theorem [Joh54] which provides a procedure for finding an optimal complete permutation for a flowshop scheduling problem with 2 machines.

## 3.3 Linked-list based PB&B

This section describes the conventional PB&B algorithm based on a linked-list. As indicated by its name, the PB&B is based on two main operations, namely *branching* and *bounding*.
In addition to these two operations, the PB&B algorithm is defined by other operations, as shown in Figure 3.2.

<figure><img src="image_2.png" /><figcaption>Figure 3.2: The conventional PB&B algorithm and its operations.</figcaption></figure>

Two of the PB&B operations are called *initializing* and *finalizing* operations. The main role of the initializing operation is to allocate memory for the data structures of the algorithm, to initialize them, and to add the empty permutation to the PB&B pool. The finalizing operation is invoked when this PB&B pool is empty. Its role is to provide the decision maker by the best found permutation, and to release the allocated data structures.

Between these two operations, the PB&B algorithm runs a large number of iterations. At each iteration, the operations of selection, costing, updating, branching, bounding, and pruning always intervene in this order. The selection operation takes a permutation from the pool. Two cases may arise: in the first case, the selected permutation is complete, and in the second case, the permutation is empty or partial.

If the first case occurs, then the costing operation calculates the cost of the permutation. If the calculated cost improves the cost of the best known complete permutation, then the updating operation saves this complete permutation as the new best known complete permutation. In the second case, costing and updating operations are not used. On the other hand, the operations of branching, bounding and pruning intervene. The role of the branching operation is to divide the selected permutation into several permutations. The bounding operation calculates the bound of each of the obtained permutations, and the pruning operation removes a certain number of permutations and puts into the pool the
un-deleted permutations. Each of the following subsections presents in more details one of the main operations of the conventional PB&B.

### 3.3.1 Branching

The branching operation works according to the political and sociological strategy *Divide and rule* (from Latin dīvide et imperā). The idea here is to divide a partial permutation *P* into a set of *E* sub-permutations. Of course, there are different ways to divide the same permutations *P* and, for each division, a different set of *E* of sub-permutations can be obtained. For example, the permutation (1, {2, 3, 4}) can be divided into two different ways:

$$ \bullet E_1 = \{ (1, 2, \{3, 4\}), (1, 3, \{2, 4\}), (1, 4, \{2, 3\}) \}. $$

$$ \bullet E_2 = \{ (1, \{3, 4\}, 2), (1, \{2, 4\}, 3), (1, \{2, 3\}, 4) \}. $$

The two divisions are made by filling respectively the first and last empty positions, using one of the three free numbers. In our work, we use only the first free and last positions. However, from a theoretical point of view, it is also possible to use the other empty positions, such as the second empty position, the third position, etc. The branching operation can use any division of a permutation *P* in a permutation set *P*₁, *P*₂, ..., *P*ₙ that meets both of the following conditions:

* $\forall (P_i, P_j), \text{Space}(P_i) \cap \text{Space}(P_j) = \emptyset$. This first condition indicates that the sub-permutations of *P* must be disjoint.

* $\text{Space}(P) = \cup_{i=1}^{n} \text{Space}(P_i)$. This second condition is that the permutation space *P* equals the union of all sub-permutation spaces of *P*.

The same permutation *P* can not be split simultaneously with different techniques. For example, *P* can not be split using both first and last free positions simultaneously. However, it is possible to use a certain technique of division of the permutation *P*, and another division technique to divide one of the sub-permutations of *P*. Therefore, it is possible to use a single division technique, but it is also possible to combine several techniques in the same PB&B algorithm. The example, Figure 3.3 shows the result of a complete division from an empty permutation of size 4. In this example, the division technique used is the one that is based on the first empty position.
<figure><img src="image_3.png" /><figcaption>Figure 3.3: Example of a complete branching tree using the first empty position division technique.</figcaption></figure>

A PB&B always starts with the division of the empty permutation of size *N*. The division of this empty permutation produces *N* partial sub-permutations of size (*N* − 1). The division of each of these partial sub-permutations produces (*N* − 2) sub-sub-permutations of size (*N* − 2). This process continues until the obtained permutations are complete and have size 0. These successive divisions make it possible to define a tree, called branching tree, where the empty permutation is the root of the tree, the complete permutations are the leaves of the tree, and partial permutations constitute the internal nodes of the tree.

The branching tree obtained is constituted of 1 empty permutation, *N* partial permu-
tations of size *N*, *N*(*N* − 1) partial permutations of size (*N* − 1), *N*(*N* − 1)(*N* − 2) partial
permutations of size (*N* − 2), ..., and finally, *N*! complete permutations. As the size of this
tree is exponential, it is impossible to explore, in a reasonable time, all permutations when
*N* exceeds a certain threshold value.

**3.3.2 Bounding**

To avoid exploring the entire branching tree, the PB&B algorithm uses another operation called the bounding. This operation receives a permutation, as input, and returns a cost, as output. If it is a minimization problem (i.e searching a permutation of minimal cost), then the bounding operation returns a lower bound. On the other hand, if it is a maximization problem, then this operation returns an upper bound. In the remainder of the manuscript, a problem of minimization is considered, but the proposed approaches are also valid for maximization problems.

As explained previously, the space of a partial permutation P, denoted *Space*(P), is equal to the set of complete permutations contained in P. For example, the permutation space of (1, {2, 3}, 4) includes the permutations (1, 2, 3, 4) and (1, 3, 2, 4). If a bounding operation receives as input the partial permutation (1, {2, 3}, 4), then this operation should return a value *Bounding*(P) such that the costs of (1, 2, 3, 4) and (1, 3, 2, 4) must be greater than *Bounding*(P). In the general case, the following condition must always be satisfied.
$$ \bullet \quad \forall P' \in \mathit{Space}(P), \mathit{Bounding}(P) \leq \mathit{Costing}(P') $$

For a permutation $P$, the higher the value of $Bounding(P)$ is, the better the bounding operation which is used. Let $P'$ be the complete permutation with the minimal possible cost in $Space(P)$. The best bounding operation is the one that returns a value of $Bounding(P)$ equals $Costing(P')$. In addition to the values of $Bounding(P)$, bounding operations can also be compared according to their computation times. The faster a bounding operation is, the better is this operation. If a permutation $P$ is complete, this manuscript assumes that $Bounding(P)$ is always equal to $Costing(P)$.

<figure><img src="image_4.png" /><figcaption>Figure 3.4: Example of a bounding tree with bound values of each node.</figcaption></figure>

### 3.3.3 Pruning

Subsection 3.3.1 shows that the branching tree is very large. To find the permutation with minimum cost, it is impossible to explore the entire tree for large problem instances. By knowing the lower bound of each partial permutation, it is possible to avoid exploring a large part of this tree. Let's assume that the algorithm knows a complete permutation $P$, which has a cost equal to $Costing(P)$, and that this algorithm encounters a partial permutation $P'$, which has a lower bound equal to $Bounding(P')$. In other words, all complete permutations in $Space(P')$ have a cost greater than $Bounding(P')$. If $Costing(P)$ is less than $Bounding(P')$, then there is no complete permutation of $Space(P')$ that can have a better cost than $P$. Therefore, it is unnecessary to further divide the partial permutation $P'$. This permutation $P'$ can be ignored and pruned from the branching tree. The pruning operation receives a cost $Costing(P)$ of complete permutation $P$ and the lower bound $Bounding(P')$ of a partial permutation $P'$ and returns a boolean value indicating whether $P'$ should be branched or pruned. Figure 3.5 gives an example of the bounded part of a branching tree. The bounded part of a tree is composed by all the nodes of the tree where the bounds are computed.
<figure><img src="image_5.png" /><figcaption>Figure 3.5: An example of the bounded part of a branching tree</figcaption></figure>

3.3.4 Selection

The role of the selection operation is to determine the order in which the permutations of
the tree are explored. In other words, its role is always to choose the next permutation, of
the pool, coded using a linked-list, to branch. In a PB&B, there are three main strategies
for selecting the next permutation of the pool to be branched.

*   **Best first:** The permutation of the pool having the smallest bound is chosen first. The advantage of this strategy is to increase the probability of finding more quickly a complete permutation of minimum cost. But its disadvantage is that the number of permutations of the pool increases rapidly.

* **Largest first:** This strategy selects the permutation of the pool with the largest size. In other words, it first selects the permutation close to the root of the branching tree. This strategy makes it possible to considerably increase the number of permutations of the pool. In parallel computing, the advantage of this strategy is to generate a pool large enough to occupy lot of computing units.

*   **Depth first:** This strategy selects the permutation with the smallest size. In other words, it first selects the furthest permutation from the root of the branching tree.

In our four PB&B approaches, the strategy used is a combination of depth and best first strategies. As already explained, a permutation is defined by its size, its bound and its last assigned number. These concepts of the last assigned number and the size of a permutation are explained in Subsection 3.2.1. Using these three values, the selection operation in our PB&B approaches combines depth and best first strategies. More exactly, the Selection operation (1) receives a permutation pool as input, (2) makes sure this pool is not empty, (3) determines the largest size S of all the permutations of the pool, (4) selects the set E of the permutations of the pool having this size S, (5) determines the best bound B of all the permutations of E, (6) selects the subset E' of all the permutations
of *E* having the bound *B*, (7) selects the permutation *P* having the lowest last number assigned, (8) removes this permutation *P* from the pool, and (9) returns the permutation *P* as the output of this Selection operation.

<figure><img src="image_6.png" /><figcaption>Figure 3.6: Example of the best found complete permutation and the content of the pool, coded using a linked-list, when using the depth best first strategy.</figcaption></figure>

For example, Figure 3.6 shows the PB&B tree during its construction and exploration. All the empty nodes of the tree are not explored, the gray permutations are already branched, the green permutations are pruned because of their bounds, the red complete permutation is the one having the minimal cost so far, and black permutations are neither pruned nor branched. In this tree, it is useless to keep in the memory the gray and green permutations. The only permutations that must be kept are those in black and red. Non-branched and non-pruned permutations are generally kept in a data structure called a linked-list, and the best complete permutation found so far can be kept in an array.

## 3.4 Related works

The related works of this section is already published in the thesis of Jan Gmys [Gmy17]. The design of parallel B&B algorithms is strongly influenced by the target architecture and the characteristics of the problem being solved [BHP05]. Therefore, and in spite of the simple and generic formulation of B&B, a large number of parallel algorithms have been proposed for different problems and architectures. [GC94] provides a complete, but over twenty year old survey of parallel B&B.

### 3.4.1 Parallel CPU B&B

Because of the simple basic formulation of B&B it is interesting to have a framework that allows users to easily customize B&B to solve their problems. Many software frameworks
have been proposed, including Bobpp [Men17, Bob], PEBBL [EHP15] and PICO [EPH01], parts of the ACRO project [ACR], ALPS/BiCePS [RLS04], BCP and SYMPHONY, which are parts of the COIN-OR project [CO].

This list includes only those frameworks which appear to be maintained at the time of writing. B&B frameworks establish an interface between the user and the parallel machine by defining abstract types for search tree nodes and solutions. As a user, one provides concrete implementations of these types as well as branching and bounding procedures, while the framework handles more generic parts of parallel B&B. The mentioned frameworks differ by the variant(s) of B&B they provide, the type of parallel model they propose and the parallel programming environment. They are implemented with these frameworks are usually designed as multi-layered class libraries, integrating additional features by building on top of existing layers. For example, BiCePS is built on top of ALPS to provide data handling capabilities required for implementing relaxation-based B&B, and PEBBL began its existence as the “core” layer of the parallel mixed integer programming (MIP) solver PICO.

The older versions of these frameworks are often based on master-worker approaches. In order to avoid that the master processor becomes a bottleneck, hierarchical organizations revealed more efficient than pure master-worker implementations [EHP15, Men17, BMT12b]. In these approaches groups of workers form clusters, cooperating locally and interacting with the master through a form of middle management. The idea is to improve locality and relieve the master process by introducing hubs, each handling several workers (master-hub-worker approach). In general, the role of hubs consists in providing work units to a set of workers and coordinating the search process locally, while limiting interaction with the master and worker processes. For the PEBBL framework near-linear speedups on over 6000 CPU cores are obtained for large B&B trees and node evaluation costs of about 10 seconds [EHP15].

Recently, [HSH<sup>+</sup>17] compared three implementations of a global optimization (GO) B&B algorithm using different levels of abstraction: the Bobpp framework, Intel Thread Building Blocks and a custom P-thread implementation. While they find the Bobpp implementation easiest to code, the authors show that the two other solutions offer better scalability for the used test case. For the optimized test functions, the authors report node processing rates of about 1 million nodes per second (Mn/s) for the sequential version of their custom implementation on a 2 GHz Sandy Bridge CPU.

[EPS09b] presents a software platform called BNB-Solver, allowing the use of serial, shared memory and distributed memory B&B. The proposed approach uses a global work pool and local work pools for each thread. Each thread stores generated nodes in its local pool during N B&B iterations. After N iterations a part of the local nodes are transferred
to the global pool. When the local pool of a thread is empty, the thread attempts to take nodes from the global pool and blocks if the global pool is empty. The algorithm terminates when the global pool is empty and all threads are blocked. The authors compare the performance of BNB-Solver with the ALPS and PEBBL frameworks and obtain results similar to [HSH$^{+}$17], in the sense that, for a knapsack problem (with a reported sequential node processing rate in the order of 1 Mn/s) BNB-Solver outperforms both frameworks.

[CMGH08] proposes two schemes for parallelizing B&B algorithms for global optimization on shared memory multi-core systems, Global and Local PAMIGO (Parallel advanced multidimensional interval analysis global optimization). Both algorithms are parallelized using POSIX threads. In Global PAMIGO, threads share a global work pool and therefore a synchronization mechanism is used for mutually exclusive accesses to the pool. For Local PAMIGO, where thread has its own pool of nodes, a dynamic load balancing mechanism is implemented. A thread stops when its local pool of nodes is empty. When the number of running threads is less than the number of available cores, and a thread has more than one node in its local pool it creates a new thread and transfers a portion of its pool to the new thread. Local PAMIGO ends when there exists no more running threads, and Global PAMIGO ends when the global pool is empty. The authors report profiling results for PAMIGO which show that memory management represents a large percentage of the computational burden. As a very large number of nodes are created in a relatively short amount of time, the kernel needs to satisfy memory allocation and deallocation requests from all threads, creating memory contention. The vast majority of parallel B&B algorithms in the literature store nodes in one or several pool(s) implemented as linked-lists (e. g. priority queues, stacks, deques).

### 3.4.2 Parallel GPU B&B

The study of [JAO$^{+}$11] provides a good overview of the challenges faced when implementing parallel backtracking on GPUs. Most of their conclusions from the investigation of GPU-based backtracking paradigm remain valid for B&B algorithm using a depth first search strategy. A fine-grained parallelization of the search space exploration and/or the node evaluation is necessary in order to make use of the GPU's massive parallel processing capabilities. This strongly depends on the nature of the problem being solved and on the choice of the parallelization model. Other critical factors include latency hiding through coalescence, saturation, and shared memory utilization [JAO$^{+}$11]. Generally speaking, the algorithmic properties of B&B, irregularity of the search space, irregular control flow and memory access patterns are at odds with the GPU programming model. Also, memory requirements for backtracking and B&B algorithms are often difficult to estimate and may
exceed the amount of memory available on GPUs. Several approaches for GPU-accelerated B&B algorithms have been proposed. These approaches correspond to different parallelization models and their design is often motivated by the nature of the problem being solved. According to the characteristics of the bounding function one may distinguish among approaches for fine-, medium- and coarse-grained problems.

The GPU B&B and backtracking algorithms for fine-grained problems proposed in [CMNLdC11, CNNdC12, FRvLP10, LLW$^+$15, RS10, ZSW11] perform massively parallel searches on the GPU, based on the parallel tree exploration model. The evaluation of a node for the *n*-Queens problem in [FRvLP10, LLW$^+$15, ZSW11] requires only a few registers of memory and only a couple of bit operations. The lower bound for the Asymmetric Traveling Salesman Problem (ATSP) used in [CMNLdC11, CNNdC12] is incrementally obtained by adding the cost of the last visited edge to the current cost and therefore has a complexity of $O(1)$. It requires an access to the distance matrix which can be stored in constant or texture memory. The size of the problems being solved is < 20 for both the ATSP and the *n*-Queens problems. These algorithms for fine-grained problems share a common approach: the search is split in two parts, an initial CPU search and a final GPU search. The upper tree of depth $d_{cutoff}$ is processed in sequential or weakly parallel manner on CPU, generating a set of active nodes at depth $d_{cutoff}$. This active set is sent to the GPU, where the lower part of the tree is processed in parallel. Each node of the active set is used as root node for an independent search, which is mapped either to a thread or a warp. This approach requires very careful tuning of the cutoff depth, which strongly influences granularity and load balancing.

Because of varying thread granularities, one of the major issues faced by such approaches is load imbalance. In all of these works the GPU search is performed without dynamic load balancing. However, as noted by [RS10], if "a job is divided into sufficiently many parts, an idle processor will be instantly fed with waiting jobs" and the "GPU's Thread Execution Manager performs that task automatically". This approach assumes two things: first, the initial CPU search is able to generate a large amount of nodes in a reasonable amount of time, and second, the work distribution among independent B&B searches is not too irregular.

For many combinatorial optimization problems the cost of the bounding operation is very high, compared to the rest of the algorithm. For instance, the most used lower bounding function for the flowshop consumes 97-99% of the sequential execution time [MCB14]. However, the cost of evaluating one node is sufficiently small to be efficiently performed by a single GPU thread. We therefore refer to this type of problem as medium-grained. For these problems, existing GPU-accelerated B&B algorithms in the literature use the GPU to evaluate large pools of nodes in parallel [CMMB13, VDM13, LEB12]. They use
conventional stacks or queues to store and manage the B&B tree on the host, offloading
the parallel evaluation of bounds to the device. Indeed, for these problems substantial
speedups can be achieved despite sequentially performing pool management on the host.
Substantial efforts have been made to port larger portions of the algorithm to the GPU
and to reduce overheads incurred by data transfers between CPU and GPU. For instance,
branching nodes on the device allows to copy only parent nodes to the GPU. Similarly,
pruning evaluated nodes on the device reduces the sequential portion and requires only
the transfer of non-pruned children nodes back to the host. Further performance im-
provements can be obtained by overlapping data transfers with GPU computations, as for
example in [VDM13]. For fine-grained problems this approach is likely to perform poorly.

For coarse-grained problems the best way to use the GPU may be as an accelerator for the bounding function itself. In [ABEB+16] a GPU-accelerated B&B algorithm for the jobshop scheduling problem is proposed. The approach also offloads nodes to the GPU but uses a block-based parallelization for each node evaluation. The number of nodes that need to be offloaded in order to saturate the GPU is therefore smaller than for medium-grained problems. A GPU-accelerated algorithm for problems with linear programming bounds is proposed in [MCA13]. Using a GPU-based LP-solver to accelerate this type of problems is very challenging. However, the authors report that for large problems above a certain density threshold their hybrid GPU-accelerated solver outperforms the sequential CLP solver of the open-source COIN-OR library.

3.4.3 Cluster B&B

There are very few works on the parallelization of B&B using multiple GPUs and CPUs in distributed heterogeneous systems. In [VDM13] a linked-list based fully distributed hybrid B&B algorithm combining multiple GPUs and CPUs is proposed. As a test case 20 jobs-on-20 machines flowshop instances are considered in their experiments using a platform composed of 20 GPUs and 128 CPUs. For load balancing a random work stealing mechanism is used. The authors propose an adaptive granularity policy to adapt the amount of stolen nodes at runtime to the normalized computing power of thief and victim. The algorithm is based on a 2-level parallel model, using GPUs for parallel evaluation of lower bounds. In order to reduce CPU-GPU communication overhead, an asynchronous implementation with overlapping host and device computations is proposed. Experimentally, near-linear mixed scalability is shown up to 20 GPUs and 128 CPUs. In [CMMT13] the combined usage of multi-core and GPU processing is investigated. An experimental comparison of concurrent and cooperative approaches shows that the cooperative approach improves the performance with respect to a GPU-only approach while the concurrent approach is not
beneficial. Among other issues, the authors identify the reduction of CPU-GPU communication overhead as a major challenge and propose overlapping communication schemes and auto tuning of the offloaded pool sizes to answer this challenge.

Some of the largest known exact resolutions of combinatorial optimization problems have been performed using the master-worker paradigm in combination with grid computing technologies (e. g. *nug30* [ABGL02]). The B&B@Grid platform [MMT07c] uses an interval encoding for work units which significantly reduces the size of messages communicated in distributed B&B. Designed for volatile computing environments, B&B@Grid is fault tolerant thanks to its checkpointing mechanism.

In [BMT12a] an adaptive multi layer hierarchical master-worker approach is applied to the B&B algorithm, using flowshop as a test case. The proposed approach evolves as new resources join the computation, and integrates three types of processes, a super master, masters and workers. Results obtained at the scale of up to 2 000 CPUs show that the multi-layered hierarchical approach clearly outperforms single-layered and classical master-worker approach in terms of efficiency for instances smaller than *Ta056*, as it minimizes bottlenecks at the level of the master and reduces idle time of the workers. In [BMT14] the authors extend their approach, proposing a fault tolerance mechanism.

## 3.5 Conclusions

This chapter presents the permutation problems, the PB&B algorithms and a state of the art of the B&B algorithm. To the best of our knowledge, permutation problems are solved with the B&B without permutation awareness. The objective of this manuscript is to present new approaches to solve permutation problems, using a B&B that takes into account the permutation aspect of these problems. The next chapters of the manuscript present four new approaches, noted PB&B@CORE, PB&B@CPU, PB&B@GPU and PB&B@CLUSTER. These approaches are developed for four hardware architectures, namely a multi-core computing node, a CPU-accelerator of cores, a GPU and a cluster of computing nodes.

* PB&B@CORE is a serial method, unlike the other three approaches which are parallel methods.

* PB&B@CPU is a method dedicated to a multi-core CPU deployment. Since several years, all commercialized CPUs contain several computing cores, and each CPU can be seen as a parallel machine.

* PB&B@GPU uses a data parallelism, unlike the other three methods based only on the parallelism of instructions. In a GPU, the same instruction can run on several data
simultaneously, according to the Single Instruction Multiple Data (SIMD) computing model.

* PB&B@CLUSTER is an approach dedicated to an architecture with distributed memory, while the other three approaches are all dedicated to an architecture with shared memory.

These PB&B@CORE, PB&B@CPU, PB&B@GPU and PB&B@CLUSTER approaches are presented respectively in the chapters 4, 5, 6 and 7.
# Single-core IVM-based permutation
B&B

## Contents

<table><tr><td>4.1 Introduction</td><td>31</td></tr><tr><td>4.2 IVM data structure</td><td>32</td></tr><tr><td>4.3 PB&B Operations</td><td>34</td></tr><tr><td>4.3.1 Selection</td><td>35</td></tr><tr><td>4.3.2 Branching</td><td>37</td></tr><tr><td>4.3.3 Bounding</td><td>37</td></tr><tr><td>4.3.4 Pruning</td><td>39</td></tr><tr><td>4.3.5 Initialization and finalization</td><td>39</td></tr><tr><td>4.4 Experiments</td><td>40</td></tr><tr><td>4.4.1 Experimental protocol</td><td>40</td></tr><tr><td>4.4.2 Obtained results</td><td>42</td></tr><tr><td>4.5 Conclusions</td><td>44</td></tr></table>

## 4.1 Introduction

Conventional B&B algorithm is designed to work on all optimization problems. Conversely, the new approach PB&B@CORE is especially dedicated to the resolution of permutation problems. In other words, the B&B algorithm is generic to all optimization problems, while our new PB&B@CORE approach is specific to permutation problems. Compared to the generic B&B algorithm, the PB&B@CORE approach manages a partial permutation pool more efficiently. This PB&B@CORE efficiency, in terms of memory and CPU usage, can be explained by our new data structure.

While B&B algorithms use a linked-list to encode a pool, our approach PB&B@CORE uses a new data structure, called Integer-Vector-Matrix (IVM). The operations of the conventional B&B, explained in Chapter 3, are designed to work on a pool encoded by a
linked-list. It is therefore necessary to rethink these operations to work effectively on an
IVM structure. This chapter, which contains three main sections, explains IVM and the
PB&B@CORE operations.

Section 4.2 describes our new IVM structure, used for the PB&B pool coding. Section 4.3 explains the PB&B operations designed to work on this IVM structure. Finally, Section 4.4 presents the experiments carried out to validate the PB&B@CORE approach.

## 4.2 IVM data structure

<figure><img src="image_7.png" /><figcaption>Figure 4.1: Example of representing a pool with a linked-list and an IVM structure.</figcaption></figure>

Figure 4.1 shows an example of the representation of a pool with, on the one hand, a linked-list and, on the other hand, our IVM structure. This figure shows the PB&B tree during its construction and exploration. In the tree, the gray permutations are already branched, green permutations are pruned because of their bounds, and black permutations are neither branched nor pruned. In the linked-list, there is no need to keep gray and green permutations. The only permutations that must be kept are those which are black. Using the depth best first strategy described in Chapter 3, the size of the linked-list can reach $N \sum_{i=1}^{N-1} i$ integers, *N* being the size of the permutation problem. The objective of the IVM structure is twofold, namely the reduction of the pool size and accelerating
the processing of this pool by its operations. The IVM structure is composed of six main sub-structures:

* **Permutation Matrix (M):** This structure is a matrix, denoted *M* in Figure 4.1, of size *N*², where *N* is the size of the permutation. In this matrix, only half of the cells are occupied. In order to encode a PB&B tree using this matrix, five rules are used:

- **Rule 1:** Each cell $M_{i,j}$ of the matrix corresponds to a permutation $Perm(M_{i,j})$ of the PB&B tree. For example, $Perm(M_{0,1})$ is equal to ({1, 3, 4, 5}, 2).

- **Rule 2:** The permutations of depth *i* of the tree are encoded at the row *i* of the matrix. For example, the depth 3 of the tree contains the permutations (5, 4, {1}, 3, 2) and (5, 1, {4}, 3, 2). These permutations are therefore encoded at the row 3 of the matrix, more exactly at cells $M_{3,0}$ and $M_{3,1}$.

- **Rule 3:** In the same row of the matrix, the permutations are ordered according to their lower bound values. For example, in the row 3, *Bounding*(5, 4, {1}, 3, 2) is smaller than *Bounding*(5, 1, {4}, 3, 2) (i.e. 21 < 23). Therefore, $Perm(M_{3,0})$ is equal to (5, 4, {1}, 3, 2), and $Perm(M_{3,1})$ is equal to (5, 1, {4}, 3, 2).

- **Rule 4:** Each cell $M_{i,j}$ contains only the last number assigned to its permutation $Perm(M_{i,j})$. For example, the cell $M_{0,1}$ has as permutation ({1, 3, 4, 5}, 2) and 2 is the last number assigned in this permutation. Therefore, the value of $M_{0,1}$ is equal to 2.

- **Rule 5:** If *k* is the last element assigned to $Perm(M_{i,j})$ and $Perm(M_{i,j})$ must be pruned, then $M_{i,j}$ is equal to $-k$. For example, since $Perm(M_{0,4})$ is equal to ({1, 2, 4, 5}, 3), the last integer assigned to this permutation is 3 and the permutation ({1, 2, 4, 5}, 3) must be pruned, then $M_{0,4}$ is equal to $-3$.

* **Selection positions (S):** This structure is a vector, denoted *S* in Figure 4.1, of size *N*, where *N* is the size of the permutation. Each position $S_i$ of this vector corresponds to the depth *i* of the PB&B tree. The value of $S_i$ is equal to the rank of the last permutation selected at the depth *i* of the tree. For example, (5, {1, 3, 4}, 2) is the last permutation selected in the depth 2 of the tree. This permutation occupies the rank 1 in this depth. Therefore, $S_2$ is equal to 1. This vector is called the selection vector because it keeps the ranks of the selected permutations at each level of the tree. As shown in Figure 4.1, each cell $S_i$ of the vector can be seen as a pointer to a cell of the row *i* of the matrix *M*.

* **Branching positions (B):** This structure is a vector, denoted *B* in Figure 4.1, of size *N*, where *N* is the size of the permutation. Each position $B_i$ of this vector
corresponds to the depth *i* of the PB&B tree. The value $B_i$ is equal to the position filled by the branching operation at the depth *i* of the tree. For example, at the depth 0 of the tree, integers are assigned to the position 4 of the permutation. Therefore $B_0$ is equal to 4. This vector is called branching vector because it keeps the filling positions at each level. Each cell $B_i$ of the vector is therefore associated with a row *i* of the matrix *M*.

* **Depth position (D):** This structure is an integer, denoted *D* in Figure 4.1, which saves the depth of the next permutation to be selected. For example, in Figure 4.1, the next permutation to select is (5, 4, {1}, 3, 2). In the PB&B tree, this permutation belongs to the depth 3. As a result, *D* is equal to 3. The value of *D* can also be viewed as a pointer to the last filled row of *M*.

* **Permutation structure (P) and its limits ($L_1$ and $L_2$):** This structure is composed of a vector, denoted *P* in Figure 4.1, and two integers, denoted $L_1$ and $L_2$. It allows to store any empty, full or partial permutation using the limits $L_1$ and $L_2$.

* **Bounding structure (C):** This structure, denoted *C* in the Figure 4.1, allows to store the lower bounds calculated by the bounding operation. The role of this structure is explained when the bounding operation is presented in Subsection 4.3.3.

## 4.3 PB&B Operations

<figure><img src="image_8.png" /><figcaption>Figure 4.2: Overview of the PB&B@CORE approach and its operations.</figcaption></figure>

The linked-list data structure is well suited for the conventional B&B algorithm. On the other hand, the IVM structure is exclusively designed for the PB&B algorithm, which is a
special case of the B&B algorithm. As another data structure is proposed, it is necessary to
adapt the operations of B&B to this structure.

4.3.1 Selection

<figure><img src="image_9.png" /><figcaption>Figure 4.3: Illustration of the IVM-based selection of PB&B.</figcaption></figure>

Like its linked-list counterpart, the selection operation returns a permutation, and writes
it in the structures *P*, *L*₁ and *L*₂. This permutation is obtained using four steps:

*   **Update of the depth integer D and the selection vector S:** This step works using five rules. These rules must be repeated until Rule 3 or Rule 5 indicates to stop the step. The five rules are given below:

- **Rule 1:** In the beginning, $S_D$ is incremented.

- **Rule 2:** If $M_{D,S_D}$ is negative, then $S_D$ is incremented. The negative value of $M_{D,S_D}$ means that $Perm(M_{D,S_D})$ is pruned.

- **Rule 3:** If $M_{D,S_D}$ is positive, then the step is stopped. In this case, a permutation is found and the pool is still not empty.

- **Rule 4:** If $S_D$ points after the last cell of the row, then the depth $D$ is decremented and $S_D$ is incremented.

- **Rule 5:** If $D$ points before the first row of the matrix, then the step is stopped. In this case, there is no permutation and the pool is empty.

In the left side of Figure 4.3, $S_D$ (i.e. $S_2$) points to $M_{2,0}$. After applying the five previous rules, the right side of Figure 4.3 is obtained. The application of these rules is done as follows:

- The step starts by incrementing $S_2$, for the first time, and points it to $M_{2,1}$.

- Because of the negative value of $M_{2,1}$, the $S_2$ is incremented, for the second time, to point to $M_{2,2}$.
- Since $M_{2,2}$ is also negative, $S_2$ is incremented, for the third time, to point to $M_{2,3}$.

- Because $M_{2,3}$ is located after the last cell of the row, $D$ is decremented and $S_D$ (i.e. $S_1$) is incremented.

- Because $S_1$ points to $M_{1,1}$, which is a positive value, the step is stopped and the right side of Figure 4.3 is obtained.

* **Copying the cells pointed by the selection vector S to the permutation vector P:** In addition to $S$, the copy of this step also takes into account the values of the branching vector $B$. In the general case, all the elements $M_{0,S_0}$, $M_{1,S_1}$, ..., $M_{D,S_D}$ in this order, and depending on whether each of the $B_0$, $B_1$, ..., $B_D$ are equal to zero or not, the elements are respectively placed at the first or the last free position of the vector $P$.

In Figure 4.3, the elements copied during this step are written in blue. In this example, there are only two copies to be made since the depth $D$ points to $S_1$:

- Since $S_0$ points to $M_{0,2}$, which is equal to 1, the integer 1 is copied to the vector $P$. Because $B_0$ is not equal to 0, the integer 1 is copied to the last free position of $P$, which is $P_4$.

- Since $S_1$ points to $M_{1,1}$, which equals 5, the element 5 is copied to the vector $P$. Because $B_1$ is equal to 0, the element 5 is copied to the first free position of $P$, which is $P_0$.

* **Copying the elements from the last filled row of M to the permutation vector P:** The integer $D$ always indicates the last filled row. On the right side of Figure 4.3, the elements of the last filled row are 4, 5, 2 and -3. Except the element 5, which is pointed by $S_D$, all these elements are still not copied into $P$. Therefore, the elements 4, 2 and 3 are copied into the three free positions of the vector $P$. The elements, copied during this step, are written in green in the Figure 4.3.

In the general case, this step copies all the absolute values of the elements $M_{D,0}$, $M_{D,1}$, ..., $M_{D,N-D}$, except the element $M_{D,S_D}$, to the free positions of the vector $P$.

* **Update of the limits $L_1$ and $L_2$:** These two pointers indicate the free positions of the permutation $P$. In other words, all $L_1$, $L_1 + 1$, ..., $L_2$ positions of $P$ are empty, and the elements $P_{L_1}$, $P_{L_1+1}$, ..., $P_{L_2}$ are free. The value of the pointer $L_1$ is equal to the number of zeros of the vector $B$, and the value of the pointer $L_2$ is equal to $N$ minus the number of non-zero values of $B$, where $N$ is the size of the permutation. In this example, $L_1$ and $L_2$ are equal to 1 and 3, respectively.
4.3.2 Branching

<figure><img src="image_10.png" /><figcaption>Figure 4.4: Example showing how branching operation works on IVM.</figcaption></figure>

The branching operation is much simpler than the selection operation. This operation is
done in three steps.

* **Step 1:** This step copies elements from the last filled row of *M* to the first empty row of *M*. During this operation, all the absolute values of the elements, except the one pointed to by *S*<sub>*D*</sub> are copied. For example, in the left part of Figure 4.4, the elements of *M*, written in red, namely 4, 2 and -3, are copied. In the right part of Figure 4.4, it is possible to see these copied elements written in blue.

* **Step 2:** The role of this step is to increment the pointer *D*. In Figure 4.4, the value of *D*, which was 1 in the left side of the figure, becomes 2 in the right side.

* **Step 3:** This step points $S_D$ to $M_{D,0}$ by writing 0 in $S_D$. As shown in Figure 4.4, $S_2$ points to $M_{2,0}$.

4.3.3 Bounding

<figure><img src="image_11.png" /><figcaption>Figure 4.5: Illustration of the IVM-based bounding of PB&B.</figcaption></figure>

The objective of the bounding operation is to update the elements of the last filled row of *M*, and the last filled position of *B*. Figure 4.5 shows an example of the state of the IVM structure before and after a bounding operation. This update is done using three steps:
* **Compute of all possible lower bounds:** In Figure 4.5, the permutations of the last filled row of the left IVM are written in red. The cells in this row correspond to those obtained after the branching of $(2, \{5, 1, 3\}, 4)$. Since the position $B_2$ is still empty, it is not known yet whether the elements $\{5, 1, 3\}$ are placed at the first empty position of the permutation $P$ or at the last empty position of this permutation. Depending on the value of $B_2$, there are two cases:

- **Case $B_2 = 0$:** This case corresponds to the case where one of the free elements is written in the first empty position of the permutation. As a result, the partial permutations of the last filled row correspond to $(2, 1, \{3, 5\}, 4)$, $(2, 3, \{1, 5\}, 4)$ and $(2, 5, \{1, 3\}, 4)$.

- **Case $B_2 = 3$:** This case corresponds to the case where one of the free elements is written in the last empty position of the permutation. In this case, the partial permutations of the last filled row are $(2, \{3, 5\}, 1, 4)$, $(2, \{1, 5\}, 3, 4)$ and $(2, \{1, 3\}, 5, 4)$.

For each of these six permutations, the obtained bound values are placed in cells $C_{1,0}$, $C_{1,1}$, $C_{1,2}$, $C_{2,0}$, $C_{2,1}$ and $C_{2,2}$. In this example, the value of these six bounds are assumed to be 50, 70, 80, 90, 85 and 75, respectively.

* **Choice of the bounding group:** The objective of this step is to choose either to place the free elements in the first free position or to place them at the last free position. In other words, the objective of this step is to determine if the value of $B_2$ is 0 or 3. In our algorithm, this choice is made by comparing the sum of the bounds of each group. Since the sum $(50 + 70 + 80 = 200)$ is smaller than $(90 + 85 + 75 = 250)$, the bounding operation chooses to place the free elements at the last empty position. This is why, in the right IVM of the figure, $B_2$ is equal to 3.

* **Sorting of the permutations according to the bound values:** In the second step, the selected partial permutations have the bounds 90, 85 and 75. These bounds correspond respectively to the filling of the free elements 5, 1 and 3 in the last empty position of the permutation $(2, \{5, 1, 3\}, 4)$. The objective of this third step is to sort these free elements 5, 1 and 3 according to the increasing values of their corresponding bounds. As shown in the right IVM of Figure 4.5, the elements are sorted according to the ascending order of their bounds.
4.3.4 Pruning

<figure><img src="image_12.png" /><figcaption>Figure 4.6: Example showing how pruning operation works on IVM.</figcaption></figure>

At the end of the bounding operation, the bounds of the permutations corresponding to
the cells of the last filled row of the IVM are calculated. The objective of the pruning
operation is to eliminate some permutation of this last filled row. As explained before, any
permutation, which has a higher bound than the cost of the best complete permutation
found so far, must be pruned.

In the left IVM of Figure 4.6, the last filled row contains three cells. These cells cor-
respond to the placement of the free elements 3, 1 and 5 in the last empty positions of
the permutation (2, {3, 1, 5}, 4). The obtained three sub-permutations are therefore respec-
tively (2, {1, 5}, 3, 4), (2, {3, 5}, 1, 4) and (2, {3, 1}, 5, 4). These permutations respectively
have the bounds 75, 85 and 90. As the best complete permutation found so far has a cost of
80, the two permutations (2, {3, 5}, 1, 4) and (2, {3, 1}, 5, 4) are pruned. In our approach,
pruning is done by multiplying cells by -1. Therefore, the elements in the last filled row
are equal to 3, -1 and -5 respectively, as shown by the right IVM of Figure 4.6.

4.3.5 Initialization and finalization

<figure><img src="image_13.png" /><figcaption>Figure 4.7: Example showing how initialization and finalization operations work on IVM.</figcaption></figure>
In addition to the four previous basic operations, any PB&B must also use an operation to initialize the IVM, and another operation to detect the end of PB&B. The left IVM of Figure 4.7 shows the status of the IVM at its initialization, and the right IVM of this figure shows the state of IVM at its end. As shown in the figure, the structure is initialized as follows:

* The value of depth *D* is equal to 0;

* The first position $P_0$ of the position vector is equal to 0;

* The first row of the matrix *M* is initialized to 1, 2, ..., *N*, where *N* is the size of the permutation;

* The limits $L_1$ and $L_2$ are initialized to the values -1 and *N*;

The right IVM of the figure shows the state of the structure at the end of the PB&B algorithm. When the IVM is in this state, the PB&B algorithm must be stopped. The value of depth *D* is equal to -1.

## 4.4 Experiments

### 4.4.1 Experimental protocol

* **Hardware testbed:** All the experiments were run on the computer Poincaré which belongs to Maison de la Simulation (Saclay, Paris). Each of Poincaré's 92 CPU nodes is composed of 2 8-core Intel Xeon Sandy Bridge E5-2670 processors running at 2.60 GHz and has 32 Gb of memory. Each of the 16 physical cores has 32 KB of L1 instruction cache, 32 KB of L1 data cache and 256 KB of L2 cache. Each of the 2 processors has 20 MB of L3 cache. The 32 GB of memory are spread across 2 NUMA nodes, one for each processor.

* **Software tools:** The serial LL-based PB&B and PB&B@CORE have been implemented in C++ and compiled using GCC 4.6 with `-O3` option. For each instance tested, the computational time spent in managing the pool of permutations is measured using the `clock_gettime` function with a nanosecond precision.

* **Problem instances:** In our experiments, we used the flowshop instances defined by Taillard [Tai93]. These standard instances are often used in the literature to evaluate the performance of methods that minimize the makespan. In the experiments of this chapter, we used the 10 instances defined with 20 jobs and 20 machines (these instances are named Ta021, Ta022, ..., and Ta030), and the 10 instances defined with 50 jobs and 10 machines (these instances are named Ta041, Ta042, ..., and
Ta050). The other instances are not used in our validation because they are either easy or difficult to solve with a serial algorithm.

* **Theoretical memory study:** It is possible to make a theoretical study of the maximum memory needed by the LL and IVM approaches to store their pool. For both approaches, the maximum size depends only on the number of jobs *N* of an instance. The size of the memory used by an IVM structure is always constant. In bytes, this size can be calculated using Equation (4.1).

$$ 
\begin{aligned}
Maximum - size(IVM) &= \left[ \frac{N(N+1)}{2} + 3N + 1 \right] \text{ bytes} \\
&= \left[ \frac{1}{2}N^2 + \frac{7}{2}N + 1 \right] \text{ bytes}
\end{aligned}
\quad (4.1) $$

Unlike IVM, the advantage of LL is to not require additional CPU calculation time for generating a permutation. In the LL approach, a conventional coding of a permutation is to write the list of scheduled jobs and the list of unscheduled jobs. By assuming that a job is encoded with 1 byte, the size of a permutation of an instance defined by *N* jobs is always equal to *N* bytes, and therefore, the size of a pool LL, which contains *X* permutations, is equal to *N* × *X* bytes. In LL, a pool reaches its maximum size when visiting the first complete permutation. At this moment, the pool contains *N* − 1 permutations with 1 job scheduled, *N* − 2 permutations with 2 jobs scheduled, ..., until *N* − (*N* − 1) (i.e. 1) permutation with *N* − 1 jobs scheduled. The maximum size of a pool LL can be calculated using Equation (4.2).

$$ 
\begin{aligned}
Maximum - size(LL) &= \left[ \sum_{i=1}^{i=N-1} (N-i) \right] \text{ permutations} \\
&= \left[ \sum_{k=1}^{k=N-1} k \right] \text{ permutations} \\
&= \left[ \frac{N(N-1)}{2} \right] \text{ permutations} \\
&= \left[ \frac{N^2(N-1)}{2} \right] \text{ bytes} \\
&= \left[ \frac{1}{2}N^3 - \frac{1}{2}N^2 \right] \text{ bytes}
\end{aligned}
\quad (4.2) $$

In terms of space, this means that the IVM structure has a maximum of $O(N^2)$ complexity, while the LL data structure has a maximum of $O(N^3)$ complexity. Therefore, IVM approach can be up to *N* better than the LL approach in terms of memory size.
Besides, Table 4.1 concretely shows the ratio between the two approaches. These ratios are computed using Equation (4.1) and Equation (4.2). The values of this table are given for different sizes of jobs (i.e. 20, 50, 100, 200 and 500) of Taillard instances.

<table><caption>Table 4.1: Comparison of PB&B@CORE and LL PB&B algorithms in terms of maximum memory.</caption><thead><tr><td>Instance size</td><td>LL max. memory size (Bytes)</td><td>IVM max. memory size (Bytes)</td><td>IVM/LL max. memory ratio</td></tr></thead><tbody><tr><td>20</td><td>3800</td><td>271</td><td>14.07</td></tr><tr><td>50</td><td>61250</td><td>1426</td><td>42.98</td></tr><tr><td>100</td><td>495000</td><td>5351</td><td>92.52</td></tr><tr><td>200</td><td>3980000</td><td>20701</td><td>192.27</td></tr><tr><td>500</td><td>62375000</td><td>126751</td><td>492.11</td></tr></tbody></table>

However, the comparison of Table 4.1, in terms of maximum memory of both approaches, is not the best indicator to get an idea about the ratio of memory sizes really used by LL and IVM approaches. It is therefore important to compare the IVM and LL structures in terms of their memory size really used during a resolution. The average size of IVM is constant and is the same than the value given in Equation (4.2), while the average size of LL can not be deduced from a theoretical study. It is therefore necessary to solve an instance to know its average memory size when the resolution uses an LL structure.

### 4.4.2 Obtained results

* **Memory evaluation:** Tables 4.2 and 4.3 give the average obtained sizes for the instances defined with 20 jobs and 50 jobs, respectively. In these two tables, the first column gives the name of the instance, the second column the average number of permutations when using LL, the third column the average size in bytes for LL, the fourth column the average size in bytes for IVM which is constant, and the last column the ratio between the sizes of LL and IVM. The last row of the table gives the average of the values of each column.

According to Table 4.1, the maximum expected ratio between the sizes of LL and IVM is equal to 20 for the instances defined by 20 jobs and 50 for the instances defined by 50 jobs. The experiments show that the average obtained ratios are respectively 9 and about 35. These results show that on average an IVM structure clearly occupies much less memory space than LL data structure.
<table><caption>Table 4.2: Comparison of the size of IVM and the average size of LL when solving the ten instances defined with 20 jobs.</caption><thead><tr><th>Instance</th><th>LL average size (Permutations)</th><th>LL average size (Bytes)</th><th>IVM size (Bytes)</th><th>LL/IVM size ratio</th></tr></thead><tbody><tr><td>Ta021</td><td>118</td><td>2360</td><td rowspan="10">271</td><td>8.74</td></tr><tr><td>Ta022</td><td>127</td><td>2540</td><td>9.41</td></tr><tr><td>Ta023</td><td>111</td><td>2220</td><td>8.22</td></tr><tr><td>Ta024</td><td>123</td><td>2460</td><td>9.11</td></tr><tr><td>Ta025</td><td>117</td><td>2340</td><td>8.67</td></tr><tr><td>Ta026</td><td>122</td><td>2440</td><td>9.04</td></tr><tr><td>Ta027</td><td>115</td><td>2300</td><td>8.52</td></tr><tr><td>Ta028</td><td>127</td><td>2540</td><td>9.41</td></tr><tr><td>Ta029</td><td>124</td><td>2480</td><td>9.19</td></tr><tr><td>Ta030</td><td>131</td><td>2620</td><td>9.70</td></tr><tr><td>Average</td><td>121.50</td><td>2430.00</td><td>271</td><td>9.00</td></tr></tbody></table>

<table><caption>Table 4.3: Comparison of the size of IVM and the average size of LL when solving the ten instances defined with 50 jobs.</caption><thead><tr><th>Instance</th><th>LL average size (Permutations)</th><th>LL average size (Bytes)</th><th>IVM size (Bytes)</th><th>LL/IVM size ratio</th></tr></thead><tbody><tr><td>Ta041</td><td>961</td><td>48050</td><td rowspan="10">1426</td><td>33.72</td></tr><tr><td>Ta042</td><td>980</td><td>49000</td><td>34.39</td></tr><tr><td>Ta043</td><td>1059</td><td>52950</td><td>37.16</td></tr><tr><td>Ta044</td><td>1085</td><td>54250</td><td>38.07</td></tr><tr><td>Ta045</td><td>824</td><td>41200</td><td>28.91</td></tr><tr><td>Ta046</td><td>1056</td><td>52800</td><td>37.05</td></tr><tr><td>Ta047</td><td>1038</td><td>51900</td><td>36.42</td></tr><tr><td>Ta048</td><td>1049</td><td>52450</td><td>36.81</td></tr><tr><td>Ta049</td><td>1094</td><td>54700</td><td>38.39</td></tr><tr><td>Ta050</td><td>993</td><td>49650</td><td>34.84</td></tr><tr><td>Average</td><td>1013.90</td><td>50695.00</td><td>1426</td><td>35.58</td></tr></tbody></table>

* **CPU Time evaluation:** As written in the previous subsection, it is clear that the IVM approach uses much less memory than the LL approach. However, unlike LL, IVM requires coding and decoding mechanisms of the permutations of the pool. A question then arises about the cost of IVM encoding and decoding mechanisms. Indeed, the gain in IVM memory should not be to the detriment of an additional computing cost to manage the pool. Therefore, the objective of this subsection is to compare the LL and IVM approaches in terms of pool management CPU time.

In our comparative study, the pool management time does not only include the CPU
time spent by reading and writing operations in the LL and IVM structures. This time
also includes the CPU time spent by the selection, pruning and branching operations.
In other words, the pool management time includes all the operations made in the
PB&B algorithm except the part of the bounding operation which computes the
bounds.

Table 4.4 gives the average time obtained for the instances defined with 20 jobs, and Table 4.5 for the instances defined with 50 jobs. In both tables, the first column gives the name of the instance, the second column the pool management CPU time of the LL approach, the third column the pool management CPU time of the IVM approach, and the last column the ratio between the pool management CPU times of LL and IVM. The last row of the table gives the average of the values of each column.

The previous paragraphs show that the IVM structure uses much less memory space
than the LL structure. As the information is encoded in IVM, unlike LL, it is intuitively
logical to expect that the management of a pool coded with IVM takes more time
than the management of a pool coded with LL. However, Table 4.4 shows that the
management of the IVM pool takes on average about 2.75 less CPU time than the
LL pool when solving instances defined by 20 jobs. In addition, Table 4.5 shows that
this pool management takes on average about 3 times less CPU time than the LL
pool when solving instances defined by 50 jobs.

These average ratios can be certainly explained by the adaptation made for the PB&B
operations. The bounding, selection, branching and pruning operations are more
optimized in the IVM approach compared to the LL approach. For example, the LL
branching operation involves dynamically creating a certain number of permutations,
and each permutation must contain a permutation almost similar to the permutation
which is branched. Unlike this LL operation, the IVM branching operation merely
copies the content of a matrix row to another row.

## 4.5 Conclusions

In this chapter, we proposed a new data structure, called IVM (Integer Vector Matrix), to
implement the pool of permutations generated by a PB&B algorithm for solving permuta-
tion optimization problems. The bounding, selection, branching and pruning operations
have been revisited to operate on this new data structure.

This new approach is validated using standard instances of the flowshop which is a
permutation problem presented in the previous chapter. This evaluation is performed in
terms of memory and CPU time usages. Experiments show that the use of IVM does not
only greatly reduce the used memory size, but also reduces the CPU time spent for the pool management. Indeed, compared to LL, these experiments show that on average the IVM structure (1) occupies 9 times less memory space for the instances defined with 20 jobs, (2) uses about 35 times less memory for the instances with 50 jobs, (3) manages the pool about 2.5 times faster for instances with 20 jobs, and (4) manages the pool about 3 times faster for instances with 50 jobs.

<table><caption>Table 4.4: Comparison of PB&B@CORE and LL-based PB&B algorithms in terms of the CPU time used for the management of the pool when solving the instances defined with 20 jobs.</caption><thead><tr><th>Instance</th><th>LL time</th><th>IVM time</th><th>LL/IVM ratio</th></tr></thead><tbody><tr><td>Ta021</td><td>411.09</td><td>149.66</td><td>2.75</td></tr><tr><td>Ta022</td><td>538.89</td><td>195.10</td><td>2.76</td></tr><tr><td>Ta024</td><td>393.28</td><td>146.01</td><td>2.69</td></tr><tr><td>Ta026</td><td>1425.14</td><td>537.28</td><td>2.65</td></tr><tr><td>Ta027</td><td>767.72</td><td>281.55</td><td>2.73</td></tr><tr><td>Ta028</td><td>133.54</td><td>47.71</td><td>2.80</td></tr><tr><td>Ta029</td><td>424.42</td><td>153.27</td><td>2.77</td></tr><tr><td>Ta030</td><td>76.45</td><td>28.10</td><td>2.72</td></tr><tr><td>Average</td><td>521.32</td><td>192.34</td><td>≈ 2.75</td></tr></tbody></table>

<table><caption>Table 4.5: Comparison of PB&B@CORE and LL-based PB&B algorithms in terms of the CPU time used for the management of the pool when solving the instances defined with 50 jobs.</caption><thead><tr><th>Instance</th><th>LL time</th><th>IVM time</th><th>LL/IVM ratio</th></tr></thead><tbody><tr><td>Ta041</td><td>7.45</td><td>2.44</td><td>3.05</td></tr><tr><td>Ta042</td><td>3235.67</td><td>1030.24</td><td>3.14</td></tr><tr><td>Ta043</td><td>690.58</td><td>228.17</td><td>3.03</td></tr><tr><td>Ta044</td><td>3.59</td><td>1.30</td><td>2.76</td></tr><tr><td>Ta045</td><td>19.60</td><td>6.08</td><td>3.22</td></tr><tr><td>Ta046</td><td>34.43</td><td>11.54</td><td>2.98</td></tr><tr><td>Ta047</td><td>153.76</td><td>51.66</td><td>2.98</td></tr><tr><td>Ta048</td><td>66.59</td><td>22.14</td><td>3.01</td></tr><tr><td>Ta049</td><td>4.74</td><td>1.78</td><td>2.66</td></tr><tr><td>Ta050</td><td>1396.59</td><td>452.31</td><td>3.09</td></tr><tr><td>Average</td><td>561.30</td><td>180.77</td><td>≈ 3.00</td></tr></tbody></table>
CHAPTER 5

# Multi-core Interval-based permutation B&B

## Contents

<table><tr><td>5.1</td><td>Introduction</td><td>47</td></tr><tr><td>5.2</td><td>Factoradic intervals</td><td>49</td></tr><tr><td>5.2.1</td><td>Number systems</td><td>49</td></tr><tr><td>5.2.2</td><td>Factorial number system</td><td>51</td></tr><tr><td>5.3</td><td>Interval-based work stealing</td><td>57</td></tr><tr><td>5.3.1</td><td>Work unit: interval of factoradics</td><td>57</td></tr><tr><td>5.3.2</td><td>Work unit communication</td><td>58</td></tr><tr><td>5.3.3</td><td>Victim selection policies</td><td>60</td></tr><tr><td>5.3.4</td><td>Granularity policies</td><td>63</td></tr><tr><td>5.4</td><td>Experiments</td><td>64</td></tr><tr><td>5.4.1</td><td>Experimental protocol</td><td>64</td></tr><tr><td>5.4.2</td><td>Pool management evaluation</td><td>66</td></tr><tr><td>5.4.3</td><td>Interval-based work stealing evaluation</td><td>67</td></tr><tr><td>5.5</td><td>Conclusions</td><td>69</td></tr></table>

## 5.1 Introduction

When running PB&B@CORE, the values of the position vector *S* are continuously updated. As described in Chapter 4, when the end of row is reached the algorithm backtracks to the previous level. Therefore, at level *D* = 0, 1, ..., *N* − 1 the value of *S*<sub>*D*</sub> is bounded by *S*<sub>*D*</sub> < *N* − *I*. The vector is equal to 00000 when the algorithm points to the first complete permutation of the PB&B tree, and equal to 43210 when it points to the last complete permutation. Between these values, the vector successively takes the values 00010, 00100, 00110, 00200, 00210, ..., 43210. For each of these values, the algorithm points to a different complete permutation. There are 120 possible values because there
are 120 complete permutations (i.e. 5!). These 120 position vector values correspond to the numbering of the 120 complete permutations using a special numbering system, called factorial number system [Knu97]. The factorial number system, also called factoradic, is a mixed radix numeral system adapted to numbering permutations. For this number system, the French term used in 1888 [Lai88b] is *numération factorielle*.

In a multi-threaded approach, it is possible to ask a thread to explore only the complete permutations with factoradic numbers in the interval $[a, b]$, where *a* and *b* are factoradic numbers. In order to explore only these numbers, the thread initializes the vector position *S* to factoradic number *a*, launches the PB&B@CORE algorithm, and stops this algorithm when the number of the vector position *S* is equal to *b*. In this approach, the work unit of a PB&B thread is therefore an interval of factoradics. When a thread finishes exploring its interval, this thread can steal from another thread a portion of its interval in order to explore this portion. In a PB&B@CPU approach, these PB&B threads can be on the same CPU, or on multiple CPUs but sharing the same memory.

Figure 5.1 gives an overview of the PB&B@CPU approach, based on units of work defined by factoradics intervals. In this approach, only one thread is started at the beginning of a resolution. The role of this thread is (1) to initialize some data structures of the algorithm, (2) to launch a number of PB&B threads to explore the entire interval $[0, N![$, and (3) to restitute the complete permutation found, when the PB&B threads have finished exploring $[0, N![$. Compared to the PB&B@CORE approach, this new PB&B@CPU approach has an additional operation to allow the work stealing between the PB&B threads.

<figure><img src="image_14.png" /><figcaption>Figure 5.1: Overview of the PB&B@CPU approach and its operations.</figcaption></figure>

This chapter is divided into three sections. Section 5.2 section gives an overview of
the factorial number system, on which the work units of our approach are based. Section 5.3 describes the work stealing of the PB&B@CPU approach. Section 5.4 presents the experiments carried out to validate the PB&B@CPU approach using the flowshop problem.

## 5.2 Factoradic intervals

As the work unit of PB&B@CPU is defined by an interval of factoradic numbers, this section recalls some concepts related to number systems and describes the factorial number system.

### 5.2.1 Number systems

A number system, called also numeral system or system of numeration, is a writing system or a mathematical notation in order to express and represent a set of numbers using symbols in a consistent manner. A number system is defined by its digits, its bases, also called radixes, and its place values. This subsection reminds these concepts before their definitions for the factorial system are given in Subsection 5.2.2.

#### Digits:

Etymologically, the word *digit* comes from ancient Latin word *digit* which means fingers. Therefore, this word is related to the decimal system where ten digits are used like the ten fingers. However, the word digit is used nowadays for all other number systems including the binary system where the word bit is more appropriate. A number is a sequence of digits which can have an arbitrary length. Each digit, in a number system, represents an integer. In the decimal and the hexadecimal systems, for instance, the digits 1 and A represent the integers one and ten, respectively. In Roman numerals, where seven symbols are used, each symbol represents also a different integer as shown in Table 5.1.

<table><caption>Table 5.1: Roman digits and their values.</caption><thead><tr><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>I</td><td>1</td></tr><tr><td>V</td><td>5</td></tr><tr><td>X</td><td>10</td></tr><tr><td>L</td><td>50</td></tr><tr><td>C</td><td>100</td></tr><tr><td>D</td><td>500</td></tr><tr><td>M</td><td>1,000</td></tr></tbody></table>
There is a particular number system defined only with one digit. This simplest numeral
system is called the unary numeral system, and can be used to represent all natural num-
bers. In order to represent any number *N*, an arbitrarily symbol, which represents the
integer 1, is simply repeated *N* times. This system is often used in tallying. For example,
using the tally mark |, the number 5 is represented as ||| |. Unlike multiplication or divi-
sion, the other operations, namely addition, subtraction and comparison, are particularly
simple to be implemented in the unary system. Compared to other numeral systems, the
unary system is not used in practice for large calculations but can be convenient for small
operations, like representing a number with fingers. Unary system is used in some data
compression algorithms such as Golomb coding [Gol66].

Place value:

Roman system, used in ancient Rome, employs combinations of letters from the Latin
alphabet in order to express numbers. For example, the first ten numbers can be expressed
as follows: I, II, III, IV, V, VI, VII, VIII, IX, X. In this system, numbers are written by combin-
ing symbols and adding or subtracting the values these symbols. For example, XIII means
thirteen by adding a ten and three ones, and IX means nine by subtracting one from ten.
Unlike the decimal system, there is no zero in Roman system and symbols do not represent
tens or hundreds according to their positions. Therefore, unlike Roman system, which is
not a positional system, decimal system is a positional numeral system.

Ancient number systems, like Roman system, were not positional, and all of the number
systems most commonly used today, like binary system, are positional systems. Place value
is a positional system of notation in which the position of a number determines its value.
In other words, the value of a number in such system is determined not just by the digits
but also by the positions of each of the digits. For example, all place values, also called
order of magnitudes, in the unary system are equal to 1, the places values of decimal are
powers of ten, like ones place, tens places, hundreds place, etc. One of the advantages of
positional notation is the use of the same symbols for different order of magnitudes. This
greatly simplifies arithmetic operations.

Radix:

Etymologically, the word radix is a Latin word for the word root, and root is a synonym
for base in the arithmetical sense. In a positional numeral system, the radix is the number
of unique digits, including zero, used to represent numbers. For example, the radix is
ten for decimal system since this system uses ten digits from 0 through 9 in order to
represent its numbers. In a positional numeral system, the number X and the radix Y are
conventionally written as $(X)_Y$. However, the radix can be implicitly assumed and not written for some systems like decimal or unary systems.

Mixed radix numeral systems are non-standard positional numeral systems. Unlike most common systems, where the base is similar to all positions, the numerical base can vary from position to position. Such representation is used when a value or a number is written using a sequence of units that are each a multiple of the next smaller unit. For example, this type of number systems can be found when expressing time. A time of 10 hours, 30 minutes, and 50 seconds might be expressed as 10 : 30 : 50 in mixed radix notation where the radix of the second and third positions is 60 and the radix of the first position is 24.

In positional fixed radix number system, where the base $R$ is fixed, each digit $a_i$ in any number $(a_{n-1}...a_0)_R$ is an integer in the range 0 to $(R-1)$ and the number is interpreted as shown in Equation (5.1).

$$ (a_{n-1}...a_0)_R = a_{n-1}R^{n-1} + ... + a_1R^1 + a_0R^0 \quad (5.1) $$

Since Equation (5.1) is a polynomial in $R$, fixed radix system can be also called polynomial system. The decimal and binary systems are both fixed radix systems, with a radix of 10 and 2, respectively. Fractional values can also be represented with the same polynomial notation.

$$ 0.a_1a_2...a_n = a_{-1}R^{-1} + a_2R^{-2} + ... + a_nR^{-n} \quad (5.2) $$

In mixed-base or radix number system, the digit $a_i$ in any number belongs to the interval 0 to $R_i$, where $R_i$ is not the same for all the values of $i$. The number is then interpreted as shown in Equation (5.3).

$$ a_{n-1}...a_0 = (...((a_{n-1} R_{n-1}) + a_{n-2})R_{n-2} + ... + a_1)R_0 + a_0 \quad (5.3) $$

### 5.2.2 Factorial number system

Factorial system, also called factoradic, is a mixed radix number system which is well adapted for numbering permutations. This system is not named like most numeral systems. For example, unary, binary and decimal are named like this because their radixes are one, two and ten, respectively. Unlike these systems, the factorial system is named according to its place value instead of its mixed radix. The term factorial number system is used the first time recently in 1998 [Knu98] while the French name *numération factorielle* is first used in 1888 [Lai88c]. The term factoradic appears to be much more recent [McC03]. General properties of mixed radix number systems also apply to the factorial system.
As explained in Table 5.2, the $i^{th}$ digit from the right has base $i$ and the place value $i!$. Therefore, the $i^{th}$ digit must be less than $i$. And in order to compute the value of a number, the value of the $i^{th}$ digit must be multiplied by $i!$. From this, it follows that the rightmost digit is always 0, the second can be 0 or 1, the third 0, 1 or 2, and so on.

<table><caption>Table 5.2: Factorial system and its radixes, place values and digits for the seven first positions.</caption><thead><tr><th>Place</th><th>...</th><th>7<sup>th</sup></th><th>6<sup>th</sup></th><th>5<sup>th</sup></th><th>4<sup>th</sup></th><th>3<sup>rd</sup></th><th>2<sup>nd</sup></th><th>1<sup>st</sup></th></tr></thead><tbody><tr><td>Radix/base</td><td>...</td><td>7</td><td>6</td><td>5</td><td>4</td><td>3</td><td>2</td><td>1</td></tr><tr><td>Place value</td><td>...</td><td>6! = 720</td><td>5! = 120</td><td>4! = 24</td><td>3! = 6</td><td>2! = 2</td><td>1! = 1</td><td>0! = 1</td></tr><tr><td>Allowed digits</td><td>...</td><td>0<br>1<br>2<br>3<br>4<br>5<br>6</td><td>0<br>1<br>2<br>3<br>4<br>5</td><td>0<br>1<br>2<br>3<br>4<br>5</td><td>0<br>1<br>2<br>3<br>4<br>5</td><td>0<br>1<br>2<br>3<br>4<br>5</td><td>0<br>1<br>2<br>3<br>4<br>5</td><td>0<br>1<br>2<br>3<br>4<br>5</td></tr></tbody></table>

It is possible to define factorial numbers without writing the rightmost digit since it is always equal to zero. In our report, a factorial number representation will be flagged by a subscript "!", so for instance $(322110)!$ stands for $(3_5!2_4!2_3!1_2!1_1!0_0!)$). In principle, the factorial system may be extended to represent fractional numbers. However, the natural extension of place values $(-1)!$, $(-2)!$, etc. are undefined. In factorial system, the symmetric choice of radix values $n = 0, 1, 2, 3, 4$, etc. after the point may be used instead. The corresponding place values are therefore $1/(0!)$, $1/(1!)$, ..., $1/(n!)$, etc. In our work, factorial fractional numbers are not used.

## Decimal to factoradic:

When converting a decimal number into its factorial representation, digits are produced from right to left. This conversion consists in repeatedly dividing the number by the radixes $1, 2, 3$, etc. After each division, the remainder should be considered as the digit. The division operation continues with the integer quotient until this quotient becomes 0. Let assume a decimal number $D = 349$ to convert into a factorial number. This conversion
is done using successive Euclidean divisions as shown in Equation (5.4).

$$ 
\begin{aligned}
349 &= 1 \times 349 + 0 \\
349 &= 2 \times 174 + 1 \\
174 &= 3 \times 58 + 0 \\
58 &= 4 \times 14 + 2 \\
14 &= 5 \times 2 + 4 \\
2 &= 6 \times 0 + 2
\end{aligned} 
\qquad (5.4) $$

Euclidean division is the operation of division of two integers, which produces a quotient and a remainder. Each line of Equation (5.4) represents an Euclidean division $Q_i = i \times Q_{i+1} + R_i$ such as:

* $Q_{i+1}$ is the quotient of the division $Q_i$

* $R_i$ is the remainder of this division

* $Q_1 = A$ is the decimal number to convert

* $Q_{n+1}$ is always equal to zero and is the last quotient

* $R_1$ is always equal to zero and is the first remainder

The factorial representation $F$ of $A$ is equal to concatenation of all reminders, as shown in Equation (5.5), is $F = (R_n...R_2R_1)! = 242010$.

$$ 
\begin{aligned}
349 &= 1 \times 349 + 0 \\
&= 1 \times (2 \times 174 + 1) + 0 \\
&= 1 \times (2 \times (3 \times 58 + 0) + 1) + 0 \\
&= 1 \times (2 \times (3 \times (4 \times 14 + 2) + 0) + 1) + 0 \\
&= 1 \times (2 \times (3 \times (4 \times (5 \times 2 + 4) + 2) + 0) + 1) + 0 \\
&= 1 \times (2 \times (3 \times (4 \times (5 \times (6 \times 0 + 2) + 4) + 2) + 0) + 1) + 0 \\
&= 2\,5! + 4\,4! + 2\,3! + 0\,2! + 1\,1! + 0\,0! \\
&= (242010)!
\end{aligned} 
\qquad (5.5) $$

Algorithm 1 explains the used method to perform this conversion.
**Algorithm 1 DECIMAL-TO-FACTORIAL(D)**

1: Place ← 1

2: **while** D ≠ 0 **do**

3: $F_{i-1} \leftarrow D \bmod i$

4: D ← D **div** i

5: i ← i + 1

6: **end while**

7: **return** F

**Factoradic to decimal:**

Converting a factorial number to a decimal number is simpler. Let $(R_{n-1}...R_1R_0)!$ a factorial number. In order to have its decimal equivalent, it suffices to calculate the value of the polynomial $\sum_{i=0}^{n-1} R_i i!$. The conversion of a factorial number to its decimal equivalent is therefore a sum of multiplications, while the conversion of a decimal number to its factorial equivalent is a concatenation of divisions. Algorithm 2 explains this sum of multiplications.

**Algorithm 2 FACTORIAL-TO-DECIMAL(F)**

1: Place-value ← 1

2: D ← 0

3: **for** i ← 1 **to** n **do**

4: Place-value ← Place-value × i

5: D ← D + $F_i$ × Place-value

6: **end for**

7: **return** D

**Addition:**

Let assume two factorial numbers $A = (A_{n-1}...A_1A_0)!$ and $B = (B_{n-1}...B_1B_0)!$ to be added and both having the size $n$. As these numbers are factorial, both conditions $\forall i, A_i \le i$ and $\forall i, B_i \le i$ are always true. The addition of $A$ and $B$ would be simple if $\forall i, A_i + B_i \le i$. In this case, the result of the addition would be $C = (C_{n-1}...C_1C_0)!$ such as $\forall i, C_i = A_i+B_i$. However, the condition $\forall i, A_i + B_i \le i$ is not always satisfied. Let assume $i$ such as
$A_i + B_i > i$. The value of $C_i$ can be calculated as explained in Equation (5.6).

$$
\begin{align*}
C_i &= A_i i! + B_i i! \\
&= [A_i + B_i]i! \\
&= [(i+1) - (i+1) + (A_i + B_i)]i! \tag{5.6} \\
&= (i+1)i! + [-(i+1) + (A_i + B_i)]i! \\
&= (i+1)! + [(A_i + B_i) - (i+1)]i!
\end{align*}
$$

Therefore, there is a carry +1 for the calculation of $C_{i+1}$, and $C_i = (A_i + B_i) - (i + 1)$ since the rule of Equation (5.7) is always true.

$$ [(A_i \le i) \text{ and } (B_i \le i) \text{ and } (A_i + B_i > i)] \Rightarrow 0 \le [(A_i + B_i) - (i+1)] \le i \quad (5.7) $$

Equation (5.8) gives the value of $C_i$ for any values of $A_i$ and $B_i$.

$$ C_i = \begin{cases} (A_i + B_i) & \text{If } (A_i + B_i) \le i \\ (A_i + B_i) - (i+1) & \text{Otherwise} \end{cases} \tag{5.8} $$

Algorithm 3 explains the method used for the addition of A and B. The algorithm proceeds from the least significant position to the most significant one, in other words, from right to left. So, it is possible to check whether there is a carry +1 when computing $C_i$. If this is the case, this carry is taken into account when computing $C_{i+1}$.

**Algorithm 3 FACTORIAL-ADDITION(A, B)**

1: **for** i ← 0 **to** (n-1) **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;C<sub>i</sub> ← C<sub>i</sub> + A<sub>i</sub> + B<sub>i</sub>

3: &nbsp;&nbsp;&nbsp;&nbsp;**if** C<sub>i</sub> > i **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;C<sub>i</sub> ← A<sub>i</sub> - (i + 1)

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;C<sub>i</sub> + 1 ← 1

6: **end if**

7: **end for**

8: **return** C

**Subtraction:**

Let assume two factorial numbers $A = (A_{n-1}...A_1A_0)!$ and $B = (B_{n-1}...B_1B_0)!$. As these numbers are factorial, both conditions $\forall i, A_i \le i$ and $\forall i, B_i \le i$ are always satisfied. The objective of this subsection is to explain how to perform a subtraction $A - B$ when
assuming $A \ge B$. To facilitate the explanation, both numbers are assumed to have the same size $n$. If the size of $B$ is smaller than $A$, then it is possible to complete $B$ with zeros at the left.

The subtraction $A - B$ would be simple if $\forall i, A_i \ge B_i$. In this case, the result of the subtraction would be $C = (C_{n-1}...C_1C_0)!$ such as $\forall i, C_i = A_i - B_i$. However, the condition $\forall i, A_i \ge B_i$ is not always satisfied. Let assume $i$ such as $A_i < B_i$. The value of $C_i$ can be computed as explained in Equation (5.9).

$$ 
\begin{aligned}
C_i &= A_i i! - B_i i! \\
&= [A_i - B_i]i! \\
&= [-(i+1) + (i+1) + (A_i - B_i)]i! \\
&= [-(i+1)]i! + [(i+1) + (A_i - B_i)]i! \\
&= [-1](i+1)! + [(i+1) + (A_i - B_i)]i!
\end{aligned} 
\qquad (5.9) $$

Therefore, there is a carry $-1$ to be taken into account when computing $C_{i+1}$, and $C_i = (A_i + B_i) - (i + 1)$ since the rule of Equation (5.10) is always true.

$$ [(A_i \le i) \text{ and } (B_i \le i) \text{ and } (A_i < B_i)] \Rightarrow 0 \le [(A_i - B_i) + (i+1)] \le i \qquad (5.10) $$

Equation (5.11) gives the value $C_i$ for any values of $A_i$ and $B_i$.

$$ C_i = \begin{cases} (A_i - B_i) & \text{If } (A_i - B_i \ge 0) \\ (A_i - B_i) + (i+1) & \text{Otherwise} \end{cases} \qquad (5.11) $$

Algorithm 4 explains the operation of the subtraction $A - B$. Like the addition, the subtraction algorithm proceeds from the least significant position to the most significant position. It is possible to check if there is a carry $-1$ when computing $C_i$. If this is the case, this carry is taken into account when computing $C_{i+1}$.

**Algorithm 4 FACTORIAL-SUBTRACTION(A, B)**

1: **for** $i \leftarrow 0$ **to** $(n-1)$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp; $C_i \leftarrow C_i + A_i - B_i$

3: &nbsp;&nbsp;&nbsp;&nbsp; **if** $C_i < 0$ **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $C_i \leftarrow C_i + (i+1)$

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $C_{i+1} \leftarrow -1$

6: &nbsp;&nbsp;&nbsp;&nbsp; **end if**

7: **end for**

8: **return** C
5.3 Interval-based work stealing

A work stealing strategy can be defined by three major components: a work unit definition, a victim selection policy and a granularity policy. The victim selection policy determines how a thief thread *R* chooses its victim *S*. The granularity policy determines the amount and which part of work thread *R* steals from thread *S*. An ideal victim selection strategy is one which (1) chooses the victim *S* with the largest amount of work, (2) and makes this choice as rapidly as possible. A good granularity policy reduces the number of work stealing operations, meaning that it allows both victim and thief to work as long as possible without initiating another work stealing event. This section presents respectively the work unit of our approach, its victim selection policy and its granularity policy.

5.3.1 Work unit: interval of factoradics

Dynamic load balancing is well suited for multi-core parallel PB&B algorithms. All threads
must share their work. To the best of our knowledge, except in rare works such as
[MMT07a], work units exchanged between threads (or processes) are sets of nodes. Since
the nodes explored in our approach are numbered, it is possible to define an interval of
node numbers as the work unit. Let's assume the interval that must be explored by a PB&B
algorithm is [00000, 43210]. It is therefore possible to have two threads T1 and T2 such as
T1 explores [00000, X] and T2 explores [X, 43210]. If T2 ends exploring its interval before
T1, then T2 sends a request to T1 to recover a portion of its interval. Therefore, T1 and
T2 can exchange their interval portions until the exploration of the whole [00000, 43210]
interval.

In order to emphasize on the difference between a conventional set of nodes and an
interval of factoradics, it should be noted that the notion of work is slightly different. In-
deed, a work unit [A, B] does not necessarily contain any nodes to branch. In that sense it
rather represents a potential amount of work. The processing of an interval of factoradics
amounts to scanning the corresponding part of the search space. In contrast, the set of
pending permutations contains only nodes that have actually been generated and not been
removed from the tree.

In [MMT07a], we have presented an original load balancing strategy where the work unit
is an interval of integers. Compared to [MMT07a], the new strategy presented in this
chapter brings two new contributions. The first contribution is that the PB&B is based on
a matrix of integers instead of a linked-list of permutations of integers, and the second
is that it is not necessary to use the fold and unfold operations defined in [MMT07a] to
transform an interval into a linked-list of nodes and vice versa. Besides, the integers of
[MMT07a] are coded using the decimal system, while in this chapter, the integers are coded using the factoradic system.

### 5.3.2 Work unit communication

In this subsection two alternative protocols for the communication of work units between PB&B@CPU threads are described. An illustration of both procedures is shown in Figure 5.2.

* **Work unit splitting at an arbitrary position:** An integer interval $[B, E[$ can be split at any integer $C$ that is convex linear combination of $B$ and $E$, i. e. $C = [(1 - \alpha)A + \alpha B]$, $0 < \alpha < 1$. The computation of a splitting point $C$ can be performed either by using decimal arithmetic operations or by implementing elementary arithmetic operations for factoradic numbers. The granularity of this procedure is controlled by the value of $\alpha$. In the example shown in Figure 5.2a, the interval of the victim is split in two parts of equal size, setting $\alpha = 0.5$.

In both cases it is necessary to have a procedure that initializes the IVM structure at an arbitrary valid position vector $V = C$. For that, it is not enough to build the matrix by iterative application of the branching operation, selecting the jobs pointed by $V$, because the information about pruned nodes is lost.

Initialization of the IVM structure at any position vector $V = C$ can be achieved as follows. Starting from $I = 0$ all nodes pointed by $V$ are expanded, bounded and pruned until the last line is reached, i. e. until $I = n$. In other words, $n$ iterations of a modified PB&B algorithm, selecting the permutations indicated by $V$, are performed. As this initialization procedure involves the bounding operation, initialization overhead can become significant.

A first observation can be made: this initialization process can actually be stopped when $V$ points to a pruned permutation. This reduces the number of initialization steps considerably. A second observation allows to further decrease the amount of time spent in initialization: Suppose that an IVM structure was used to explore an interval $[B, E[$ and that it has reached the end of this task, that is $B = E$. Let $l$ be the level at which the exploration stopped. Now we want to initialize this IVM at a new position $V = \tilde{B}$. If $B[i] = \tilde{B}[i]$ for $i = 0, 1, ..., k$ with $k < l$, then these first $k$ lines of the matrix $M$ are already correctly initialized. The initialization process described before can thus begin at $I = k$.
<figure><img src="image_15.png" /><figcaption>(a) Arithmetic interval splitting</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Subtree-based interval splitting</figcaption></figure>

<caption>Figure 5.2: Illustration of work unit splitting.</caption>

* **Subtree-based work unit splitting:** In [Ler15b] a method for exchanging work units between PB&B@CPU threads without initialization is proposed. This procedure is based on splitting the interval of the victim IVM directly in its factoradic form, without converting it to decimals. In addition to the new position and end vectors, an initialized matrix is transferred from the victim IVM to the thief IVM. Figure 5.2b illustrates this procedure. The transfer of a work unit from IVM *S* (Sender) to IVM *R* (Receiver) can be performed as follows.

1. Let $[B_S, E_S]$ be the interval to split. Let *l* be the smallest index such that $B_S[l] \neq E_S[l]$.

2. For $i = 0, 1, ..., l-1$, copy row *i* of IVM from *S* to *R* (position, end, direction vectors and matrix).

3. Split the tree at level *l* by choosing *C* such that $B_S[l] < C \le E_S[l]$. For the receiving IVM *R*, set $B_R[l] = C$ and $E_R[l] = E_S[l]$. For the sending IVM *S*, $E_S[l] = C - 1$.
4. For $i = l + 1, \dots, n - 1$, set $E_S[i] = n - i - 1$ and $B_R[i] = 0$.

5. The receiving IVM $R$ is now initialized and exploration can start at level $l$. Therefore $I_R$ is set to $I_R = l$.

This initialization method proposed by [Ler15b] requires no additional computation of bounds. Therefore it reduces the overhead induced by work stealing operations, compared to the previously introduced initialization procedure. Another important advantage of this method is that it avoids redundant computations. Indeed, if intervals are split at arbitrary points, some permutations may be branched redundantly. In this example, thread T1 explores interval $[0000, 1100[$ and thread T2 explores $[1110, 2110[$, which may cause redundant computation of bounds along the frontier $11XX$.

However, there are also disadvantages. Especially in a distributed memory setting the size of data transfers should be kept low. Indeed, this is the primary motivation for using interval encoded work units [MMT07c]. Also, in the subtree-based interval splitting, the granularity is controlled in a coarser way because full subtrees are communicated. This method is similar to stack splitting strategies where nodes (and implicitly the subtrees rooted in these nodes) are transferred from non-empty to empty pools. Dividing intervals by the first method allows a finer control of granularity.

We use the second (subtree-based) work splitting method for communicating work units via shared memory and the first method in distributed memory contexts.

### 5.3.3 Victim selection policies

In this Subsection, four victim selection policies are described. Two of them, the random and ring policies have a low computational complexity and require no access to shared data structures or knowledge about the global workload repartition. The two other selection policies, namely the largest and the honest policies, use simple heuristics which aim at selecting a stealing victim holding a large or difficult piece of work. These victim selection policies use the available information about the workers' activity and require some additional computation and protected accesses to shared data structures.

* Ring victim selection policy: In this deterministic policy, threads are connected to each other with an unidirectional ring. A thread $R$ always steals its precedent thread $R'$. If the thread $R$ is different from the thread 1, then the thread $R'$ is equal to the thread $R - 1$. Otherwise, the thread $R'$ is equal to the thread $T$, where $T$ is the number of threads. In this policy, the work stealing operation of a thread $R$ is a blocking event when the thread $R'$ has no work. In this case, the work stealing operation will be satisfied when the thread $R'$ will receive work. This policy is
also used, for instance in [KRR88]. If the thread numbering is matched with the
underlying architecture the deterministic nature of this policy can be used to reduce
communication costs.

The average distance to another thread is $T/2$. For a large number of threads this distance should be decreased and each thread should have more than one neighbor. Moreover the possibility of two threads selecting each other mutually must be excluded. For instance, the lifeline scheme [SKK<sup>+</sup>11], based on cyclic hypercubes, satisfies these properties and has been shown to be scalable to clusters of several thousand cores. The ring strategy corresponds indeed to the simplest form of the lifeline scheme, to which it could be extended when a large number of threads is used. As shown in function *choose-ring*, the cost of the victim selection function is very low. The main issue of this strategy is that work units may not propagate fast enough through the ring. No locking or access to shared data structure is required by this strategy: the only information an idle worker needs to select a victim is its own thread/process ID.

*   **Random victim selection policy:** The random selection policy is provably efficient [BL99] and the most frequently used in the literature. In this policy, a thread victim $R'$ is randomly selected when a thread $R$ initiates a work stealing operation. Unlike the ring policy, this work stealing operation is not a blocking event. In other words, the thread $R$ continues to choose other threads randomly as long as it does not find a thread with a non-empty interval or linked-list. Besides its own thread/process ID a thief only needs to access a variable that indicates the state of the randomly selected victim.

*   **Largest victim selection policy:** In a PB&B algorithm, it is often impossible to determine the hardness of a work item. This policy is based on a simple heuristic to choose the thread with the most difficult work to finish. Indeed, the largest policy assumes that probably the larger the size of a work is, the more difficult this work will be. Therefore, this policy computes the amount of work of each thread, chooses the thread with the biggest size, and returns the rank $R'$ of this thread. In the linked-list based approach, the size of a linked-list is equal to the number of its nodes, and in the interval-based approach, the size of an interval $[A, B]$ is equal to $B - A$. As shown in function *choose-largest*, this policy has a higher computational complexity than the three other victim selection policies. In particular a thief requires locked accesses the *length* or *size* variable of each busy worker. Although this polling may compromise the scalability of this strategy, good results for this policy are reported in [ASW<sup>+</sup>14].
<caption>**Algorithm 5** Pseudocode of the victim selection policies.</caption>

1: **function** CHOOSE-THREAD(R, strategy)
2:     **switch** strategy **do**
3:         **case** RING:
4:             **return** CHOOSE-RING(R)
5:         **case** RANDOM:
6:             **return** CHOOSE-RANDOM(R)
7:         **case** LARGEST:
8:             **return** CHOOSE-LARGEST(R)
9:         **case** HONEST:
10:         **return** CHOOSE-HONEST(R)
11: **end function**
12: **function** CHOOSE-RING(R)
13:     **return** (R − 1)%T
14: **end function**
15: **function** CHOOSE-RANDOM(R)
16:     **while** true **do**
17:         R' ← random(1,T)
18:         **if** (has-work(R') AND (R'≠R)) **then**
19:             **return** R'
20:         **end if**
21:     **end while**
22: **end function**
23: **function** CHOOSE-LARGEST(R)
24:     max-size ← 0
25:     **for all** R" ∈ {1, 2, ..., T} AND (R"≠R) **do**
26:         **if** (size(R") > max-size) **then**
27:             R' ← R"
28:             max-size ← size(R")
29:         **end if**
30:     **end for**
31:     **return** R'
32: **end function**
33: **function** CHOOSE-HONEST(R)
34:     remove(rank-threads,R)
35:     **while** not-empty(rank-threads) **do**
36:         R' ← pop-front(rank-threads)
37:         **if** (has-work(R') **then**
38:             push-back(rank-threads,R)
39:             **return** R'
40:         **end if**
41:     **end while**
42: **end function**
*   **Honest victim selection policy:** This strategy is based on another heuristic to determine the thread with the most difficult work to finish. The heuristic assumes that if a thread $R_1$ has stolen work less recently than a thread $R_2$, then the thread $R_1$ has probably a work which is more difficult than the work of the thread $R_2$. Therefore, the thread $R$ steals the work from the thread victim $R'$ which is the least recent thief. As shown in function *choose-honest*, this policy has a higher computational complexity than the ring and random policies but a smaller computational complexity than the largest policy. In the largest victim selection policy, it is important to compute the amount of work of each pool and computing the size of any pool is a blocking operation for the thread which owns this pool. In the honest victim selection policy, one operation of removing is performed on the rank threads list, and this operation is non-blocking. The operations on the global rank threads list must be protected by locks.

### 5.3.4 Granularity policies

When a thread $R'$ is contacted by a thread $R$, the thread $R$ must determine the amount of work of its thread victim $R'$ to steal. In this chapter, two granularity policies are used.

*   **Steal half policy:** This policy indicates that the thread $R$ steals the second half of the work of the thread $R'$ and leaves the other half for the thread $R'$. In the linked-list based approach, the work of a thread $R'$ is constituted by a set of $N$ nodes. The thread $R$ steals the last $N/2$ nodes and leaves the other nodes for the thread $R'$. Nodes are always stolen from the tail, i.e. from the end which is opposite to the working end of the private deque. While in the interval-based approach, the work of a thread $R'$ is constituted by an interval $[A, B[$. The thread $R$ steals the interval $[(A + B)/2, B[$ and leaves the interval $[A, (A + B)/2[$ for the thread $R'$. Leaving the first half of the interval $[A, B[$ avoids the thread $R'$ to initialize its matrix and vectors.

*   **Steal $T^{th}$ policy:** Theoretically, steal half policy may not be appropriate for certain victim selection policies. Suppose, for instance, four threads where thread 1 has a certain amount of work $W$, and threads 2, 3 and 4 have completed their work. The amount of work $W$ may be the number of nodes or the size of the interval. In a ring selection, the threads 2, 3 and 4 steal work from the threads 1, 2 and 3, respectively. Using the steal half policy and the ring selection, the amounts of work $W/2$, $W/4$, $W/8$ and $W/16$ are allocated to the threads 1, 2, 3 and 4, respectively. Steal $T^{th}$ policy indicates that the thread $R$ leaves $W/T$ of the work to its thread victim $R'$, where $T$ is the number of threads, and steals $(T - 1)W/T$ of the work. In the previous
example, using steal $T^{th}$ policy and the ring selection allocate the amount of works $W/4$, $3W/16$, $9W/64$ and $27W/64$ to the threads 1, 2, 3 and 4, respectively. For this example, steal $T^{th}$ policy gives a better granularity policy than the steal half policy. In our experiments, steal $T^{th}$ policy is tested only for the ring selection. Indeed, steal half policy seems to be theoretically appropriate for the other victim selection policies.

## 5.4 Experiments

In this section, we present an experimental study to evaluate the performance of the different work stealing strategies for our PB&B@CPU approach on a multi-core system. These strategies are compared with the multi-threaded linked-list based PB&B (CPU LL PB&B) approach when solving some hard flowshop scheduling problem instances.

### 5.4.1 Experimental protocol

* **Hardware testbed:** All the experiments are run on a computer composed of 2 8-core Sandy Bridge E5-2650 processors and 32 GB of memory.

* **Software testbed:** The operating system installed is a CentOS 6.5 Linux distribution. The UNIX *time* command is used to measure the elapsed execution time for each flowshop instance and time measures for specific parts of the algorithms are obtained using the *clock_gettime* function with nanosecond precision.

* **Problem instances:** In our experiments, we used only the 10 instances where the number of machines and the number of jobs are equal to 20. Instances where the number of machines is equal to 5 or 10 are easy to solve. For these instances, the used bounding operation gives such good lower bounds that it is possible to solve them in few seconds using a multi-core PB&B. Instances where the number of jobs is equal to 50, 100, 200, or 500, and the number of machines is equal to 20 are very hard to solve.

* **PB&B initialization:** When an instance is solved twice using a PB&B performing a parallel tree exploration, the number of explored permutations is often different between the two resolutions, because the nodes are expanded in a different order. Moreover, the parallel algorithm is likely to explore a (much) larger or smaller number of nodes than the serial version of the algorithm. This may lead to speed-up anomalies, i.e. speed-ups greater than *P* on *P* cores.
To compare the performance of two PB&B algorithms, the explored search space should be exactly the same between the different tests. In order to study the performance of the parallel algorithm in the absence of speed-up anomalies, we choose to always initialize our PB&B by the optimal complete permutation of the instance to be solved. With this initialization it is ensured that both, the sequential and parallel algorithms explore exactly the critical subtree. The critical subtree is composed of all nodes for which the bounding operation gives a lower bound smaller than the cost of the optimal complete permutation. It corresponds to the part of the search space which a PB&B algorithm must explore to prove the optimality of the optimal complete permutation.

It is not unusual to use this technique in the study of parallel PB&B algorithms (see, for instance [KRR88]). Although this situation is not realistic in practice, an initialization to $\infty$ (or $-\infty$ in the case of a maximization problem) is neither realistic in practice. Indeed, the number of branched nodes that do not belong to the critical subtree can be reduced by improving the selection strategy or by initializing the algorithm with a good, suboptimal complete permutation produced by a metaheuristic. However, even with a suboptimal initialization of the algorithm the size of the explored search space still varies between two parallel resolutions. The initialization of the algorithm to the optimal cost ensures that the size and share of the search tree do not depend on the decrease of the best complete permutation found so far and that the number of explored permutations is the same between the two sequential or parallel resolutions. This allows a fair comparison between the two versions.

* **Sequential PB&B:** Table 6.1 shows the number of nodes branched during the resolution of instances *Ta021-Ta030* initialized with the optimal complete permutation. This number represents the total amount of work to be done and ranges from 1.6 for the smallest to 140.8 million nodes for the largest instance. It also shows the elapsed execution times for the resolution of these instances with sequential linked-list PB&B and PB&B@CORE algorithms. It should be noted that the resolution time is not directly proportional to the number of branched nodes. The node processing speed (in nodes/sec), also shown in Table 6.1, varies from one instance to another, as a result of the trees irregularity (the cost of the bounding operation depends on the depth of a permutation). As the flowshop bounding operation consumes between 97% and 99% of the elapsed time, the choice of the data structure, LL or IVM, has no significant impact on the elapsed execution time. In order to speed up the exploration process which lasts on average for more than 6 hours, parallel processing is necessary.
<table><caption>Table 5.3: Number of branched permutations during the resolution of Taillard's instances *Ta021–Ta030* initialized with the optimal cost (in millions of permutations) and sequential execution times using LL-based and IVM-based pool management.</caption><thead><tr><th rowspan="2">Inst.</th><th colspan="2">Branched nodes</th><th colspan="2">Elapsed time (sec)</th><th colspan="2">Nodes/sec</th><th rowspan="2">Rate</th></tr><tr><th>×10<sup>6</sup></th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>CPU LL PB&B</th><th>PB&B@CPU</th></tr></thead><tbody><tr><td>Ta021</td><td>41.4</td><td>24489</td><td>23889</td><td>1691</td><td>1734</td><td>1.03</td></tr><tr><td>Ta022</td><td>22.1</td><td>11758</td><td>11450</td><td>1877</td><td>1927</td><td>1.03</td></tr><tr><td>Ta023</td><td>140.8</td><td>79322</td><td>77298</td><td>1776</td><td>1822</td><td>1.03</td></tr><tr><td>Ta024</td><td>40.1</td><td>19753</td><td>19367</td><td>2028</td><td>2069</td><td>1.02</td></tr><tr><td>Ta025</td><td>41.4</td><td>25332</td><td>24661</td><td>1636</td><td>1680</td><td>1.03</td></tr><tr><td>Ta026</td><td>71.4</td><td>34562</td><td>33722</td><td>2065</td><td>2117</td><td>1.02</td></tr><tr><td>Ta027</td><td>57.1</td><td>28295</td><td>27535</td><td>2018</td><td>2074</td><td>1.03</td></tr><tr><td>Ta028</td><td>8.1</td><td>4569</td><td>4440</td><td>1770</td><td>1822</td><td>1.03</td></tr><tr><td>Ta029</td><td>6.8</td><td>3674</td><td>3583</td><td>1845</td><td>1892</td><td>1.03</td></tr><tr><td>Ta030</td><td>1.6</td><td>898</td><td>873</td><td>1835</td><td>1888</td><td>1.03</td></tr><tr><td>Avg</td><td>43.1</td><td>23265</td><td>22682</td><td>1854</td><td>1902</td><td>1.03</td></tr></tbody></table>

### 5.4.2 Pool management evaluation

<table><caption>Table 5.4: Comparison of CPU LL PB&B and PB&B@CPU in terms of time spent managing the pool of permutations (in seconds).</caption><thead><tr><th rowspan="2">Inst.</th><th colspan="3">Random 1/2</th><th colspan="3">Honest 1/2</th><th colspan="3">Largest 1/2</th><th colspan="3">Ring 1/2</th><th colspan="3">Ring 1/T</th></tr><tr><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>LL</th><th>IVM</th><th>Rate</th></tr></thead><tbody><tr><td>Ta021</td><td>1275</td><td>121</td><td>10.6</td><td>1304</td><td>121</td><td>10.8</td><td>1276</td><td>121</td><td>10.6</td><td>1284</td><td>127</td><td>0.7</td><td>3189</td><td>121</td><td>26.3</td></tr><tr><td>Ta022</td><td>704</td><td>67</td><td>10.5</td><td>706</td><td>67</td><td>10.5</td><td>682</td><td>67</td><td>10.1</td><td>657</td><td>1675</td><td>0.4</td><td>1393</td><td>68</td><td>20.5</td></tr><tr><td>Ta023</td><td>4700</td><td>399</td><td>11.8</td><td>4442</td><td>400</td><td>11.1</td><td>4620</td><td>400</td><td>11.5</td><td>4645</td><td>10117</td><td>0.5</td><td>9009</td><td>401</td><td>22.5</td></tr><tr><td>Ta024</td><td>1248</td><td>111</td><td>11.2</td><td>1256</td><td>111</td><td>11.2</td><td>1290</td><td>111</td><td>11.6</td><td>1247</td><td>1714</td><td>0.4</td><td>2529</td><td>112</td><td>22.6</td></tr><tr><td>Ta025</td><td>1197</td><td>115</td><td>10.5</td><td>1176</td><td>115</td><td>10.3</td><td>1192</td><td>115</td><td>10.4</td><td>1216</td><td>6419</td><td>0.2</td><td>1965</td><td>116</td><td>16.9</td></tr><tr><td>Ta026</td><td>2177</td><td>203</td><td>10.7</td><td>2175</td><td>203</td><td>10.7</td><td>2169</td><td>203</td><td>10.7</td><td>2160</td><td>2431</td><td>0.9</td><td>4498</td><td>204</td><td>22.1</td></tr><tr><td>Ta027</td><td>1894</td><td>156</td><td>12.1</td><td>1934</td><td>156</td><td>12.4</td><td>1953</td><td>156</td><td>12.5</td><td>1917</td><td>5926</td><td>0.3</td><td>3556</td><td>156</td><td>22.7</td></tr><tr><td>Ta028</td><td>262</td><td>24</td><td>11.1</td><td>261</td><td>24</td><td>11.0</td><td>263</td><td>24</td><td>11.1</td><td>261</td><td>1556</td><td>0.2</td><td>600</td><td>25</td><td>24.5</td></tr><tr><td>Ta029</td><td>216</td><td>20</td><td>11.0</td><td>219</td><td>20</td><td>11.2</td><td>220</td><td>19</td><td>11.3</td><td>225</td><td>1415</td><td>0.2</td><td>396</td><td>21</td><td>19.0</td></tr><tr><td>Ta030</td><td>54</td><td>5</td><td>11.1</td><td>55</td><td>5</td><td>10.8</td><td>53</td><td>5</td><td>10.9</td><td>54</td><td>374</td><td>0.1</td><td>118</td><td>6</td><td>20.0</td></tr><tr><td>Avg</td><td>1373</td><td>122</td><td>11.3</td><td>1353</td><td>122</td><td>11.1</td><td>1372</td><td>122</td><td>11.2</td><td>1367</td><td>3487</td><td>0.4</td><td>2725</td><td>123</td><td>22.2</td></tr></tbody></table>

<table><caption>Table 5.5: Comparison of multi-core IVM-based and linked-list based parallel PB&B algorithms in terms of speed-up over their sequential counterparts.</caption><thead><tr><th rowspan="2">Inst.</th><th colspan="3">Random 1/2</th><th colspan="3">Honest 1/2</th><th colspan="3">Largest 1/2</th><th colspan="3">Ring 1/2</th><th colspan="3">Ring 1/T</th></tr><tr><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>LL</th><th>IVM</th><th>Rate</th></tr></thead><tbody><tr><td>Ta021</td><td>14.0</td><td>14.3</td><td>1.02</td><td>14.0</td><td>14.3</td><td>1.02</td><td>14.0</td><td>14.3</td><td>1.02</td><td>14.0</td><td>13.0</td><td>0.93</td><td>11.0</td><td>14.3</td><td>1.30</td></tr><tr><td>Ta022</td><td>13.7</td><td>14.0</td><td>1.02</td><td>13.7</td><td>14.1</td><td>1.03</td><td>13.8</td><td>14.1</td><td>1.03</td><td>13.8</td><td>11.7</td><td>0.85</td><td>11.4</td><td>14.1</td><td>1.21</td></tr><tr><td>Ta023</td><td>13.8</td><td>14.1</td><td>1.02</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.8</td><td>12.1</td><td>0.88</td><td>11.6</td><td>14.1</td><td>1.22</td></tr><tr><td>Ta024</td><td>13.6</td><td>13.9</td><td>1.03</td><td>13.6</td><td>13.9</td><td>1.03</td><td>13.6</td><td>13.9</td><td>1.02</td><td>13.6</td><td>11.7</td><td>0.86</td><td>11.0</td><td>13.9</td><td>1.27</td></tr><tr><td>Ta025</td><td>14.2</td><td>14.4</td><td>1.01</td><td>14.2</td><td>14.4</td><td>1.02</td><td>14.2</td><td>14.4</td><td>1.02</td><td>14.1</td><td>10.9</td><td>0.77</td><td>12.9</td><td>14.4</td><td>1.12</td></tr><tr><td>Ta026</td><td>13.7</td><td>14.0</td><td>1.02</td><td>13.7</td><td>14.0</td><td>1.02</td><td>13.7</td><td>14.0</td><td>1.02</td><td>13.7</td><td>12.9</td><td>0.94</td><td>11.1</td><td>14.0</td><td>1.27</td></tr><tr><td>Ta027</td><td>13.6</td><td>13.7</td><td>1.01</td><td>13.5</td><td>13.8</td><td>1.02</td><td>13.5</td><td>13.7</td><td>1.02</td><td>13.5</td><td>11.0</td><td>0.82</td><td>11.2</td><td>13.7</td><td>1.22</td></tr><tr><td>Ta028</td><td>13.9</td><td>14.2</td><td>1.02</td><td>13.9</td><td>14.1</td><td>1.02</td><td>13.9</td><td>14.2</td><td>1.02</td><td>13.8</td><td>10.0</td><td>0.72</td><td>11.1</td><td>14.0</td><td>1.26</td></tr><tr><td>Ta029</td><td>13.8</td><td>14.1</td><td>1.02</td><td>13.8</td><td>14.1</td><td>1.02</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.7</td><td>9.7</td><td>0.70</td><td>12.0</td><td>14.0</td><td>1.17</td></tr><tr><td>Ta030</td><td>13.6</td><td>13.8</td><td>1.02</td><td>13.6</td><td>13.6</td><td>1.00</td><td>13.7</td><td>13.5</td><td>0.98</td><td>13.7</td><td>9.7</td><td>0.71</td><td>11.2</td><td>13.4</td><td>1.19</td></tr><tr><td>Avg</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.8</td><td>14.0</td><td>1.02</td><td>13.8</td><td>11.3</td><td>0.82</td><td>11.4</td><td>14.0</td><td>1.2</td></tr></tbody></table>

In our experiments five work stealing strategies are evaluated, *Ring-1/T*, *Ring-1/2*, *Random-1/2*, *Largest-1/2* and *Honest-1/2*. As explained in Subsection 5.3.4 the steal-T<sup>th</sup> policy is theoretically not appropriate for the random, largest and honest victim selection
policies. For each of the ten instances *Ta021–Ta030* Table 5.4 shows the time spent managing the pool of permutations. The time spent managing the pool is the cumulated time spent outside of bound computation. For each work stealing strategy, the table shows the obtained results for the LL-based approach, the IVM-based approach, and the ratio between these two approaches. For example when solving instance *Ta021* using the *Random-1/2* strategy the LL-based approach spends 1275 seconds managing the linked-lists of permutations, while the IVM-based approach spends 121 seconds managing the IVM structure, which makes the IVM-based pool management 10.6 times faster than the LL-based pool management.

Except for the *Ring-1/2* strategy the IVM-based approach performs on average at least 11 times better than its LL-based counterpart. While the *Ring-1/2* is clearly inappropriate for the PB&B@CPU, this strategy shows good performances, comparable with *Random*, for the LL-based approach. Conversely, the *steal-T*^th^ granularity policy leads to good results for the IVM-based approach while it is not well-suited to the LL-based approach. Table 5.5 shows the speed-up over the respective serial algorithm for the different work stealing strategies and both approaches. Our PB&B@CPU algorithm reaches a reasonable average speed-up of up to about 14× using 16 threads. For work stealing strategies *Largest*, *Honest*, *Random* the PB&B@CPU solves the ten flowshop instances on average 14.0 times faster than the serial PB&B@CORE, while the LL-based multi-core algorithm yields a slightly inferior average speed-up of 13.8 over its serial counterpart, so the PB&B@CPU shows a slightly better scalability. Although the use of the IVM structure allows a significant reduction of the pool management time (for most strategies) this has only a modest impact on the overall execution time. For all strategies the total execution time for the IVM-based approach is only 4 – 5% faster than its LL-based counterpart (comparing *Ring-1/T* for IVM with *Ring-1/2* for LL). This is due to the predominant cost of the bounding operation.

5.4.3 Interval-based work stealing evaluation

In this subsection we report more detailed experimental results comparing the five pro-
posed work stealing strategies. For each strategy and both approaches, CPU LL PB&B and
PB&B@CPU, Table 5.6 shows the time spent by threads waiting for new work units. The
waiting time is the cumulated time, for all threads, that elapses between the moment an
idle thread starts searching for a victim and the reception of a non-empty interval or set of
nodes, i.e. a successful steal attempt. These measured time intervals completely cover the
time a victim thread is answering the request. The double of this waiting time is therefore
an upper bound for the time spent outside of useful work, i.e. waiting and communication
time.
For example, using the *Random* strategy and solving instance *Ta021*, LL-based threads wait 4.67 seconds for new nodes, which is 15.9 times longer than IVM-based threads, which wait 0.29 seconds for new intervals. For all instances and both approaches, the *Largest* strategy results in the lowest waiting time, between 0.10 and 0.17 seconds. The *Random* and *Honest* strategies result in longer waiting times which, however, remain below < 0.5 seconds for all tested instances. For all these three strategies the average waiting time for the LL approach is at least 10 times higher. When using the ring topology, if one wants to keep the waiting time low, then the granularity policy must be chosen in accordance with the data structure used for the pool management. For instance, using the *Ring-1/2* strategy in combination with the IVM approach threads wait on average for 3360 seconds, which corresponds to ≈ 15% of the average serial execution time. Similarly, using the *Ring-1/T* strategy in a LL-based PB&B algorithm results in an average waiting time of 1176 seconds and over 2 million work stealing operations on average. These important waiting times explain the high pool management times reported for these strategies in Table 5.4. Choosing the unsuitable granularity policy for the ring topology, i.e. Steal-1/2 for IVM and Steal-Tth for LL, limits the average speed-up over the serial algorithm to 11.3, respectively 11.4.

It can be seen from Table 5.6 that the waiting time for the IVM-based approach is much less sensitive to the instance's size and shape than for the LL-based approach. To illustrate this, Figure 5.3 shows the number of performed work stealing operations and the time spent waiting for new work units as a function of the instance size. In order to analyze the behavior of the proposed work stealing strategies according to varying tree sizes, the figure shows the averaged values for groups of instances in which the number of branched nodes is similar: *Ta028-Ta030* are considered as a group of small instances, *Ta021*, *Ta022*, *Ta024* and *Ta025* as medium, *Ta026* and *Ta027* as large, and finally *Ta023* as a very large instance. Figure 5.3 does not show the unsuitable strategies, i.e. *Ring-1/2* using IVM and *Ring-1/T* using LL. In this figure the experimental results for the IVM-based approach using strategies *Random*, *Honest* and *Largest* are barely distinguishable because these values are very small compared to the LL-based approach for large instances.

Figure 5.3 illustrates an interesting feature of the IVM-based work stealing strategies. The resolution of the largest instance *Ta023* requires less work stealing operations and induces less waiting time than the smallest instance *Ta030*, even though its exploration lasts 80 times as long. While the LL-based PB&B@CPU algorithm shows a strong variation in the number of work stealing operations according to the size of the instance being solved, this correlation is apparently very weak or even absent in the IVM-based algorithm.

The largest policy is also the most costly selection policy in terms of time spent per victim selection (including the time spent in synchronization on shared data). Selecting
a victim with the largest policy takes more than twice as much time as with the honest, and almost 10 times as much as with the random selection policy. However, the total time measured for victim selection with the largest policy is 13 ms against 3 ms for the random policy, which corresponds to about 10% of the waiting time for *Largest*, respectively 1% of the waiting time for *random*. As a percentage of the total elapsed time this cost for the victim selection is certainly negligible – however, if the thread-count is increased the cost of the *Largest* and *Honest* strategies may become important, due to contention on the shared data structures accessed by the victim selection function.

<table><caption>Table 5.6: Comparison of multi-core IVM and LL-based parallel PB&B algorithms in terms of time spent waiting for work units; cumulated waiting time for 16 threads (in seconds)</caption><thead><tr><th rowspan="2">Inst.</th><th colspan="3">Random 1/2</th><th colspan="3">Honest 1/2</th><th colspan="3">Largest 1/2</th><th colspan="3">Ring 1/2</th><th colspan="3">Ring 1/T</th></tr><tr><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th><th>CPU LL PB&B</th><th>PB&B@CPU</th><th>Rate</th></tr></thead><tbody><tr><td>Ta021</td><td>4.67</td><td>0.29</td><td>15.9</td><td>8.95</td><td>0.35</td><td>25.5</td><td>2.84</td><td>0.12</td><td>24.4</td><td>3.98</td><td>1610</td><td>&lt;0.01</td><td>1677</td><td>1.11</td><td>&gt;100</td></tr><tr><td>Ta022</td><td>1.89</td><td>0.29</td><td>6.8</td><td>3.69</td><td>0.32</td><td>11.4</td><td>1.07</td><td>0.17</td><td>6.4</td><td>1.49</td><td>1609</td><td>&lt;0.01</td><td>616</td><td>1.20</td><td>&gt;100</td></tr><tr><td>Ta023</td><td>7.87</td><td>0.30</td><td>26.0</td><td>15.78</td><td>0.44</td><td>35.5</td><td>6.09</td><td>0.13</td><td>46.1</td><td>6.33</td><td>9699</td><td>&lt;0.01</td><td>3850</td><td>1.63</td><td>&gt;100</td></tr><tr><td>Ta024</td><td>2.68</td><td>0.25</td><td>10.7</td><td>3.83</td><td>0.28</td><td>13.7</td><td>0.97</td><td>0.16</td><td>6.2</td><td>3.04</td><td>3114</td><td>&lt;0.01</td><td>1076</td><td>0.70</td><td>&gt;100</td></tr><tr><td>Ta025</td><td>3.39</td><td>0.33</td><td>10.4</td><td>3.83</td><td>0.34</td><td>11.3</td><td>1.18</td><td>0.13</td><td>9.1</td><td>3.01</td><td>6288</td><td>&lt;0.01</td><td>670</td><td>1.69</td><td>&gt;100</td></tr><tr><td>Ta026</td><td>3.08</td><td>0.26</td><td>11.9</td><td>5.77</td><td>0.30</td><td>19.0</td><td>1.81</td><td>0.10</td><td>17.3</td><td>3.20</td><td>2234</td><td>&lt;0.01</td><td>1983</td><td>1.34</td><td>&gt;100</td></tr><tr><td>Ta027</td><td>2.68</td><td>0.31</td><td>8.7</td><td>4.25</td><td>0.37</td><td>11.4</td><td>1.76</td><td>0.11</td><td>15.8</td><td>2.33</td><td>5758</td><td>&lt;0.01</td><td>1365</td><td>0.83</td><td>&gt;100</td></tr><tr><td>Ta028</td><td>1.34</td><td>0.22</td><td>6.2</td><td>2.87</td><td>0.40</td><td>7.1</td><td>0.88</td><td>0.12</td><td>7.2</td><td>1.49</td><td>1529</td><td>&lt;0.01</td><td>303</td><td>1.02</td><td>&gt;100</td></tr><tr><td>Ta029</td><td>0.85</td><td>0.25</td><td>3.4</td><td>1.60</td><td>0.27</td><td>5.9</td><td>0.48</td><td>0.13</td><td>3.8</td><td>1.51</td><td>1392</td><td>&lt;0.01</td><td>157</td><td>1.56</td><td>&gt;100</td></tr><tr><td>Ta030</td><td>0.56</td><td>0.25</td><td>2.3</td><td>1.42</td><td>0.39</td><td>3.6</td><td>0.42</td><td>0.13</td><td>3.2</td><td>1.51</td><td>368</td><td>&lt;0.01</td><td>57</td><td>1.25</td><td>45.8</td></tr><tr><td>Aver.</td><td>2.90</td><td>0.27</td><td>10.6</td><td>5.20</td><td>0.35</td><td>14.9</td><td>1.75</td><td>0.13</td><td>13.9</td><td>2.79</td><td>3360</td><td>&lt;0.01</td><td>1176</td><td>1.23</td><td>&gt;100</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Figure 5.3: Number of work stealing operations and waiting time for the CPU LL PB&B and PB&B@CPU using different work stealing strategies.</figcaption></figure>

## 5.5 Conclusions

In this chapter, we propose a new approach to manage the pools of permutations of a parallel PB&B algorithm using work stealing parallel strategies. Our new PB&B approach
aims at reducing the number of work stealing operations and the total execution time of
the algorithm.

To the best of our knowledge, all other parallel PB&B algorithms published in the literature are based on conventional structures, like linked-lists for the depth first strategy, to store the permutations of each thread-private pool. Our new approach is based on the factorial number system, which is a mixed radix numeral system adapted to numbering permutations, and the use of an IVM structure to manage the permutations of the private pools.

A work stealing strategy can be defined by its victim selection policy and its granularity policy. The victim selection policy indicates the thread victim that a thread can steal, and the granularity policy determines the amount of work stolen by the thread. In this chapter, five work stealing strategies are defined for the CPU LL PB&B and the PB&B@CPU algorithm. These strategies are *Ring 1/T*, *Ring 1/2*, *Random 1/2*, *Largest 1/2*, and *Honest 1/2*.

The obtained results, on 10 flowshop instances, demonstrate that our new IVM-based approach outperforms the conventional linked-list based approach in terms of total elapsed CPU time and performed work stealing operations. Also, compared to the linked-list based approach, the IVM-based approach allows to reduce the idle time of workers. We have shown that, especially for large instances, the IVM-based work stealing performs much less work stealing operations than the LL-based approach.

Moreover, the IVM-based approach, contrary to the linked-list based approach, has constant and controllable memory requirements. The memory usage pattern is a huge advantage for many-core architectures, such as GPUs, where the performance depends a lot on memory usage and memory access patterns. For 20-job instances, more than 100 IVM structures fit easily into 48KB of shared memory available per multiprocessor on most devices. Although the implementation of a deque, a stack, etc. is theoretically feasible in CUDA, this type of data structure would necessarily reside in global memory, which is much slower than shared memory. The introduction of Unified Memory since CUDA 6 is likely to make the use of linked-lists on the GPU more practical. However, fitting into the fastest part of the hierarchical GPU memory, the IVM is clearly the better adapted data structure for the GPU architecture. The next chapter introduces the PB&B@GPU approach based on the use of IVM and interval concepts on GPU accelerators.
CHAPTER 6

# GPU-centric permutation
Branch-and-Bound

## Contents

<table><tr><td>6.1 Introduction</td><td>71</td></tr><tr><td>6.2 IVM-based Data structures</td><td>73</td></tr><tr><td>6.2.1 Memory requirements</td><td>73</td></tr><tr><td>6.2.2 Memory constraint: Number of IVMs</td><td>74</td></tr><tr><td>6.3 Operations</td><td>75</td></tr><tr><td>6.3.1 Bound mapping operation</td><td>75</td></tr><tr><td>6.3.2 IVM mapping operation</td><td>79</td></tr><tr><td>6.3.3 Work stealing operation on GPU</td><td>80</td></tr><tr><td>6.4 Experiments</td><td>83</td></tr><tr><td>6.4.1 Experimental protocol</td><td>83</td></tr><tr><td>6.4.2 Evaluation of bound mapping schemes</td><td>84</td></tr><tr><td>6.4.3 Evaluation of IVM mapping schemes</td><td>86</td></tr><tr><td>6.4.4 Evaluation of work stealing schemes</td><td>88</td></tr><tr><td>6.4.5 Comparison of PB&B@GPU and GPU LL PB&B</td><td>91</td></tr><tr><td>6.5 Conclusions</td><td>92</td></tr></table>

## 6.1 Introduction

In this chapter, the PB&B@CPU algorithm is revisited for GPUs. In contrast to the other existing GPU-accelerated B&B algorithms, this chapter deals with a GPU-centric implementation of PB&B, meaning that it completely bypasses the CPU. To the best of our knowledge, our PB&B@GPU algorithm is the first to perform all PB&B operations on the GPU device. The key idea to achieving this is to use the IVM structure for pool management, which is much better suited for GPUs than dynamic LL data structure. In addition, the use of intervals, as work units, makes work stealing strategies on GPUs more efficient. As shown
in Figure 6.1, the PB&B@GPU approach is very similar to the PB&B@CPU approach. How-
ever, PB&B@GPU differs from PB&B@CPU in several aspects:

* The number of PB&B threads of PB&B@GPU is much higher than the number of PB&B threads of PB&B@CPU.

<figure><img src="image_18.png" /><figcaption>Figure 6.1: Overview of the PB&B@GPU approach and its operations.</figcaption></figure>

* The work stealing operation of PB&B@GPU is different from the work stealing operation of PB&B@CPU; Subsection 6.3.3 describes the work stealing operation of the PB&B@GPU.

* In addition to the other operations, the PB&B@GPU approach also uses a supplementary operation called a bound mapping and IVM mapping operations; Subsection 6.3.1 and Subsection 6.3.2 present respectively these bound mapping and IVM mapping operations.

* Unlike the bounding operation of PB&B@CPU, where this operation calculates the set of bounds of all the sub-permutations of a partial permutation, the bounding operation of PB&B@GPU calculates the bound of a single sub-permutation. Therefore, for each sub-permutation, it is necessary to launch up to *N* GPU threads, where *N* is the size of the permutation problem. Each of these threads calculates the bound
of a single sub-permutation. In a PB&B@GPU approach with *T* PB&B threads, it is possible to run up to *N* × *T* threads for a single bounding operation.

* In PB&B@CPU, the operations of the different PB&B threads are run asynchronously. In other words, in this approach, a PB&B thread may run the branching operation at the same time when another PB&B thread runs the bounding operation. On the other hand, in the PB&B@GPU approach, the operations of different PB&B threads are executed synchronously. It means that, if a PB&B thread of PB&B@GPU executes a certain operation, then all the other threads perform the same operation. For example, if a PB&B thread runs the bounding operation, then all the other threads run this operation. The synchronization between all the operations is represented in Figure 6.1 by the synchronization points, drawn in green.

This chapter is divided into three main sections. In addition to Section 6.2 and Section 6.3, which respectively describe the data structures and the operations of PB&B@GPU, Section 6.4 presents the experiments performed to validate PB&B@GPU approach.

# 6.2 IVM-based Data structures

This section focuses on the data structure of our PB&B@GPU based on IVM. The memory requirements of the IVM structure are very advantageous for a GPU implementation of the PB&B algorithm. The required amount of memory and possible data placements in the hierarchical device memory are discussed in Subsection 6.2.1. The amount of used memory depends on the number of IVM structures used by the algorithm. This number also has a direct impact on the degree of parallelism which is analyzed in Subsection 6.2.2.

## 6.2.1 Memory requirements

Compared to a conventional LL-based approach, the IVM structure allows to reduce the CPU time and memory required for the storage and management of the pool of permutations [MLMT14b]. Contrary to LL, the IVM structure is well adapted to the GPU memory model. Instead of using a variable length queue that requires dynamic memory allocations and tends to be scattered in memory, the IVM structures are constant in size and need only one allocation of contiguous memory. For a problem instance with *N* elements, the storage of the matrix *M* requires *N*² bytes of memory (for *N* < 127, using 1-byte integers). Moreover 3*N* bytes are needed to store the position, end and branching vectors, 1 byte to store the depth, and *N* bytes to store permutations before calling the bounding operation. In total, the IVM structure requires a constant amount of 1 + 4*N* + *N*² bytes of memory, i.e.
481 bytes per IVM for a 20-element instance. It is possible to store only the upper triangular part of M, requiring $1 + 4N + \frac{N(N+1)}{2}$ bytes per IVM, i.e. 291 bytes when $N = 20$. For $N = 20$ it is therefore possible to fit $\ge 100$ IVM structures into 48KB of shared memory. In our PB&B@GPU, only the upper triangular part of M is stored.

From a programming perspective the IVM structures are easy to handle. The components of all IVMs are merged into single one dimensional arrays. For instance, solving a $N$-element instance using $T$ IVM structures, the matrices are stored in a one dimensional array matrices of size $T \times \frac{N(N+1)}{2}$ allocated in global device memory. The element $M(i, j)$ of the $k^{th}$ IVM is accessed by matrices[indexM(i, j, k)], where indexM is a wrapper function defined as in Equation (6.1) if M is stored as a square and as in Equation (6.2) if the upper triangular part of M is stored.

$$indexM(i, j, k) = k \times N \times N + i \times N + j \qquad (6.1)$$

$$indexM(i, j, k) = k \times \frac{N(N+1)}{2} + i \times N - \frac{i(i-1)}{2} + j \qquad (6.2)$$

For the flowshop problem, the data needed for the computation of the lower bounds is mostly read-only and require 34.5KB of memory. These data are stored in the constant memory space, residing in global device memory but accessed through a cache on each Streaming Multiprocessor (SMX). Some of the data structures used for the bounding may be loaded to shared memory during the computation of the lower bounds. Concerning the use of shared memory for those data structures, this chapter follows the recommendations made in [Cha13], where this difficult choice is examined.

### 6.2.2 Memory constraint: Number of IVMs

Often, the bounding operation is by far the most time-consuming part of a PB&B algorithm. In the case of flowshop it amounts for about 97 - 99% [MCB14] of the total execution time for a sequential PB&B. It is therefore crucial for the performance of our PB&B@GPU that the parallel bounding operation makes the best use of the GPU resources. The choice of the number of PB&B processes (=IVMs) to use is therefore guided by its impact on the performance of the bounding kernel. On the one hand, if too few IVMs participate in the exploration process, the bounding kernel under-utilizes the GPU. On the other hand, if too many IVMs are used, then the number of generated permutations per iteration exceeds the maximum occupancy of the device and the computation of bounds is partially serialized. The number of permutations generated per IVM per iteration is variable and unpredictable in shape and size. However, the workload for the bounding kernel can be roughly estimated. For flowshop instances of 20 jobs, the bulk of permutations is situated at
depth 10, leading to approximately 20 bound evaluations per IVM and iteration. Supposing that the number of empty IVMs is low thanks to dynamic load balancing, and given the approximative number of 20,000 concurrent threads at full occupancy, the number *T* of used IVMs should be around *T* = 1000.

## 6.3 Operations

This section introduces the two mapping operations and the work stealing operation.

### 6.3.1 Bound mapping operation

The shape of the tree explored by a PB&B algorithm is highly irregular and unpredictable, resulting in an irregular workload, irregular control flow and irregular memory access pattern. If not addressed properly, these irregularities may cause a low occupancy of the device, serialized execution of instructions and poor bandwidth usage due to uncoalesced memory accesses. Both, the application's memory access pattern and the divergent behavior of threads depend strongly on the chosen mapping of threads onto the data. When a GPU application runs, each streaming multiprocessor is assigned one or more thread block(s) to execute. Those threads are partitioned into groups of 32 threads<fnref n="1" />, called *warps*, which are scheduled for execution. CUDA's single instruction multiple-thread (SIMT) execution model assumes that a warp executes one common instruction at a time. Consequently, full efficiency is realized when all 32 threads of a warp agree on their execution path. However, if threads of a warp diverge due to a data dependent conditional branch, the warp serially executes each branch path taken. Threads that are not on that path are disabled, and when all paths complete, the threads converge back to the same execution path. This phenomenon is called *thread divergence* and often causes serious performance degradations. In a very similar way, if the threads in a warp agree on the location of a requested piece of data, it may be fetched in single cycle, otherwise serialization of the data accesses occurs. In this subsection the focus is put on reducing thread divergence and increasing warp execution efficiency by making judicious mapping choices.

#### Static mapping method

The most straightforward approach probably consists in mapping each thread onto a child permutation directly from its `threadId`. This naive and static approach is shown in Algorithm 6. For instance, launching $2 \times N \times \#IVM$ threads (line 2), the first $N \times \#IVM$

<footnote n="1">¹We assume using the GK110 model</footnote>
threads place unscheduled elements in the beginning, the second $N \times \#IVM$ threads in the end. Regardless of the IVM's state or current depth in the tree, $2 \times N$ threads are reserved for each IVM. Each thread is assigned an IVM to work on and an element to schedule, like shown in lines 4-6 of Algorithm 6. The approach of Algorithm 6 has several disadvantages. The *if*-conditionals in line mask many of the launched threads, precisely $2 \times k$ threads per father permutation of depth $k$, plus $2N$ threads per empty IVM. Moreover, different lanes in the same warp work on different IVMs, thus thread divergence occurs due to different values of $L_1$ and $L_2$. If $T \times N$ is a multiple of warp-size, then the *if-else* conditional (lines 10 and 14) does not cause any thread divergence.

<caption>**Algorithm 6 Kernel STATIC-BOUND-MAPPING**</caption>

IN: Fathers (father, $L_1, L_2$)
OUT: Lower bounds begin, lower bounds end, sums of lower bounds

1: **kernel** *NAIVE-BOUND*
2:     <<<< 2 \* \#elements \* \#IVM threads>>>
3:     thId←block Id.x\*blockDim.x + threadId.x
4:     **if** (state[ivm] == not-empty) **then**
5:         **if** ($L_1$[ivm] < element < $L_2$[ivm]) **then**
6:             **if** (dir == 0) **then**
7:                 swap( schedule, $L_1$[ivm] + 1, element )
8:                 LB-begin[ivm][element]←computeLB( schedule )
9:                 sum-begin[ivm] += LB-begin[ivm][element]
10:             **else if** (dir == 1) **then**
11:                 swap( schedule, $L_2$[ivm] - 1, element )
12:                 LB-end[ivm][element]←computeLB( schedule )
13:                 sum-end[ivm] += LB-end[ivm][element]
14:             **end if**
15:         **end if**
16:     **end if**
17: **end kernel**

Remapping method (First step)

Algorithm 7 describes how to build the maps *ivm-map* and *element-map* sequentially. However, sequential execution of this procedure on the device has prohibitive cost, exceeding 25% of the total execution time. The remapping should therefore be built in parallel. The parallelization of the outer *for*-loop (Algorithm 7, line 3) is not straightforward, because
it is unknown at which location the data for each IVM is to be written to. Computing the prefix-sum of a vector containing the number of elements to be scheduled per IVM allows its parallelization.

<caption>**Algorithm 7 Build mapping (serial)**</caption>

1: **kernel** SERIAL PREPARE BOUND
2:   running-index ← 0
3:   **for** (ivm = 0 → T)) **do**
4:     **if** (state[ivm] = not-empty) **then**
5:       **for** (element = L₁[ivm] + 1 → L₂[ivm]) **do**
6:         ivm-map[running-index] ← ivm
7:         element-map[running-index] ← element
8:         running-index++
9:     **end for**
10:   **end if**
11:   **end for**
12:   todo ← running-index
13: **end kernel**

<caption>**Algorithm 8 STEP1-REMAPPED-BOUND-MAPPING**</caption>

1: **for all** (non-empty ivm) **do**
2:   todo-per[ivm] ← (L₂[ivm] - L₁[ivm] - 1)
3:   **end for**
4:   Aux ← parallel-prefix-sumerti
5:   prepare-bound<<< #IVM × #ELEMENTS >>>
6:   **kernel** [KERNEL] PREPARE-BOUND
7:   thId ← blockIdx.x*blockDim.x + threadIdx.x
8:   ivm ← thId / N
9:   thPos ← thId % N
10:  **if** (thPos < todo-per[ivm]) **then**
11:   ivm-map[Aux[ivm] + thPos] ← ivm
12:   element-map[Aux[ivm] + thPos] ← L₁[ivm] + 1 + thPos
13:   **end if**
14:   todo ← Aux[#IVM] + todo-per[#IVM]
15: **end kernel**
The operation *prefix-sum* is defined as

$$ \mathit{prefix-sum}: [a_0 \ a_1 \ a_2 \ \dots \ a_n] \mapsto [0 \ a_0 \ (a_0+a_1) \ (a_0+a_1+a_2) \ \dots \ \sum_{i=0}^{n-1} a_i]. $$

Efficient parallel CUDA implementations for this operation have been proposed in the literature [HSO07]. It is also available in the CUDA Thrust library. However, for relatively small vectors it may be preferable to reimplement the operation, in order to avoid casting the input data to a `thrust::device_ptr`.

A first building step consists in filling an array `todo-per-IVM` with $L_2 - L_1 - 1$ for each IVM. The element $R$ of `prefix-sum'(todo-per-IVM)` indicates at which position of `ivm-map` and `element-map` the data of an IVM $R$ starts to be written. The complete parallelized building of the mapping is shown in Algorithm 8. The building of the mapping ranges over several kernels.

## Remapping method (Second step)

**Algorithm 9 Kernel STEP2-REMAPPED-BOUND-MAPPING**

in:fathers (father, $L_1, L_2$), ivm-map, element-map

out:lower bounds begin, lower bounds end, sums of lower bounds

1: **kernel REMAPPED-BOUND**

2:     `<<< 2 * todo threads >>>`

3:     thId ← blockIdx.x *atisDim.x + threadIdx.x

4:     dir ← thId mod 2

5:     ivm ← ivm-map[thId/2]

6:     element ← element-map[thId/2]

7:     schedule ← fathers[ivm]

8:     toSwap ← (1-dir) * (L₁[ivm] + 1) + dir * (L₂[ivm] - 1)

9:     swap( schedule, toSwap, element )

10:    LB[dir][ivm][element] ← computeLB( schedule )

11:    sum[dir][ivm] += LB[dir][ivm][element]

12: **end kernel**

The goal of the remapping procedure which prepares the bounding is to build two maps `ivm-map` and `element-map` which contain, for `todo` threads, the information which IVM to work on and which element to swap. Using an even/odd pattern these maps provide sufficient information for both groups of threads. After building these maps, the bounding kernel (as shown in Algorithm 9) is called with $2 \times todo$ threads, where:
* threads 0 and 1 work on IVM *ivm-map[0]*, swapping element *element-map[0]* respectively to begin/end,

* threads 2 and 3 work on IVM *ivm-map[1]*, swapping element *element-map[1]* respectively to begin/end,

* ...

* threads $2 \times \text{todo} - 2$ and $2 \times \text{todo} - 1$ work on IVM *ivm-map[\text{todo}-1]*,...

The remapped bounding kernel is launched at each iteration with a kernel configuration of $(2 * \text{todo}/\text{blockDim}) + 1$ blocks (simplified in Algorithm 9) which is adapted to the workload. The proposed approach is known as *stream compaction* in the literature. It reduces the number of idle lanes per warp as well as the number of threads launched per kernel invocation. However, any thread divergence resulting from the begin-end distinction should also be avoided, as this involves a serialization of the costly *computeLB* procedure. To achieve this, the bodies of the *if-else* conditional (Algorithm 6) can be merged into a single one (Algorithm 9, lines 8–11). Two different arguments of the same type, occurring on the right-hand side of a statement can often be re-factored into a single one, like in Algorithm 9, line 8. The different arrays on the left hand side are merged into larger ones. This allows to merge the statements of lines 8,9 and 12,13 of Algorithm 6 into single statements (Algorithm 9, lines 10,11). The separation of data within these merged arrays is assured by indexing with the variable *dir*, which evaluates differently for even/odd threads.

### 6.3.2 IVM mapping operation

The IVM management kernels (i.e work stealing, branching, selection and pruning) require a single thread per IVM. The naive approach consists in launching *T* threads and mapping thread *k* on IVM *k*, for *k* = 0, 1, ..., *T* − 1 (see Algorithm 10). Given the high number of conditional instructions in the IVM management kernels it is very unlikely that all 32 threads in a warp follow the same execution path if this mapping is used. Indeed, in these kernels control flow divergence results from different IVM states, different numbers of scheduled elements at both ends of the active permutation and from the search for the next node which requires an unknown number of iterations. An alternative mapping, shown in Algorithm 11, can solve this issue. An entire warp is assigned to each IVM, so all threads belonging to the same warp follow the same execution path. This strategy goes in the opposite direction of the stream compaction approach proposed for the bounding kernel. As only one thread per IVM is needed, all lanes in a warp except this first are
masked. Thus, the kernels are launched with 32× as many threads as necessary (i.e. 32× T). Using this mapping, the overhead induced by thread divergence completely disappears (although technically, the disabled threads are diverging at line 5 of Algorithm 11). The drawback is obviously the launching of 31T idle threads. However, in Subsection 6.2.2 we argued that T should be chosen around T = 1000, which is small compared to #SM×(max. threads per SM). This, and the fact that the control flow irregularity is very high, justifies the approach of using 1 warp per IVM. Moreover, using only 4-8 IVM structures per block allows to store them into shared memory without limiting the theoretical occupancy of the device. The loading of data from global to shared memory can be done very efficiently, using the additional threads which are not used for computation.

<caption>**Algorithm 10 Kernel THREAD-IVM-MAPPING**</caption>

1: **kernel**<<<#IVM threads>>>
2: ivm←blockIdx.x*blockDim.x + threadIdx.x
3: do-something-with[ivm]

<caption>**Algorithm 11 Kernel WARP-IVM-MAPPING**</caption>

1: **kernel**<<<warblock×#IVM threads>>>
2: thId←blockIdx.x*blockDim.x + threadIdx.x
3: ivm←thId/32
4: thPos←thId%32
5: **if** (thPos == 0) **then**
6:     do-something-with[ivm]
7: **end if**

### 6.3.3 Work stealing operation on GPU

Work stealing is well adapted for irregular applications. Like threads of a PB&B@CPU, the IVM structures must share their work units. In a multi-core CPU environment, a thread that runs out of work becomes a thief that attempts to steal a portion of work from a victim thread which is selected according to a victim selection strategy. The same principle can be applied to the PB&B@GPU. The proposed load balancing strategy is conceptually different in the sense that an IVM-based PB&B process does not necessarily correspond to any particular thread but only to a segment of data. Secondly, compared to PB&B@CPU work stealing operation, the work stealing PB&B@GPU operations between IVMs are lock-free and performed synchronously. The kernel implements the 1D-Ring work stealing strategy. Algorithm 12 shows the pseudo-code of this procedure. Although
designed for PB&B@CPU, the 1D-Ring strategy suits the synchronous execution mode of the GPU. The *T* IVM structures are numbered *R* = 0, 1, ..., *T* − 1 and are arranged as an oriented ring, i.e. such that IVM 0 is IVM (*T* − 1)’s successor. Each empty IVM *R* tries to steal work from its predecessor (*R* − 1)%*T*. This operation can be performed in parallel, as the mapping of empty IVMs onto their respective victims is one-to-one. If the selected victim has a non-empty interval, then all but 1/*T*th of its interval is stolen. The function computeNewPos (line 6) receives the victim’s interval [*A*, *B*] as input and returns a point *M* = (1 − 1/*T*)*A* + 1/*T* *B*. The division of intervals can be performed directly on the factoradic numbers without explicitly converting them to decimals. The IVM which got stolen continues the exploration of the remaining interval [*A*, *M*[, while the stealing IVM needs to initialize its matrix at the new position vector *M* before starting the exploration of [*M*, *B*]. Its state variable is therefore set to *initializing* (line 10). Each IVM cycles through three distinct states, from *exploring* to *empty* to *initializing* and back to *exploring*. An IVM can be in one of these three states at any given stage of the algorithm. Depending on the state of an IVM, different actions are performed during an iteration. In this kernel one thread per IVM is required. More parallelism can hardly be exposed. However, more threads can eventually be used to assign vectors in one parallel operation (lines 7 − 9).

<caption>**Algorithm 12 Kernel WORK-STEALING**</caption>

1: **kernel** WORK-STEALING
2:     thId ← block Id.x *block Dim.x + thread Id.x
3:     ivm ← map(thId)
4:     victim ← (ivm-1)%T
5:     **if** (state[ivm]=empty .and. state[victim]=exploring) **then**
6:         new-pos ← computeNewPos(pos[victim], end[victim])
7:         pos[ivm] ← new-pos
8:         end[ivm] ← end[victim]
9:         end[victim] ← new-pos
10:        state[ivm] ← initializing
11:    **end if**
12: **end kernel**

The topology used in the work stealing strategy, already described, is a unidirectional 1-dimensional ring (1D-Ring). The maximal distance between two IVMs in the 1D-Ring is *T*. Work units propagate through the ring as they are passed downstream from exploring to empty IVMs. As most of the explored PB&B nodes are actually contained in a relatively small interval, the workload tends to be concentrated in some part of the ring. Thus,
workers situated far away from the source are only kept busy if the overall workload is large enough. With an increasing number *T* of IVM structures it becomes more likely that no work is dripping down to some of the workers. A topology that reduces the maximum distance between two workers should therefore improve the scaling with *T*.

The 1D-Ring can be easily generalized to a 2D-Ring, or torus, topology. Instead of using a single ring, IVMs are arranged in *R* rings of ring-size *C* = *T*/*R*. In a first step each empty IVM attempts to steal from its left neighbor within the same ring. A second step connects the rings between each other: each empty IVM selects the IVM with the corresponding ID in the preceding ring (with ring *R* − 1 being connected to ring 0). The roles played by both directions are symmetric. Ideally, the number *R* is therefore such that *R* = √*T*, which is only possible if *T* is square. In that case the 2D-Ring reduces the maximum distance between two IVMs to 2√*T*. If *C* ≠ *R*, then the diameter of the 2D-Ring is (*C* + *R*).

The 2D-Ring topology is implemented by two subsequent calls of kernel work stealing (Algorithm 12), where only line 4 of the algorithm needs to be modified. In particular, line 4 of Algorithm 12 is replaced by the following.

$$
\text{In Step 1 IVM } i \text{ selects } victim(i) = \begin{cases} i - 1, & \text{if } i \bmod C \neq 0 \\ i + (C - 1), & \text{otherwise} \end{cases}
$$

$$
\text{In Step 2 IVM } i \text{ selects } victim(i) = \begin{cases} i - C, & \text{if } i > (C - 1) \\ (R - 1)C + i, & \text{otherwise} \end{cases}
$$

Figure 6.2 illustrates the 2D-Ring topology in the form of a 2D grid. A torus, used in the 2D-Ring work stealing strategy is obtained by connecting the upper with the lower and the leftmost with the rightmost cells. Similarly, the topology can be extended to a hypercube, which is used for instance in [SKK$^{+}$11] for unbalanced tree search.

<table><caption>Figure 6.2: Illustration of 2D-Ring topology for *T* IVMs using *R* rings of ring size *C*.</caption>
  <tr>
    <td>0</td>
    <td>1</td>
    <td>2</td>
    <td>...</td>
    <td>(C-1)</td>
  </tr>
  <tr>
    <td>C</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
  </tr>
  <tr>
    <td>2C</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
  </tr>
  <tr>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
  </tr>
  <tr>
    <td>(R-1)C</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>RxC-1<br>=T-1</td>
  </tr>
</table>
## 6.4 Experiments

In this section the performance of PB&B@GPU is analyzed for different mapping choices, a varying number of IVMs and different work stealing strategies. In Subsection 6.4.2, the two strategies for the bound mapping operation are evaluated. Subsection 6.4.3 compares the two mapping strategies for the IVM mapping operation. The algorithm's scalability and load balancing issues are examined in Subsection 6.4.4. Finally, our PB&B@GPU algorithm is compared to the GPU LL PB&B algorithm presented in [CM13].

### 6.4.1 Experimental protocol

* **Hardware testbed:** All the experiments are run on a computer equipped with an NVIDIA Tesla K20m GPU based on the GK110 architecture. The device is composed of 2496 CUDA cores (clock speed 705MHz). Its maximum power consumption is 225W. The CPU is a 8-core Sandy Bridge E5-2650 processor.

* **Software tools:** Version 6.5.14 of the CUDA Toolkit is used. The operating system installed is a CentOS 6.5 Linux distribution. For the evaluation of the elapsed execution time the UNIX *time* command is used. The duration of each CUDA kernel and profiling of the kernels is done with the *nvprof* command line profiler. In order to reduce the profiling time, sample data was collected every 100 iterations of the algorithm.

* **Configuration choices:** The chosen size for the thread blocks is 128. The configurable size of the device's shared memory/L1 cache is set to 48/16KB for kernels except bound, where the opposite configuration 16/48KB is used. For the comparison of the mapping strategies the number of used IVM structures is set to *T* = 768, according to preliminary experiments. The best mapping found in Subsection 6.4.2 and Subsection 6.4.3 is used to determine an optimal value for *T* and the better work stealing strategy in Subsection 6.4.4.

* **Problem instances:** In our experiments, the validation is performed using the 10 instances where the number of machines and the number of jobs are equal to 20 which belong to the group 20x20. When an instance is solved twice using a PB&B performing a parallel tree exploration, the number of explored permutations is often different between the two resolutions, because the order of exploration varies. To compare the performance of two PB&B algorithms, the number of explored permutations should be exactly the same between the different tests. Therefore, we choose to always initialize our PB&B by the optimal complete permutation of the instance
to be solved. This initialization ensures that the tree-shape does not depend on the decrease of the best complete permutation found so far and that the number of explored permutations is the same between the two resolutions. Table 6.1 shows the number of branched nodes during the resolution of instances *Ta021-Ta030* initialized with the optimal complete permutation. This number represents the total amount of work to be done and ranges from 1.6 (for the smallest) to 140.8 million nodes (for the largest) instance.

<table><caption>Table 6.1: Number of branched permutations during the resolution of Tail-lard's instances *Ta021-Ta030* initialized with the optimal cost (in millions of nodes).</caption><thead><tr><th>Instance</th><th>Ta021</th><th>Ta022</th><th>Ta023</th><th>Ta024</th><th>Ta025</th><th>Ta026</th><th>Ta027</th><th>Ta028</th><th>Ta029</th><th>Ta030</th><th>Average</th></tr></thead><tbody><tr><td>#Nodes (in millions)</td><td>41.4</td><td>22.1</td><td>140.8</td><td>40.1</td><td>41.4</td><td>71.4</td><td>57.1</td><td>8.1</td><td>6.8</td><td>1.6</td><td>43.1</td></tr></tbody></table>

### 6.4.2 Evaluation of bound mapping schemes

<figure><img src="image_19.png" /><figcaption>Figure 6.3: Execution time for instances *Ta021-Ta030* for thread data mappings *static* (S) and *remap* (M) for the kernel bound.</figcaption></figure>

In this subsection the two mapping schemes for the bounding kernel, presented in Subsection 6.3.1, are compared to each other in terms of elapsed execution time of the algorithm. The first, using the remapping shown in Algorithm 9 is referred to as *remap*, the second, using the static mapping of Algorithm 6, as *static*. Figure 6.3 shows the total elapsed time for solving instances *Ta021-Ta030*. For both mappings and for each instance it shows the portion of time spent in the kernel bounding, in the IVM management kernels (i.e work stealing, branching, selection and pruning) as well as in the remapping operation (for