new

[book_gpu.git] / BookGPU / Chapters / chapter9 / ch9.tex

\chapterauthor{Malika Mehdi and Ahc\`{e}ne Bendjoudi}{CERIST Research Center, DTISI, 3 rue des frères Aissou, 16030 Ben-Aknoun, Algiers, Algeria}

\chapterauthor{Lakhdar Loukil}{University of Oran, Algeria}

%\chapterauthor{Ahc\`{e}ne Bendjoudi}{CERIST Research Center, DTISI, 3 rue des frères Aissou, 16030 Ben-Aknoun, Algiers, Algeria}

\chapterauthor{Nouredine Melab}{Université Lille 1, LIFL/UMR CNRS 8022, 59655-Villeneuve d'Ascq cedex, France}

\chapter{Parallel GPU-accelerated metaheuristics}
\label{chapter9}
\section{Introduction}
This chapter presents GPU-based parallel
metaheuristics\index{metaheuristics!parallel metaheuristics},
challenges, and issues related to the particularities of the GPU
architecture and a synthesis on the different implementation
strategies used in the literature. The implementation of parallel
metaheuristics on
GPUs is not straightforward. The traditional models used in CPUs
must be rethought to meet the new requirements of GPU
architectures. This chapter is organized as follows. Combinatorial
optimization\index{combinatorial optimization} and  resolution
methods are introduced in Section~\ref{ch8:sec:optim}. The main
traditional parallel models used for metaheuristics are recalled
in Section~\ref{ch8:sec:paraMeta}.
Section~\ref{ch8:sec:challenges} highlights the  main challenges
related to the GPU implementation of metaheuristics. A
state-of-the-art of GPU-based parallel
metaheuristics is
summarized in Section~\ref{ch8:sec:state}. In Section~\ref{ch8:sec:frameworks}, the main developed GPU-based
frameworks for metaheuristics are described. Finally, a case study
is presented in Section~\ref{ch8:sec:case} and some concluding
remarks are given in Section~\ref{ch8:conclusion}

\section{Combinatorial optimization}
\label{ch8:sec:optim}

Combinatorial optimization (CO) is a branch of applied and discrete mathematics.
It consists in finding optimal configuration(s) among a finite set of possible configurations
(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})


%\begin{minipage}{0.5\linewidth}
A combinatorial problem $P=(S,f)$ can be defined by
\begin{itemize}
\item a set of decision variables $X$,
\item an objective function $f$ to optimize (minimize or maximize) over the set $S$,
\item subject to constraints on the decision variables.
\end{itemize}
%\end{minipage}

COPs are generally formulated as mixed integer programming
problems (MIPS) and most of them are NP-hard~\cite{garey}.
According to the quality level of solutions and deadlines required
for solving an optimization problem, two classes of optimization
methods can be distinguished: \emph{exact methods} and
\emph{approximate methods}. Exact methods allow one to reach
optimal solution(s) of the handled optimization problem with a
proof of its or their optimality. The known methods of this
class are the \emph{branch and bound technique}, \emph{dynamic
programming}, \emph{constraint programming}, and \emph{A*
algorithm}. However, optimization problems, whether practical or
academic, are often complex and NP-hard. Moreover, a large number
of real-life optimization problems encountered in science,
engineering, economics, and business are usually large-sized
problems for which the size of the potential solution domain
increases dramatically with the size of the problem instance. Such
problems cannot be tackled using exact methods due to the
excessive computation time needed by these methods to find optimal
solution(s). In such a situation, approximate methods (or
\textit{metaheuristics}) offer an alternative approach to solve
these problems. Indeed, these methods allow one to reach good
quality solutions in reasonable computation time compared to exact
methods but with no guarantee to find optimal or even bounded
solutions. This is due to the nature of the search process adopted
by these approaches which consists of performing a
search in a subset of the whole search space.

Regarding the number of solutions considered at each iteration in
the search process, two classes of metaheuristics can be
distinguished~\cite{talbi2009mfdti}: \textit{solution-based} and
\textit{population-based} metaheuristics. In the rest of this
chapter, the term \emph{s-metaheuristic} refers to solution-based
metaheuristic and the term \emph{p-metaheuristic} refers to
population-based metaheuristic. In s-metaheuristics, the search
process starts with a single solution (generally set at random)
and iteratively improves it by exploring its neighborhood in the
search space. The most known methods in this class are local
search methods that include \emph{simulated
annealing}\index{metaheuristics!simulated annealing}~\cite{Kirkpatrick1983SA},
\emph{tabu search}~\cite{Glover1989TS}, \emph{iterated local
search\index{metaheuristics!iterated local
search}}~\cite{stutzle2006ILSforQAP}, and \emph{variable
neighborhood search}~\cite{HansenMladenovic1997VNS}.

Unlike s-metaheuristics, p-metaheuristics start with a population
of solutions and implement an iterative process that evolves the
current population towards a new population of better quality
solutions. The process is repeated until a stopping criterion is
satisfied. \emph{Evolutionary algorithms}, \emph{swarm
optimization}, and \emph{ant colonies} fall into this class.


\section{Parallel models for metaheuristics}\label{ch8:sec:paraMeta}
Optimization problems, whether real-life or academic, are more
often NP-hard and CPU time and/or memory consuming. Metaheuristics
allow the significant reduction of the computational time of the search
process but remain time-consuming particularly when it comes
dealing with large-sized problems. 

The use of parallelism in the design of metaheuristics is a relevant
approach that is widely adopted by the combinatorial optimization
community for various reasons:

\begin{itemize}
  \item One of the main goals of parallelism is to reduce the search
time. This will allow the design of high performance optimization
methods and the solving of large-sized optimization problems.

  \item Sequential processor architectures have reached their
physical limit which prevents creating faster processors. The
current trend of microprocessor manufacturers consists of placing
multicores on a single chip. Nowadays, laptops and
workstations are multicore processors. In addition, the evolution
of network technologies and the proliferation of broadband
networks have made possible the emergence of clusters of
workstations (COWs), networks of workstations (NOWS), and
large-scale networks of machines (grids) as platforms for
high-performance computing.
\end{itemize}

From the granularity of parallelism point of view, three major parallel
models for metaheuristics can be distinguished~\cite{talbi2009mfdti}: \emph{algorithmic-level}\index{metaheuristics!algorithmic-level parallelism},
\emph{iteration-level}\index{metaheuristics!iteration-level parallelism}, and \emph{solution-level} as illustrated in Figure~\ref{ch8:fig:paraMeta}. 

\begin{figure}[h!]
\centerline{\includegraphics[width=0.6\textwidth]{Chapters/chapter9/figures/paraMeta.pdf}}
\caption{Parallel models for metaheuristics.}
\label{ch8:fig:paraMeta}
\end{figure}

\begin{itemize}

\item{In the
algorithmic-level parallel model, several self-contained metaheuristics
are launched in parallel. The parallel
metaheuristics\index{metaheuristics!parallel metaheuristics} may
start with identical or different solutions (s-metaheuristics
case) or populations (p-metaheuristics case). Their parameter
settings such as the size of tabu list for tabu
search\index{metaheuristics!tabu search}, transition probabilities
for ant colonies, mutation and crossover probabilities for
evolutionary
algorithm\index{metaheuristics!evolutionary algorithm}s may be the
same or different. The parallel processes may evolve independently
or in a cooperative manner. In cooperative parallel models, the
algorithms exchange information related to the search during
evolution in order to find better and more robust solutions.}

\item{In the iteration-level parallel model, the focus is on the parallelization
of each iteration of the metaheuristic. Indeed, metaheuristics are
generally iterative search processes. Moreover, the most
resource-consuming part of a metaheuristic is the evaluation of
the generated solutions at each iteration. For s-metaheuristics
(e.g., tabu search\index{metaheuristics!tabu search}, simulated
annealing, variable neighborhood search), the evaluation and
generation of the neighborhood is the most time-consuming step of
the algorithm particularly when it comes to dealing with large
neighborhood sets. In this parallel model, the neighborhood is
decomposed into partitions, and each partition is evaluated in a
parallel way. For p-metaheuristics (e.g., evolutionary
algorithms, ant
colonies, swarm optimization), the
iteration-level
parallel model arises naturally since these metaheuristics deal
with a population of independent solutions. In evolutionary
algorithms, for
instance, the iteration-level model consists of decomposing the whole population
into several partitions where each partition is evaluated in
parallel.}

\item{In the
solution-level
parallel model, the focus is on the parallelization of the
evaluation of a single solution. This model is useful when the
objective function and/or the constraints are time and/or memory
consuming. Unlike the two previous parallel models, the
solution-level\index{metaheuristics!solution-level parallelism}
parallel model is problem-dependent.}
\end{itemize}
\clearpage
\section{Challenges for the design of GPU-based metaheuristics}
\label{ch8:sec:challenges}

Developing GPU-based parallel
metaheuristics\index{metaheuristics!parallel metaheuristics} is
not straightforward. The parallel models have to be rethought to
meet the new requirements of the GPU architecture. Several major
issues have to be taken into account both at design and
implementation levels to develop efficient metaheuristics on GPU.
These issues are mainly related to the size and latency of the GPU
memories, thread synchronization and divergence, the distribution
of tasks, and data transfer between the CPU and
GPU~\cite{luongMultiStart}.

\subsection{Data placement on a hierarchical memory}
\index{GPU!data placement} During the execution of
metaheuristics on GPU, the different threads may access multiple
data structures from multiple memory spaces. These memories have
different sizes and access latencies. Nevertheless, faster
memories (registers, shared and constant memories) are generally
very small in size, and the larger memories (global memory) are
relatively slower. However, p-metaheuristics require the
exploration of a large amount of individuals to diversify the
search. Moreover, an efficient execution of s-metaheuristics
requires exploring large neighborhoods. Thus, programmers have to
take into account this point to efficiently place the different
data structures of the metaheuristic on the different memories to
benefit from both the faster memories and the larger ones. A deep
study has to be performed on both the metaheuristic data
structures (size and access frequency) and the GPU memories (size
and access latency) to identify which data will be placed on which
memory. Generally, the most accessed ones should be put on faster
memories (registers, shared memory) and larger ones on the larger
memories (global memory). Also, an efficient mapping between
threads and the corresponding metaheuristic elements (one neighbor
per thread, one individual per thread, single population per
threads block, one ant per thread, etc.) must be defined to ensure
a maximum occupancy of the GPU and
to cover CPU/GPU communication\index{GPU!CPU/GPU communication} and memory access times.

According to the used metaheuristic and to the handled problem, the data
values may have different types and different ranges of their values. The data
types should then be chosen carefully and the ranges of the data values should
be analyzed to minimize the amount of occupied memory space.

In addition to the size and latency of GPU memory issues, the
memory access pattern is another important issue to be dealt with
to speedup GPU-based metaheuristics. Indeed, the different
memories have been designed to achieve specific features that
programmers must take into account to optimize their codes and
then to benefit from these features. For instance, the global
memory is optimized for coalesced accesses. The texture and the
constant memory are read-only memories. The texture is optimized
for uncoalesced accesses and the constant memory is optimized for
simultaneously accesses of all threads to the same
location~\cite{luongMultiStart}. Therefore, to improve the
performance of the kernel execution, the programmers should put
coalesced data on global memory, uncoalesced read-only data (e.g.,
metaheuristic instance data) on the texture, concurrent read-only
data (e.g., data for fitness evaluation that all threads
concurrently access) on the constant memory, and the most accessed
data structures (e.g., population of individuals for a CUDA thread
block) on the shared memory.

\subsection{Threads synchronization}
\index{GPU!threads synchronization} The thread
synchronization issue is caused by both the GPU architecture and
the synchronization requirements of the implemented method.
Indeed, GPUs are based on a multicore architecture organized into
several multi-processors (Streaming Multiprocessor SM) supporting
the SPMD model (Single Program Multiple Data). Each SM contains
several cores executing the same instruction of different threads
following the SIMD model (Single Instruction Multiple Data). These
threads belong to a warp (a group of 32 threads) and handle
different data elements. An efficient execution of applications on
GPU is achieved when launching a large amount of threads
(thousands of threads)~\cite{CUDA}. However, the execution order
of these thousands of threads is unknown by the programmer which
makes the prediction of their execution order a challenging issue.
Plus, the developer has to control explicitly the
threads through the insertion of barrier synchronizations in the
codes to avoid concurrent accesses to data structures and to meet
some requirements related to data-dependent synchronizations.

\subsection{Thread divergence}

Thread divergence\index{GPU!thread divergence} is
another challenging issue in GPU-based
metaheuristics~\cite{cecilia, pugace, audreyANT}. Generally,
metaheuristics contain irregular loops and conditional
instructions when generating and evaluating the neighborhood
(s-metaheuristics), and the population (p-metaheuristics) in the
same block. In addition,  the decision to apply a crossover or a
mutation on an individual in a genetic
algorithm and the
exploration of different paths using an ant
colony\index{metaheuristics!ant colony optimization} are random
operations. Threads of the same warp have to execute
instructions simultaneously  leading to different branches whereas
in an SIMD model the threads of a same warp execute the same
instruction. Consequently, the different branches of a conditional
instruction which is data-dependent lead to a serial execution of
the different threads degrading the performance of the application
in terms of execution time. The challenge here is then to revisit
the traditional irregular metaheuristic codes to eliminate these
divergences.
\clearpage
\subsection{Task distribution and CPU/GPU communication}

The performance of GPU-based metaheuristics in terms of execution
time could be improved by choosing the most appropriate parallel
model (algorithmic-level, instruction-level,
solution-level).
Moreover, an efficient decomposition of the metaheuristic and an
efficient assignment of code portions between the CPU and GPU
should be adopted. The objective is to take benefit from the GPU
computing power without affecting the efficiency and the behavior
of the metaheuristic and without losing performance in CPU/GPU
communication\index{GPU!CPU/GPU communication} and
memory accesses. In order to decide which part of the
metaheuristic will be executed on which component, one should
perform a careful analysis on the serial code of the
metaheuristic, identify the compute-intensive tasks (e.g.,
exploration of the neighborhood, evaluation of individuals), and
then offload them to the GPU, while the remaining
tasks still run on the CPU in a serial way. 

The CPU/GPU communication is done through the global
memory which is a slow memory making the memory transfer between
the CPU and GPU time-consuming which can significantly degrade the
performance of the application. Accesses to this memory should be
optimized by minimizing the amount of transferred data between the
two components in order to reduce the communication time and,
therefore, the whole execution time of the metaheuristic.

\section{State-of-the-art parallel metaheuristics on GPUs}
\label{ch8:sec:state}
After more than two decades of research by the combinatorial optimisation
community devoted to developing adequate parallel metaheuristics for different types of
parallel architectures (clusters, supercomputers and grids), the actual developement
of General Perpose GPU (GPGPU) brings new challenges for parallel metaheuristics on SIMD architectures.

The first works on metaheuristic algorithms implemented on GPUs
started on old graphics cards before the appearance of modern GPUs
equipped with high-level programming interfaces such as CUDA and
OpenCL. Among these pioneering works we cite the work of Wong et al.~\cite{wongOldGPU2006} dealing with the
implementation
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
is implemented on old GPU architectures. Yu et al.~\cite{yu2005} and
Li et al.~\cite{li2007} proposed a full parallelization of genetic
algorithms  on old GPU architectures using
shader libraries based on Direct3D and OpenGL.

Such architectures are based on preconfigured pipelined stages
used to accelerate the transformation of 3D geometric primitives
into pixels. Implementing a general-purpose algorithm on such
preconfigured architectures is very hard and requires the
tailoring of the algorithm's data and instructions to fit the
pipelined stages of the GPU. Since then, GPU architectures have
evolved to become programmable using high-level programming
interfaces. In this section we will focus only on recent
state-of-the-art works dealing with metaheuristics implementation
on modern programmable GPUs. In this review two classes are
considered: (1) s-metaheuristics on GPUs and (2) p-metaheuristics
on GPUs. A comparative study is done of the main works and a
classification of the different existing strategies is proposed in
Section~\ref{ch8:sec:synthesis}.

\subsection{Implementing solution-based metaheuristics on GPUs}


A very basic local search algorithm starts with an initial
solution generated either at random or by the mean of a specific
heuristic and is based on two elementary components: the
generation of neighborhood structures using an elementary move
function and a selection function that determines which solution
in the current neighborhood will replace the actual search point.
The search continues as long as there is improvement and stops
when there is no better solution in the current neighborhood. The
exploration (or evaluation) of the different moves of a given
neighborhood is done independently for each move. Thus, the
easiest way to accelerate a local search algorithm is to
parallelize the evaluation of the neighborhood (instruction-level
parallelism). This is  by far the most used scheme in the
literature for parallelizing local search algorithms on GPUs.
Nevertheless, small neighborhoods may lead to  nonoptimal
occupation of the CUDA threads which may lead in turn to an
overhead due to the communication and memory latencies. Therefore,
large neighborhoods are necessary for efficient implementation of
local searches on GPUs.

Luong et al.~\cite{luong2010large} proposed efficient
mappings for large neighborhood structures on GPUs. In this work,
three different neighborhoods are studied and mapped to the
hierarchical GPU for binary problems. The three neighborhoods are
based on the Hamming distance. The move operators used in
the three neighborhoods consider Hamming distances of 1, 
2, and 3 (this consists on flipping the binary value of one, two,
or three positions at a time in the candidate binary vector).
In~\cite{luong2010large}, each thread is associated to a unique
solution in the neighborhood. The addressed issue is how to
efficiently map the different neighborhoods on the device memory,
more explicitly, how to calculate the memory index of each
solution associated to each CUDA thread's \textit{id}.
%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.
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.\\

In the same context, Deevacq et al.~\cite{audreyANT}
proposed two parallelization strategies inspired by the multiwalk
parallelization strategy, of a 3-opt iterated local
search algorithm (ILS)
over a CPU/GPU architecture. In the first strategy, each Local
Search (LS) is associated to a unique CUDA thread and improves a
unique solution by generating its neighborhood. The second
strategy associates each solution to a CUDA block and the
neighborhood exploration is parallelized on the block threads. In
the first strategy, since several LS are executed on different
solutions on each Multi Processor (MP), the data structures should
be stored on the global memory, while the exploration of a single
solution at a time in the second strategy allows the use of the
shared memory to store the related data structures. The two
strategies have been experimented on standard benchmarks of
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$.

The same strategy is used by Luong et al.
in~\cite{luongMultiStart} to implement multistart parallel local
search algorithms (a special case of the
algorithmic-level parallel model where several homogeneous LS
algorithms are used). The multistart model is combined with
iteration-level
parallelism: several LS algorithms are managed by the CPU and the
neighborhood evaluation step of each algorithm is parallelized on
the GPU (each GPU thread is associated with one neighbor and
executes the same evaluation function kernel). The advantage of
such a model is that it allows a  high occupancy of the GPU
threads. Nevertheless, memory management causes new issues due to the
quantity of data to store and to communicate between CPU     and
GPU. A second proposition for implementing the same model on GPU
consists of implementing the whole LS processes on GPU with each
GPU thread being associated to a unique LS algorithm. This solves
the communication issue encountered in the first model. In
addition, a memory management strategy is proposed to improve the
efficiency of the
algorithmic-level model: texture memory is used to avoid memory latency
due to uncoalesced memory accesses. The proposed approaches are
implemented on the quadratic assignment problem (QAP) using CUDA.
The acceleration rates obtained for the
algorithmic-level with usage of texture memory rise from $7.8\times$ to
$12\times$ (for different
QAP benchmark sizes). 

Janiak et al.~\cite{Janiak_et_al_2008} implemented two
algorithms for TSP and the flow-shop scheduling problem (FSP).
These algorithms are based on a multistart tabu
search model. Both of the 
algorithms exploit multicore CPU and GPU. A full parallelization
on GPU is adopted using shader libraries where each thread is
mapped with one tabu search.
However, even though their experiments report that the use of GPU
speedups the serial execution almost $16 \times$, the mapping of
one thread with one tabu search
requires a large number of local search algorithms to
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.

Luong et al.~\cite{luong2012ppsn} proposed a GPU-based
implementation of hybrid metaheuristics on heterogeneous parallel
architectures (multicore CPU  coupled to one GPU).  The challenge
of using such a heterogeneous architecture is how to distribute
tasks between the CPU cores and the GPU in such a way to have
optimal performances. Among the three traditional parallel models
(solution-level,
iteration-level
and algorithmic-level), the authors point out that the most convenient
model for the considered heterogeneous architecture is a hybrid
model combining
iteration-level
and algorithmic-level models. Several CPU threads execute several instances
of the same S-metaheuristic in parallel while the GPU device is
associated to one CPU core and  used to accelerate the
neighborhood calculation of several S-metaheuristics at the same
time.
 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
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}. 

All the previously noted works  exploit the same parallel models
used on CPUs based on the task parallelism. A different
implementation approach is proposed by Paul in~\cite{gerald2012}
to implement a simulated
annealing (SA) algorithm
for the QAP on GPUs. Indeed, the author used a preinitialized
matrix \emph{delta} in which the incremental evaluation of simple
swap moves are calculated and stored relative to the initial
permutation $p$. For the GPU implementation, the authors used the
parallel implementation of neighborhood exploration. The
time-consuming tasks in the SA-matrix are identified by the
authors as updating the matrix and passing through it to select
the next accepted move. To initialize the delta matrix, several
threads from different blocks explore different segments of the
matrix (different moves) at the same time.  In order to select the
next accepted swap, several threads are also used. Starting from
the last move, a group of threads explores different subsets of
the delta matrix. The shared memory is used to preload all the
necessary elements for a given group of threads responsible for
the updating of the delta matrix. The main difference in this work
compared to the previous works resides in the introduction of a
data parallelism using the precalculated delta matrix. The use of
this matrix allows the increase in the number of threads
involved in the evaluation of a single move. Experimentations are
done on standard QAP instances from the
QAPLIB~\cite{burkard1991qaplib}. Speedups up to $10 \times$ are
achieved by the GPU implementation compared
to the same sequential implementation on CPU using SA-matrix.

\subsection[Implementing population-based metaheuristics\hfill\break on GPUs]{Implementing population-based metaheuristics on GPUs}

State-of-the-art works dealing with the implementation of
p-metaheuristics on GPUs generally rely on parallel models and
research efforts done for parallelizing different classes of
p-metaheuristics over different types of parallel architectures:
supercomputers, clusters, and computational grids. Three main
classes of p-metaheuristics are considered in this section:
evolutionary
algorithms (EAs), ant
colony optimization
(ACO), and particle
swarm optimization\index{metaheuristics!particle swarm
optimization} (PSO).
\clearpage
\subsubsection*{Evolutionary algorithms}

Traditional parallel models for EAs are classified in three main
classes: coarse-grain models using several subpopulations
(islands), master-slave models used for the parallelization of CPU
intensive steps (evaluation and transformation), and cellular
models based on the use of one population disposed (generally) on
a toroidal grid.

The three traditional models have been implemented on GPUs by several researchers for different
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.

In~\cite{kannan}, Kannan and Ganji present a CUDA implementation
of the drug discovery application Autodock (molecular docking
application). Autodock uses a genetic
algorithm to find optimal
docking positions of a ligand to a protein. The most
time-consuming task in Autodock is the fitness function
evaluation. The fitness function used for a docking problem
consists of calculating the energy  of the ligand-protein complex
(sum of intermolecular energies). The authors explore two
different approaches to evaluate the fitness function on GPU. In
the first approach, each GPU thread calculates the energy function
of a single individual. This approach requires the use of
large-sized populations to saturate the GPU (thousands of
individuals). In the second approach each individual is associated
with one thread block. The evaluation of the energy function is
performed by the threads of the same block. This allows the use of
medium population sizes (hundreds of individuals) and the
acceleration of a single fitness evaluation. Another great
advantage of the per block approach resides in the use of shared
memory instead of global memory to store all the information
related to each individual. The obtained speedups range from  $10
\times$ to $47 \times$ for population sizes
ranging from $50$ to $10000$. 

Maitre et al.~\cite{maitre2009} also exploited the
master-slave parallelism of EAs on GPUs using EASEA. EASEA is a
C-like metalanguage for easy development of EAs. The user writes a
description  of the problem-specific components (fitness function,
problem representation, etc) in EASEA. The code  is then compiled
to obtain a ready-to-use evolutionary
algorithm. The EASEA
compiler  uses genetic
algorithm LIB and EO
Libraries to produce C++ or JAVA written EA codes.
In~\cite{maitre2009}, the authors  proposed an extension of EASEA
to produce CUDA code from the EASEA files. This extension has been
used  to generate a master-slave parallel EA in which the
sequential algorithm is maintained on CPU and the population is
sent to GPU for evaluation. Two problems have been considered
during the experiments: a benchmark mathematical function and a
real problem (molecular structure prediction). In order to
maximize the GPU occupation, very large populations are used (from
$2000$ to $20000$). Even though transferring such large
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.

The master-slave  model is efficient when the fitness function is
highly time intensive. Nevertheless, it requires the use of
large-sized populations in order to saturate the GPU unless the
per-block is used (one individual perblock) when the acceleration
of the fitness function itself is possible. The use of many
subpopulations of medium sizes is another way to obtain a maximum
occupation of the GPU. This is coarse-grained parallelism (island
model).

The coarse-grained model is  used by Pospichal et al.
in~\cite{pospichal10} to implement a parallel genetic
algorithm on GPU. In this
work the entire genetic
algorithm is implemented
on GPU. This choice is motivated by the overhead engendered by the
CPU/GPU communication
when only population evaluation is performed on GPU. Each
population island is mapped with a CUDA thread block and each
thread is responsible for a unique individual. Subpopulations are
stored on shared memory of each block. Nevertheless, because
interblock communications are not possible on the CUDA
architecture, the islands evolve independently in each block, and
migrations are performed
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.

The same strategy is also adopted by Tsutsui and Fujimoto
in~\cite{tsutsuiGAQAP}  to implement a coarse-grained genetic
algorithm on GPU to solve
the QAP. Initially, several subpopulations are created on CPU and
transferred to the global memory. The subpopulations are organized
in the global memory into blocks of $8$ individuals in such a way
to allow coalesced memory access by the threads of the same thread
block. Each sub-population is allocated to a single thread block
in the GPU and transfered to the shared memory to start evolution.
Population evaluation and transformation are done in parallel by
the different threads of a given block. Migration is also done
through the global memory. Experiments are performed on standard
QAP benchmarks from the QAPLIB~\cite{burkard1991qaplib}. The GPU
implementation reached speedups of $2.9\times$ to $12.6 \times$
compared to a single core implementation of a coarse-grained
genetic algorithm on a
Intel i7 processor.

Nowotniak and Kucharski~\cite{nowotniak} proposed a GPU-based
implementation of a Quantum Inspired Genetic Algorithm (QIGA). The 
parallel model used is a hierarchical model based on two levels: each
thread in a block transforms a unique individual and a different
population is assigned to each block. The algorithm is run
entirely on GPU. A different instance of the QIGA is run on each
block and the computations have been shared between 8 GPUs. This
approach is very convenient to speed up the experimental process
on metaheuristics that require several independent runs of the
same algorithm in order to assess statistical efficiency. The
populations are stored in the shared memory, the data matrix used
for fitness evaluation is placed in read only constant memory, and
finally seeds for random numbers generated on the GPU are stored
in the global memory.

In coarse-grained parallelism, the use of the per-block approach
to implement the islands (one subpopulation per thread block) is
almost natural  and it allows the use of shared memory to store
the subpopulations. Fine-grained parallel models for EAs (cellular
EAs) have been used by many authors to implement parallel EAs on
GPUs. Indeed, the fine-grained parallelism of EAs fits perfectly
to the SIMD architecture of the GPU. 

Pinel et al. in~\cite{pinel2012JPDC} developed a highly
parallel synchronous cellular genetic
algorithm (CGA), called
GraphCell, to solve the independent task scheduling problem on GPU
architectures. In CGAs, the population is arranged into a
two-dimensional toroidal grid where only neighboring solutions are
allowed to interact with each other  during the recombination
step. In GraphCell, two recombination operators for CGA are
especially designed to run efficiently on GPU. Indeed, instead of
assigning a single thread to each solution of the population to
perform the recombination, in GraphCell, a single thread is
assigned to each task  of each solution. Offsprings are created by
independently modifying the assignment of some tasks in the
current solution. Mainly, each thread chooses one neighboring
solution in the grid as second parent using different selection
strategies and assigns one  task of  the solution (first parent)
to the machine for which the same task is assigned in the second
parent.  This way, the number of threads used for the
recombination step is equal to population size $\times$  size of
the solution (number of tasks). This leads to a high number of
threads used to accelerate the recombination operators especially
when dealing with large instances of the problem. In addition to
the recombination operators, the rest of the CGA steps are  also
parallelized on GPU (fitness evaluation, mutation, and
replacement).

A similar work is proposed by Vidal and Alba  in~\cite{albaCGAGPU}
where a  CGA using a toroidal grid is completely implemented on
GPU. A direct mapping between the population and the GPU threads
is done. At each step, several threads execute the same operations
on each individual independently: initialization, computing the
neighborhood, selection, crossover, mutation, and evaluation.  A
synchronization is done for all threads to perform the replacement
stage and form the next generation. All the data of the algorithm
is placed on the global memory. Several experiments have been
performed on a set of standard benchmark functions with different
grid sizes ranging from $32^2$ to $512^2$. The speedups reached by
the GPU implementation against the CPU version range from
$5\times$ to $24\times$ and increase as the size of the population
increases. However, the CPU implementation runs faster than the GPU
version for all the tested  benchmarks when the size of the
population is set to $32^2$. When the size of the population is
too small, there are not enough computations to cover the overhead
created by the call of kernel functions, CPU/GPU
communications,
synchronization, and access to global memory. Finally, an
interesting review on GPU parallel computation in bio-inspired
algorithms is proposed by Arenas et al.
in~\cite{arenas2011}. 

\subsubsection*{Ant colony optimization}

Ant colony optimization
(ACO) is another
p-metaheuristic subject to parallelization on GPUs.
State-of-the-art works on parallelizing
ACO focus on
accelerating the tour construction step performed by each ant by
taking a task-based parallelism approach, with pheromone
deposition on the CPU.

In~\cite{cecilia}, Cecilia et al.  present a GPU-based
implementation of
ACO for TSP where
the two steps (tour construction and pheromone update) are
parallelized on the GPU. A data parallelism approach is used to
enhance the performance of the tour  construction step. The
authors use two categories of artificial ants: queen ants
associated with CUDA thread-blocks and worker ants associated with
CUDA threads. A queen ant represents a simulated ant and worker
ants collaborate with the queen ant to accelerate the decision
about the  next city to visit. The tour construction step of each
queen ant is accelerated. Each worker ant  maintains a history of
the search in a tabu list containing a chronological ordering of
the already visited cities. This memory is used to determine the
feasible neighborhoods. After all queen ants have constructed
their tours, the pheromone trails are updated. For pheromone
update, several GPU techniques are also used to increase the data
bandwidth of the application mainly by the use of precalculated
matrices that are easily updated by several threads (one thread
per matrix entry). The achieved speedups are $21 \times$ for tour
construction and  $20 \times$ for pheromone updates.

In another work, Tsutsui and Fujimoto~\cite{tsutsui} propose  a
hybrid algorithm combining
ACO metaheuristic
and Tabu Search (TS) implemented
on GPU to solve the QAP. A solution of QAP is represented as a
permutation of ${1,2,..,n}$ with $n$ being the size of the
problem. The TS algorithm is based on the 2-opt neighborhood
(swapping of two elements $(i,j)$ in the permutation). The authors
point out that the move cost of each neighbor depends on the
couple $(i,j)$. Two groups of moves are formed according to the
move cost. In order to avoid thread divergence\index{GPU!thread divergence} within the same warp, the
neighborhood evaluation is parallelized in such a way to assign
only moves of the same cost to each thread warp. This strategy is
called MATA for Move-cost Adjusted Thread Assignment. Concerning
the memory management\index{GPU!memory management}, all
the data of the ACO\index{metaheuristics!ant colony optimization}
(population, pheromone matrix), QAP matrices, and tabu list are
placed on the global memory of the GPU. Shared memory is used only
for working data common to all threads in a given block.
 All the
steps of the hybrid algorithm
ACO-TS
(ACO initialization,
pheromone update, construct solutions, applying TS) are
implemented as kernel functions on the GPU. The GPU/CPU
communications are only used to transfer the best-so-far solution
in order to verify termination conditions. The experiments
performed on standard QAP benchmarks showed that the GPU
implementation using MATA obtained a  speedup of $19 \times$
compared to the CPU implementation,  compared with a speedup of
only $5 \times$
when MATA is not used. 

\subsubsection*{Particle swarm optimization}
In~\cite{zhou2009} Zhou and Tan propose a full GPU implementation
of a standard PSO algorithm. All the data is stored in global
memory (velocities, positions, swarm population, etc). Only
working data is copied to shared memory by each thread. The four
steps of the PSO have been parallelized on GPU: fitness evaluation
of the swarm, update of local best and global best  of each
particle, and update of velocity and position of each particle.
The same strategy is used to parallelize the first and last
steps: the evaluation of fitness functions is performed in
parallel for each dimension  by several threads. It is the case
for the update of position and velocities of each particle: one
dimension at a time is updated for the whole swarm by several
threads. Random numbers needed for  updating the velocities and
positions for the whole PSO processes are generated on CPU at the
starting of the algorithm and transferred to the GPU global
memory. For the steps 2 and 3 (update of local best and global
best of each particle), the GPU threads are mapped to the \emph{N}
particles of the swarm one to one. Experiments done on four
benchmark functions show speedups ranging from $3.7 \times$ to
$9.0 \times$ for swarm sizes
ranging from $400$ to $2800$.

\subsection{Synthesis of the implementation strategies}
\label{ch8:sec:synthesis} After reviewing some works dealing with
GPU-based implementation of metaheuristics, in this section we
will try to come up with a classification of the different
strategies used in the literature for the implementation of
metaheuristics on GPUs. One may distinguish between the parallel
models adopted in each metaheuristic (design level)  and the way
they are exploited on GPU architectures (implementation level).
Indeed, even though the parallelization models used in most works
for GPUs are derived from the traditional parallel models of each
metaheuristic (on CPU), their implementation could take a
different way and sometimes it may result in new parallel models
customized for GPUs.

Traditional parallel models for metaheuristics are based on an
intuitive task parallelism: the independent tasks of the
algorithms are simply parallelized. For example, in the case of
s-metaheuristics, the evaluation of large neighborhoods could be
done in parallel since there is no synchronization at this step of
the algorithm. That is the case of EAs  when it comes to applying
the evaluation step. Nevertheless, because of the particularity of
the GPU architecture, some authors have used new implementation
techniques to enhance the data parallelism in the sequential
algorithms in order to increase the data throughput of the
application.

From this observation we propose the following classification
based on 2 levels: design level and implementation level as
illustrated in Figure~\ref{ch8:fig:classification}. The design
level regroups the three classes of parallel models used in
metaheuristics
(solution-level,
iteration-level,
algorithmic-level) with examples for s-metaheuristics, EAs,
ACO and PCO. This
classification  is principally built from the reviewed
state-of-the-art works  in the previous section. The
implementation level refers to the way these parallel design
models are implemented on GPU. This classification focuses only on
the mapping strategies between the GPU threads and the
parallelized tasks (neighborhood evaluation, solution
construction, and so on). The different implementation strategies
are explained in the following sections.

\begin{figure}[h]
\centerline{\includegraphics[width=1\textwidth]{Chapters/chapter9/figures/classification.pdf}}
\caption{A two level classification of state-of-the-art GPU-based parallel metaheuristics.}
%[Comparison of number of cores in a CPU and in a GPU]
\label{ch8:fig:classification}
\end{figure}


\subsubsection*{GPU thread mapping for solution-level parallelism}
\index{GPU!thread mapping} Parallel
models at solution level consist of parallelizing a time intensive
atomic task of the algorithm. Generally, it consists of the
fitness evaluation~\cite{kannan}. Nevertheless, crossover
operators have also been parallelized by some
authors~\cite{pinel2012JPDC}. These kinds of models are not always
possible as they are problem-dependent. The GPU implementation of
solution-level models uses the perblock mapping: each solution is
associated to a block of threads. A second level of parallelism is
used inside each block to parallelize the fitness evaluation of a
single solution. This kind of mapping allows the use of shared
memory to store the data structures of the solution and does not
require the use of very large neighborhoods or populations.
%data parallelism in SA-matrix to parallelize

\subsubsection*{GPU thread mapping for iteration-level parallelism}
\index{GPU!thread mapping}
Iteration-level parallelism consists of parallelizing the tasks
performed independently on different solutions. Different mapping
strategies are adopted in the reviewed works to implement these
models.

In Figure \ref{ch8:fig:classification}, the first example of
iteration-level
parallelism is the parallel evaluation of neighborhoods in
s-metaheuristics. In most of the reviewed works,  a per-thread
mapping approach is used: each solution of the neighborhood is
mapped to a unique thread in the GPU for
evaluation~\cite{luong2010large, audreyANT}. The same mapping is
used for master-slave parallel EAs to accelerate the evaluation of
large populations. This kind of mapping is only efficient for very
large neighborhoods and very large populations (to saturate the
GPU). Many authors have pointed out that the use of such large
populations (or neighborhoods) may lead to an overhead due to the
communication costs between the CPU and the GPU (if the sequential
part of the algorithm is placed on CPU). In order to circumvent
this issue, many authors have implemented the entire algorithm on
GPU~\cite{pospichal10}. On the other hand, as several solutions
may be mapped with the same thread block in the GPU, shared memory
could not be used and all the data should be placed on global
memory.
%pheromone update data parallelism matrices

\subsubsection*{GPU thread mapping  for algorithmic-level parallelism}
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).

The implementation of the multistart model is based on two
different mapping strategies~\cite{luongMultiStart, audreyANT}:
(1) each Local Search (LS) is associated to a unique thread and
(2) each solution (from multiple neighborhoods) is associated to a
unique thread. The first mapping strategy (one thread per LS)
presents a big drawback: the number of LS to use should be very
large to saturate the GPU and cover the memory access latency. On
the other hand, the evaluation of the neighborhood of a single LS
by one thread is time intensive. Furthermore, shared memory could
not be used to store the huge data generated by the different
neighborhoods. In the second mapping strategy, the LS algorithms
are placed on CPU and the neighborhood evaluation of each LS is
parallelized on GPU using per-thread mapping strategy (one thread
per solution). This consists of a hierarchical parallel model
combining algorithmic-level parallelism (multistart) with
iteration-level
parallelism (master-worker).

In the island model, the same mapping is used in all the reviewed
works~\cite{tsutsuiGAQAP, nowotniak, maitre2009}: each
subpopulation (island) is associated to one thread block in the
GPU. A second level of parallelism is used inside the block to
parallelize the evaluation step of the local population.
Migrations are always performed through the global memory as
interblock communications  are impossible in CUDA. The first
advantage of this hierarchical implementation is that it allows
the occupation of a large number of threads even for medium
population sizes. The second advantage consists of using shared
memory to store subpopulations to benefit from the low latency of
shared memory.

\section{Frameworks for metaheuristics on GPUs}
\label{ch8:sec:frameworks}

After the first pioneering works of metaheuristics on graphics
processing units, the next challenge is to provide easy-to-use
frameworks and libraries for rapid development of metaheuristics
on GPUs. Although the works on this subject are not yet mature and
do not cover the main metaheuristic algorithms, we will present
the only three works to our knowledge, which propose open source
frameworks for developing
metaheuristics on GPUs.

The three frameworks reviewed in this section are
PUGACE\index{GPU-based frameworks!PUGACE}~\cite{pugace},
ParadisEO-MO-GPU\index{GPU-based
frameworks!ParadisEO-MO-GPU}~\cite{paradiseoGPU}, and
libCUDAOptimize\index{GPU-based
frameworks!libCUDAOptimize}~\cite{libcuda}. PUGACE\index{GPU-based
frameworks!PUGACE} is a framework for implementing EAs on GPUs.
ParadisEO-MO-GPU is an extension of the framework
ParadisEO~\cite{paradiseo} for implementing s-metaheuristics on
GPU. Finally, libCUDAOptimize\index{GPU-based
frameworks!libCUDAOptimize} is a  library intended for the
implementation of p-metaheuristics on GPU. The three frameworks
are presented in more detail in the following.

\subsection{PUGACE: framework for implementing evolutionary computation on GPUs}
PUGACE is a generic  framework
for easy implementation of cellular evolutionary algorithms on
GPUs implemented using C and CUDA. It is based on the frameworks
MALLBA and JCell (a framework for cellular algorithms). The
authors justified the choice of cellular evolutionary
algorithm by the good
feedback found in the literature concerning its efficient
implementation on GPUs compared to other parallel models for EAs
(island, master-slave). The main standard evolutionary operators
are already implemented in PUGACE: different selection strategies, standard
crossover, and mutation operators (PMX, swap, 2-exchange,
etc.). Different problem encoding is also supported. The framework
is organized as a set of modules in which the different
functionalities  are separated as much as possible in order to
facilitate the extension of the framework. All
the functions and procedures that execute on GPU are implemented in the same file kernel.cu. 

The implementation strategy adopted on the GPU is as follows.
Population initialization is done on the CPU side and the
population is transferred to the GPU. On the GPU side, each
individual is associated to a unique CUDA thread. The function
evaluation and mutation are done on the GPU while selection and
replacement are maintained on the CPU. In order to avoid thread
divergence\index{GPU!thread divergence} appearing in
the same CUDA thread block at the crossover step (because of the
probability of application which may give different results from
one thread to the other), the decision of whether to apply a
crossover is taken at the block level and applied to all the
individuals within the block. It is the decision on the choose
of the cutting point for the crossover.

The framework is validated on standard benchmarks of the QAP. Speedups of $15.44 \times$  to $18.41
\times$ are achieved compared to a CPU implementation of a cEA
using population sizes rising from 512 to 1024 and from 1024 to
2048.

% The authors noticed an
% improvement on the speedup when the size of the population increases. This is best explained
% by the linear increasing of execution time on CPU side while a sublinear increasing is
% occuring on GPU side (due to a better exploitation of the GPU threads). Neverethless,
% dispite the positive impact of large populations on th parallel efficieny of the implementation
% on GPU, a study on the impact on the final solution quality when using such large populations
% is to be done. Indeed, too large populations on EAs may lead to high diversification of the
% search and poor intensification capabilities, which leads to poor results.
% For the moment, plateforms whith more than one CPU and one GPU are not yet supported in PUGACE\index{GPU-based frameworks!PUGACE}.
% No details on how to get PUGACE\index{GPU-based frameworks!PUGACE} code source are given in the paper.

\subsection{ParadisEO-MO-GPU}

Melab et al.~\cite{paradiseoGPU} developed a reusable
framework ParadisEO-MO-GPU\index{GPU-based
frameworks!ParadisEO-MO-GPU} for parallel local search
metaheuristics (s-metaheuristics) on GPUs. It focuses on the
iteration-level
parallel model of s-metaheuristics which consists of exploring in
parallel the neighborhood of a problem solution on GPU. The
framework, implemented using C++ and CUDA, is an extension of the
ParadisEO~\cite{paradiseo} framework previously developed by the
same team for parallel and distributed metaheuristics on both
dedicated parallel hardware platforms and computational grids. The
objective of this framework is to facilitate the development of
GPU-based metaheuristics providing a transparent use of GPUs to
users who are unfamiliar with advanced features of all
parallelization techniques and deployment on GPUs. The framework
allows one to efficiently manage the hierarchical organization of
the memories (latencies and sizes) of the GPU and its
communication with the CPU as well as the minimizing of the user
involvement in its management.

\begin{figure}[h]
\centerline{\includegraphics[width=0.8\textwidth]{Chapters/chapter9/figures/paradiseoGPU.pdf}}
\caption{The skeleton of the ParadisEO-MO-GPU\index{GPU-based
frameworks!ParadisEO-MO-GPU}.}
%[Comparison of number of cores in a CPU and in a GPU]
\label{ch1:fig:paradiseoGPU}
\end{figure}

The framework is based on a master-worker model where the master
is the CPU and the workers are threads executed by the processing
cores of the GPU. The CPU executes the serial part of the
metaheuristic and sends only the current solution to the GPU to
minimize the transfer cost. The GPU, on its side, generates the
neighboring of the received solution and evaluates them at each
iteration. All the threads execute the same kernel and according
to a static mapping table between the threads and the neighbors
where each thread is associated with exactly one neighbor
evaluation. After all the neighborhood is generated and evaluated,
it is sent back to the CPU which selects the best solution (See
Figure~\ref{ch1:fig:paradiseoGPU}).

\subsection{libCUDAOptimize: an open source library of GPU-based metaheuristics}
\index{GPU-based frameworks!libCUDAOptimize}
LibCUDAOptimize~\cite{libcuda} is a C++/CUDA open source library
for the design and implementation of metaheuristics on GPUs. Until
now, the metaheuristics supported by LibCUDAOptimize are: scatter
search, differential evolution, and particle swarm
optimization\index{metaheuristics!particle swarm optimization}.
Nevertheless, the library is designed in such a way to allow
further extensions for other metaheuristics and it is still in
development phase by the authors. The parallelization strategy
adopted on GPU is principally based on fitness evaluation. The
sequential structure of the optimization algorithms is maintained
on CPU.

\section{Case study: accelerating large neighborhood LS method on GPUs for solving the Q3AP}
\label{ch8:sec:case}

In this case study, a large neighborhood GPU-based local
search\index{GPU!based local search}
method for solving the Quadratic 3-dimensional Assignment Problem
(Q3AP) will be presented. The local search method is an Iterated
Local Search\index{metaheuristics!iterated local search}
(ILS)~\cite{stutzle2006ILSforQAP} using an embedded
TS algorithm. The ILS principle
consists of executing iteratively the embedded local search, each
iteration which starts from a disrupted local optima reached by
the previous local search process. The disruption heuristic is a
performance parameter of an ILS algorithm and should be
judiciously defined. A template of an
ILS algorithm is given by the Algorithm~\ref{ch8:ils_algorithm_template}.

\begin{algorithm}[H]
    \SetAlgoLined
    %\KwData{this text}
    %\KwResult{how to write algorithm with \LaTeX2e }
   $s_{0}$=GenerateInitSol()\;
   $s^*$=LocalSearch($s_0$)\;
    \Repeat{a maximum number of iterations is reached}{
     $s'$=Perturbate($s^*,history$)\;
     $s^{*'}$=LocalSearch($s'$)\;
     $s^*$=AcceptationCriteria($s^*,s^{*'},history$)\;

    }
\caption{iterated local search template}
\label{ch8:ils_algorithm_template}
\end{algorithm}

%
%   \begin{algorithm}
%  % \label{ch8:ils_algorithm_template}
%   \caption{iterated local search\index{Metaheuristics!iterated local search} template}
%   \begin{algorithmic}[1]
%
%   \STATE $s_{0}$=\texttt{GenerateInitSol}();
%
%   \STATE $s^*$=\texttt{LocalSearch}($s_0$);
%
%   \REPEAT
%     \STATE $s'$=\texttt{Perturbate}($s^*,history$);
%     \STATE $s^{*'}$=\texttt{LocalSearch}($s'$);
%     \STATE $s^*$=\texttt{AcceptationCriteria}($s^*,s^{*'},history$);
%
%   \UNTIL {a maximum number of iterations is reached}
%
%   \end{algorithmic}
%
%   \end{algorithm}

\subsection{The quadratic 3-dimensional assignment problem}

The Q3AP is an extension of the well-known NP-hard QAP. The latter
is one of the most studied problems by the combinatorial
optimization community due to its wide range of practical
applications (site planning, schedule problem, computer-aided
design, etc.) and to its computational challenges since it is
considered as one of the most computationally difficult
optimization problems.

The Q3AP was first introduced by William P. Pierskalla in
1967~\cite{Pierskalla1967Q3AP} and, unlike the QAP, the Q3AP is a
less studied problem. Indeed, the Q3AP was revisited only this
past year and has recently been used to model some advanced
assignment problems such as the symbol-mapping problem encountered
in wireless communication systems and described
in~\cite{Hahn2008q3ap}. The Q3AP optimization problem can be
mathematically expressed as follows:

% debut -----------------------------------------------------------------------------
\begin{equation*} min \left \{
\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}
\quad+ \right . \end{equation*}

\begin{equation} \label{Q3APFormulation} \left .
\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\sum_{l=0}^{n-1}b_{ijl}x_{ijl}
\right \} \end{equation}

where

\begin{eqnarray}
  X=(x_{ijl})\in I \cap J \cap L, \label{Q3APConstraints1}\\
  x_{ijl}\in \{0,1\}, \quad i,j,l=0,1,...,n-1 \label{Q3APConstraints2}
\end{eqnarray} $I$, $J$, and $L$ sets are defined as follows:
\begin{displaymath} I=\lbrace X=(x_{ijl}):
\sum_{j=0}^{n-1}\sum_{l=0}^{n-1}x_{ijl}=1, \quad
i=0,1,...,n-1\rbrace \mathrm{;} \end{displaymath}
\begin{displaymath} J=\lbrace X=(x_{ijl}):
\sum_{i=0}^{n-1}\sum_{l=0}^{n-1}x_{ijl}=1, \quad
j=0,1,...,n-1\rbrace \mathrm{;} \end{displaymath}
\begin{displaymath} L=\lbrace X=(x_{ijl}):
\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}x_{ijl}=1, \quad
l=0,1,...,n-1\rbrace  \end{displaymath}
% fin------------------------------------------------------------------------------

Other equivalent formulations of the Q3AP can be found in the literature. A particularly useful one is the
\textit{permutation-based formulation} wherein the Q3AP can be expressed as follows:

 \begin{equation*}
 min\left\{ f(\pi_1,\pi_2) =
 \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
 + \right.
\end{equation*}

\begin{equation}
 \left. \sum_{i=0}^{n-1}b_{i\pi_{1(i)}\pi_{2(i)}}\right
 \} \label{Q3APPermutationBasedFormulation}
\end{equation}
where $\pi_1$ and $\pi_2$ are permutations over the set
$\left\lbrace 0,1,\ldots,n-1\right\rbrace$. According to this
formulation, minimizing the Q3AP consists of finding a double
permutation $(\pi_1,\pi_2)$ which minimizes the objective function
(\ref{Q3APPermutationBasedFormulation}).

The Q3AP is proven to be an NP-hard problem since it is an
extension of the quadratic assignment problem and of the axial
 3-dimensional assignment problem which are both NP-hard problems. It is particularly computationally difficult since the
number of feasible solutions of an instance of size $n$ is $n! \times n!$.

\subsection{Iterated tabu search algorithm for the Q3AP}

To tackle large-sized instances of the Q3AP and speed up the
search process, a parallel ILS algorithm has been designed. The local search embedded in
the ILS is a TS. A TS
procedure~\cite{Glover1989TS} starts from an initial feasible
solution and tries, at each step, to move to a neighboring
solution that minimizes the fitness (for a minimization case). If
no such move exists, the neighbor solution that less degrades the
fitness is chosen as a next move. This enables the TS process to
escape local optima. However, this strategy may generate cycles,
i.e., previous moves can be selected again. To avoid cycles, the
TS manages a short-term memory that contains the moves that have
been recently performed. A TS
template is given by Algorithm \ref{TS_pseudo_code}.
%
%  \begin{algorithm}
%  \label{TS_pseudo_code}
%  \caption{Tabu search template}
%  \begin{algorithmic}[1]
%  \STATE GenerateInitSol($s_0$);
%  \STATE TabuList=$\phi$;
%  \STATE $t=0$;
%  \REPEAT
%          \STATE $m(t)$ = SelectBestMove($s(t)$);
%          \STATE $s(t+1)$ = ApplyMove($m(t),s(t)$);
%          \STATE TabuList = TabuList $\bigcup \{m(t)\}$;
%          \STATE $t = t + 1$;
%  \UNTIL{a maximum number of iterations is reached}
%
%  \end{algorithmic}
%  \end{algorithm}

\begin{algorithm}[H]
    \SetAlgoLined
    %\KwData{this text}
    %\KwResult{how to write algorithm with \LaTeX2e }
    GenerateInitSol($s_0$)\;
    TabuList=$\phi$\;
    $t=0$\;
    \Repeat{a maximum number of iterations is reached}{
          $m(t)$ = SelectBestMove($s(t)$)\;
          $s(t+1)$ = ApplyMove($m(t),s(t)$)\;
          TabuList = TabuList $\bigcup \{m(t)\}$\;
          $t = t + 1$\;

    }
\caption{tabu search template}
\label{TS_pseudo_code}
\end{algorithm}


\subsection{Neighborhood functions for the Q3AP}

The neighborhood function is a crucial parameter in any local
search algorithm. Indeed, if the neighborhood function is not
adequate to the problem and/or does not consider the targeted
computing framework, any local search algorithm will fail to reach
good quality solutions of the search space.

Regarding the Q3AP, many neighborhood structures can be
considered. A basic neighborhood was proposed and investigated in
different works of the
literature~\cite{Hahn2008q3ap,DBLP:journals/ijfcs/LoukilMMTB12}
and consists of the set of all solutions (double-permutations)
generated from the current one by an exchange of two positions in
either the first ($\pi_1$) or the second ($\pi_2$) permutation.
This neighborhood that we denote by $N_b$ can be formalized as
follows:

\begin{equation}
\begin{array}{rl}
N_b(\pi_1,\pi_2)=\{&(\pi_1^{'},\pi_2^{'}): \pi_1^{'}[k]=\pi_1[l], \pi_1^{'}[l]=\pi_1[k]\\
                       & 0\leq k \neq l < n;\\
                & \pi_1^{'}[i]=\pi_1[i], \quad 0\leq i\neq k,l < n;\\
                & \pi_2^{'}[j]=\pi_2[j], \quad 0\leq j<n\\ \} \\
\bigcup \\
\{& (\pi_1^{'},\pi_2^{'}): \pi_2^{'}[k]=\pi_2[l], \pi_2^{'}[l]=\pi_2[k]\\
                & 0\leq k \neq l < n;\\
                & \pi_2^{'}[i]=\pi_2[i], \quad 0\leq i < n;\\
                & \pi_1^{'}[j]=\pi_1[j], \quad 0\leq j\neq k,l < n\\
\}& \end{array}
\end{equation}
where $(\pi_1,\pi_2)$ is the current solution and $n$ is the size
of the Q3AP instance. The size $|N_b|$ of such neighborhood is
equal to

\begin{equation}
 |N_b| = n \times (n-1)
\end{equation}


In our GPU-based implementations, a large-sized neighborhood
structure has been used for experimentation. In fact, theoretical
and experimental studies have shown that the use of large
neighborhood structures may improve the effectiveness of LS
algorithms~\cite{Ahuja:2007:VLN:1528422.1528438}. However, for
such a neighborhood, the generation/evaluation step of an LS
becomes a time-consuming task and may dramatically increase the
computational time of the LS process. This
justifies the use of intensive data-parallelism provided by GPUs where all neighboring solutions may be concurrently evaluated. 

The considered large-sized neighborhood consists of swapping two
positions in both permutations $\pi_1$ and $\pi_2$. This
neighborhood structure, $N(\pi_1,\pi_2)$, can be formalized as
follows:

\begin{equation} \begin{array}{rl}
N(\pi_1,\pi_2)=\{&(\pi_1^{'},\pi_2^{'}): \pi_1^{'}[k]=\pi_1[l], \pi_1^{'}[l]=\pi_1[k],\\
                & \pi_2^{'}[r]=\pi_2[s], \pi_2^{'}[s]=\pi_2[r]\\
                & 0\leq k \neq l < n,0\leq r \neq s < n;\\
                & \pi_1^{'}[i]=\pi_1[i], \quad 0\leq i\neq k,l < n;\\
                & \pi_2^{'}[j]=\pi_2^{'}[j], \quad 0\leq j\neq r,s < n\\

\}& \end{array}
\end{equation}

So, for a given Q3AP instance of size $n$, the size $|N|$ of the advanced neighborhood set can be expressed by the following formula:

\begin{equation}
  |N|=\left(\frac{n \times (n-1)}{2}\right )^2
\end{equation}

\subsection{Design and implementation of a GPU-based iterated tabu search algorithm for the Q3AP} \label{ch8:ITS-Q3APSection}

The use of GPU devices to speed up the search process of local
search methods is not a straightforward task. It requires one to
consider, at the same time, the characteristics and underlying
issues of the GPU architecture and the parallel models of LS
methods. The main challenges that must be faced when designing a
local search algorithm are the efficient distribution of the
search process between the CPU and the GPU minimizing the data
transfer between them, the hierarchical memory
management\index{GPU!memory management} and the
capacity constraints of GPU memories, and the thread
synchronization. All these issues must be regarded when designing
parallel LS models to allow
solving of large scale optimization problems on GPU architectures.

To go back to our problem (i.e., Q3AP), we propose in
Algorithm~\ref{ch8:algoITS} an iterated tabu
search on GPU (GPU-ITS). The
parallel model is in agreement with the
iteration-level
parallel model of LS methods presented in Section
\ref{ch8:sec:paraMeta} (Fig. \ref{ch8:fig:paraMeta}). This
algorithm can be seen as a cooperative model between the CPU and
the GPU. Indeed, the GPU is used as a coprocessor in a synchronous
manner. The resource-consuming part, i.e., the generation and
evaluation kernel, is calculated by the GPU and the rest of the LS
process is handled by the CPU. First of all, at initialization
stage, memory allocations on GPU are made; the input matrices
(distance and flow matrices) and the initial candidate solution of
the Q3AP must be allocated (lines 4 and 5). Since GPUs require
massive computations with predictable memory accesses, a structure
has to be allocated for storing all the neighborhood fitnesses at
different addresses (line 6). Second, the matrices and the initial
candidate solution have to be copied on the GPU (lines 7 and 8).
We notice that the input matrices are read-only structures and do
not change during all the execution of the LS algorithm.
Therefore, their associated memory is copied only once during all
the execution. Third, comes the parallel
iteration-level,
in which each neighboring solution is generated, evaluated, and
copied into the neighborhood fitnesses structure (from lines 10 to
14). Fourth, since the order in which candidate neighbors are
evaluated is undefined, the neighborhood fitnesses structure has
to be copied to the host CPU (line 15). Then the selection
strategy is applied to this structure (line 16): the exploration
of the neighborhood fitnesses structure is done by the CPU.
Finally, after a new candidate has been selected, this information
is copied to the GPU (line 18). The process is repeated until a
given number of iterations has been reached.

% \begin{algorithm}
% \caption{Template of an iterated tabu search\index{Metaheuristics!tabu~search} on GPU for solving the Q3AP} \label{LSGPU}
%  \begin{algorithmic}[1]
%  \STATE Choose an initial solution
%  \STATE Evaluate the solution
%  \STATE Initialize the tabu list
%  \STATE Allocate the Q3AP input data on GPU device memory
%  \STATE Allocate a solution on GPU device memory
%  \STATE Allocate a neighborhood fitnesses structure on GPU device memory
%  \STATE Copy the Q3AP input data on GPU device memory
%  \STATE Copy the solution on GPU device memory
%  \REPEAT
%      \FOR{each generated neighbor on GPU}
%          \STATE Incremental evaluation of the candidate solution
%          \STATE Insert the resulting fitness into the neighborhood fitnesses structure
%      \ENDFOR
%      \STATE Copy the neighborhood fitnesses structure on CPU host memory
%      \STATE Select the best admissible neighboring solution
%      \STATE Update the tabu list
%      \STATE Copy the chosen solution on GPU device memory
%  \UNTIL{a maximum number of iterations reached}
%  \end{algorithmic}
%  \end{algorithm}

\begin{algorithm}[H]
    \SetAlgoLined
    %\KwData{this text}
    %\KwResult{how to write algorithm with \LaTeX2e }
    Choose an initial solution\;
    Evaluate the solution\;
    Initialize the tabu list\;
    Allocate the Q3AP input data on GPU device memory\;
    Allocate a solution on GPU device memory\;
    Allocate a neighborhood fitnesses structure on GPU device memory\;
    Copy the Q3AP input data on GPU device memory\;
    Copy the solution on GPU device memory\;
    $t=0$\;
    \Repeat{a maximum number of iterations is reached}{
        \For{each generated neighbor on GPU}{
          Incremental evaluation of the candidate solution\;
          Insert the resulting fitness into the neighborhood fitnesses structure\;
        }
        Copy the neighborhood fitnesses structure on CPU host memory\;
        Select the best admissible neighboring solution\;
        Update the tabu list\;
        Copy the chosen solution on GPU device memory\;
    }
\caption{template of an iterated tabu search on GPU for solving the Q3AP}
\label{ch8:algoITS}
\end{algorithm}


\subsection {Experimental results}

In this section, some experimental results related to the approach
presented in Section \ref{ch8:ITS-Q3APSection} are reported. We
recall that the approach is a GPU-based iterated tabu
search (GPU-ITS) method
consisting in an iterated local search (ILS) embedding a tabu
search (TS) and where the
generation/evaluation step of the TS process is executed on GPU.
The ILS is used to improve the quality of successive local optima
provided by TS methods. This is achieved by perturbing the local
optima reached by the current TS process and reconsidering it as
initial solution of the following TS process. Regarding our
algorithm, the applied perturbation is a random number $\mu $ of
swaps in either the first or the second permutation where $\mu \in
[2:n]$ ($n$ is the instance size).

Experiments have been carried out on a node of the Chirloute
cluster of the Lille site. This is one of the 10 sites that
currently make up Grid5000~\cite{grid5000}, the French
experimental computational grid. A Chirloute cluster node consists
of an Intel Xeon E5620 CPU and a NVIDIA Tesla Fermi M2050 (448
cores) GPU type. The number of ILS iterations and the number of TS
iterations were set respectively to 20 and 500. The tabu list size
has been initalized to $\frac{m}{4}$, $m$ being the size of the
neighborhood set.

Table \ref{ch8:ITSQ3APResults} reports the obtained results for
the GPU-ITS using our large-sized neighborhood structure. The
method has been tested on 5 Q3AP instances derived from the QAP
Nugent instances in QAPLIB~\cite{burkard1991qaplib}. The average
time measurement for 20 executions is reported in seconds and
acceleration factors compared to a standalone CPU are also
considered. The algorithm is stopped when a maximum number of
iterations has been reached or when the optimal (or best known)
value has been discovered. Average and max values of the
evaluation function of the parallel GPU version have been
measured. The number of successful tries (hits) and the average
number of ILS iterations to converge to the optimal/best known
value are also represented. The associated standard deviation for
each average measurement is shown in small type characters.

\begin{table}
%\tiny
%\footnotesize
\small
\begin{tabular}{|l|r|r|r|r|r|r|r|r|} \hline
\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} \\
& \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
Nug12-d & $580$   & $615.4$    & $744$   & $35$\% & $87.7$ & $2.5$ & $34.7 \times$& $16$ \\
    &   & \tiny{$41.7$}    & && \tiny{$40.9$} & \tiny{$0.9$} & & \tiny{$6.3$} \\ \hline
Nug13-d & $1912$  & $1985.4$   & $2100$  & $20$\% & $209.2$ & $3.3$ & $63.5 \times$ & $17$ \\
&  & \tiny{$51.0$}   &  & & \tiny{$1.3$} & \tiny{$1.0$} &  & \tiny{$5.6$} \\ \hline
Nug15-d & $2230$  & $2418.1$  & $2580$  & $30$\% & $305.5$ & $5.2$ & $58.8 \times$ & $17$ \\
 &  & \tiny{$135.3$}  & & & \tiny{$164.5$} & \tiny{$1.3$} & & \tiny{$5.0$} \\ \hline
Nug18-d & $17836$ & $18110.2$ & $18506$ & $10$\% & $1375.9$ & $12.8$ & $107.4 \times$ & $19$ \\
& & \tiny{$157.8$} & & & \tiny{$123.5$} & \tiny{$2.6$} & & \tiny{$4.2$} \\ \hline
Nug22-d & $42476$ & $43282.1$ & $44140$ & $25$\% & $4506.5$ & $32.7$ & $137.8 \times$ & $18$ \\
& & \tiny{$529.6$} & & & \tiny{$341.1$} & \tiny{$6.6$} & & \tiny{$4.0$} \\ \hline
\end{tabular}
\caption{Results of the GPU-based iterated tabu search for
different Q3AP instances.} \label{ch8:ITSQ3APResults} % \center %Shashi
\end{table}

%\begin{table}
%\tiny
%%\footnotesize
%\caption{Results of the GPU-based iterated tabu search for different Q3AP instances}
%\label{ch8:ITSQ3APResults} \center
%\begin{tabular}{|c|c|c|c|c|c|c|c|c|} \hline
%Instance & Opt./BKV\footnotemark & Average val. & Max val. & Hits & CPU time & GPU time & Accel. & ITS iter.\\ \hline
%Nug12-d & $580$   & $615.4_{41.7}$    & $744$   & $35$\% & $87.7_{40.9}$ & $2.5_{0.9}$ & $\times 34.7 $& $16_{6.3}$ \\ \hline
%Nug13-d & $1912$  & $1985.4_{51.0}$   & $2100$  & $20$\% & $209.2_{1.3}$ & $3.3_{1.0}$ & $\times 63.5$ & $17_{5.6}$ \\ \hline
%Nug15-d & $2230$  & $2418.1_{135.3}$  & $2580$  & $30$\% & $305.5_{164.5}$ & $5.2_{1.3}$ & $\times 58.8 $ & $17_{5.0}$ \\ \hline
%Nug18-d & $17836$ & $18110.2_{157.8}$ & $18506$ & $10$\% & $1375.9_{123.5}$ & $12.8_{2.6}$ & $\times 107.4$ & $19_{4.2}$ \\ \hline
%Nug22-d & $42476$ & $43282.1_{529.6}$ & $44140$ & $25$\% & $4506.5_{341.1}$ & $32.7_{6.6}$ & $\times 137.8$ & $18_{4.0}$ \\ \hline
%\end{tabular}
%\end{table}

\footnotetext{Best known value.}

Regarding the execution time, the generation and evaluation of the
neighborhood in parallel on GPU provides an efficient way to
speedup the search process in comparison with a single CPU. In
fact, for the smaller instance Nug12-d, the GPU version starts to
be faster than the CPU one (acceleration factor of $34.7~\times$).
As long as the problem size increases, the speedup grows
significantly (up to $137.8 \times$ for the Nug22-d instance).
This means that the use of GPU provides an efficient way to deal
with large neighborhoods. Indeed, the proposed neighborhoods were
unpractical in terms of single CPU computational resources for
large Q3AP instances such as Nug22-d (estimated to around $2$
hours per run). So, implementing this algorithm on GPU has allowed
the exploitation of  parallelism in such neighborhood and improved
the robustness/quality of provided solutions.

\section{Conclusion}
\label{ch8:conclusion}

This chapter has presented state-of-the-art GPU-based parallel
metaheuristics and
a case study on implementing large neighborhood local search
methods on GPUs for solving large benchmarks of the quadratic
three-dimensional assignment problem (Q3AP). 

Implementing parallel
metaheuristics on
GPU architectures poses new issues and challenges such as memory
management;  finding
efficient mapping strategies between tasks to parallelize; and the
GPU threads, thread divergence, and synchronization. Actually, most
of metaheuristics have been implemented on GPU using different
implementation strategies. In this chapter, a two-level
classification of the reviewed works has been proposed: design
level and implementation level. Design level regroups traditional
parallel models used for metaheuristics while implementation level
refers to the way those models are mapped to the GPU architecture.
This classification focuses mainly on  the mapping between the
metaheuristic tasks to parallelize and the GPU threads. Indeed,
the choice of a given mapping strategy strongly influences the
other challenges (memory usage, communication, thread
divergence).

\putbib[Chapters/chapter9/biblio9]