[book_gpu.git] / BookGPU / Chapters / chapter8 / ch8.tex

\chapterauthor{Imen Chakroun}{Universit\'e Lille 1 CNRS/LIFL, INRIA Lille Nord Europe, Cit\'e scientifique - 59655, Villeneuve d'Ascq cedex, France\\}
\chapterauthor{Nouredine Melab}{Universit\'e Lille 1 CNRS/LIFL, INRIA Lille Nord Europe, Cit\'e scientifique - 59655, Villeneuve d'Ascq cedex, France\\}

\chapter{GPU-accelerated Tree-based Exact Optimization Methods}

\section{Introduction}
\label{ch8:introduction}

In practice, a wide range of problems can be modeled as NP-hard combinatorial optimization problems (COPs). Those problems consist in choosing the best combination out of a large finite set of possible combinations and are known to be large in size and difficult to solve to optimality. One of the most popular methods for solving exactly a COP (finding a solution having the optimal cost), is the Branch-and-Bound (B\&B) algorithm. This algorithm is based on an implicit enumeration of all the feasible solutions of the tackled problem. Enumerating the solutions of a problem consists in building a dynamically generated search tree whose nodes are subsets of solutions of the considered problem. The construction of such tree and its exploration is performed using four operators: branching, bounding, selection and pruning. Due to the exponentially increasing number of potential solutions, the B\&B algorithm explores only promising nodes of the search tree using an estimated optimal solution called ``lower bound'' of the associated sub-problem. 

\vspace{0.3cm}

Although this bounding mechanism allows to considerably reduce the exploration time, often only small or moderately-sized instances of COPs can be practically solved. For this reason, over the last decades, parallel computing has been revealed as an attractive way to deal with larger instances of COPs. However, while many contributions have been proposed for parallel B\&B methods using Massively Parallel Processors \cite{ch8:Allen_1997}, Networks or Clusters of Workstations \cite{ch8:Quinn_1990} and SMP machines \cite{ch8:Casadoa_2008}, very few contributions have been proposed for redesigning B\&B algorithms on Graphical Processing Units (GPUs) \cite{ch8:Carneiro_2011}. For years, the use of GPU accelerators was limited to graphics and video applications. Driven by the demand for high-definition 3D graphics on personal computers, GPUs have evolved into a highly parallel, multi-threaded and many-core environment. Their utilization has recently been extended to other application domains such as scientific computing \cite{ch8:Kurzak_2010}. 

\vspace{0.3cm}

In this work, we rethink the design and implementation of irregular tree-based algorithms such as B\&B algorithm on top of GPUs. During the execution of the B\&B algorithm, the number of new generated nodes and the number of not yet explored but promising nodes are variable and depend on the level of the tree being explored and on the best solution found so far. Therefore, due to such unstructured and unpredictable nature of its search tree, designing efficient B\&B on top of GPUs is not straightforward. We investigate two different approaches for designing GPU-based B\&B starting from the parallel models for B\&B identified in \cite{ch8:MelabHDR_2005}. The first one is based on the ``parallel tree exploration'' paradigm. This approach consists in exploring in parallel different sub-spaces of the tree. The second approach is based on the ``parallel evaluation of bounds'' approach. The two approaches have been applied to the permutation \index{Flowshop Scheduling Problem} (FSP)(see Section~\ref{ch8:BB-FSP}) which is an NP-hard combinatorial optimization problem. The lower bound function used in this work for FSP is the one proposed in~\cite{ch8:Johnson_1954} for two machines and generalized in~\cite{ch8:Lenstra_1978} to more than two machines.

\vspace{0.3cm}

When rethinking those two parallel model for GPU's architectures, our main focus was on the lower bound function. Indeed, preliminary experiments we carried out on some Taillard's problem instances \cite{ch8:Taillard_1993} show that computing the lower bounds takes on average between 98\% and 99\% of the total execution time of the B\&B. The GPU-based lower bound's implementation raises mainly two challenges. On the one hand, having in mind that the execution model of GPUs is SIMD, irregular computations (containing loops and conditional instructions) contained in the lower bound function may lead to a very challenging issue: the thread or branch divergence. This problem drops down the performance and arises when threads of a same warp (the smallest executable unit of parallelism on the GPU) execute different data-dependent instructions. On the other hand, the lower bound computation usually uses large in size and frequently accessed data structures. Since GPU is a many-core co-processor device that provides a hierarchy of memories having different sizes and access latencies, the placement and sharing of these data sets become challenging.

\vspace{0.3cm}

The scope of this chapter is to design parallel B\&B algorithms on GPU accelerators to allow highly efficient solving of permutation-based COPs. To do so, our contributions consist in: (1) rethinking two approaches for parallel B\&B on top of GPUs, discussing the performances of each and identifying which best suits the GPU accelerators. (2) proposing a new approach for thread/branch divergence reduction through a thorough analysis of the different loops and conditional instructions of the bounding function. (2) defining an optimal mapping of the data structures of the bounding function on the hierarchy of memories provided in the GPU device through a careful analysis of both the data structures (size and access frequency) and the GPU memories (size and access latency). 

\vspace{0.3cm}

The chapter is organized in seven main sections. Section \ref{ch8:BB} presents the B\&B algorithm. Section \ref{ch8:Parallel-BB} introduces the different models used to parallelize B\&B algorithms. Section \ref{ch8:BB-FSP} briefly describes the Flowshop Scheduling permutation Problem. In Section~\ref{ch8:approach1}, we describe the GPU-accelerated B\&B based on the parallel tree exploration. In Section~\ref{ch8:approach2},  details about the second approach, the GPU-accelerated B\&B based on the parallel evaluation of lower bounds, are given. In Section \ref{ch8:ThreadDivergence}, the thread divergence issue related to the location of nodes in the B\&B tree and to the control flow instructions within the bounding operator is described. In Section \ref{ch8:DataAccessOpt}, the memory access optimization challenge is addressed and an overview of the GPU memory hierarchy and the used memory access pattern is given. In Section~\ref{ch8:Experiments}, we report experimental results showing the performances of each of two studied approaches compared to a sequential CPU-based execution of the B\&B and demonstrating the efficiency of the proposed optimizations.

\section{\index{Branch-and-Bound} algorithm}
\label{ch8:BB}

Branch-and-bound algorithms are by far the most widely used methods for exactly solving large scale NP-hard combinatorial optimization problems. Indeed, they allow to find the optimal solution of a problem with proof of optimality. 

\vspace{0.3cm}

The basic idea of the B\&B algorithm consists in implicitly enumerating all the solutions of the original problem by only examining a subset of feasible solutions and eliminating the others when they are not likely to lead to a feasible or an optimal solution. Enumerating the solutions of a problem consists in building a dynamically generated search tree whose nodes are subsets of solutions of the considered problem. The construction of such tree and its exploration are performed using four operators: branching, bounding, selection and pruning.

\vspace{0.3cm}

The algorithm proceeds in several iterations during which the best solution found so far is progressively improved. During the exploration process, the search space is described by a pool of unexplored nodes and the best solution found so far. The generated and not yet examined (pending) nodes are kept into a list initialized with the original problem. At each iteration of the algorithm, the following steps are performed:

\begin{itemize}
 \item The {\it selection operator} chooses one node to process among the pending nodes according to a defined strategy. If the selection is based on the depth of the sub-problem in the B\&B tree, we speak about a depth-first exploration strategy. A selection based on the breadth of the sub-problem is called a breadth-first exploration. A best-first selection strategy could also be used. It is based on the presumed capacity of the node to yield good solutions.
 \item The {\it branching operator} subdivides a solution space into two or more disjointed sub-spaces to be investigated in a subsequent iteration.
 \item The {\it bounding operator} computes a bound value of the optimal solution of each generated sub-problem.
 \item Each sub-problem having a greater bound than the upper-bound, i.e. the cost of the best solution found so far, is eliminated using the {\it pruning operator}.
\end{itemize}

Algorithm \ref{ch8:algoBB} gives the general template of the Branch-and-Bound method.

\begin{algorithm}[H]

\SetAlgoLined  

\vspace{0.2cm}

Create the initial problem; \\
Inset the initial problem into the tree; \\
Set the Upper\_Bound to $\propto$;  \\
Set the Best\_Solution to $\emptyset$; \\

\While{ not\_empty\_tree() }
{
    \vspace{0.2cm}      

    Sub\_Problem = Take\_sub\_problem();  

    \If{ Is\_leaf ( Sub\_Problem ) } 
    {
            Upper\_Bound = Cost\_Of( Sub\_Problem );\\
            Best\_Solution = Sub\_Problem;
    }
    \Else
    {
       Lower\_Bound = compute\_lower\_bound(Sub\_Problem);
      
       \If{ Lower\_Bound $\leq$ Upper\_Bound } 
       {        
          Branch(Sub\_Problem); \\
          Insert child sub problems into the tree; 
       }        
       \Else
       {
          Prune (Sub\_Problem); 
       }
    }
}

\caption{General template of the Branch-and-Bound Algorithm.}
\label{ch8:algoBB}
\end{algorithm}

\section{Parallel Branch-and-Bound algorithms}
\label{ch8:Parallel-BB}

Thanks to the bounding operator, B\&B allows to significantly reduce the computing time needed to explore the whole solution space. However, finding an optimal solution for large instances remains unpractical using a sequential B\&B. Therefore, parallel processing of these algorithms has been widely studied in the literature. In \cite{ch8:MelabHDR_2005}, a taxonomy of the various existing parallel paradigm used to parallelize the B\&B algorithm is presented.

\vspace{0.2cm}

This taxonomy based on the classification proposed in \cite{ch8:Gendron_1994} identified several models to accelerate the B\&B search. The first model we consider in this chapter is called ``parallel tree exploration model'' and belongs to the ``Tree-based'' strategies that aim to build and explore the B\&B tree in parallel. The second model called ``parallel evaluation of bounds model'' (evaluation of bounds in parallel) belong to the parallelization approach called ``Node-based''. This strategy aims to accelerate the execution of a particular operation at the node level.

\vspace{0.2cm}

\subsection{The parallel tree exploration model}
\label{ch8:para_tree}

Tree-based strategies consist in building and/or exploring the solution tree in parallel by performing operations on several sub-problems simultaneously. This coarse-grained type of parallelism affects the general structure of the B\&B algorithm and makes it highly irregular.\\

The \index{parallel tree exploration} model, illustrated in Figure \ref{ch8:parallel_tree}, consists in visiting in parallel different paths of the same tree. The search tree is explored in parallel by performing the branching, selection, bounding and elimination operators on several sub-problems simultaneously.\\

\begin{figure}
  \begin{center}
\includegraphics[scale=0.5]{Chapters/chapter8/figures/parallel_exploration.eps}%

\caption{Illustration of the parallel tree exploration model}
\label{ch8:parallel_tree}
  \end{center}
\end{figure}

\subsection{The parallel evaluation of bounds model}
\label{ch8:Node_parallel}

Node-based strategies introduce parallelism when performing the operations on a single problem. For instance, they consist in executing the bounding operation in parallel for each sub-problem to accelerate the execution. This type of parallelism has no influence on the general structure of the B\&B algorithm and is particular to the problem being solved.\\

The \index{parallel evaluation of bounds} model, as shown in Figure \ref{ch8:bounds_parallel}, allows the parallelization of the bounding of sub-problems generated by the branching operator. This model is used in the case where the bounding operator is performed several times after the branching operator. The model does not change the order and the number of explored sub-problems in the parallel B\&B algorithm compared to the sequential B\&B.

\begin{figure}
  \begin{center}
\includegraphics[scale=0.5]{Chapters/chapter8/figures/parallel_bounding.eps}%

\caption{Illustration of the parallel evaluation of bounds model}
\label{ch8:bounds_parallel}
  \end{center}
\end{figure}

\section{The Flowshop Scheduling Problem}
\label{ch8:BB-FSP}

\subsection{Definition of the Flowshop Scheduling Problem} 
\label{ch8:LB-FSP}

As a case study for our GPU-based Branch-and-Bound, we considered the NP-hard and well-known problem in the scheduling theory: the "Permutation Flow-shop Scheduling Problem" (FSP). 
In this work, the mono-objective case is considered. The FSP aims to find the optimal schedule of n jobs on m machines so that the overall completion time of all jobs, called {\it makespan}, is minimized. 

\vspace{0.3cm}

Let us suppose the set of jobs is represented by J = {$j_1$, $j_2$, . . . $j_n$} and the set of machines is represented by M = {$m_1$,$m_2$, . . .$m_m$} organized 
in the line. Each job $j_i$ is a sequence of operations ji = { $oi_1$, $oi_2$, . . . $oi_m$ } where oim is the duration required for the job ji on the machine m. 
A feasible solution of the flowshop permutation should satisfy these constraints:

\begin{itemize}
 \item A machine can not start processing a job if all the machines, which are located upstream, did not finish their treatment. Thus, the operation $oi_j$ cannot be processed by the machine $m_j$ if it is not completed on $m_j$ - 1. 
\item An operation can not be interrupted, and the machines are critical resources, because a machine processes one job at a time.
\item The sequence of jobs should be the same on every machine, e.g. if j3 is treated in position 2 on the first machine, j3 is also executed in position 2 on all machines.
\end{itemize}

Figure~\ref{flow-shop} illustrates a solution of a flow-shop problem instance defined by 6 jobs and 3 machines. 

\begin{figure}[h!]
\centering
\includegraphics[height=1.7cm,width=6.8cm]{Chapters/chapter8/figures/FlowShop.eps}
\caption{Flow-shop problem instance with 3 jobs and 6 machines.}
\label{flow-shop}
\end{figure}

\vspace{0.3cm}

\subsection{\index{Lower Bound} for the Flowshop Scheduling Problem}
\label{ch8:LB-FSP}

The lower bounding technique provides a lower bound (LB) for each sub-problem generated by the branching operator. The more the bound is accurate, the more it allows to eliminate not promising nodes from the search tree. Therefore, the efficiency of a B\&B algorithm depends strongly on the quality of its lower bound function. In this chapter, we use the lower bound proposed by Lenstra {\it et al.}~\cite{ch8:Lenstra_1978} for FSP, based on the Johnson's algorithm~\cite{ch8:Johnson_1954}.

\vspace{0.2cm}

The Johnson's algorithm allows to solve optimally FSP with two machines ($m=2$) using the following transitive rule $\preceq$:

$$J_i \preceq J_j \Leftrightarrow \min(p_{i,1}\ ;\ p_{j,2}) \leq
\min(p_{i,2}\ ;\ p_{j,1})$$

We recall that $p_{k,l}$ designates the processing time of the job $J_k$ on the machine $M_l$. From the above rule, it follows the Johnson's theorem: \\

\textbf{Jonhson's theorem} \emph{Given $P$ an FSP with $m=2$, if $J_i\preceq J_j$ there exists an optimal schedule for $P$ in which the job $J_i$ precedes the job $J_j$.}\\

According to Johnson's theorem, FSP with $m=2$ is solved with a time complexity of $O(n.log n)$. The optimal solution is obtained by first sorting in increasing order the jobs having a
processing time shorter on the first machine than on the second one~; Second, sorting in decreasing order the jobs having a shorter processing time on the second machine. 

\vspace{0.2cm}

In~\cite{ch8:JRJackson_1956} and~\cite{ch8:LGMitten_1959}, the Johnson's rule has been extended by Jackson and Mitten with lags which allowed further Lenstra {\it et al.} to propose a lower bound for FSP with $m \geq 3$. A lag~$l_j$ designates the minimum duration between the starting time of the job $J_j$ on the second machine and its finishing time on the first machine. Jackson and Mitten demonstrated that the optimal solution for FSP with $m=2$ can be obtained using the following transitive rule $\preceq$:

$$J_i \preceq J_j \Leftrightarrow \min(p_{i,1}+l_i\ ;\ l_j+p_{j,2})
\leq \min(l_i+p_{i,2}\ ;\ p_{j,1}+l_j)$$

Based on this rule, Lenstra {\it et al.}~\cite{ch8:Lenstra_1978} have proposed the following lower bound for a sub-problem associated to a partial schedule where a set {\Large $\jmath$} of jobs have to be scheduled on $m$ machines. $P_{Ja}^*(\jmath,M_k,M_l)$ represents the Jackson-Mitten optimal solution for the sub-problem that consists in scheduling the set {\Large $\jmath$} of jobs on the two machines $M_k$ and~$M_l$. The term  $r_{i,k} = \sum_{l<k} p_{i,l}$ designates the starting time of the job $J_i$ on the machine $M_k$. The other term $q_{j,l} = \sum_{k>l} p_{j,k}$ refers to the latency between the finishing time of $J_j$ on $M_l$ and the finishing time of the schedule.

$$LB(\jmath)=\max\limits_{1 \leq k < l \leq m}\{P_{Ja}^*(\jmath,M_k,M_l)+\min\limits_{(i,j)\in \jmath^2, i \neq
j}(r_{i,k}+q_{j,l}) \}$$

According to this $LB$ expression, the lower bound for the scheduling of a subset {\Large $\jmath$} of jobs is calculated by applying the Johnson's rule with lags considering all the couples
$(k,l)$ for $1 \leq k,l \leq m$ and $k<l$. As illustrated in Figure~\ref{LagKLExample}, the lag $l_j$ of a job $J_j$ for a couple $(k,l)$ of machines is the sum of the processing times of the job on
all the machines between~$k$~and~$l$.

\begin{figure}
  \begin{center}
\includegraphics[width=8cm]{Chapters/chapter8/figures/johnson_with_lags.eps}%
\caption{The lag $l_j$ of a job $J_j$ for a couple $(k,l)$ of machines is the sum of the processing times of the job on all the machines between~$k$~and~$l$.}
\label{LagKLExample}
  \end{center}
\end{figure}

\section{GPU-accelerated B\&B based on the parallel tree exploration (GPU-PTE-BB)}
\label{ch8:approach1} 

The first approach we investigate for designing B\&B on GPUs consists in exploring in parallel the generated search tree. The idea is to divide the global search space into disjoint sub-spaces that are explored in parallel by the GPU threads. As explained in Section \ref{ch8:BB}, during the execution of a B\&B, the search space is described by a list of unexplored (pending) nodes and the best solution found so far. In the considered GPU-based scheme, a set of parent nodes is selected from this list according to their depth: deepest pending nodes are the first selected. The selected pool of nodes is off loaded to the GPU where each thread builds its own local search tree by applying the {\it branching}, {\it bounding} and {\it pruning} operators to the assigned node.

\vspace{0.2cm}

\begin{figure}[h!]
\centering
\includegraphics[height=8cm, width=8.1cm]{Chapters/chapter8/figures/Diagram1.eps}
\caption{The overall architecture of the parallel tree exploration-based GPU-accelerated Branch-and-Bound algorithm.} 
\label{tree_approach}
\end{figure}

\vspace{0.2cm}

According to the CUDA threading model, each thread has a unique identifier used to determine its assigned role, assigns specific input and output positions and selects work to perform. Therefore, each node (problem) from the pending list is mapped to a thread to ensure that each sub-space of the solution space is evaluated concurrently and is disjoint from others. Figure \ref{tree_approach} illustrates the scheme of the parallel tree exploration-based GPU-accelerated B\&B.

\section{GPU-accelerated B\&B based on the parallel evaluation of bounds (GPU-PEB-BB) }
\label{ch8:approach2} 

In the GPU-accelerated B\&B based on the parallel evaluation of bounds, illustrated in Figure~\ref{ch8:approach}, the generation of the sub-problems (elimination, selection and branching operations) to be solved is performed on CPU and the evaluation of their lower bounds (bounding operation) is executed on the GPU device. The pool of sub-problems generated on CPU is off-loaded to the GPU device to be evaluated by a pool of threads partitioned into blocks. Each thread applies the lower bound function to one sub-problem. Once the evaluation is completed, the lower bound values corresponding to the different sub-problems is returned to the CPU to be used by the elimination operator to decide either to be pruned or to be decomposed. The process is iterated until the exploration is completed and the optimal solution is found. 

\vspace{0.2cm}

\begin{figure}[h!]
  \begin{center}
\includegraphics[scale=0.3]{Chapters/chapter8/figures/approach.eps}%
\caption{The overall architecture of the GPU-accelerated Branch-and-Bound algorithm based on the parallel evaluation of bounds.}
\label{ch8:approach}
  \end{center}
\end{figure}

\vspace{0.2cm}

In both considered approaches, GPU-PEB-BB and GPU-PTE-BB, the GPU-based lower bound's implementation raises mainly two challenges. The first one is related to the ``single instruction multiple data'' (SIMD) model of the GPU and to the implementation of the LB. Indeed, although typically every GPU thread will run the identical lower bound function, the body of the lower bound can contains conditions on thread identifiers and data. This implies that different instructions are executed in some threads. In SIMD architectures like GPUs this behavior leads to the thread or branch divergence issue. This problem arises when threads of a same warp execute different data-dependent instructions. It might causes serious performance declining since computation occurs in parallel only when the same instructions are being performed. The second challenge consists in adjusting the pattern of accesses to the GPU device memory. Good placement of data over the different memory hierarchy grants programmers to further improve the throughput of many high-performance CUDA applications. For B\&B applied to FSP, threads of the same block perform concurrent accesses to the six data structures of the problem when they execute the lower bound function. These data structures have different sizes and access frequencies and should be wisely placed on the different memories of the GPUs that also have different sizes and latencies. 

\vspace{0.2cm}

In the following, we present how we dealt with the thread/branch divergence issue and maps the different data structures on the memory hierarchy of the GPU device taking into account the characteristics of the data structures and those of the different GPU memories.

\vspace{-0.4cm}

\section{\index{Thread divergence}}
\label{ch8:ThreadDivergence}

\subsection{The thread divergence issue}

During the execution of an application on GPU, to each GPU multiprocessor is assigned one or more thread block(s) to execute. Those threads are partitioned into warps that get scheduled for execution. For each  instruction of the flow, the multiprocessor selects a warp that is ready to be run. A warp executes one common instruction at a time, so full efficiency is realized when all threads of a warp agree on their execution path. In this chapter, the G80 model, in which a warp is a pool of 32 threads, is used. If threads of a warp diverge via a data-dependent conditional branch, the warp serially executes each branch path taken. Threads that are not on the taken path are disabled, and when all paths complete, the threads converge back to the same execution path. This phenomenon is called thread/branch divergence and often causes serious performance degradations. Branch divergence occurs only within a warp; different warps execute independently regardless of whether they are executing common or disjointed code paths.

\vspace{0.2cm}

This section discusses thread divergence issue encountered when computing the bounds by GPU. The thread divergence occurs for two main reasons, namely the locations of nodes in the search tree and the control flow instructions within the bounding operator. 

\vspace{0.3cm}

\textbf{Divergence related to the location of nodes}

\vspace{0.3cm}

This divergence is related to the positions of the nodes in the B\&B search tree. Below is given an example from the source code of the used LB showing that the execution flow depends on the position of the node in the search tree. In the following piece of code, three methods are used {\it is\_leaf()}, {\it makespan()} and {\it lower\_bound()}. {\it is\_leaf()} tests if the node {\it \_node} is a leaf or an internal node. If {\it \_node} is a leaf, {\it makespan()} computes the cost of its makespan. Otherwise, {\it \_node} is an internal node and {\it lower\_bound()} computes the value of its lower bound. 

\begin{verbatim}
        if (_node.is_leaf())
           return _node.makespan();
        else
           return _node.lower_bound();
\end{verbatim}

\vspace{0.3cm}

\textbf{Divergence related to the control flow instructions}

\vspace{0.2cm}

Control flow refers to the order in which the instructions, statements or function calls are executed in a program. This flow is determined by instructions such as {\it if-then-else}, {\it for}, {\it while-do}, {\it switch-case}, etc. There are a dozen of such instructions in the implementation of our bounding operator. The source code examples given below show two scenarios in which this kind of instructions is used.

\begin{itemize}
\item Example 1:\\ \vspace{-0.4cm}
\begin{verbatim}
   if( pool[thread_idx].begin != 0 ) 
      time = TimeMachines[1] ;
   else
      time = TimeArrival[1] ;
\end{verbatim}

\item Example 2:\\ \vspace{-0.4cm} 
\begin{verbatim}
   for(int k = 0 ; k < pool[thread_idx].begin; k++)
      jobTime = jobEnd[k] ;  
\end{verbatim}

\end{itemize}

In these two examples, {\it thread\_idx} is the index associated to the current thread. Let suppose that the code of Example 1 is executed by $32$ threads, {\it pool[thread\_idx].begin} is equal to $0$ for the first thread, and {\it pool[thread\_idx].begin} is not equal to $0$ for the other $31$ threads. When the first thread executes the statement {\it ``time = TimeArrival[1];''}, 
all the other $31$ threads remain idle. Therefore, the GPU cores on which these $31$ threads are executed remain idle and can not be used during the execution of the statement {\it ``time = TimeArrival[1];``}. 

\vspace{0.2cm}

The same scenario occurs during the execution of Example 2. Let us suppose that the instruction is executed by $32$ threads, {\it pool[thread\_idx].begin} is equal to $100$ for the first thread, and {\it pool[thread\_idx].begin} is equal to $0$ for the other $31$ threads. When the first thread executes the loop $for$, all the other $31$ threads remain idle. 

\vspace{0.2cm}

Existing techniques for handling branch divergence either demand hardware support \cite{ch8:Fung} or require host-GPU interaction \cite{ch8:Zhang}, which incurs overhead. Some other works such as \cite{ch8:Han} intervene at the code level. They expose a branch distribution method that aims to reduce the divergent portion of a branch by factoring out structurally similar code from the branch paths. In our work, we have also opted for software-based optimizations like \cite{ch8:Han}. In fact, we figure out how to literally rewrite the branching instructions into basic ones in order to make thread execution paths uniform. We also demonstrate that we could ameliorate performances only by judiciously reordering data being assigned to each thread.

\subsection{Mechanisms for reducing branch divergence}

\vspace{0.3cm}

 \textbf{Thread-data reordering}

\vspace{0.2cm}

At each iteration of our GPU-accelerated B\&B approach, several thousands of sub-problems are sent to the GPU. The GPU groups the received sub-problems into several warps according to their reception order. The first 32 sub-problems belong to the first warp, the following 32 sub-problems belong to the second warp, etc. Therefore, thread-data reordering technique sorts sub-problems before sending them to the GPU. These sub-problems are sorted according to their position in the B\&B tree. This sort of sub-problems allows to have warps containing more homogeneous sub-problems, and reduces the number of thread divergences. 

\vspace{0.2cm}

 \textbf{Branch refactoring}

\vspace{0.2cm}

As quoted above, thread or branch divergence occurs when the kernel includes conditional instructions and loops that make the threads performing different control flows leading to their serial execution. In this chapter, we investigate the branch refactoring approach to deal with thread divergence. Branch refactoring consists in rewriting the conditional instructions so that threads of the same warp execute an uniform code avoiding their divergence. To do that, two major ``if" scenarios are studied and some optimizations are proposed accordingly. These two scenarios correspond to the conditional instructions contained in the $LB$ kernel code. In the first scenario, the conditional expression is a comparison of the content of a variable to 0. For instance, the following example extracted from the pseudo-code of the lower bound $LB$ illustrates such scenario.

\vspace{0.3cm}

\begin{tabular}{l}
\\
\small
\textsf{ if ( pool[thread\_idx].limit1 $\neq$ 0 ) tmp = MM[1];  }\\
\small
\textsf{ else  tmp = RM[1] ; }\\  \\
\end{tabular}

\vspace{0.2cm}

The refactoring idea is to replace the conditional expression by two functions namely $f$ and $g$ as shown in Equation~\ref{ch8:Eq1}.

\vspace{0.2cm}

The behavior of $f$ and $g$ fits the cosine trigonometric function. These functions return values between $0$ and $1$. An integer variable is used to store the result of the cosine function. Its value is $0$ or $1$ since it is rounded to $0$ if it is not equal to~$1$. In order to increase the performance the CUDA runtime math operations are used: $sinf(x)$, $expf(x)$ and so forth. Those functions are mapped directly to the hardware level~\cite{ch8:cuda}. They are faster but provide lower accuracy which does not matter in our case because the results are rounded to $int$. 

\begin{equation}
\begin{array}{lllllllll}
\small
    &\multicolumn{8}{l}{ if (x \neq 0) ~ a = b[1]; ~~~~~~ if (x \neq 0) ~ a = b[1] + 0 \times c[1];} \\
    & \multicolumn{2}{l}{} & ~~~~~~\Rightarrow & \multicolumn{2}{l}{} \\
    &\multicolumn{2}{l}{\emph{else}}    $a = c[1];$ &     &\multicolumn{2}{l}{\emph{else}} $a = 0 \times b[1] + c[1];$ \\\\
    & \multicolumn{6}{l}{\Rightarrow a = f(x) \times b[1] + g(x) \times c[1];}\\ \\
    &\multicolumn{6}{l}{\emph{where:}}\\
    &&\multicolumn{6}{l}{ f(x)=\left\{
                    \begin{array}{lll}
                        f(x) = 0    & if &x = 0\\
                        1           &else\\
                    \end{array}
                \right.}\\
     &\multicolumn{6}{l}{\emph{and}}\\
    &&\multicolumn{6}{l}{g(x)=\left\{
                    \begin{array}{lll}
                        g(x) = 1    & if &x = 0\\
                        0           &else & \\
                    \end{array}
                  \right.}
\end{array}
\label{ch8:Eq1}
\end{equation}

\vspace{0.3cm}

The throughput of $sinf(x)$, $cosf(x)$, $expf(x)$ is one operation per clock cycle~\cite{ch8:cuda}. The refactoring result for the ``if" pseudo-code given above is the following:

\vspace{0.3cm}

\begin{tabular}{l}
\\
\small
\textsf{int coeff = \_\_cosf (pool[thread\_idx].limit1);}\\
\small
\textsf{tmp = (1 - coeff) $\times$ MM[1] +  coeff $\times$ RM[1];}\\ \\
\end{tabular}

\vspace{0.3cm}

The second "if" scenario considered in our study compares two values between themselves as shown in Equation~\ref{ch8:Eq2}.

\vspace{0.2cm}

\begin{equation}
\begin{array}{lllllllll}
\small
    &\multicolumn{8}{l}{if (x > y)  a = b[1];~~~ \Rightarrow if (x - y \geq 1)   a = b[1];}\\\\
    \Rightarrow &\multicolumn{6}{l}{if (x - y - 1 \geq 0)}& a = b[1];~~~~~ (x, y) \in N\\\\
    \Rightarrow &\multicolumn{4}{l}{a = f(x, y) \times b[1] + g(x,y) \times a;}\\\\
    &\multicolumn{6}{l}{\emph{where:}}\\
    &&\multicolumn{6}{l}{ f(x,y)=\left\{
                    \begin{array}{lll}
                        1    & if &x - y - 1 \geq 0\\
                        0    &if&x - y - 1 < 0\\
                    \end{array}
                \right.}\\
    &\multicolumn{6}{l}{\emph{and}}\\
    &&\multicolumn{6}{l}{g(x,y)=\left\{
                    \begin{array}{llll}
                        0    & if &x - y - 1 \geq 0\\
                        1    &if&x - y - 1 < 0 & \\
                    \end{array}
                \right.}\\\\
\end{array}
\label{ch8:Eq2}
\end{equation}

\vspace{0.3cm}

For instance, the following example extracted from the pseudo-code of the lower bound $LB$ illustrates such scenario.

\vspace{0.3cm}

\footnotesize
\begin{tabular}{ll}
\\
\multicolumn{2}{l}{\textsf{if(RM[1]] $>$ MIN )}\{}  \textsf{Best\_idx = Current\_idx;}  \textsf{\}}\\\\
\end{tabular}
\normalsize

\vspace{0.3cm}

The same transformations as those applied for the first scenario are applied here using the exponential function. Recall that the exponential is a positive function which is equal to $1$ when applied to $0$. Thus, if $x$ is greater than $y$ then $expf(x-y-1)$ returns a value between $0$ and $1$. If the result is rounded to an integer value $0$ will be obtained. Now, if $x$ is less than $y$ then $expf(x-y-1)$ returns a value greater than $1$ and since the minimum between $1$ and the exponential is get, the returned result would be $1$. Such behavior satisfies exactly our prerequisites. The above ``if" instruction pseudo-code is now equivalent to:

\vspace{0.3cm}

\small
\begin{tabular}{l}
\\
\textsf{int coeff = min(1, \_\_expf(RM[1] - MIN - 1)); }\\
\textsf{Best\_idx = coeff $\times$ Current\_idx + ( 1 - coeff ) $\times$ Best\_idx ;}\\
\end{tabular}
\normalsize

\section{Memory access optimization}
\label{ch8:DataAccessOpt}

\index{Memory access optimizations} are by far the most studied area for improving GPU-based application performances. Indeed, adjusting the pattern of accesses to the GPU device memory grants programmers to further improve the throughput of many high-performance CUDA applications. The goal of memory access optimizations is generally to use as much fast memory and as little slow-access memory as possible. This section discusses how best to set up data LB items on the various kinds of memory on the device.   

\vspace{0.2cm} 

CUDA enabled devices use several memory spaces, which have different characteristics in term of sizes and access latencies. These memory spaces include global memory, local memory , shared memory, texture memory , and registers. Devices of compute capability 2.0 have also an L1 $/$ L2 cache hierarchy that is used to cache local and global memory accesses.

\begin{itemize}
\item At the thread-level, each thread has its own allocated registers and a private local memory. CUDA uses this local memory for thread-private variables that do not fit in the threads registers, as well as for stack frames and register spilling. \item At the thread block-level, each thread block has a shared memory visible to all its associated threads. \item At the grid-level, all threads have access to the same global memory. Texture and constant cached memories are two other memories accessible by all threads. 
\end{itemize}

The data access optimization challenge is to find the best mapping of the data structures of the application at hand (different sizes and access frequencies) and the GPU hierarchy of memories (different sizes and access latencies). For instance, of these different memory spaces, global memory is the most plentiful but the one with the highest access latency. On the contrary, shared memory is smaller in size but has much higher bandwidth and lower latency than the global memory.

\subsection{Complexity analysis of the memory usage of the Lower Bound }
\label{ch8:MemComplex}

In this section, the characteristics of the data structures used by the lower bound function are studied in terms of sizes and access frequencies. For an efficient implementation of the LB, six data structures are required: the  matrix $PTM$ of the processing times of the jobs, the matrix of lags $LM$, the Johnson's matrix $JM$, the matrix $RM$ of the earliest starting times of jobs, the matrix $QM$ of their lowest latency times and the matrix $MM$ containing the couples of machines. The complexities of the different data structures are summarized in Table~\ref{ch8:tabMemComplex} where the columns represent respectively the name of the data structure, its size and the number of times it is accessed.

\vspace{0.2cm} 

In the $LB$ expression, the computation of the term $P_{Ja}^*(\jmath,M_k,M_l)$ requires the calculation of the lag of each remaining job to be scheduled on the couple $(M_k,M_l)$ of machines using its processing times on these machines (Johnson's rule with lags). Such computation is repeated for each couple $(M_k,M_l)$ of machines with $1 \leq k,l \leq m$ and $k<l$. To avoid the repetitive computation of the lags, they are computed once at the beginning of the algorithm and stored in the matrix $LM$. The dimension of $LM$ is $n \times \frac{m\times (m-1)}{2}$, where $n$ and $m$ are respectively the number of jobs to be scheduled and $m$ the number of machines. $LM$ is accessed $n' \times \frac{m \times (m-1)}{2}$ times, $n'$ being the number of remaining jobs to be scheduled in the sub-problem for which the lower bound is being calculated. The processing times of all the jobs on all the machines are stored in the matrix $PTM$. This matrix has a dimension of $n \times m$ and is accessed $n' \times m \times (m-1)$ times.

\vspace{0.2cm} 

In addition, in order to avoid relaunching the Johnson's algorithm for each couple of machines and each subset of jobs, the Johnson's algorithm is computed once to find the optimal solutions on the couples of machines. These optimal solutions are then stored in the Johnson's matrix $JM$. This matrix has the same dimension as $LM$ and is accessed $n \times \frac{m \times (m-1)}{2}$ times during the computation of the lower bound. Finally, the $MM$ matrix that contains all the couples of machines has a dimension and access frequency of $m \times (m-1)$.  

\vspace{0.2cm} 

To reduce the computation time cost of the term $\min\limits_{(i,j)\in \jmath^2, i \neq j}(r_{i,k}+q_{j,l})$ in the $LB$ expression, two matrices are defined, namely $RM$ and $QM$. They are used to store respectively the lowest starting and latency times of all the jobs on each machine. Their dimension is $m$ and are accessed $ m \times (m-1)$ times and $ \frac{m \times (m-1)}{2}$ times respectively.

\begin{table}
  \centering
\begin{tabular}{|c|c|c|}
\hline
  \textbf{Matrix} & \textbf{Size} & \textbf{Number of accesses} \\
 \hline
 \hline
   PTM &  $n \times m$ & $n' \times m \times (m-1)$ \\
 \hline
   LM & $n \times \frac{m \times (m-1)}{2}$ & $n' \times \frac{m \times (m-1)}{2}$ \\
 \hline
   JM & $n \times \frac{m \times (m-1)}{2}$ & $n \times \frac{m \times (m-1)}{2}$ \\
 \hline
   RM &  $m$ & $m \times (m-1)$ \\
 \hline
   QM &  $m$ & $\frac{m \times (m-1)}{2}$ \\
 \hline
   MM &  $m \times (m-1)$ & $m \times (m-1)$ \\
 \hline
\end{tabular}
\vspace{0.5cm}
 \caption{The different data structures of the $LB$ algorithm and their associated complexities in memory size and numbers of accesses. The parameters $n$, $m$ and $n'$ designate respectively the total number of jobs, the total number of machines and the number of remaining jobs to be scheduled for the sub-problems the lower bound is being computed.}
\label{ch8:tabMemComplex}
\end{table}

\subsection{Data placement pattern of the Lower Bound on GPU}
\label{ch8:MemComplex}

This section discusses how best to map the six data structures identified above on the various kinds of memories of the GPU device.

\vspace{0.2cm} 

The focus is put on the shared memory which is a key enabler for many high-performance CUDA applications. Indeed, because it is on-chip, shared memory has much higher bandwidth and lower latency than local and global memory. However, for large problem instances (large $n$ and $m$) the data structures especially JM and LM (see Table \ref{ch8:tabMemSizes}), do not fit in the shared memory for some GPU configurations. 

\vspace{0.2cm} 

In order to achieve further performances, we also take care of adequately use the global memory by judiciously configuring the L1 cache which greatly enables improving performance over direct access to global memory. Indeed, the GPU device we are using in our experiments is based on the NVIDIA Fermi architecture which introduced two new hierarchies of memories (L1 $/$ L2 cache)
compared to older architectures. 

\begin{table*}
  \centering
  \footnotesize
  \begin{tabular}{|r|r|r|r|r|r|}
    \hline
Prob. instance & JM & LM & PTM & RM, QM & MM \\
    \hline
    \hline
$200 \times 20$ & 38.000 (38KB) & 38.000 (76KB) & 4.000 (4KB) & 20 (0.04KB) & 380 (0.76KB) \\
    \hline
$100 \times 20$ & 19.000 (19KB) & 19.000 (38KB) & 2.000 (2KB) & 20 (0.04KB) & 380 (0.76KB) \\
    \hline
$50 \times 20$ & 9.500 (9.5KB) & 9.500 (19KB) & 1.000 (1KB) & 20 (0.04KB) & 380 (0.76KB) \\
    \hline
$20 \times 20$ & 3.800 (3.8KB) & 3.800 (7.6KB) & 400 (0.4KB) & 20 (0.04KB) & 380 (0.76KB) \\
   \hline
  \end{tabular}
\vspace{0.5cm}
\caption{The sizes of each data structure for the different experimented problem instances. The sizes are given in number of elements and in bytes (between brackets).}
\label{ch8:tabMemSizes}
\end{table*}

\vspace{0.2cm} 

Taking into consideration the sizes of each data structure presented in Table \ref{ch8:tabMemSizes}, our challenge is to find which data structure has to be mapped on which memory and in some cases how to split the data  structures on different memories and efficiently manage their accesses. The sizes in bytes reported in Table \ref{ch8:tabMemSizes}, are computed knowing that in our implementation the elements of $JM$ and $PTM$ are unsigned chars (one byte) and that the elements of $LM$, $RM$, $QM$ and $MM$ are unsigned short ints (2 bytes). It is important here to highlight that the types of the data of the used matrices impact the size of each matrix. For instance, a matrix of $100$ integers has a size of $400$ octets while the same matrix with $100$ unsigned chars has a size of $100$ octets. In order to minimize the size of each of the used matrices, we analyzed the ranges of their values and defined their data types accordingly. For instance, in PTM all the processing times have positive values varying between $0$ and $100$. Therefore, we defined PTM as a matrix of \verb|unsigned char| having values in the range $[0, 255]$. Using the \verb|unsigned char| type instead of the integer type allows us to reduce by $4$ times the memory space occupied by PTM.

\vspace{0.2cm} 

According to the Table \ref{ch8:tabMemSizes} :

\begin{itemize}
 \item The data structures $RM$, $QM$ and $MM$ are small sized matrices. Therefore, their impact on the performances is not significant whatever is the memory to which they are off-loaded. In particular, preliminary experiments proves that putting them on the shared memory would allows a very poor performance improvement. 
\item The $LM$ data structure is the double of the $JM$ in memory size but with a much lower access frequency. It is thus better to map $JM$ on the shared memory.
\item The $PTM$ has almost the same access frequency than $JM$ but requires less memory space.
\end{itemize}
  
\vspace{0.2cm} 

Consequently, the focus is put on the study of the performance impact of the placement of $JM$ and $PTM$ on the shared memory. Three placement scenarios of $JM$ and $PTM$ are experimented and studied: (1) Only $PTM$ is stored in shared memory and all others are placed in global memory~; (2) Only $JM$ is stored in shared memory and all others are placed on global memory~; (3) $PTM$ and $JM$ are stored together in shared memory and all others are placed on global memory. 

\vspace{0.2cm} 

Taking profit from the configurable storage space provided in the new Fermi-based devices, the $64$ KB of local storage was spitted between the shared memory and the L1 cache according to the experimented scenario.

\begin{itemize}
\item For the scenario were the data structures are put on the shared memory the $64$ KB of available storage are split on $48$ KB for shared memory and $16$ KB for L1 cache.
\item For the scenario where the data sets are put on global memory we used $16$ KB for shared memory and $48$ KB for L1 cache.
\end{itemize}

\section{Experiments}
\label{ch8:Experiments}

In the following, we present the experimental study  we have performed with the aim to evaluate the performance impact of the GPU-accelerated bounding, the techniques for reducing the thread divergence and the proposed approach for data placement on the GPU memories. 

\subsection{Parameters settings} 

In our experiments, we used the flow-shop instances defined by Taillard \cite{ch8:Taillard_1993}. These standard instances are often used in the literature to evaluate the performance of methods that minimize the makespan. Optimal solutions of some of these instances are still not known. These instances are divided into groups of $10$ instances. In each group, the $10$ instances  are defined by the same number of jobs and the same number of machines. The groups of 10 instances have different numbers of jobs, namely $20$, $50$, $10$, $200$ and $500$, and different numbers of machines, namely $5$, $10$ and $20$. For example, there are $10$ instances with $200$ jobs and $20$ machines belonging to the same group of instances.

\vspace{0.2cm} 

In this work, we used only the instances where the number of machines is equal to $20$. Indeed, instances where the number of machines is equal to $5$ or $10$ are easy to solve. For these instances, the used bounding operator gives so good lower bounds that it is possible to solve them in few minutes using a sequential B\&B. Therefore, these instances do not require the use of a GPU.

\vspace{0.2cm} 

Our approach has been implemented using C-CUDA 4.0. The experiments have been carried out using a an Intel Xeon E5520 bi-processor coupled with a GPU device. The bi-processor is 64-bit, quad-core and has a clock speed of 2.27GHz. The GPU device is an Nvidia Tesla C2050 with 448 CUDA cores (14 multiprocessors with 32 cores each), a clock speed of 1.15GHz, a 2.8GB global memory, a 49.15KB shared memory, and a warp size of 32.

\subsection{Experimental protocol: computing the speed up}
\label{ch8:Protocol}

We need to compute the speed up of our approach to evaluate its performances. This speed up is obtained by comparing our GPU B\&B version to a sequential B\&B version deployed on one CPU core. However, all the instances used in our experiments are extremely hard to solve. Indeed, the resolution of each of these instances requires several months of computation on one CPU core. For example, the optimal solution of one of these instances defined by $50$ jobs and $20$ machines is obtained after $25$ days of computation using an average of $328$ CPU cores \cite{ch8:Mezmaz_2007}. 

\vspace{0.2cm} 

Using the approach defined in \cite{ch8:Mezmaz_2007}, it is possible to obtain a random list $L$ of sub-problems such as the resolution of $L$ lasts $T$ minutes with a sequential B\&B. So by initializing the pool of our sequential B\&B with the sub-problems of this list $L$, we are sure that the resolution of the sequential B\&B will last $T{cpu}$ minutes such as $T{cpu}$ will be approximately equal to $T$. Therefore, it will be possible to initialize the pool of our GPU B\&B with the same list $L$ of sub-problems in order to compute the speed up. Let suppose that the resolution of the GPU B\&B will last $T{gpu}$ minutes. So the speed up of our GPU algorithm will be equal to $Tcpu/Tgpu$. With this experimental protocol, the sub-problems explored by the GPU and CPU B\&B versions will be exactly the same. So to find the speed up associated to an instance, we:

\begin{itemize}
\item compute, using the approach defined in \cite{ch8:Mezmaz_2007}, a list $L$ of sub-problems such as the resolution of $L$ lasts $T$ minutes with a sequential B\&B, 
\item initialize the pool of our sequential B\&B with the sub-problems of this list $L$, 
\item solve the sub-problems of this pool with our sequential B\&B ,
\item get the sequential resolution time $T{cpu}$ and the number of explored sub-problems $N{cpu}$, 
\item check that $T{cpu}$ is approximately equal to $T$, 
\item initialize the pool of our GPU B\&B with the sub-problems of the list $L$, 
\item solve the sub-problems of this pool with our GPU B\&B,
\item get the GPU resolution time $T{gpu}$ and the number of explored sub-problems $N{gpu}$, 
\item check that $N{gpu}$ is exactly equal to $N{cpu}$, 
\item and finally compute the speed up associated to this instance by dividing $T{cpu}$ on $T{gpu}$ (i.e. $Tcpu/Tgpu$).
\end{itemize}

\vspace{0.2cm} 

Table \ref{ch8:instance_time} gives, for each instance according to its number of jobs and its number of machines, the used resolution time with a sequential B\&B. For example, the sequential resolution time of each instance defined with $20$ jobs and $20$ machines is approximately 10 minutes. Of course, the computation time of the lower bound of a sub-problem defined with $20$ jobs and $20$ machines is on average greater than the computation time of the lower bound of a sub-problem defined with $50$ jobs and $20$ machines. Therefore, as shown in this table, the sequential resolution time increases with the size of the instance in order to be sure that the number of sub-problems explored is significant for all instances.  

\begin{table*}
\setlength{\tabcolsep}{0.2cm}
\renewcommand{\arraystretch}{1.2}
\centering
\footnotesize
\begin{tabular}{|r|r|r|r|r|}
\hline
Instance (No. of jobs x No. of machines) & 20$\times$20 & 50$\times$20  & 100$\times$20 & 200$\times$20 \\
\hline
Sequential resolution time (minutes) & 10 & 50  & 150 & 300 \\
\hline
\end{tabular}
\vspace{0.3cm}
\caption{The sequential resolution time of each instance according to its number of jobs and machines}
\label{ch8:instance_time}
\end{table*}

\subsection{Performance impact of GPU-based parallelism}

The objective of the experimental study presented in this section is to compared the performances of both proposed approaches for designing B\&B on top of GPUs. 

Table \ref{ch8:ParaGPU1} and Table~\ref{ch8:ParaGPU2} report respectively the speedups obtained with the GPU-PTE-BB and GPU-PEB-BB approaches for different problem instances. The first part of both tables gives the size of the pool generated and evaluated on the GPU. The second part of the tables gives the average speedup for each group of instances and for each pool size. Each line corresponds to a group of $10$ instances defined by the same number of jobs and the same number of machines. 

The results obtained with the GPU-PTE-BB approach (see Table \ref{ch8:ParaGPU1}) show that exploring in parallel the tree search allows to speedup the execution of the B\&B compared to a CPU-based execution. Indeed, an acceleration factor up to 40.50 is obtained for the 20 $\times$ 20 problem instances using a pool of 262144 sub-problems. 

The results show also that the parallel efficiency decreases with the size of the problem instance. For a fixed number of machines (here 20 machines) and a fixed pool size, the obtained speedup decline accordingly with the number of jobs. For instance for a pool size of 262144, the acceleration factor obtained with 200 jobs (13.4) while it is (40.50) for the instances with 20 jobs. This behavior is mainly due to the overhead induced by the transfer of the pool of resulting sub-problems between the CPU and the GPU. For example, for the instances with 200 jobs the size of the pool to exchange between the CPU and the GPU is ten times bigger than the size of the pool for the instances with 20 jobs.

\begin{table*}
\setlength{\tabcolsep}{0.2cm}
\renewcommand{\arraystretch}{1.2}
  \centering
  \footnotesize
 \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines) & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
$200 \times $20 & 1.12 & 2.89 & 3.57 & 4.23 & 6.442 & 8.32 & 13.4\\
    \hline
$100 \times $20 & 1.33 & 1.88 & 3.45 & 6.45 & 12.38 & 20.40 & 28.76 \\
    \hline
$50 \times $20 & 2.70 & 3.80 & 6.82 & 13.04 & 23.53 & 30.94 & 37.66\\
    \hline
$20 \times $20 & 6.43 & 11.43 & 20.14 & 27.78 & 30.12 & 35.74 & 40.50\\
    \hline
    \hline
% Total average speedup & 2.895 & 5.0 & 8.495 & 14.625 & 22.61 & 30.6 & 41.65\\   
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
  \caption{Speedups for different problem instances and pool sizes with the GPU-PTE-BB approach.}
\label{ch8:ParaGPU1}
\end{table*}

The results obtained with the GPU-PEB-BB approach (see Table \ref{ch8:ParaGPU2}) show that evaluating in parallel the bounds of a selected pool, allow to significantly speedup the execution of the B\&B. Indeed, an acceleration factor up to 71.69 is obtained for the 200 $\times$ 20 problem instances using a pool of 262144 sub-problems. The results show also that the parallel efficiency grows with the size of the problem instance. For a fixed number of machines (here 20 machines) and a fixed pool size, the obtained speedup grows accordingly with the number of jobs. For instance for a pool size of 262144, the acceleration factor obtained with 200 jobs (71.69) is almost the double of the one obtained with 20 jobs (38.40). 

As far the pool size tuning is considered, we could notice that this parameter depends strongly on the problem instance being solved. Indeed, while the best acceleration is obtained with a pool size of 8192 sub-problems for the instances 50 $\times$ 20 and 20 $\times$ 20, the best speedups are obtained with a pool size of 262144 sub-problems with the instances 200 $\times$ 20 and 100 $\times$ 20.\\

\begin{table*}
\setlength{\tabcolsep}{0.2cm}
\renewcommand{\arraystretch}{1.2}
  \centering
  \footnotesize
 \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines) & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
$200 \times $20 & 42.83 & 56.23 & 57.68 & 61.21 & 66.75 & 68.30 & \textbf{71.69}\\
    \hline
$100 \times $20 & 42.59 & 56.18 & 57.53 & 60.95 & 65.52 & 65.70 & \textbf{65.97}\\
    \hline
$50 \times $20 & 42.57 & \textbf{56.15} & 55.69 & 55.49 & 55.39 & 55.27 & 55.14\\
    \hline
$20 \times $20 & 38.74 & \textbf{46.47} & 45.37 & 41.92 & 39.55 & 38.90 & 38.40\\
    \hline
    \hline
% Total average speedup & 41.68 & 53.76 & 54.07 & 54.89 & 56.80 & 57.04 & 57.80\\
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
  \caption{Speedups for different problem instances and pool sizes with the GPU-PEB-BB approach.}
\label{ch8:ParaGPU2}
\end{table*}

Compared to the parallel tree exploration-based GPU-accelerated B\&B approach, the parallel evaluation of bounds approach is by far much more efficient wherever the instance is. For example, while the GPU-PEB-BB approach reaches speedup of $\times$71.69 for the instance with 200 jobs on 20 machines, a speedup of a $\times$13.4 is measured with the parallel tree exploration-based approach which corresponds to an acceleration of $\times$5.56 . Moreover, on the contrary to the GPU-PEB-BB approach, in the GPU-PTE-BB the speedups decrease when the problem instance becomes higher. Remember here that while in the GPU-PEB-BB approach all threads evaluate only one node each whatever the permutation size is. In the GPU-PTE-BB, each thread branches all the children of its assigned parent node. Therefore, the bigger the size of the permutation is, the bigger the amount of work performed by each thread is and the bigger the difference between the workload is. Indeed, let us suppose that for the instance with $200$ jobs, the thread $0$ handles a node from the level $2$ of the tree and the thread $100$ handles a node from the level $170$ of the tree. In this case, the thread $0$ generates and evaluates $198$ nodes while the thread $100$ decomposes and bounds only $30$ nodes. The problem in this example is that the kernel execution would last until the thread $0$ finishes its work while the other threads might have ended their works and stayed idle.

\subsection{Thread divergence reduction}

The objective of this section is to demonstrate that the thread divergence reduction mechanisms we propose has an impact on the performance of the GPU accelerated B\&B and to evaluate how this impact is significant. 
In the following, the reported results are obtained with the GPU-accelerated B\&B based on the parallel evaluation of bounds.

\begin{table*}[!h]
\setlength{\tabcolsep}{0.2cm}
\renewcommand{\arraystretch}{1.2}
  \centering
  \footnotesize
 \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines)  & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
    \hline
$200 \times $20 & 46.63 & 60.88 & 63.80 & 67.51 & 73.47 & 75.94 & \textbf{77.46}\\
    \hline
$100 \times $20 & 45.35 & 58.49 & 60.15 & 62.75 & 66.49 & 66.64 & \textbf{67.01}\\
    \hline
$50 \times $20 & 44.39 & \textbf{58.30} & 57.72 & 57.68 & 57.37 & 57.01 & 56.42\\
    \hline
$20 \times $20 & 41.71 & \textbf{50.28} & 49.19 & 45.90 & 42.03 & 41.80 & 41.65\\
    \hline
    \hline
% Total average speedup & 44.52 & 56.99 & 57.72 & 58.46 & 59.84 & 60.35 & 60.64\\
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
  \caption{Speedups for different instances and pool sizes using thread divergence management.}
\label{ch8:ParaDivergence}
\end{table*}

Table~\ref{ch8:ParaDivergence} shows the experimental results obtained using the sorting process and the refactoring approach presented in Section \ref{ch8:ThreadDivergence}. Results show that the proposed optimizations emphasize the GPU acceleration reported in Table~\ref{ch8:ParaGPU2} and obtained without thread divergence reduction. For example, for the instances of 200 jobs over 20 machines and a pool size of 262144, the average reported speedup is 77.46 while the average acceleration factor obtained without thread divergence management for the same instances and the same pool size is 71.69 which corresponds to an improvement of 7.68\%. Such considerable but not outstanding improvement is predictable, as claimed in \cite{ch8:Han}, since the factorized part of the branches in the FSP lower bound is very small. 

\subsection{Data access optimization}

The objective of the experimental study presented in this section is to find the best mapping of the six data structures of the lower bound LB kernel on the memories of the GPU device. In the following, the reported results are obtained with the GPU-accelerated B\&B based on the parallel evaluation of bounds.

Table~\ref{ch8:PTM-on-SM} reports the speedups obtained for the first experimented scenario where only the matrix $PTM$ is put on the shared memory. Results show that the speedup grows on average with the growing of the pool size in the same way as in Table~\ref{ch8:ParaDivergence}. For the largest problem instance and pool size, putting the PTM matrix on the shared memory improves the speedups up to ($14\%$) compared to those obtained when $PTM$ is on global memory reaching an acceleration of $\times 90.51$ for the problem instances $200 \times 20$ and a pool size of $262144$ sub-problems .  

\begin{table*}
  \centering
  \footnotesize
  \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines)  & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
    \hline
$200 \times $20 & 54.03 & 67.75 & 68.43 & 72.17 & 82.01 & 88.35 & \textbf{90.51}\\
    \hline
$100 \times $20 & 52.92 & 66.57 & 66.25 & 71.21 & 76.63 & 79.76 & \textbf{83.01}\\
    \hline
$50 \times $20 & 49.85 & \textbf{65.68} & 64.40 & 59.91 & 58.57 & 57.36 & 55.09\\
    \hline
$20 \times $20 & 41.94 & \textbf{60.10} & 48.28 & 39.86 & 39.61 & 38.93 & 37.79 \\
     \hline
    \hline
% Average Speedup& 49.69 & 65.03 & 61.84 & 60.79 & 64.21 & 66.10 & 66.60 \\                                             
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
 \caption{Speedup for different FSP instances and pool sizes obtained with data access optimization. $PTM$ is placed in shared memory and all others are placed in global memory.}
\label{ch8:PTM-on-SM}
\end{table*}

Table~\ref{ch8:JM-on-SM} reports the behavior of the speedup averaged on the different problem instances (sizes) as a function of the pool size for the scenario where the Johnson's matrix is put on the shared memory. Results show that putting the $JM$ matrix on the shared matrix improves more the performances comparing to the first scenario where $PTM$ is put on the shared memory. Indeed, according to Table~\ref{ch8:tabMemComplex}, matrix $JM$ is accessed more frequently than matrix $PTM$. Putting $JM$ matrix on the shared memory allows accelerations up to $\times 97.83$ for the problem instances $200 \times 20$.

\begin{table*}
  \centering
  \footnotesize
  \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines) & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
    \hline
$200 \times $20 & 63.01 & 79.40 & 81.40 & 84.02 & 93.61 & 96.56 & \textbf{97.83}\\
    \hline
$100 \times $20 & 61.70 & 77.79 & 79.32 & 81.25 & 86.73 & 87.81 & \textbf{88.69}\\
    \hline
$50 \times $20 & 59.79 & \textbf{75.32} & 72.20 & 71.04 & 70.12 & 68.74 & 68.07 \\
    \hline
$20 \times $20 & 49.00 & \textbf{60.25} & 55.50 & 45.88 & 44.47 & 43.11 & 42.82 \\
     \hline
    \hline
% Average Speedup& 58.37 & 73.19 & 72.11 & 70.55 & 73.73 & 74.06 & 74.35 \\ 
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
 \caption{Speedup for different FSP instances and pool sizes obtained with data access optimization. 
$JM$ is placed in shared memory and all others are placed in global memory.}
\label{ch8:JM-on-SM}
\end{table*}            

Table~\ref{ch8:JM-PTM-on-SM} reports the behavior of the average speedup for the different problem instances (sizes) with $20$ machines for the data placement scenario where both $PTM$ and $JM$ are put on shared memory. According to the underlying Table, the scenarios~(3) ($JM$ together or without $PTM$ in shared memory) is clearly better than the scenarii~(1)and~(2) (respectively $PTM$ in shared memory and $JM$ in shared memory) whatever is the problem instance (size). 

\begin{table*}
  \centering
  \footnotesize
  \begin{tabular}{|r|r|r|r|r|r|r|r|}
    \hline
Pool size & 4096 & 8192 & 16384 & 32768 & 65536 & 131072 & 262144\\
    \hline
    \hline
(NJobs $\times$ NMachines)  & \multicolumn{7}{|c|}{Average speedup for each group of 10 instances}\\
    \hline
    \hline
$200 \times $20 & 66.13 & 87.34 & 88.861 & 95.23 & 98.83 & 99.89 & \textbf{100.48}\\
    \hline
$100 \times $20 & 65.85 & 86.33 & 87.60 & 89.18 & 91.41 & 92.02 & \textbf{92.39}\\
    \hline
$50 \times $20 & 64.91 & \textbf{81.50} & 78.02 & 74.16 & 73.83 & 73.25 & 72.71\\
    \hline
$20 \times $20 & 53.64 & \textbf{61.47} & 59.55 & 51.39 & 47.40 & 46.53 & 46.37\\
     \hline
    \hline
% Average Speedup & 62.63 & 79.16 & 78.51 & 77.49 & 77.87 & 77.92 & 77.99\\
%     \hline
%     \hline
  \end{tabular}
\vspace{0.3cm}
 \caption{Speedup for different FSP instances and pool sizes obtained with data access optimization. 
$PTM$ and $JM$ are placed together in shared memory and all others are placed in global memory.}
\label{ch8:JM-PTM-on-SM}
\end{table*}    

By carefully analyzing each of the scenarii of data placement on the memory hierarchies of the GPU, the recommendation is to put in the shared memory the Johnson's and the processing time matrices ($JM$ and $PTM$) if they fit in together. Otherwise, the whole or a part of the Johnson's matrix has to be put in priority in the shared memory. The other data structures are mapped to the global memory.

\section{Conclusion and Future Work}
\label{ch8:Conclusion}

In this chapter, we have revisited the design of parallel B\&B algorithms on GPU accelerators to allow highly efficient solving of permutation-based COPs. To do so, our contributions consist in: (1) rethinking two approaches for parallel B\&B on top of GPUs, discussing the performances of each and identifying which best suits the GPU accelerators. (2) proposing a new approach for thread/branch divergence reduction through a thorough analysis of the different loops and conditional instructions of the bounding function. (3) defining an optimal mapping of the data structures of the bounding function on the hierarchy of memories provided in the GPU device through a careful analysis of both the data structures (size and access frequency) and the GPU memories (size and access latency). 

In the first parallel tree-exploration-based B\&B, a set of pending nodes is selected from this list according to their depth and off-loaded to the GPU where each thread builds its own local search tree by applying
the branching, bounding and pruning operators to the assigned node. In the GPU-accelerated B\&B based on the parallel evaluation of bounds, the generation of the sub-problems (branching, selection and pruning operations) is performed on CPU and the evaluation of their lower bounds (bounding operation) is executed on the GPU device. Pools of sub-problems are off-loaded from CPU to GPU to be evaluated by blocks of threads. After evaluation, the lower bounds are returned to the CPU.

In both considered approaches, our focus is on the GPU-based lower bound's implementation and the associated thread divergence and data placement challenges. The proposed mechanisms for reducing the thread divergence issue are based on a thorough analysis of the different loops and conditional instructions of the lower bound function. On the one hand, the sorting process aims to homogenize the data of the sub-problems off-loaded to the GPU to minimize the number of threads that diverge on loop instructions. On the other hand, the technique of branch refactoring rewrite the conditional instructions into uniform instructions so that threads of the same warp execute a same code. The proposed data access optimization is based on a preliminary analysis of the lower bound function. Such analysis allowed us to identify six data structures for which we have proposed a complexity analysis in terms of memory size and access frequency. Due to the limited size of the shared memory the matrices do not fit in all together. According to the complexity study, the recommendation is to put in the shared memory the Johnson's and the processing time matrices ($JM$ and $PTM$) if they fit in together. Otherwise, the whole or a part of the Johnson's matrix has to be put in priority in the shared memory. The other data structures are mapped to the global memory. Such recommendation has been confirmed through extensive experiments using a recent C2050 Tesla GPU card. 

The Flowshop Scheduling Problem has been considered as a case study. The proposed approaches have been experimented using a Tesla C2050 GPU card on different classes of FSP instances. The experimental results show that the parallel evaluation of bounds is the parallelization paradigm that performs better on top of GPU accelerators. Compared to the parallel tree-exploration model, accelerations up to $\times$5.56 are achieved.

Experiments show also that the proposed refactoring approach improves the parallel efficiency whatever the FSP instance and the pool size are. However, the improvement was not significant because the factorized part of the branches in the FSP lower bound is very small. The optimizations obtained with the proposed thread reduction mechanisms allowed us to achieve accelerations up to $\times$77.46 compared to a sequential B\&B. The data access optimizations grant accelerations up to $\times 100$ compared to a single CPU-based B\&B.

In the near future, we plan to extend this work to a cluster of GPU-accelerated multi-core processors. From the application point of view, the objective is to optimally solve challenging and unsolved Flow-Shop instances as we did it for one 50$\times$20 problem instance with grid computing \cite{ch8:Mezmaz_2007}. Finally, we plan to investigate other lower bound functions to deal with other combinatorial optimization problems. 

\putbib[Chapters/chapter8/biblio8]