1 \chapterauthor{Malika Mehdi and Ahc\`{e}ne Bendjoudi}{CERIST Research Center, DTISI, 3 rue des frères Aissou, 16030 Ben-Aknoun, Algiers, Algeria}
3 \chapterauthor{Lakhdar Loukil}{University of Oran, Algeria}
5 %\chapterauthor{Ahc\`{e}ne Bendjoudi}{CERIST Research Center, DTISI, 3 rue des frères Aissou, 16030 Ben-Aknoun, Algiers, Algeria}
7 \chapterauthor{Nouredine Melab}{Université Lille 1, LIFL/UMR CNRS 8022, 59655-Villeneuve d'Ascq cedex, France}
9 \chapter{Parallel GPU-accelerated metaheuristics}
11 \section{Introduction}
12 This chapter presents GPU-based parallel
13 metaheuristics\index{metaheuristics!parallel metaheuristics},
14 challenges, and issues related to the particularities of the GPU
15 architecture and a synthesis on the different implementation
16 strategies used in the literature. The implementation of parallel
18 GPUs is not straightforward. The traditional models used in CPUs
19 must be rethought to meet the new requirements of GPU
20 architectures. This chapter is organized as follows. Combinatorial
21 optimization\index{combinatorial optimization} and resolution
22 methods are introduced in Section~\ref{ch8:sec:optim}. The main
23 traditional parallel models used for metaheuristics are recalled
24 in Section~\ref{ch8:sec:paraMeta}.
25 Section~\ref{ch8:sec:challenges} highlights the main challenges
26 related to the GPU implementation of metaheuristics. A
27 state-of-the-art of GPU-based parallel
29 summarized in Section~\ref{ch8:sec:state}. In Section~\ref{ch8:sec:frameworks}, the main developed GPU-based
30 frameworks for metaheuristics are described. Finally, a case study
31 is presented in Section~\ref{ch8:sec:case} and some concluding
32 remarks are given in Section~\ref{ch8:conclusion}
34 \section{Combinatorial optimization}
37 Combinatorial optimization (CO) is a branch of applied and discrete mathematics.
38 It consists in finding optimal configuration(s) among a finite set of possible configurations
39 (or solutions) of a given combinatorial optimization problem (COP). The set of all possible solutions noted $S$ is called solution space or search space. Each solution in $S$ is defined by its real cost calculated by an objective function. COPs are generally defined as follows~\cite{blumMeta}:\\ %(Definition~\ref{def:cops})
42 %\begin{minipage}{0.5\linewidth}
43 A combinatorial problem $P=(S,f)$ can be defined by
45 \item a set of decision variables $X$,
46 \item an objective function $f$ to optimize (minimize or maximize) over the set $S$,
47 \item subject to constraints on the decision variables.\\
51 COPs are generally formulated as mixed integer programming
52 problems (MIPS) and most of them are NP-hard~\cite{garey}.
53 According to the quality level of solutions and deadlines required
54 for solving an optimization problem, two classes of optimization
55 methods can be distinguished: \emph{exact methods} and
56 \emph{approximate methods}. Exact methods allow one to reach
57 optimal solution(s) of the handled optimization problem with a
58 proof of its or their optimality. The known methods of this
59 class are the \emph{branch and bound technique}, \emph{dynamic
60 programming}, \emph{constraint programming}, and \emph{A*
61 algorithm}. However, optimization problems, whether practical or
62 academic, are often complex and NP-hard. Moreover, a large number
63 of real-life optimization problems encountered in science,
64 engineering, economics, and business are usually large-sized
65 problems for which the size of the potential solution domain
66 increases dramatically with the size of the problem instance. Such
67 problems cannot be tackled using exact methods due to the
68 excessive computation time needed by these methods to find optimal
69 solution(s). In such a situation, approximate methods (or
70 \textit{metaheuristics}) offer an alternative approach to solve
71 these problems. Indeed, these methods allow one to reach good
72 quality solutions in reasonable computation time compared to exact
73 methods but with no guarantee to find optimal or even bounded
74 solutions. This is due to the nature of the search process adopted
75 by these approaches which consists of performing a
76 search in a subset of the whole search space.\\
78 Regarding the number of solutions considered at each iteration in
79 the search process, two classes of metaheuristics can be
80 distinguished~\cite{talbi2009mfdti}: \textit{solution-based} and
81 \textit{population-based} metaheuristics. In the rest of this
82 chapter, the term \emph{s-metaheuristic} refers to solution-based
83 metaheuristic and the term \emph{p-metaheuristic} refers to
84 population-based metaheuristic. In s-metaheuristics, the search
85 process starts with a single solution (generally set at random)
86 and iteratively improves it by exploring its neighborhood in the
87 search space. The most known methods in this class are local
88 search methods that include \emph{simulated
89 annealing}\index{metaheuristics!simulated annealing}~\cite{Kirkpatrick1983SA},
90 \emph{tabu search}~\cite{Glover1989TS}, \emph{iterated local
91 search\index{metaheuristics!iterated local
92 search}}~\cite{stutzle2006ILSforQAP}, and \emph{variable
93 neighborhood search}~\cite{HansenMladenovic1997VNS}.\\
95 Unlike s-metaheuristics, p-metaheuristics start with a population
96 of solutions and implement an iterative process that evolves the
97 current population towards a new population of better quality
98 solutions. The process is repeated until a stopping criterion is
99 satisfied. \emph{Evolutionary algorithms}, \emph{swarm
100 optimization}, and \emph{ant colonies} fall into this class.
103 \section{Parallel models for metaheuristics}\label{ch8:sec:paraMeta}
104 Optimization problems, whether real-life or academic, are more
105 often NP-hard and CPU time and/or memory consuming. Metaheuristics
106 allow the significant reduction of the computational time of the search
107 process but remain time-consuming particularly when it comes
108 dealing with large-sized problems. \\
110 The use of parallelism in the design of metaheuristics is a relevant
111 approach that is widely adopted by the combinatorial optimization
112 community for various reasons:
115 \item One of the main goals of parallelism is to reduce the search
116 time. This will allow the design of high performance optimization
117 methods and the solving of large-sized optimization problems.
119 \item Sequential processor architectures have reached their
120 physical limit which prevents creating faster processors. The
121 current trend of microprocessor manufacturers consists of placing
122 multicores on a single chip. Nowadays, laptops and
123 workstations are multicore processors. In addition, the evolution
124 of network technologies and the proliferation of broadband
125 networks have made possible the emergence of clusters of
126 workstations (COWs), networks of workstations (NOWS), and
127 large-scale networks of machines (grids) as platforms for
128 high-performance computing.
131 From the granularity of parallelism point of view, three major parallel
132 models for metaheuristics can be distinguished~\cite{talbi2009mfdti}: \emph{algorithmic-level}\index{metaheuristics!algorithmic-level parallelism},
133 \emph{iteration-level}\index{metaheuristics!iteration-level parallelism}, and \emph{solution-level} as illustrated in Figure~\ref{ch8:fig:paraMeta}. \\
136 \centerline{\includegraphics[width=0.6\textwidth]{Chapters/chapter9/figures/paraMeta.pdf}}
137 \caption{Parallel models for metaheuristics.}
138 \label{ch8:fig:paraMeta}
144 algorithmic-level parallel model, several self-contained metaheuristics
145 are launched in parallel. The parallel
146 metaheuristics\index{metaheuristics!parallel metaheuristics} may
147 start with identical or different solutions (s-metaheuristics
148 case) or populations (p-metaheuristics case). Their parameter
149 settings such as the size of tabu list for tabu
150 search\index{metaheuristics!tabu search}, transition probabilities
151 for ant colonies, mutation and crossover probabilities for
153 algorithm\index{metaheuristics!evolutionary algorithm}s may be the
154 same or different. The parallel processes may evolve independently
155 or in a cooperative manner. In cooperative parallel models, the
156 algorithms exchange information related to the search during
157 evolution in order to find better and more robust solutions.}
159 \item{In the iteration-level parallel model, the focus is on the parallelization
160 of each iteration of the metaheuristic. Indeed, metaheuristics are
161 generally iterative search processes. Moreover, the most
162 resource-consuming part of a metaheuristic is the evaluation of
163 the generated solutions at each iteration. For s-metaheuristics
164 (e.g., tabu search\index{metaheuristics!tabu search}, simulated
165 annealing, variable neighborhood search), the evaluation and
166 generation of the neighborhood is the most time-consuming step of
167 the algorithm particularly when it comes to dealing with large
168 neighborhood sets. In this parallel model, the neighborhood is
169 decomposed into partitions, and each partition is evaluated in a
170 parallel way. For p-metaheuristics (e.g., evolutionary
172 colonies, swarm optimization), the
174 parallel model arises naturally since these metaheuristics deal
175 with a population of independent solutions. In evolutionary
177 instance, the iteration-level model consists of decomposing the whole population
178 into several partitions where each partition is evaluated in
183 parallel model, the focus is on the parallelization of the
184 evaluation of a single solution. This model is useful when the
185 objective function and/or the constraints are time and/or memory
186 consuming. Unlike the two previous parallel models, the
187 solution-level\index{metaheuristics!solution-level parallelism}
188 parallel model is problem-dependent.}
191 \section{Challenges for the design of GPU-based metaheuristics}
192 \label{ch8:sec:challenges}
194 Developing GPU-based parallel
195 metaheuristics\index{metaheuristics!parallel metaheuristics} is
196 not straightforward. The parallel models have to be rethought to
197 meet the new requirements of the GPU architecture. Several major
198 issues have to be taken into account both at design and
199 implementation levels to develop efficient metaheuristics on GPU.
200 These issues are mainly related to the size and latency of the GPU
201 memories, thread synchronization and divergence, the distribution
202 of tasks, and data transfer between the CPU and
203 GPU~\cite{luongMultiStart}.
205 \subsection{Data placement on a hierarchical memory}
206 \index{GPU!data placement} During the execution of
207 metaheuristics on GPU, the different threads may access multiple
208 data structures from multiple memory spaces. These memories have
209 different sizes and access latencies. Nevertheless, faster
210 memories (registers, shared and constant memories) are generally
211 very small in size, and the larger memories (global memory) are
212 relatively slower. However, p-metaheuristics require the
213 exploration of a large amount of individuals to diversify the
214 search. Moreover, an efficient execution of s-metaheuristics
215 requires exploring large neighborhoods. Thus, programmers have to
216 take into account this point to efficiently place the different
217 data structures of the metaheuristic on the different memories to
218 benefit from both the faster memories and the larger ones. A deep
219 study has to be performed on both the metaheuristic data
220 structures (size and access frequency) and the GPU memories (size
221 and access latency) to identify which data will be placed on which
222 memory. Generally, the most accessed ones should be put on faster
223 memories (registers, shared memory) and larger ones on the larger
224 memories (global memory). Also, an efficient mapping between
225 threads and the corresponding metaheuristic elements (one neighbor
226 per thread, one individual per thread, single population per
227 threads block, one ant per thread, etc.) must be defined to ensure
228 a maximum occupancy of the GPU and
229 to cover CPU/GPU communication\index{GPU!CPU/GPU communication} and memory access times.\\
231 According to the used metaheuristic and to the handled problem, the data
232 values may have different types and different ranges of their values. The data
233 types should then be chosen carefully and the ranges of the data values should
234 be analyzed to minimize the amount of occupied memory space.\\
236 In addition to the size and latency of GPU memory issues, the
237 memory access pattern is another important issue to be dealt with
238 to speedup GPU-based metaheuristics. Indeed, the different
239 memories have been designed to achieve specific features that
240 programmers must take into account to optimize their codes and
241 then to benefit from these features. For instance, the global
242 memory is optimized for coalesced accesses. The texture and the
243 constant memory are read-only memories. The texture is optimized
244 for uncoalesced accesses and the constant memory is optimized for
245 simultaneously accesses of all threads to the same
246 location~\cite{luongMultiStart}. Therefore, to improve the
247 performance of the kernel execution, the programmers should put
248 coalesced data on global memory, uncoalesced read-only data (e.g.,
249 metaheuristic instance data) on the texture, concurrent read-only
250 data (e.g., data for fitness evaluation that all threads
251 concurrently access) on the constant memory, and the most accessed
252 data structures (e.g., population of individuals for a CUDA thread
253 block) on the shared memory.
255 \subsection{Threads synchronization}
256 \index{GPU!threads synchronization} The thread
257 synchronization issue is caused by both the GPU architecture and
258 the synchronization requirements of the implemented method.
259 Indeed, GPUs are based on a multicore architecture organized into
260 several multi-processors (Streaming Multiprocessor SM) supporting
261 the SPMD model (Single Program Multiple Data). Each SM contains
262 several cores executing the same instruction of different threads
263 following the SIMD model (Single Instruction Multiple Data). These
264 threads belong to a warp (a group of 32 threads) and handle
265 different data elements. An efficient execution of applications on
266 GPU is achieved when launching a large amount of threads
267 (thousands of threads)~\cite{CUDA}. However, the execution order
268 of these thousands of threads is unknown by the programmer which
269 makes the prediction of their execution order a challenging issue.
270 Plus, the developer has to control explicitly the
271 threads through the insertion of barrier synchronizations in the
272 codes to avoid concurrent accesses to data structures and to meet
273 some requirements related to data-dependent synchronizations.
275 \subsection{Thread divergence}
277 Thread divergence\index{GPU!thread divergence} is
278 another challenging issue in GPU-based
279 metaheuristics~\cite{cecilia, pugace, audreyANT}. Generally,
280 metaheuristics contain irregular loops and conditional
281 instructions when generating and evaluating the neighborhood
282 (s-metaheuristics), and the population (p-metaheuristics) in the
283 same block. In addition, the decision to apply a crossover or a
284 mutation on an individual in a genetic
286 exploration of different paths using an ant
287 colony\index{metaheuristics!ant colony optimization} are random
288 operations. Threads of the same warp have to execute
289 instructions simultaneously leading to different branches whereas
290 in an SIMD model the threads of a same warp execute the same
291 instruction. Consequently, the different branches of a conditional
292 instruction which is data-dependent lead to a serial execution of
293 the different threads degrading the performance of the application
294 in terms of execution time. The challenge here is then to revisit
295 the traditional irregular metaheuristic codes to eliminate these
298 \subsection{Task distribution and CPU/GPU communication}
300 The performance of GPU-based metaheuristics in terms of execution
301 time could be improved by choosing the most appropriate parallel
302 model (algorithmic-level, instruction-level,
304 Moreover, an efficient decomposition of the metaheuristic and an
305 efficient assignment of code portions between the CPU and GPU
306 should be adopted. The objective is to take benefit from the GPU
307 computing power without affecting the efficiency and the behavior
308 of the metaheuristic and without losing performance in CPU/GPU
309 communication\index{GPU!CPU/GPU communication} and
310 memory accesses. In order to decide which part of the
311 metaheuristic will be executed on which component, one should
312 perform a careful analysis on the serial code of the
313 metaheuristic, identify the compute-intensive tasks (e.g.,
314 exploration of the neighborhood, evaluation of individuals), and
315 then offload them to the GPU, while the remaining
316 tasks still run on the CPU in a serial way. \\
318 The CPU/GPU communication is done through the global
319 memory which is a slow memory making the memory transfer between
320 the CPU and GPU time-consuming which can significantly degrade the
321 performance of the application. Accesses to this memory should be
322 optimized by minimizing the amount of transferred data between the
323 two components in order to reduce the communication time and,
324 therefore, the whole execution time of the metaheuristic.
326 \section{State-of-the-art parallel metaheuristics on GPUs}
327 \label{ch8:sec:state}
328 After more than two decades of research by the combinatorial optimisation
329 community devoted to developing adequate parallel metaheuristics for different types of
330 parallel architectures (clusters, supercomputers and grids), the actual developement
331 of General Perpose GPU (GPGPU) brings new challenges for parallel metaheuristics on SIMD architectures.\\
333 The first works on metaheuristic algorithms implemented on GPUs
334 started on old graphics cards before the appearance of modern GPUs
335 equipped with high-level programming interfaces such as CUDA and
336 OpenCL. Among these pioneering works we cite the work of Wong et al.~\cite{wongOldGPU2006} dealing with the
338 of EAs on graphics processing cards and the work by Catala et al. in~\cite{catala2007} where the ACO\index{metaheuristics!ant colony optimization} algorithm
339 is implemented on old GPU architectures. Yu et al.~\cite{yu2005} and
340 Li et al.~\cite{li2007} proposed a full parallelization of genetic
341 algorithms on old GPU architectures using
342 shader libraries based on Direct3D and OpenGL.\\
344 Such architectures are based on preconfigured pipelined stages
345 used to accelerate the transformation of 3D geometric primitives
346 into pixels. Implementing a general-purpose algorithm on such
347 preconfigured architectures is very hard and requires the
348 tailoring of the algorithm's data and instructions to fit the
349 pipelined stages of the GPU. Since then, GPU architectures have
350 evolved to become programmable using high-level programming
351 interfaces. In this section we will focus only on recent
352 state-of-the-art works dealing with metaheuristics implementation
353 on modern programmable GPUs. In this review two classes are
354 considered: (1) s-metaheuristics on GPUs and (2) p-metaheuristics
355 on GPUs. A comparative study is done of the main works and a
356 classification of the different existing strategies is proposed in
357 Section~\ref{ch8:sec:synthesis}.
359 \subsection{Implementing solution-based metaheuristics on GPUs}
362 A very basic local search algorithm starts with an initial
363 solution generated either at random or by the mean of a specific
364 heuristic and is based on two elementary components: the
365 generation of neighborhood structures using an elementary move
366 function and a selection function that determines which solution
367 in the current neighborhood will replace the actual search point.
368 The search continues as long as there is improvement and stops
369 when there is no better solution in the current neighborhood. The
370 exploration (or evaluation) of the different moves of a given
371 neighborhood is done independently for each move. Thus, the
372 easiest way to accelerate a local search algorithm is to
373 parallelize the evaluation of the neighborhood (instruction-level
374 parallelism). This is by far the most used scheme in the
375 literature for parallelizing local search algorithms on GPUs.
376 Nevertheless, small neighborhoods may lead to nonoptimal
377 occupation of the CUDA threads which may lead in turn to an
378 overhead due to the communication and memory latencies. Therefore,
379 large neighborhoods are necessary for efficient implementation of
380 local searches on GPUs.\\
382 Luong et al.~\cite{luong2010large} proposed efficient
383 mappings for large neighborhood structures on GPUs. In this work,
384 three different neighborhoods are studied and mapped to the
385 hierarchical GPU for binary problems. The three neighborhoods are
386 based on the Hamming distance. The move operators used in
387 the three neighborhoods consider Hamming distances of 1,
388 2, and 3 (this consists on flipping the binary value of one, two,
389 or three positions at a time in the candidate binary vector).
390 In~\cite{luong2010large}, each thread is associated to a unique
391 solution in the neighborhood. The addressed issue is how to
392 efficiently map the different neighborhoods on the device memory,
393 more explicitly, how to calculate the memory index of each
394 solution associated to each CUDA thread's \textit{id}.
395 %For 1-Hamming neighborhoods, as there is exactly n solutions in the neighborhood, the mapping of this neighborhood to CUDA threads is obvious: the CPU host offloads to GPU exactly $n$ threads, and each thread id is associated to one index in the binary vector. In the case of 2-Hamming and 3-Hamming neighborhoods, each thread id should be mapped respectively to two and three indexes in the candidate vector.
396 The three neighborhoods are implemented and experimented on the Permuted Perceptron Problem (PPP) using a tabu search\index{metaheuristics!tabu search} algorithm (TS). Accelerations from $9.9 \times$ to $18.5 \times$ are obtained on different problem sizes.\\ % The experiments are performed on an Intel Xeon 8 cores 3GHz coupled with an NVIDIA GTX 280 card.\\
398 In the same context, Deevacq et al.~\cite{audreyANT}
399 proposed two parallelization strategies inspired by the multi-walk
400 parallelization strategy, of a 3-opt iterated local
401 search algorithm (ILS)
402 over a CPU/GPU architecture. In the first strategy, each Local
403 Search (LS) is associated to a unique CUDA thread and improves a
404 unique solution by generating its neighborhood. The second
405 strategy associates each solution to a CUDA block and the
406 neighborhood exploration is parallelized on the block threads. In
407 the first strategy, since several LS are executed on different
408 solutions on each Multi Processor (MP), the data structures should
409 be stored on the global memory, while the exploration of a single
410 solution at a time in the second strategy allows the use of the
411 shared memory to store the related data structures. The two
412 strategies have been experimented on standard benchmarks of
413 the Traveling Salesman Problem (TSP) with sizes varying from $100$ to $3038$ cities. The results indicate that increasing the number of solutions to be explored simultaneously improves the speedup in the two strategies, but at the same time it decreases final solution quality. The greater speedup factor reached by the second strategy is $6 \times$.\\
415 The same strategy is used by Luong et al.
416 in~\cite{luongMultiStart} to implement multistart parallel local
417 search algorithms (a special case of the
418 algorithmic-level parallel model where several homogeneous LS
419 algorithms are used). The multistart model is combined with
421 parallelism: several LS algorithms are managed by the CPU and the
422 neighborhood evaluation step of each algorithm is parallelized on
423 the GPU (each GPU thread is associated with one neighbor and
424 executes the same evaluation function kernel). The advantage of
425 such a model is that it allows a high occupancy of the GPU
426 threads. Nevertheless, memory management causes new issues due to the
427 quantity of data to store and to communicate between CPU and
428 GPU. A second proposition for implementing the same model on GPU
429 consists of implementing the whole LS processes on GPU with each
430 GPU thread being associated to a unique LS algorithm. This solves
431 the communication issue encountered in the first model. In
432 addition, a memory management strategy is proposed to improve the
434 algorithmic-level model: texture memory is used to avoid memory latency
435 due to uncoalesced memory accesses. The proposed approaches are
436 implemented on the quadratic assignment problem (QAP) using CUDA.
437 The acceleration rates obtained for the
438 algorithmic-level with usage of texture memory rise from $7.8\times$ to
439 $12\times$ (for different
440 QAP benchmark sizes). \\
442 Janiak et al.~\cite{Janiak_et_al_2008} implemented two
443 algorithms for TSP and the flow-shop scheduling problem (FSP).
444 These algorithms are based on a multistart tabu
445 search model. Both of the
446 algorithms exploit multicore CPU and GPU. A full parallelization
447 on GPU is adopted using shader libraries where each thread is
448 mapped with one tabu search.
449 However, even though their experiments report that the use of GPU
450 speedups the serial execution almost $16 \times$, the mapping of
451 one thread with one tabu search
452 requires a large number of local search algorithms to
453 cover the memory access latency. The same mapping policy is adopted by Zhu et al. in~\cite{zhu_et_al_2008} (one thread is associated to one local search) solving the quadratic assignment problem but using the CUDA toolkit instead of shader libraries.\\
455 Luong et al.~\cite{luong2012ppsn} proposed a GPU-based
456 implementation of hybrid metaheuristics on heterogeneous parallel
457 architectures (multicore CPU coupled to one GPU). The challenge
458 of using such a heterogeneous architecture is how to distribute
459 tasks between the CPU cores and the GPU in such a way to have
460 optimal performances. Among the three traditional parallel models
463 and algorithmic-level), the authors point out that the most convenient
464 model for the considered heterogeneous architecture is a hybrid
467 and algorithmic-level models. Several CPU threads execute several instances
468 of the same S-metaheuristic in parallel while the GPU device is
469 associated to one CPU core and used to accelerate the
470 neighborhood calculation of several S-metaheuristics at the same
472 In order to efficiently exploit the remaining CPU cores, a load-balancing heuristic is also proposed in order to decide on the number of additional
473 S-metaheuristics to launch on the remaining CPU cores relative to the efficiency of the GPU calculations. The proposed approach is applied to the QAP using several instances of the Fast Ant Colony Algorithm (FANT)~\cite{taillardFant}. \\
475 All the previously noted works exploit the same parallel models
476 used on CPUs based on the task parallelism. A different
477 implementation approach is proposed by Paul in~\cite{gerald2012}
478 to implement a simulated
479 annealing (SA) algorithm
480 for the QAP on GPUs. Indeed, the author used a preinitialized
481 matrix \emph{delta} in which the incremental evaluation of simple
482 swap moves are calculated and stored relative to the initial
483 permutation $p$. For the GPU implementation, the authors used the
484 parallel implementation of neighborhood exploration. The
485 time-consuming tasks in the SA-matrix are identified by the
486 authors as updating the matrix and passing through it to select
487 the next accepted move. To initialize the delta matrix, several
488 threads from different blocks explore different segments of the
489 matrix (different moves) at the same time. In order to select the
490 next accepted swap, several threads are also used. Starting from
491 the last move, a group of threads explores different subsets of
492 the delta matrix. The shared memory is used to preload all the
493 necessary elements for a given group of threads responsible for
494 the updating of the delta matrix. The main difference in this work
495 compared to the previous works resides in the introduction of a
496 data parallelism using the precalculated delta matrix. The use of
497 this matrix allows the increase in the number of threads
498 involved in the evaluation of a single move. Experimentations are
499 done on standard QAP instances from the
500 QAPLIB~\cite{burkard1991qaplib}. Speedups up to $10 \times$ are
501 achieved by the GPU implementation compared
502 to the same sequential implementation on CPU using SA-matrix.\\
504 \subsection{Implementing population-based metaheuristics on GPUs}
506 State-of-the-art works dealing with the implementation of
507 p-metaheuristics on GPUs generally rely on parallel models and
508 research efforts done for parallelizing different classes of
509 p-metaheuristics over different types of parallel architectures:
510 supercomputers, clusters, and computational grids. Three main
511 classes of p-metaheuristics are considered in this section:
513 algorithms (EAs), ant
516 swarm optimization\index{metaheuristics!particle swarm
519 \subsubsection*{Evolutionary algorithms}
521 Traditional parallel models for EAs are classified in three main
522 classes: coarse-grain models using several subpopulations
523 (islands), master-slave models used for the parallelization of CPU
524 intensive steps (evaluation and transformation), and cellular
525 models based on the use of one population disposed (generally) on
528 The three traditional models have been implemented on GPUs by several researchers for different
529 optimization problems. The main chalenges to be raised when implementing the traditional models on GPUs concern (1) the saturation of the GPU in order to cover memory latency by calculations, and (2) efficent usage of the hierarchical GPU memories.\\
531 In~\cite{kannan}, Kannan and Ganji present a CUDA implementation
532 of the drug discovery application Autodock (molecular docking
533 application). Autodock uses a genetic
534 algorithm to find optimal
535 docking positions of a ligand to a protein. The most
536 time-consuming task in Autodock is the fitness function
537 evaluation. The fitness function used for a docking problem
538 consists of calculating the energy of the ligand-protein complex
539 (sum of intermolecular energies). The authors explore two
540 different approaches to evaluate the fitness function on GPU. In
541 the first approach, each GPU thread calculates the energy function
542 of a single individual. This approach requires the use of
543 large-sized populations to saturate the GPU (thousands of
544 individuals). In the second approach each individual is associated
545 with one thread block. The evaluation of the energy function is
546 performed by the threads of the same block. This allows the use of
547 medium population sizes (hundreds of individuals) and the
548 acceleration of a single fitness evaluation. Another great
549 advantage of the per block approach resides in the use of shared
550 memory instead of global memory to store all the information
551 related to each individual. The obtained speedups range from $10
552 \times$ to $47 \times$ for population sizes
553 ranging from $50$ to $10000$. \\
555 Maitre et al.~\cite{maitre2009} also exploited the
556 master-slave parallelism of EAs on GPUs using EASEA. EASEA is a
557 C-like metalanguage for easy development of EAs. The user writes a
558 description of the problem-specific components (fitness function,
559 problem representation, etc) in EASEA. The code is then compiled
560 to obtain a ready-to-use evolutionary
562 compiler uses genetic
564 Libraries to produce C++ or JAVA written EA codes.
565 In~\cite{maitre2009}, the authors proposed an extension of EASEA
566 to produce CUDA code from the EASEA files. This extension has been
567 used to generate a master-slave parallel EA in which the
568 sequential algorithm is maintained on CPU and the population is
569 sent to GPU for evaluation. Two problems have been considered
570 during the experiments: a benchmark mathematical function and a
571 real problem (molecular structure prediction). In order to
572 maximize the GPU occupation, very large populations are used (from
573 $2000$ to $20000$). Even though transferring such large
574 populations from the CPU to the GPU device memory at every generation is very costly, the authors report important speedups on the two problems on a GTX260 card: $105 \times$ is reported for the benchmark function while for the real problem the reported speedup is $60 \times$. This may be best explained by the complexity of the fitness functions. Nevertheless, there is no indication in the paper about the memory management of the populations on GPU.\\
576 The master-slave model is efficient when the fitness function is
577 highly time intensive. Nevertheless, it requires the use of
578 large-sized populations in order to saturate the GPU unless the
579 per-block is used (one individual perblock) when the acceleration
580 of the fitness function itself is possible. The use of many
581 subpopulations of medium sizes is another way to obtain a maximum
582 occupation of the GPU. This is coarse-grained parallelism (island
585 The coarse-grained model is used by Pospichal et al.
586 in~\cite{pospichal10} to implement a parallel genetic
587 algorithm on GPU. In this
588 work the entire genetic
589 algorithm is implemented
590 on GPU. This choice is motivated by the overhead engendered by the
591 CPU/GPU communication
592 when only population evaluation is performed on GPU. Each
593 population island is mapped with a CUDA thread block and each
594 thread is responsible for a unique individual. Subpopulations are
595 stored on shared memory of each block. Nevertheless, because
596 interblock communications are not possible on the CUDA
597 architecture, the islands evolve independently in each block, and
598 migrations are performed
599 asynchronously through the global memory. That is, after a given number of generations, selected individuals for migration from each island are copied to the GPU global memory part of the neighbor island and then to its shared memory to replace the worst individuals in the local population. The experiments are performed on three benchmark mathematical functions. During these experiments, the island sizes are varied from $2$ to $256$ individuals and island numbers from $1$ to $1024$. The maximum performance is achieved for high number of islands and increasing population sizes.\\
601 The same strategy is also adopted by Tsutsui and Fujimoto
602 in~\cite{tsutsuiGAQAP} to implement a coarse-grained genetic
603 algorithm on GPU to solve
604 the QAP. Initially, several subpopulations are created on CPU and
605 transferred to the global memory. The subpopulations are organized
606 in the global memory into blocks of $8$ individuals in such a way
607 to allow coalesced memory access by the threads of the same thread
608 block. Each sub-population is allocated to a single thread block
609 in the GPU and transfered to the shared memory to start evolution.
610 Population evaluation and transformation are done in parallel by
611 the different threads of a given block. Migration is also done
612 through the global memory. Experiments are performed on standard
613 QAP benchmarks from the QAPLIB~\cite{burkard1991qaplib}. The GPU
614 implementation reached speedups of $2.9\times$ to $12.6 \times$
615 compared to a single core implementation of a coarse-grained
616 genetic algorithm on a
617 Intel i7 processor.\\
619 Nowotniak and Kucharski~\cite{nowotniak} proposed a GPU-based
620 implementation of a Quantum Inspired Genetic Algorithm (QIGA). The
621 parallel model used is a hierarchical model based on two levels: each
622 thread in a block transforms a unique individual and a different
623 population is assigned to each block. The algorithm is run
624 entirely on GPU. A different instance of the QIGA is run on each
625 block and the computations have been shared between 8 GPUs. This
626 approach is very convenient to speed up the experimental process
627 on metaheuristics that require several independent runs of the
628 same algorithm in order to assess statistical efficiency. The
629 populations are stored in the shared memory, the data matrix used
630 for fitness evaluation is placed in read only constant memory, and
631 finally seeds for random numbers generated on the GPU are stored
632 in the global memory.\\
634 In coarse-grained parallelism, the use of the per-block approach
635 to implement the islands (one subpopulation per thread block) is
636 almost natural and it allows the use of shared memory to store
637 the subpopulations. Fine-grained parallel models for EAs (cellular
638 EAs) have been used by many authors to implement parallel EAs on
639 GPUs. Indeed, the fine-grained parallelism of EAs fits perfectly
640 to the SIMD architecture of the GPU. \\
642 Pinel et al. in~\cite{pinel2012JPDC} developed a highly
643 parallel synchronous cellular genetic
644 algorithm (CGA), called
645 GraphCell, to solve the independent task scheduling problem on GPU
646 architectures. In CGAs, the population is arranged into a
647 two-dimensional toroidal grid where only neighboring solutions are
648 allowed to interact with each other during the recombination
649 step. In GraphCell, two recombination operators for CGA are
650 especially designed to run efficiently on GPU. Indeed, instead of
651 assigning a single thread to each solution of the population to
652 perform the recombination, in GraphCell, a single thread is
653 assigned to each task of each solution. Offsprings are created by
654 independently modifying the assignment of some tasks in the
655 current solution. Mainly, each thread chooses one neighboring
656 solution in the grid as second parent using different selection
657 strategies and assigns one task of the solution (first parent)
658 to the machine for which the same task is assigned in the second
659 parent. This way, the number of threads used for the
660 recombination step is equal to population size $\times$ size of
661 the solution (number of tasks). This leads to a high number of
662 threads used to accelerate the recombination operators especially
663 when dealing with large instances of the problem. In addition to
664 the recombination operators, the rest of the CGA steps are also
665 parallelized on GPU (fitness evaluation, mutation, and
668 A similar work is proposed by Vidal and Alba in~\cite{albaCGAGPU}
669 where a CGA using a toroidal grid is completely implemented on
670 GPU. A direct mapping between the population and the GPU threads
671 is done. At each step, several threads execute the same operations
672 on each individual independently: initialization, computing the
673 neighborhood, selection, crossover, mutation, and evaluation. A
674 synchronization is done for all threads to perform the replacement
675 stage and form the next generation. All the data of the algorithm
676 is placed on the global memory. Several experiments have been
677 performed on a set of standard benchmark functions with different
678 grid sizes ranging from $32^2$ to $512^2$. The speedups reached by
679 the GPU implementation against the CPU version range from
680 $5\times$ to $24\times$ and increase as the size of the population
681 increases. However, the CPU implementation runs faster than the GPU
682 version for all the tested benchmarks when the size of the
683 population is set to $32^2$. When the size of the population is
684 too small, there are not enough computations to cover the overhead
685 created by the call of kernel functions, CPU/GPU
687 synchronization, and access to global memory. Finally, an
688 interesting review on GPU parallel computation in bio-inspired
689 algorithms is proposed by Arenas et al.
690 in~\cite{arenas2011}. \\
692 \subsubsection*{Ant colony optimization}
694 Ant colony optimization
696 p-metaheuristic subject to parallelization on GPUs.
697 State-of-the-art works on parallelizing
699 accelerating the tour construction step performed by each ant by
700 taking a task-based parallelism approach, with pheromone
701 deposition on the CPU.\\
703 In~\cite{cecilia}, Cecilia et al. present a GPU-based
706 the two steps (tour construction and pheromone update) are
707 parallelized on the GPU. A data parallelism approach is used to
708 enhance the performance of the tour construction step. The
709 authors use two categories of artificial ants: queen ants
710 associated with CUDA thread-blocks and worker ants associated with
711 CUDA threads. A queen ant represents a simulated ant and worker
712 ants collaborate with the queen ant to accelerate the decision
713 about the next city to visit. The tour construction step of each
714 queen ant is accelerated. Each worker ant maintains a history of
715 the search in a tabu list containing a chronological ordering of
716 the already visited cities. This memory is used to determine the
717 feasible neighborhoods. After all queen ants have constructed
718 their tours, the pheromone trails are updated. For pheromone
719 update, several GPU techniques are also used to increase the data
720 bandwidth of the application mainly by the use of precalculated
721 matrices that are easily updated by several threads (one thread
722 per matrix entry). The achieved speedups are $21 \times$ for tour
723 construction and $20 \times$ for pheromone updates.\\
725 In another work, Tsutsui and Fujimoto~\cite{tsutsui} propose a
726 hybrid algorithm combining
728 and Tabu Search (TS) implemented
729 on GPU to solve the QAP. A solution of QAP is represented as a
730 permutation of ${1,2,..,n}$ with $n$ being the size of the
731 problem. The TS algorithm is based on the 2-opt neighborhood
732 (swapping of two elements $(i,j)$ in the permutation). The authors
733 point out that the move cost of each neighbor depends on the
734 couple $(i,j)$. Two groups of moves are formed according to the
735 move cost. In order to avoid thread divergence\index{GPU!thread divergence} within the same warp, the
736 neighborhood evaluation is parallelized in such a way to assign
737 only moves of the same cost to each thread warp. This strategy is
738 called MATA for Move-cost Adjusted Thread Assignment. Concerning
739 the memory management\index{GPU!memory management}, all
740 the data of the ACO\index{metaheuristics!ant colony optimization}
741 (population, pheromone matrix), QAP matrices, and tabu list are
742 placed on the global memory of the GPU. Shared memory is used only
743 for working data common to all threads in a given block.
745 steps of the hybrid algorithm
748 pheromone update, construct solutions, applying TS) are
749 implemented as kernel functions on the GPU. The GPU/CPU
750 communications are only used to transfer the best-so-far solution
751 in order to verify termination conditions. The experiments
752 performed on standard QAP benchmarks showed that the GPU
753 implementation using MATA obtained a speedup of $19 \times$
754 compared to the CPU implementation, compared with a speedup of
756 when MATA is not used. \\
758 \subsubsection*{Particle swarm optimization}
759 In~\cite{zhou2009} Zhou and Tan propose a full GPU implementation
760 of a standard PSO algorithm. All the data is stored in global
761 memory (velocities, positions, swarm population, etc). Only
762 working data is copied to shared memory by each thread. The four
763 steps of the PSO have been parallelized on GPU: fitness evaluation
764 of the swarm, update of local best and global best of each
765 particle, and update of velocity and position of each particle.
766 The same strategy is used to parallelize the first and last
767 steps: the evaluation of fitness functions is performed in
768 parallel for each dimension by several threads. It is the case
769 for the update of position and velocities of each particle: one
770 dimension at a time is updated for the whole swarm by several
771 threads. Random numbers needed for updating the velocities and
772 positions for the whole PSO processes are generated on CPU at the
773 starting of the algorithm and transferred to the GPU global
774 memory. For the steps 2 and 3 (update of local best and global
775 best of each particle), the GPU threads are mapped to the \emph{N}
776 particles of the swarm one to one. Experiments done on four
777 benchmark functions show speedups ranging from $3.7 \times$ to
778 $9.0 \times$ for swarm sizes
779 ranging from $400$ to $2800$.\\
781 \subsection{Synthesis of the implementation strategies}
782 \label{ch8:sec:synthesis} After reviewing some works dealing with
783 GPU-based implementation of metaheuristics, in this section we
784 will try to come up with a classification of the different
785 strategies used in the literature for the implementation of
786 metaheuristics on GPUs. One may distinguish between the parallel
787 models adopted in each metaheuristic (design level) and the way
788 they are exploited on GPU architectures (implementation level).
789 Indeed, even though the parallelization models used in most works
790 for GPUs are derived from the traditional parallel models of each
791 metaheuristic (on CPU), their implementation could take a
792 different way and sometimes it may result in new parallel models
793 customized for GPUs.\\
795 Traditional parallel models for metaheuristics are based on an
796 intuitive task parallelism: the independent tasks of the
797 algorithms are simply parallelized. For example, in the case of
798 s-metaheuristics, the evaluation of large neighborhoods could be
799 done in parallel since there is no synchronization at this step of
800 the algorithm. That is the case of EAs when it comes to applying
801 the evaluation step. Nevertheless, because of the particularity of
802 the GPU architecture, some authors have used new implementation
803 techniques to enhance the data parallelism in the sequential
804 algorithms in order to increase the data throughput of the
807 From this observation we propose the following classification
808 based on 2 levels: design level and implementation level as
809 illustrated in Figure~\ref{ch8:fig:classification}. The design
810 level regroups the three classes of parallel models used in
814 algorithmic-level) with examples for s-metaheuristics, EAs,
816 classification is principally built from the reviewed
817 state-of-the-art works in the previous section. The
818 implementation level refers to the way these parallel design
819 models are implemented on GPU. This classification focuses only on
820 the mapping strategies between the GPU threads and the
821 parallelized tasks (neighborhood evaluation, solution
822 construction, and so on). The different implementation strategies
823 are explained in the following sections.\\
826 \centerline{\includegraphics[width=1\textwidth]{Chapters/chapter9/figures/classification.pdf}}
827 \caption{A two level classification of state-of-the-art GPU-based parallel metaheuristics.}
828 %[Comparison of number of cores in a CPU and in a GPU]
829 \label{ch8:fig:classification}
833 \subsubsection*{GPU thread mapping for solution-level parallelism}
834 \index{GPU!thread mapping} Parallel
835 models at solution level consist of parallelizing a time intensive
836 atomic task of the algorithm. Generally, it consists of the
837 fitness evaluation~\cite{kannan}. Nevertheless, crossover
838 operators have also been parallelized by some
839 authors~\cite{pinel2012JPDC}. These kinds of models are not always
840 possible as they are problem-dependent. The GPU implementation of
841 solution-level models uses the perblock mapping: each solution is
842 associated to a block of threads. A second level of parallelism is
843 used inside each block to parallelize the fitness evaluation of a
844 single solution. This kind of mapping allows the use of shared
845 memory to store the data structures of the solution and does not
846 require the use of very large neighborhoods or populations.\\
847 %data parallelism in SA-matrix to parallelize
849 \subsubsection*{GPU thread mapping for iteration-level parallelism}
850 \index{GPU!thread mapping}
851 Iteration-level parallelism consists of parallelizing the tasks
852 performed independently on different solutions. Different mapping
853 strategies are adopted in the reviewed works to implement these
856 In Figure \ref{ch8:fig:classification}, the first example of
858 parallelism is the parallel evaluation of neighborhoods in
859 s-metaheuristics. In most of the reviewed works, a per-thread
860 mapping approach is used: each solution of the neighborhood is
861 mapped to a unique thread in the GPU for
862 evaluation~\cite{luong2010large, audreyANT}. The same mapping is
863 used for master-slave parallel EAs to accelerate the evaluation of
864 large populations. This kind of mapping is only efficient for very
865 large neighborhoods and very large populations (to saturate the
866 GPU). Many authors have pointed out that the use of such large
867 populations (or neighborhoods) may lead to an overhead due to the
868 communication costs between the CPU and the GPU (if the sequential
869 part of the algorithm is placed on CPU). In order to circumvent
870 this issue, many authors have implemented the entire algorithm on
871 GPU~\cite{pospichal10}. On the other hand, as several solutions
872 may be mapped with the same thread block in the GPU, shared memory
873 could not be used and all the data should be placed on global
875 %pheromone update data parallelism matrices
877 \subsubsection*{GPU thread mapping for algorithmic-level parallelism}
878 Algorithmic-level parallelism consists of launching several self-contained algorithms in parallel. In the previously reviewed works two algorithmic-level models have been used: the multistart model and the island model (parallel EAs).\\
880 The implementation of the multistart model is based on two
881 different mapping strategies~\cite{luongMultiStart, audreyANT}:
882 (1) each Local Search (LS) is associated to a unique thread and
883 (2) each solution (from multiple neighborhoods) is associated to a
884 unique thread. The first mapping strategy (one thread per LS)
885 presents a big drawback: the number of LS to use should be very
886 large to saturate the GPU and cover the memory access latency. On
887 the other hand, the evaluation of the neighborhood of a single LS
888 by one thread is time intensive. Furthermore, shared memory could
889 not be used to store the huge data generated by the different
890 neighborhoods. In the second mapping strategy, the LS algorithms
891 are placed on CPU and the neighborhood evaluation of each LS is
892 parallelized on GPU using per-thread mapping strategy (one thread
893 per solution). This consists of a hierarchical parallel model
894 combining algorithmic-level parallelism (multistart) with
896 parallelism (master-worker).\\
898 In the island model, the same mapping is used in all the reviewed
899 works~\cite{tsutsuiGAQAP, nowotniak, maitre2009}: each
900 subpopulation (island) is associated to one thread block in the
901 GPU. A second level of parallelism is used inside the block to
902 parallelize the evaluation step of the local population.
903 Migrations are always performed through the global memory as
904 interblock communications are impossible in CUDA. The first
905 advantage of this hierarchical implementation is that it allows
906 the occupation of a large number of threads even for medium
907 population sizes. The second advantage consists of using shared
908 memory to store subpopulations to benefit from the low latency of
911 \section{Frameworks for metaheuristics on GPUs}
912 \label{ch8:sec:frameworks}
914 After the first pioneering works of metaheuristics on graphics
915 processing units, the next challenge is to provide easy-to-use
916 frameworks and libraries for rapid development of metaheuristics
917 on GPUs. Although the works on this subject are not yet mature and
918 do not cover the main metaheuristic algorithms, we will present
919 the only three works to our knowledge, which propose open source
920 frameworks for developing
921 metaheuristics on GPUs.\\
923 The three frameworks reviewed in this section are
924 PUGACE\index{GPU-based frameworks!PUGACE}~\cite{pugace},
925 ParadisEO-MO-GPU\index{GPU-based
926 frameworks!ParadisEO-MO-GPU}~\cite{paradiseoGPU}, and
927 libCUDAOptimize\index{GPU-based
928 frameworks!libCUDAOptimize}~\cite{libcuda}. PUGACE\index{GPU-based
929 frameworks!PUGACE} is a framework for implementing EAs on GPUs.
930 ParadisEO-MO-GPU is an extension of the framework
931 ParadisEO~\cite{paradiseo} for implementing s-metaheuristics on
932 GPU. Finally, libCUDAOptimize\index{GPU-based
933 frameworks!libCUDAOptimize} is a library intended for the
934 implementation of p-metaheuristics on GPU. The three frameworks
935 are presented in more detail in the following.
937 \subsection{PUGACE: framework for implementing evolutionary computation on GPUs}
938 PUGACE is a generic framework
939 for easy implementation of cellular evolutionary algorithms on
940 GPUs implemented using C and CUDA. It is based on the frameworks
941 MALLBA and JCell (a framework for cellular algorithms). The
942 authors justified the choice of cellular evolutionary
943 algorithm by the good
944 feedback found in the literature concerning its efficient
945 implementation on GPUs compared to other parallel models for EAs
946 (island, master-slave). The main standard evolutionary operators
947 are already implemented in PUGACE: different selection strategies, standard
948 crossover, and mutation operators (PMX, swap, 2-exchange,
949 etc.). Different problem encoding is also supported. The framework
950 is organized as a set of modules in which the different
951 functionalities are separated as much as possible in order to
952 facilitate the extension of the framework. All
953 the functions and procedures that execute on GPU are implemented in the same file kernel.cu. \\
955 The implementation strategy adopted on the GPU is as follows.
956 Population initialization is done on the CPU side and the
957 population is transferred to the GPU. On the GPU side, each
958 individual is associated to a unique CUDA thread. The function
959 evaluation and mutation are done on the GPU while selection and
960 replacement are maintained on the CPU. In order to avoid thread
961 divergence\index{GPU!thread divergence} appearing in
962 the same CUDA thread block at the crossover step (because of the
963 probability of application which may give different results from
964 one thread to the other), the decision of whether to apply a
965 crossover is taken at the block level and applied to all the
966 individuals within the block. It is the decision on the choose
967 of the cutting point for the crossover.\\
969 The framework is validated on standard benchmarks of the QAP. Speedups of $15.44 \times$ to $18.41
970 \times$ are achieved compared to a CPU implementation of a cEA
971 using population sizes rising from 512 to 1024 and from 1024 to
974 % The authors noticed an
975 % improvement on the speedup when the size of the population increases. This is best explained
976 % by the linear increasing of execution time on CPU side while a sublinear increasing is
977 % occuring on GPU side (due to a better exploitation of the GPU threads). Neverethless,
978 % dispite the positive impact of large populations on th parallel efficieny of the implementation
979 % on GPU, a study on the impact on the final solution quality when using such large populations
980 % is to be done. Indeed, too large populations on EAs may lead to high diversification of the
981 % search and poor intensification capabilities, which leads to poor results.
982 % For the moment, plateforms whith more than one CPU and one GPU are not yet supported in PUGACE\index{GPU-based frameworks!PUGACE}.
983 % No details on how to get PUGACE\index{GPU-based frameworks!PUGACE} code source are given in the paper.
985 \subsection{ParadisEO-MO-GPU}
987 Melab et al.~\cite{paradiseoGPU} developed a reusable
988 framework ParadisEO-MO-GPU\index{GPU-based
989 frameworks!ParadisEO-MO-GPU} for parallel local search
990 metaheuristics (s-metaheuristics) on GPUs. It focuses on the
992 parallel model of s-metaheuristics which consists of exploring in
993 parallel the neighborhood of a problem solution on GPU. The
994 framework, implemented using C++ and CUDA, is an extension of the
995 ParadisEO~\cite{paradiseo} framework previously developed by the
996 same team for parallel and distributed metaheuristics on both
997 dedicated parallel hardware platforms and computational grids. The
998 objective of this framework is to facilitate the development of
999 GPU-based metaheuristics providing a transparent use of GPUs to
1000 users who are unfamiliar with advanced features of all
1001 parallelization techniques and deployment on GPUs. The framework
1002 allows one to efficiently manage the hierarchical organization of
1003 the memories (latencies and sizes) of the GPU and its
1004 communication with the CPU as well as the minimizing of the user
1005 involvement in its management.\\
1008 \centerline{\includegraphics[width=0.8\textwidth]{Chapters/chapter9/figures/paradiseoGPU.pdf}}
1009 \caption{The skeleton of the ParadisEO-MO-GPU\index{GPU-based
1010 frameworks!ParadisEO-MO-GPU}.}
1011 %[Comparison of number of cores in a CPU and in a GPU]
1012 \label{ch1:fig:paradiseoGPU}
1015 The framework is based on a master-worker model where the master
1016 is the CPU and the workers are threads executed by the processing
1017 cores of the GPU. The CPU executes the serial part of the
1018 metaheuristic and sends only the current solution to the GPU to
1019 minimize the transfer cost. The GPU, on its side, generates the
1020 neighboring of the received solution and evaluates them at each
1021 iteration. All the threads execute the same kernel and according
1022 to a static mapping table between the threads and the neighbors
1023 where each thread is associated with exactly one neighbor
1024 evaluation. After all the neighborhood is generated and evaluated,
1025 it is sent back to the CPU which selects the best solution (See
1026 Figure~\ref{ch1:fig:paradiseoGPU}).
1028 \subsection{libCUDAOptimize: an open source library of GPU-based metaheuristics}
1029 \index{GPU-based frameworks!libCudaOptimize}
1030 LibCUDAOptimize~\cite{libcuda} is a C++/CUDA open source library
1031 for the design and implementation of metaheuristics on GPUs. Until
1032 now, the metaheuristics supported by LibCUDAOptimize are: scatter
1033 search, differential evolution, and particle swarm
1034 optimization\index{metaheuristics!particle swarm optimization}.
1035 Nevertheless, the library is designed in such a way to allow
1036 further extensions for other metaheuristics and it is still in
1037 development phase by the authors. The parallelization strategy
1038 adopted on GPU is principally based on fitness evaluation. The
1039 sequential structure of the optimization algorithms is maintained
1042 \section{Case study: accelerating large neighborhood LS method on GPUs for solving the Q3AP}
1043 \label{ch8:sec:case}
1045 In this case study, a large neighborhood GPU-based local
1046 search\index{GPU!based local search}
1047 method for solving the Quadratic 3-dimensional Assignment Problem
1048 (Q3AP) will be presented. The local search method is an Iterated
1049 Local Search\index{metaheuristics!iterated local search}
1050 (ILS)~\cite{stutzle2006ILSforQAP} using an embedded
1051 TS algorithm. The ILS principle
1052 consists of executing iteratively the embedded local search, each
1053 iteration which starts from a disrupted local optima reached by
1054 the previous local search process. The disruption heuristic is a
1055 performance parameter of an ILS algorithm and should be
1056 judiciously defined. A template of an
1057 ILS algorithm is given by the Algorithm~\ref{ch8:ils_algorithm_template}.\\
1059 \begin{algorithm}[H]
1062 %\KwResult{how to write algorithm with \LaTeX2e }
1063 $s_{0}$=GenerateInitSol()\;
1064 $s^*$=LocalSearch($s_0$)\;
1065 \Repeat{a maximum number of iterations is reached}{
1066 $s'$=Perturbate($s^*,history$)\;
1067 $s^{*'}$=LocalSearch($s'$)\;
1068 $s^*$=AcceptationCriteria($s^*,s^{*'},history$)\;
1071 \caption{iterated local search template}
1072 \label{ch8:ils_algorithm_template}
1077 % % \label{ch8:ils_algorithm_template}
1078 % \caption{iterated local search\index{Metaheuristics!iterated local search} template}
1079 % \begin{algorithmic}[1]
1081 % \STATE $s_{0}$=\texttt{GenerateInitSol}();
1083 % \STATE $s^*$=\texttt{LocalSearch}($s_0$);
1086 % \STATE $s'$=\texttt{Perturbate}($s^*,history$);
1087 % \STATE $s^{*'}$=\texttt{LocalSearch}($s'$);
1088 % \STATE $s^*$=\texttt{AcceptationCriteria}($s^*,s^{*'},history$);
1090 % \UNTIL {a maximum number of iterations is reached}
1096 \subsection{The quadratic 3-dimensional assignment problem}
1098 The Q3AP is an extension of the well-known NP-hard QAP. The latter
1099 is one of the most studied problems by the combinatorial
1100 optimization community due to its wide range of practical
1101 applications (site planning, schedule problem, computer-aided
1102 design, etc.) and to its computational challenges since it is
1103 considered as one of the most computationally difficult
1104 optimization problems.\\
1106 The Q3AP was first introduced by William P. Pierskalla in
1107 1967~\cite{Pierskalla1967Q3AP} and, unlike the QAP, the Q3AP is a
1108 less studied problem. Indeed, the Q3AP was revisited only this
1109 past year and has recently been used to model some advanced
1110 assignment problems such as the symbol-mapping problem encountered
1111 in wireless communication systems and described
1112 in~\cite{Hahn2008q3ap}. The Q3AP optimization problem can be
1113 mathematically expressed as follows:
1115 % debut -----------------------------------------------------------------------------
1116 \begin{equation*} min \left \{
1117 \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\sum_{l=0}^{n-1}\sum_{k=0}^{n-1}\sum_{s=0}^{n-1}\sum_{r=0}^{n-1}c_{ijlksr}x_{ijl}x_{ksr}
1118 \quad+ \right . \end{equation*}
1120 \begin{equation} \label{Q3APFormulation} \left .
1121 \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\sum_{l=0}^{n-1}b_{ijl}x_{ijl}
1122 \right \} \end{equation}
1127 X=(x_{ijl})\in I \cap J \cap L, \label{Q3APConstraints1}\\
1128 x_{ijl}\in \{0,1\}, \quad i,j,l=0,1,...,n-1 \label{Q3APConstraints2}
1129 \end{eqnarray} $I$, $J$, and $L$ sets are defined as follows:
1130 \begin{displaymath} I=\lbrace X=(x_{ijl}):
1131 \sum_{j=0}^{n-1}\sum_{l=0}^{n-1}x_{ijl}=1, \quad
1132 i=0,1,...,n-1\rbrace \mathrm{;} \end{displaymath}
1133 \begin{displaymath} J=\lbrace X=(x_{ijl}):
1134 \sum_{i=0}^{n-1}\sum_{l=0}^{n-1}x_{ijl}=1, \quad
1135 j=0,1,...,n-1\rbrace \mathrm{;} \end{displaymath}
1136 \begin{displaymath} L=\lbrace X=(x_{ijl}):
1137 \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}x_{ijl}=1, \quad
1138 l=0,1,...,n-1\rbrace \end{displaymath}
1139 % fin------------------------------------------------------------------------------
1141 Other equivalent formulations of the Q3AP can be found in the literature. A particularly useful one is the
1142 \textit{permutation-based formulation} wherein the Q3AP can be expressed as follows:
1145 min\left\{ f(\pi_1,\pi_2) =
1146 \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i\pi_{1(i)}\pi_{2(i)}j\pi_{1(j)}\pi_{2(j)}} \quad
1151 \left. \sum_{i=0}^{n-1}b_{i\pi_{1(i)}\pi_{2(i)}}\right
1152 \} \label{Q3APPermutationBasedFormulation}
1154 where $\pi_1$ and $\pi_2$ are permutations over the set
1155 $\left\lbrace 0,1,\ldots,n-1\right\rbrace$. According to this
1156 formulation, minimizing the Q3AP consists of finding a double
1157 permutation $(\pi_1,\pi_2)$ which minimizes the objective function
1158 (\ref{Q3APPermutationBasedFormulation}).\\
1160 The Q3AP is proven to be an NP-hard problem since it is an
1161 extension of the quadratic assignment problem and of the axial
1162 3-dimensional assignment problem which are both NP-hard problems. It is particularly computationally difficult since the
1163 number of feasible solutions of an instance of size $n$ is $n! \times n!$.
1165 \subsection{Iterated tabu search algorithm for the Q3AP}
1167 To tackle large-sized instances of the Q3AP and speed up the
1168 search process, a parallel ILS algorithm has been designed. The local search embedded in
1169 the ILS is a TS. A TS
1170 procedure~\cite{Glover1989TS} starts from an initial feasible
1171 solution and tries, at each step, to move to a neighboring
1172 solution that minimizes the fitness (for a minimization case). If
1173 no such move exists, the neighbor solution that less degrades the
1174 fitness is chosen as a next move. This enables the TS process to
1175 escape local optima. However, this strategy may generate cycles,
1176 i.e., previous moves can be selected again. To avoid cycles, the
1177 TS manages a short-term memory that contains the moves that have
1178 been recently performed. A TS
1179 template is given by Algorithm \ref{TS_pseudo_code}.\\
1182 % \label{TS_pseudo_code}
1183 % \caption{Tabu search template}
1184 % \begin{algorithmic}[1]
1185 % \STATE GenerateInitSol($s_0$);
1186 % \STATE TabuList=$\phi$;
1189 % \STATE $m(t)$ = SelectBestMove($s(t)$);
1190 % \STATE $s(t+1)$ = ApplyMove($m(t),s(t)$);
1191 % \STATE TabuList = TabuList $\bigcup \{m(t)\}$;
1192 % \STATE $t = t + 1$;
1193 % \UNTIL{a maximum number of iterations is reached}
1198 \begin{algorithm}[H]
1201 %\KwResult{how to write algorithm with \LaTeX2e }
1202 GenerateInitSol($s_0$)\;
1205 \Repeat{a maximum number of iterations is reached}{
1206 $m(t)$ = SelectBestMove($s(t)$)\;
1207 $s(t+1)$ = ApplyMove($m(t),s(t)$)\;
1208 TabuList = TabuList $\bigcup \{m(t)\}$\;
1212 \caption{tabu search template}
1213 \label{TS_pseudo_code}
1217 \subsection{Neighborhood functions for the Q3AP}
1219 The neighborhood function is a crucial parameter in any local
1220 search algorithm. Indeed, if the neighborhood function is not
1221 adequate to the problem and/or does not consider the targeted
1222 computing framework, any local search algorithm will fail to reach
1223 good quality solutions of the search space.\\
1225 Regarding the Q3AP, many neighborhood structures can be
1226 considered. A basic neighborhood was proposed and investigated in
1227 different works of the
1228 literature~\cite{Hahn2008q3ap,DBLP:journals/ijfcs/LoukilMMTB12}
1229 and consists of the set of all solutions (double-permutations)
1230 generated from the current one by an exchange of two positions in
1231 either the first ($\pi_1$) or the second ($\pi_2$) permutation.
1232 This neighborhood that we denote by $N_b$ can be formalized as
1237 N_b(\pi_1,\pi_2)=\{&(\pi_1^{'},\pi_2^{'}): \pi_1^{'}[k]=\pi_1[l], \pi_1^{'}[l]=\pi_1[k]\\
1238 & 0\leq k \neq l < n;\\
1239 & \pi_1^{'}[i]=\pi_1[i], \quad 0\leq i\neq k,l < n;\\
1240 & \pi_2^{'}[j]=\pi_2[j], \quad 0\leq j<n\\ \} \\
1242 \{& (\pi_1^{'},\pi_2^{'}): \pi_2^{'}[k]=\pi_2[l], \pi_2^{'}[l]=\pi_2[k]\\
1243 & 0\leq k \neq l < n;\\
1244 & \pi_2^{'}[i]=\pi_2[i], \quad 0\leq i < n;\\
1245 & \pi_1^{'}[j]=\pi_1[j], \quad 0\leq j\neq k,l < n\\
1248 where $(\pi_1,\pi_2)$ is the current solution and $n$ is the size
1249 of the Q3AP instance. The size $|N_b|$ of such neighborhood is
1253 |N_b| = n \times (n-1)
1257 In our GPU-based implementations, a large-sized neighborhood
1258 structure has been used for experimentation. In fact, theoretical
1259 and experimental studies have shown that the use of large
1260 neighborhood structures may improve the effectiveness of LS
1261 algorithms~\cite{Ahuja:2007:VLN:1528422.1528438}. However, for
1262 such a neighborhood, the generation/evaluation step of an LS
1263 becomes a time-consuming task and may dramatically increase the
1264 computational time of the LS process. This
1265 justifies the use of intensive data-parallelism provided by GPUs where all neighboring solutions may be concurrently evaluated. \\
1267 The considered large-sized neighborhood consists of swapping two
1268 positions in both permutations $\pi_1$ and $\pi_2$. This
1269 neighborhood structure, $N(\pi_1,\pi_2)$, can be formalized as
1272 \begin{equation} \begin{array}{rl}
1273 N(\pi_1,\pi_2)=\{&(\pi_1^{'},\pi_2^{'}): \pi_1^{'}[k]=\pi_1[l], \pi_1^{'}[l]=\pi_1[k],\\
1274 & \pi_2^{'}[r]=\pi_2[s], \pi_2^{'}[s]=\pi_2[r]\\
1275 & 0\leq k \neq l < n,0\leq r \neq s < n;\\
1276 & \pi_1^{'}[i]=\pi_1[i], \quad 0\leq i\neq k,l < n;\\
1277 & \pi_2^{'}[j]=\pi_2^{'}[j], \quad 0\leq j\neq r,s < n\\
1282 So, for a given Q3AP instance of size $n$, the size $|N|$ of the advanced neighborhood set can be expressed by the following formula:
1285 |N|=\left(\frac{n \times (n-1)}{2}\right )^2
1288 \subsection{Design and implementation of a GPU-based iterated tabu search algorithm for the Q3AP} \label{ch8:ITS-Q3APSection}
1290 The use of GPU devices to speed up the search process of local
1291 search methods is not a straightforward task. It requires one to
1292 consider, at the same time, the characteristics and underlying
1293 issues of the GPU architecture and the parallel models of LS
1294 methods. The main challenges that must be faced when designing a
1295 local search algorithm are the efficient distribution of the
1296 search process between the CPU and the GPU minimizing the data
1297 transfer between them, the hierarchical memory
1298 management\index{GPU!memory management} and the
1299 capacity constraints of GPU memories, and the thread
1300 synchronization. All these issues must be regarded when designing
1301 parallel LS models to allow
1302 solving of large scale optimization problems on GPU architectures.\\
1304 To go back to our problem (i.e., Q3AP), we propose in
1305 Algorithm~\ref{ch8:algoITS} an iterated tabu
1306 search on GPU (GPU-ITS). The
1307 parallel model is in agreement with the
1309 parallel model of LS methods presented in Section
1310 \ref{ch8:sec:paraMeta} (Fig. \ref{ch8:fig:paraMeta}). This
1311 algorithm can be seen as a cooperative model between the CPU and
1312 the GPU. Indeed, the GPU is used as a coprocessor in a synchronous
1313 manner. The resource-consuming part, i.e., the generation and
1314 evaluation kernel, is calculated by the GPU and the rest of the LS
1315 process is handled by the CPU. First of all, at initialization
1316 stage, memory allocations on GPU are made; the input matrices
1317 (distance and flow matrices) and the initial candidate solution of
1318 the Q3AP must be allocated (lines 4 and 5). Since GPUs require
1319 massive computations with predictable memory accesses, a structure
1320 has to be allocated for storing all the neighborhood fitnesses at
1321 different addresses (line 6). Second, the matrices and the initial
1322 candidate solution have to be copied on the GPU (lines 7 and 8).
1323 We notice that the input matrices are read-only structures and do
1324 not change during all the execution of the LS algorithm.
1325 Therefore, their associated memory is copied only once during all
1326 the execution. Third, comes the parallel
1328 in which each neighboring solution is generated, evaluated, and
1329 copied into the neighborhood fitnesses structure (from lines 10 to
1330 14). Fourth, since the order in which candidate neighbors are
1331 evaluated is undefined, the neighborhood fitnesses structure has
1332 to be copied to the host CPU (line 15). Then the selection
1333 strategy is applied to this structure (line 16): the exploration
1334 of the neighborhood fitnesses structure is done by the CPU.
1335 Finally, after a new candidate has been selected, this information
1336 is copied to the GPU (line 18). The process is repeated until a
1337 given number of iterations has been reached.
1340 % \caption{Template of an iterated tabu search\index{Metaheuristics!tabu~search} on GPU for solving the Q3AP} \label{LSGPU}
1341 % \begin{algorithmic}[1]
1342 % \STATE Choose an initial solution
1343 % \STATE Evaluate the solution
1344 % \STATE Initialize the tabu list
1345 % \STATE Allocate the Q3AP input data on GPU device memory
1346 % \STATE Allocate a solution on GPU device memory
1347 % \STATE Allocate a neighborhood fitnesses structure on GPU device memory
1348 % \STATE Copy the Q3AP input data on GPU device memory
1349 % \STATE Copy the solution on GPU device memory
1351 % \FOR{each generated neighbor on GPU}
1352 % \STATE Incremental evaluation of the candidate solution
1353 % \STATE Insert the resulting fitness into the neighborhood fitnesses structure
1355 % \STATE Copy the neighborhood fitnesses structure on CPU host memory
1356 % \STATE Select the best admissible neighboring solution
1357 % \STATE Update the tabu list
1358 % \STATE Copy the chosen solution on GPU device memory
1359 % \UNTIL{a maximum number of iterations reached}
1363 \begin{algorithm}[H]
1366 %\KwResult{how to write algorithm with \LaTeX2e }
1367 Choose an initial solution\;
1368 Evaluate the solution\;
1369 Initialize the tabu list\;
1370 Allocate the Q3AP input data on GPU device memory\;
1371 Allocate a solution on GPU device memory\;
1372 Allocate a neighborhood fitnesses structure on GPU device memory\;
1373 Copy the Q3AP input data on GPU device memory\;
1374 Copy the solution on GPU device memory\;
1376 \Repeat{a maximum number of iterations is reached}{
1377 \For{each generated neighbor on GPU}{
1378 Incremental evaluation of the candidate solution\;
1379 Insert the resulting fitness into the neighborhood fitnesses structure\;
1381 Copy the neighborhood fitnesses structure on CPU host memory\;
1382 Select the best admissible neighboring solution\;
1383 Update the tabu list\;
1384 Copy the chosen solution on GPU device memory\;
1386 \caption{template of an iterated tabu search on GPU for solving the Q3AP}
1391 \subsection {Experimental results}
1393 In this section, some experimental results related to the approach
1394 presented in Section \ref{ch8:ITS-Q3APSection} are reported. We
1395 recall that the approach is a GPU-based iterated tabu
1396 search (GPU-ITS) method
1397 consisting in an iterated local search (ILS) embedding a tabu
1398 search (TS) and where the
1399 generation/evaluation step of the TS process is executed on GPU.
1400 The ILS is used to improve the quality of successive local optima
1401 provided by TS methods. This is achieved by perturbing the local
1402 optima reached by the current TS process and reconsidering it as
1403 initial solution of the following TS process. Regarding our
1404 algorithm, the applied perturbation is a random number $\mu $ of
1405 swaps in either the first or the second permutation where $\mu \in
1406 [2:n]$ ($n$ is the instance size).\\
1408 Experiments have been carried out on a node of the Chirloute
1409 cluster of the Lille site. This is one of the 10 sites that
1410 currently make up Grid5000~\cite{grid5000}, the French
1411 experimental computational grid. A Chirloute cluster node consists
1412 of an Intel Xeon E5620 CPU and a NVIDIA Tesla Fermi M2050 (448
1413 cores) GPU type. The number of ILS iterations and the number of TS
1414 iterations were set respectively to 20 and 500. The tabu list size
1415 has been initalized to $\frac{m}{4}$, $m$ being the size of the
1418 Table \ref{ch8:ITSQ3APResults} reports the obtained results for
1419 the GPU-ITS using our large-sized neighborhood structure. The
1420 method has been tested on 5 Q3AP instances derived from the QAP
1421 Nugent instances in QAPLIB~\cite{burkard1991qaplib}. The average
1422 time measurement for 20 executions is reported in seconds and
1423 acceleration factors compared to a standalone CPU are also
1424 considered. The algorithm is stopped when a maximum number of
1425 iterations has been reached or when the optimal (or best known)
1426 value has been discovered. Average and max values of the
1427 evaluation function of the parallel GPU version have been
1428 measured. The number of successful tries (hits) and the average
1429 number of ILS iterations to converge to the optimal/best known
1430 value are also represented. The associated standard deviation for
1431 each average measurement is shown in small type characters.\\
1437 \begin{tabular}{|l|r|r|r|r|r|r|r|r|} \hline
1438 \multicolumn{1}{|c|}{Instance }& \multicolumn{1}{|c|}{Optimal}& \multicolumn{1}{|c|}{Average} & \multicolumn{1}{|c|}{Maximal} & \multicolumn{1}{|c|}{Hits} & \multicolumn{1}{|c|}{CPU} & \multicolumn{1}{|c|}{GPU} & \multicolumn{1}{|c|}{Speed-} & \multicolumn{1}{|c|}{ITS} \\
1439 & \multicolumn{1}{|c|}{/BKV\footnotemark }& \multicolumn{1}{|c|}{value }& \multicolumn{1}{|c|}{value }& \multicolumn{1}{|c|}{} & \multicolumn{1}{|c|}{time }& \multicolumn{1}{|c|}{time }& \multicolumn{1}{|c|}{up}& \multicolumn{1}{|c|}{iter.}\\ \hline
1440 Nug12-d & $580$ & $615.4$ & $744$ & $35$\% & $87.7$ & $2.5$ & $34.7 \times$& $16$ \\
1441 & & \tiny{$41.7$} & && \tiny{$40.9$} & \tiny{$0.9$} & & \tiny{$6.3$} \\ \hline
1442 Nug13-d & $1912$ & $1985.4$ & $2100$ & $20$\% & $209.2$ & $3.3$ & $63.5 \times$ & $17$ \\
1443 & & \tiny{$51.0$} & & & \tiny{$1.3$} & \tiny{$1.0$} & & \tiny{$5.6$} \\ \hline
1444 Nug15-d & $2230$ & $2418.1$ & $2580$ & $30$\% & $305.5$ & $5.2$ & $58.8 \times$ & $17$ \\
1445 & & \tiny{$135.3$} & & & \tiny{$164.5$} & \tiny{$1.3$} & & \tiny{$5.0$} \\ \hline
1446 Nug18-d & $17836$ & $18110.2$ & $18506$ & $10$\% & $1375.9$ & $12.8$ & $107.4 \times$ & $19$ \\
1447 & & \tiny{$157.8$} & & & \tiny{$123.5$} & \tiny{$2.6$} & & \tiny{$4.2$} \\ \hline
1448 Nug22-d & $42476$ & $43282.1$ & $44140$ & $25$\% & $4506.5$ & $32.7$ & $137.8 \times$ & $18$ \\
1449 & & \tiny{$529.6$} & & & \tiny{$341.1$} & \tiny{$6.6$} & & \tiny{$4.0$} \\ \hline
1451 \caption{Results of the GPU-based iterated tabu search for
1452 different Q3AP instances.} \label{ch8:ITSQ3APResults} % \center %Shashi
1458 %\caption{Results of the GPU-based iterated tabu search for different Q3AP instances}
1459 %\label{ch8:ITSQ3APResults} \center
1460 %\begin{tabular}{|c|c|c|c|c|c|c|c|c|} \hline
1461 %Instance & Opt./BKV\footnotemark & Average val. & Max val. & Hits & CPU time & GPU time & Accel. & ITS iter.\\ \hline
1462 %Nug12-d & $580$ & $615.4_{41.7}$ & $744$ & $35$\% & $87.7_{40.9}$ & $2.5_{0.9}$ & $\times 34.7 $& $16_{6.3}$ \\ \hline
1463 %Nug13-d & $1912$ & $1985.4_{51.0}$ & $2100$ & $20$\% & $209.2_{1.3}$ & $3.3_{1.0}$ & $\times 63.5$ & $17_{5.6}$ \\ \hline
1464 %Nug15-d & $2230$ & $2418.1_{135.3}$ & $2580$ & $30$\% & $305.5_{164.5}$ & $5.2_{1.3}$ & $\times 58.8 $ & $17_{5.0}$ \\ \hline
1465 %Nug18-d & $17836$ & $18110.2_{157.8}$ & $18506$ & $10$\% & $1375.9_{123.5}$ & $12.8_{2.6}$ & $\times 107.4$ & $19_{4.2}$ \\ \hline
1466 %Nug22-d & $42476$ & $43282.1_{529.6}$ & $44140$ & $25$\% & $4506.5_{341.1}$ & $32.7_{6.6}$ & $\times 137.8$ & $18_{4.0}$ \\ \hline
1470 \footnotetext{Best known value.}
1472 Regarding the execution time, the generation and evaluation of the
1473 neighborhood in parallel on GPU provides an efficient way to
1474 speedup the search process in comparison with a single CPU. In
1475 fact, for the smaller instance Nug12-d, the GPU version starts to
1476 be faster than the CPU one (acceleration factor of $34.7~\times$).
1477 As long as the problem size increases, the speedup grows
1478 significantly (up to $137.8 \times$ for the Nug22-d instance).
1479 This means that the use of GPU provides an efficient way to deal
1480 with large neighborhoods. Indeed, the proposed neighborhoods were
1481 unpractical in terms of single CPU computational resources for
1482 large Q3AP instances such as Nug22-d (estimated to around $2$
1483 hours per run). So, implementing this algorithm on GPU has allowed
1484 the exploitation of parallelism in such neighborhood and improved
1485 the robustness/quality of provided solutions.
1487 \section{Conclusion}
1488 \label{ch8:conclusion}
1490 This chapter has presented state-of-the-art GPU-based parallel
1492 a case study on implementing large neighborhood local search
1493 methods on GPUs for solving large benchmarks of the quadratic
1494 three-dimensional assignment problem (Q3AP). \\
1496 Implementing parallel
1498 GPU architectures poses new issues and challenges such as memory
1500 efficient mapping strategies between tasks to parallelize; and the
1501 GPU threads, thread divergence, and synchronization. Actually, most
1502 of metaheuristics have been implemented on GPU using different
1503 implementation strategies. In this chapter, a two-level
1504 classification of the reviewed works has been proposed: design
1505 level and implementation level. Design level regroups traditional
1506 parallel models used for metaheuristics while implementation level
1507 refers to the way those models are mapped to the GPU architecture.
1508 This classification focuses mainly on the mapping between the
1509 metaheuristic tasks to parallelize and the GPU threads. Indeed,
1510 the choice of a given mapping strategy strongly influences the
1511 other challenges (memory usage, communication, thread
1514 \putbib[Chapters/chapter9/biblio9]