[book_gpu.git] / BookGPU / Chapters / chapter15 / ch15.tex

\chapterauthor{Pierre Fortin}{Laboratoire d'Informatique de Paris 6, University Paris 6}
\chapterauthor{Rachid Habel}{T\'el\'ecom SudParis}
\chapterauthor{Fabienne J\'ez\'equel}{Laboratoire d'Informatique de Paris 6, University Paris 6}
\chapterauthor{Jean-Luc Lamotte}{Laboratoire d'Informatique de Paris 6, University Paris 6}
\chapterauthor{Stan Scott}{School of Electronics, Electrical Engineering \& Computer Science,
The Queen's University of Belfast}

\newcommand{\fixme}[1]{{\bf #1}}

\chapter[Numerical validation and performance optimization on GPUs in atomic physics]{Numerical validation and performance optimization on GPUs of
an application in atomic physics} 
\label{chapter15}

\section{Introduction}\label{ch15:intro}
As described in Chapter~\ref{chapter1}, GPUs are characterized by hundreds 
of cores and theoretically perform one order of magnitude better than CPUs.  
An important factor to consider when programming on GPUs
 is the cost of
data transfers between CPU memory and GPU memory. Thus, to have good
performance on
GPUs, applications should be coarse-grained and have a high arithmetic
intensity 
($i.e.$ the ratio of arithmetic operations to memory operations). 
Another important aspect of GPU programming is that floating-point
operations are preferably performed in single precision, if the
 validity of results is not impacted by that format.   
The GPU compute power for floating-point operations is indeed greater in
single precision than in double precision.  
The peak performance ratio between single precision and double
precision varies for example for NVIDIA GPUs from $12$ for the first Tesla
GPUs (C1060), 
to $2$ for the Fermi GPUs (C2050 and C2070)  
and to $3$ for the latest Kepler architecture (K20/K20X).  
As far as AMD GPUs are concerned, the latest AMD GPU (Tahiti HD 7970)
presents a ratio of $4$.  
Moreover, GPU internal memory accesses and CPU-GPU data transfers are
faster in single precision than in double precision,  
because of the different format lengths. 

This chapter describes the deployment on GPUs of PROP, a program of the
2DRMP~\cite{FARM_2DRMP,2DRMP} suite which models electron collisions
with H-like atoms and ions at intermediate energies. 2DRMP operates successfully on serial
computers, high performance clusters and supercomputers.  The primary
purpose of the PROP program is to propagate a global
R-matrix~\cite{Burke_1987}, $\Re$, in the two-electron configuration
space.
The propagation needs to be performed for all collision energies, 
for instance hundreds of energies,
which are independent.
Propagation equations are dominated by matrix multiplications involving sub-matrices of $\Re$.
However, the matrix multiplications are not
straightforward in the sense that $\Re$ dynamically changes the designation of its rows and
columns and increases in size as the propagation proceeds \cite{VECPAR}.

In a preliminary investigation PROP was selected by GENCI\footnote{GENCI: Grand Equipement National
  de Calcul Intensif, \url{www.genci.fr}} and
CAPS\footnote{CAPS is a software company providing products and solutions
  for manycore application programming and deployment,
  \url{www.caps-entreprise.com}},
following their first call for projects in 2009-2010 
aimed at
deploying applications on hybrid systems based on GPUs.
First CAPS  
recast the propagation equations with larger matrices.  
For matrix products the GPU performance gain over CPU increases indeed 
with the matrix size, since the 
CPU-GPU transfer overhead becomes less significant and since CPUs are
still more efficient for fine computation grains. 
Then, using HMPP\footnote{
HMPP or {\em CAPS compiler}, see: \url{www.caps-entreprise.com/hmpp.html}},
a commercial 
hybrid and parallel compiler, CAPS 
developed a version of PROP, in
which matrix multiplications are performed on
the GPU or the CPU, depending on the matrix size.
Unfortunately this partial GPU implementation of PROP does not offer
significant acceleration. 

The work described in this chapter, which is based on a study presented in \cite{PF_PDSEC2011}, aims at 
 improving PROP performance on 
GPUs by exploring two directions. First, because the original version of PROP is written
in double precision, 
we study the numerical stability of PROP in single precision. 
Second, we deploy  the whole
computation code of PROP on
GPUs to avoid the overhead generated by
data transfers 
and we propose successive improvements
(including a specific one to the Fermi architecture) 
in order to optimize the GPU code.




\section{2DRMP and the PROP program}
\label{s:2DRMP_PROP}
\subsection{Principles of R-matrix propagation}
2DRMP~\cite{FARM_2DRMP,2DRMP} is part of the CPC library\footnote{CPC:
Computer Physics Communications, 
\url{http://cpc.cs.qub.ac.uk/}}.  
It is a suite of seven 
programs aimed at creating virtual experiments on high performance and grid
architectures to enable the study of electron scattering from H-like
atoms and ions at intermediate energies.  The 2DRMP suite  uses the
two-dimensional $R$-matrix propagation approach~\cite{Burke_1987}. 
In 2DRMP the two-electron configuration  space  ($r_1$,$r_2$) is
divided into sectors. 
Figure~\ref{prop} shows the division of the two-electron configuration
space ($r_1$,$r_2$) into 4 vertical $strips$ representing 10 $sectors$. 
The key computation in 2DRMP, performed by the PROP program, is the
propagation of a global 
$R$-matrix, $\Re$, from sector to sector across the internal region, as shown in Fig.~\ref{prop}. 

\begin{figure}[h]
\begin{center}
\includegraphics*[width=0.65\linewidth]{Chapters/chapter15/figures/prop.pdf} 
\caption{\label{prop} Subdivision of the configuration space 
($r_1$,$r_2$) into a set of connected sectors.}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics*[width=0.8\linewidth]{Chapters/chapter15/figures/Domain.pdf} 
\caption{\label{domain} Propagation of the R-matrix from domain D to domain D'.}
\end{center}
\end{figure}

We consider the general situation in
Fig.~\ref{domain} where we assume that we already know
the global $R$-matrix, $\Re^{I}$, associated with the boundary defined
by edges 5, 2, 1 and 6 
in domain $D$ and we wish to
evaluate the new global $R$-matrix, $\Re^{O}$, associated with edges 5, 3, 4 and 6 
in domain $D'$ following propagation across subregion $d$.
Input edges are denoted by I (edges 1 and~2), output edges by O (edges 3 and 4) and
common edges by X (edges 5 and~6).
Because of symmetry, only the lower half of domains $D$ and $D'$ has to be considered. 
The global $R$-matrices, $\Re^{I}$ in domain $D$ and $\Re^{O}$ in
domain $D'$, can be written as:
\begin{equation}
\Re^{I} = \left(\begin{array}{cc}
      \Re_{II}^{I} & \Re_{IX}^{I}\\
      \Re_{XI}^{I} & \Re_{XX}^{I}
    \end{array}\right) 
\ \
\Re^{O} = \left(\begin{array}{cc}
      \Re_{OO}^{O} & \Re_{OX}^{O}\\
      \Re_{XO}^{O} & \Re_{XX}^{O}
    \end{array}\right).
\label{eq:RI_RO}
\end{equation}



From the set of local $R$-matrices, $\mathbf{R}_{ij}$ ($i,j\in \{1,2,3,4\}$)
associated 
with subregion $d$, we can define
\begin{subequations}
\begin{eqnarray}
    \mathbf{r}_{II} = \left(\begin{array}{cc}
      \mathbf{R}_{11} & \mathbf{R}_{12}\\
      \mathbf{R}_{21} & \mathbf{R}_{22}
    \end{array}\right), \label{eqaa} & 
    \mathbf{r}_{IO} = \left(\begin{array}{cc}
      \mathbf{R}_{13} & \mathbf{R}_{14}\\
      \mathbf{R}_{23} & \mathbf{R}_{24}
    \end{array}\right), \label{eqbb}\\
    \mathbf{r}_{OI} = \left(\begin{array}{cc}
      \mathbf{R}_{31} & \mathbf{R}_{32}\\
      \mathbf{R}_{41} & \mathbf{R}_{42}
    \end{array}\right), \label{eqcc} &
   \mathbf{r}_{OO} = \left(\begin{array}{cc}
      \mathbf{R}_{33} & \mathbf{R}_{34}\\
      \mathbf{R}_{43} & \mathbf{R}_{44}
    \end{array}\right),\label{eqdd}
\end{eqnarray}
\end{subequations}
where $I$ represents the input edges 1 and 2, and $O$ represents
the output edges 3 and 4 (see Fig.~\ref{domain}). 
The propagation across each sector is characterized by equations~(\ref{eq1}) to (\ref{eq4}).
\begin{subequations}
\begin{eqnarray}
\Re^{O}_{OO} &=& \mathbf{r}_{OO} - \mathbf{r}_{IO}^T (r_{II} + \Re^{I}_{II})^{-1}\mathbf{r}_{IO}, \label{eq1} \\
\Re^{O}_{OX} &=& \mathbf{r}_{IO}^T (\mathbf{r}_{II} + \Re^{I}_{II})^{-1}\Re^{I}_{IX},  \label{eq2} \\
   \Re^{O}_{XO} &=& \Re^{I}_{XI}(\mathbf{r}_{II} + \Re^{I}_{II})^{-1}\mathbf{r}_{IO},  \label{eq3} \\
   \Re^{O}_{XX} &=& \Re^{I}_{XX} - \Re^{I}_{XI}(\mathbf{r}_{II} +\Re^{I}_{II})^{-1}\Re^{I}_{IX}. \label{eq4}
\end{eqnarray}
\end{subequations}

The matrix inversions are not explicitly performed. To compute
$(r_{II} + \Re^{I}_{II})^{-1}\mathbf{r}_{IO}$ and $(\mathbf{r}_{II} + \Re^{I}_{II})^{-1}\Re^{I}_{IX}$,
two linear systems are solved.


\medskip

While equations (\ref{eq1})-(\ref{eq4}) can be applied to the
propagation across a general subregion two special situations should be
noted: propagation across a diagonal subregion and propagation across
a subregion bounded by the $r_{1}$-axis at the beginning of a new
strip.

In the case of a diagonal subregion, from symmetry considerations,
edge 2 is identical to edge 1 and edge 3 is identical to edge~4. 
Accordingly, with only one input edge and one output edge equations
(\ref{eqaa})-(\ref{eqdd}) become: 
\begin{subequations}
\begin{eqnarray}
\mathbf{r}_{II} = 2\mathbf{R}_{11}, \ 
\mathbf{r}_{IO} = 2\mathbf{R}_{14}, \label{eq4b}\\
\mathbf{r}_{OI} = 2\mathbf{R}_{41}, \ 
\mathbf{r}_{OO} = 2\mathbf{R}_{44}. \label{eq4d}
\end{eqnarray}
\end{subequations}
In the case of a subregion bounded by the $r_1$-axis at the beginning
of a new strip, we note that the input boundary $I$ consists of only
one edge. When propagating across the first subregion in the second
strip there is no common boundary $X$: in this case only equation
(\ref{eq1}) needs to be solved.

\medskip

Having obtained the global $R$-matrix $\Re$ on the boundary of the
innermost subregion (labeled $0$ in Fig.~\ref{prop}), $\Re$ is propagated across
each subregion in the order indicated in Fig.~\ref{prop},
working systematically from the
$r_1$-axis at the bottom of each strip across all subregions to the
diagonal. 



\subsection{Description of the PROP program}
\label{sec:PROP}

\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
 \hline
 \multirow{2}{0.09\linewidth}{\centering Data set} &
 \multirow{2}{0.15\linewidth}{\centering
  Local $R$-\\matrix size} &
 \multirow{2}{0.07\linewidth}{\centering Strips} &
 \multirow{2}{0.09\linewidth}{\centering Sectors} &
 \multirow{2}{0.19\linewidth}{\centering Final global \\$R$-matrix size} &
 \multirow{2}{0.15\linewidth}{\centering Scattering\\energies} \\
  & & & & & \\
  \hline
  Small & 90x90  & 4 & 10 & 360x360 & 6\\
 \hline
  Medium  & 90x90  & 4 & 10 & 360x360 & 64\\
 \hline
  Large  & 383x383  & 20 & 210 &  7660x7660 & 6\\
 \hline
  Huge  & 383x383  & 20 & 210 &  7660x7660 & 64\\ \hline
\end{tabular}
\caption{\label{data-sets}Characteristics of four data sets}
\end{center}
\end{table}

The PROP program computes the propagation of the $R$-matrix across the sectors of the internal region.  
Table~\ref{data-sets} shows four different 
data sets used in this study and highlights the principal parameters of the PROP program. 
 PROP execution can be described by Algorithm~\ref{prop-algo}. 
First, amplitude arrays and
correction data are read from data files generated by the preceding
program of the 2DRMP suite.  
Then, local $R$-matrices are constructed from amplitude arrays. 
Correction data is used to compute correction vectors added to the diagonal of the local 
$R$-matrices. The local $R$-matrices, together with the input $R$-matrix, 
 $\Re^{I}$,
computed on the previous sector, are used to compute the current output
$R$-matrix, 
$\Re^{O}$. 
 At the end of a sector evaluation,
the output $R$-matrix becomes  the input $R$-matrix 
for the next evaluation.  

\begin{algorithm}
\caption{\label{prop-algo}PROP algorithm}
\begin{algorithmic}
\FOR{all scattering energies}
 \FOR{all sectors}
 \STATE Read amplitude arrays
 \STATE Read correction data
\STATE Construct local $R$-matrices
\STATE From $\Re^{I}$ and local $R$-matrices, compute $\Re^{O}$
\STATE $\Re^{O}$ becomes $\Re^{I}$ for the next sector 
 \ENDFOR
 \STATE Compute physical $R$-Matrix 
\ENDFOR
\end{algorithmic}
\end{algorithm}


On the first sector, there is no input $R$-matrix yet. To bootstrap
the propagation, the first output $R$-matrix is constructed using only
one local $R$-matrix.  On the last sector, that is, on the boundary of
the inner region, a physical $R$-matrix corresponding to the output
$R$-matrix is computed and stored into an output file.

In the PROP program, sectors are characterized into four types,
depending on the computation performed: 
\begin{itemize}
\item the starting sector (labeled 0 in Fig.~\ref{prop})
\item the axis sectors (labeled 1, 3 and 6 in Fig.~\ref{prop})
\item the diagonal sectors (labeled 2, 5 and 9 in Fig.~\ref{prop})
\item the off-diagonal sectors (labeled 4, 7 and 8 in Fig.~\ref{prop}).
\end{itemize}


The serial version of PROP is implemented in Fortran~90 and uses
for linear algebra operations BLAS and LAPACK routines
which are fully optimized for x86 architecture.
This 
program 
serially propagates the $R$-matrix for
all scattering energies. 
Since the propagations for these different 
energies are independent, there also 
exists an embarrassingly parallel version of 
PROP 
that spreads the computations of 
several energies
among multiple CPU nodes via  
MPI. 



\subsection{CAPS implementation}
\label{caps}

In order to handle larger matrices, and thus obtain better GPU  speedup, CAPS  
recast equations (\ref{eq1}) to (\ref{eq4}) into one equation.
The output $R$-matrix $\Re^{O}$ defined by equation~(\ref{eq:RI_RO}) is now computed as follows. 
\begin{equation}\label{eq_CAPS_1}
\Re^{O} = \Re^{O^{\ \prime}} + U A^{-1} V, 
\end{equation}
\begin{equation}\label{eq_CAPS_2}
{\rm with} \
\Re^{O^{\ \prime}}= \left(\begin{array}{cc}
      \mathbf{r}_{OO} & 0\\
      0 & \Re^I_{XX}    \end{array}\right), \
U= \left(\begin{array}{c}
      \mathbf-{r}_{IO}^{T}\\
      \Re^I_{XI}    \end{array}\right), 
\end{equation}
\begin{equation}\label{eq_CAPS_3}
A= \mathbf{r}_{II} + \Re^I_{II} \ {\rm and}  \
V= (\mathbf{r}_{IO}\ \ \ -\Re^I_{IX}).
\end{equation}

To compute $W=A^{-1}V$, no matrix inversion is performed. The matrix
system $AW=V$ is solved. 
This reimplementation of PROP reduces the number of equations to be
solved and the number of matrix copies for evaluating each sector.
For instance, for an off-diagonal sector, 
copies fall from 22 to 5, matrix multiplications from 4 to~1 and calls
to a linear equation solver from 2 to 1. 

To implement this version, CAPS 
used HMPP, a 
commercial hybrid and parallel compiler, 
based on compiler directives like the new OpenACC standard\footnote{See: \url{www.openacc-standard.org}}.  
If the matrices are large enough (the limit sizes are experimental parameters), 
they are multiplied on the GPU, otherwise on the CPU. 
CAPS 
 used the MKL BLAS implementation on an Intel Xeon
x5560 quad core CPU (2.8 GHz) 
and the CUBLAS library (CUDA 2.2) on one Tesla C1060 GPU. 
On the large data set (see Table~\ref{data-sets}), CAPS 
 obtained a speedup of 1.15 for the GPU 
version over the CPU one (with multi-threaded MKL calls on the four
CPU cores). This limited gain in performance is mainly
due to the use of double precision computation 
and to the small or medium sizes of most matrices.
For these matrices, the computation gain on  
the GPU is indeed
strongly affected by the overhead 
generated by transferring these matrices from  
the CPU memory to the GPU memory to perform each matrix multiplication and then
transferring the result back to the CPU memory. 

Our goal is to speedup PROP more significantly by porting the whole
code to the GPU and therefore avoiding 
the 
intermediate data transfers between
the host (CPU) and the GPU. We will also study the
stability of PROP in single precision because 
single precision computation is faster on the GPU  
and CPU-GPU data transfers are twice as fast as those performed in
double precision. 



\section{Numerical validation of PROP in single precision}
\label{single-precision}

\begin{comment}
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|}
  \hline
   relative error interval & \# occurrences \\
  \hline
   [0, 1.E-8) & 18 \\
  \hline
   [1.E-8, 1.E-6) & 1241 \\
  \hline
   [1.E-6, 1.E-4) & 48728 \\
  \hline
   [1.E-4, 1.E-2) & 184065 \\
  \hline
   [1.E-2, 1) & 27723 \\
  \hline
   [1, 100) & 304 \\
  \hline
   [100, $+\infty$) & 1 \\
  \hline
\end{tabular}
\end{center}
\caption{\label{sp-distrib}Error distribution for medium case in single precision}
\end{table}
\end{comment}


Floating-point input data, computation and output data of PROP are 
originally in double precision format.
PROP produces a standard $R$-matrix H-file \cite{FARM_2DRMP}
 and a collection of Rmat00X files (where X
ranges from 0 to the number of scattering energies - 1)
holding the physical R-matrix for each 
energy.
The  H-file  and the Rmat00X files are binary input files of the FARM program \cite{FARM_2DRMP}
(last program of the 2DRMP suite).
Their text equivalent are the  prop.out 
and the prop00X.out files. 
To study the validity of PROP results in  single precision,
first,
reference results are 
 generated by running the serial version of PROP in double precision.
Data used in the most costly computation parts are read from input files in
double precision format and then 
cast to single precision format.
PROP results  (input of FARM) are computed in single precision and  written
into files in double precision. 

\subsection{Medium case study}
\begin{figure}[h]
\begin{center}
\includegraphics*[width=0.9\linewidth]{Chapters/chapter15/figures/error.pdf} 
\caption{\label{fig:sp-distrib} Error distribution for medium case in single precision}
\end{center}
\end{figure}

The physical $R$-matrices, in
the prop00X.out files, are compared to the
reference ones for the medium case (see Table~\ref{data-sets}). 
 The relative
error distribution is
given in Fig.~\ref{fig:sp-distrib}. 
We focus on the largest errors. 
\begin{itemize}
\item Errors greater than 100: the only impacted value is of order 1.E-6
and is negligible compared to the other ones 
in the same prop00X.out file.

\item Errors between 1 and 100: the values corresponding to the
  largest errors are of order 1.E-3 and are negligible compared to
  the majority of the other values which range between 1.E-2 and
  1.E-1.

\item Errors between 1.E-2 and 1: the largest errors ($\ge$ 6\%)
  impact values the order of magnitude of which is at most 1.E-1.
  These values are negligible. 
  Relative errors of approximately 5\% impact values the order of
  magnitude of which is at most 1.E2. 
  For instance, the value 164 produced by the reference version of
  PROP becomes 172 in the single precision version.

\end{itemize}

To study the impact of the single precision version of PROP on the
FARM program, the cross-section
results files corresponding to 
transitions 
{1s1s}, 
{1s2s}, {1s2p}, {1s3s}, {1s3p}, {1s3d},
{1s4s}, {2p4d} are compared to the reference ones.  
Table~\ref{sp-farm} shows that all cross-section files are impacted by
errors.  Indeed in the  {2p4d} file,  four relative errors are 
greater than one and the maximum relative error is 1.60. 
However the largest errors impact negligible values. For example, the maximum
error (1.60) impacts a reference value which is 4.5E-4.  The largest 
values are impacted by low errors. For instance, the maximum value
(1.16) is impacted by a relative error of the order 1.E-3. 

\begin{table}[t] 
\begin{center}
\begin{tabular}{|c|c||c|c|} \hline
  file & largest relative error & file & largest relative error\\ \hline
 {1s1s} & 0.02& {1s3p} & 0.11  \\ \hline
 {1s2s} & 0.06 &  {1s3d} &  0.22 \\ \hline
 {1s2p} &  0.08 & {1s4s} &  0.20  \\ \hline
 {1s3s} &  0.17 &2p4d & 1.60  \\ \hline
\end{tabular}
\caption{\label{sp-farm}Impact  on FARM  of the single precision version of PROP}
\end{center}
\end{table}

To examine in more detail the impact of PROP on FARM, 
cross sections above the ionization threshold (1 Ryd)
are compared in single and
double precision  for 
transitions amongst the 1s, \dots 4s, 2p, \dots 4p, 3d, 4d target states.  
This comparison is carried out by generating 45 plots.  In all the
 plots, results in single and double precision match except for few
 scattering energies which are very close to pseudo-state thresholds. 
For example Fig.~\ref{1s2p} and \ref{1s4d} present the scattering energies corresponding to the
{1s2p} and {1s4d} cross-sections computed in single and double precision.   For some cross-sections, 
increasing a threshold parameter from 1.E-4 to 1.E-3 in the FARM
program
results in energies close to threshold being avoided
 and therefore
the cross-sections in double and single precision match more
accurately. 
This is the case for instance for cross-section 1s2p (see Fig.~\ref{1s2p3}). 
However for other cross-sections (such as 1s4d) some problematic energies remain even if the 
threshold parameter in  the FARM
program is increased to 1.E-3 (see Fig.~\ref{1s4d3}).  A higher 
threshold parameter would be required for such cross-sections. 

\begin{figure}[t]
\centering
\subfigure[threshold = 1.E-4]{ 
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s2p.pdf}
   \label{1s2p}
 }
\subfigure[threshold = 1.E-3]{
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s2p3.pdf}
 \label{1s2p3}
 }
\label{fig:1s2p_10sectors}
\caption{1s2p cross-section, 10 sectors}
\end{figure}

\begin{figure}[t]
\centering
\subfigure[threshold = 1.E-4]{
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s4d.pdf}
   \label{1s4d}
 }
\subfigure[threshold = 1.E-3]{
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s4d3.pdf}
 \label{1s4d3}
 }
\label{fig:1s4d_10sectors}
\caption{1s4d cross-section, 10 sectors}
\end{figure}

As a conclusion, the medium case study shows that the execution of
PROP in single precision leads to a few inexact scattering energies to
be computed by the FARM program for some cross-sections.
Thanks to a suitable threshold parameter in the FARM program these problematic energies may possibly 
be skipped. 
Instead of investigating deeper the choice of such a parameter for the medium case, we analyze the 
single precision computation  in a more
realistic case in Sect.~\ref{huge}. 
\begin{comment}
The conclusion of the medium case study is that running PROP in single
precision gives relatively stable results provided that suitable
parameter values are used in the FARM program in order to skip the
problematic energies that are too close to the pseudo-state
thresholds.  To verify if this conclusion is still valid with a larger
data set, the single precision computation is analyzed in a more
realistic case in Sect.~\ref{huge}.
\end{comment}

\subsection{Huge case study}\label{huge}


\begin{figure}[t] 
  \centering
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s2pHT.pdf}
\caption{\label{1s2pHT}1s2p cross-section, threshold = 1.E-4, 210 sectors}
\end{figure}

\begin{figure}[t] 
  \centering
\includegraphics*[width=.76\linewidth]{Chapters/chapter15/figures/1s2pHT.pdf}
\caption{\label{1s4dHT}1s4d cross-section, threshold = 1.E-4, 210 sectors}
\end{figure}

We study here the impact on FARM of the PROP program run in
single precision for the huge case (see Table~\ref{data-sets}).
The cross-sections
corresponding to all
atomic target states 1s \dots 7i are explored, which
leads to 
406 comparison plots. 
It should be noted that in this case, over the same energy range above the ionization threshold, the density of pseudo-state thresholds is significantly increased compared to the medium case.
As expected, all the plots exhibit large differences between single and double
precision cross-sections. 
For example Fig.~\ref{1s2pHT} and  \ref{1s4dHT} present the 1s2p and 1s4d cross-sections computed in
single and double precision for the huge case.  
We can conclude that PROP in single precision gives invalid results 
for realistic simulation cases above the ionization threshold.
Therefore the  deployment of PROP on GPU, described in Sect.~\ref{gpu-implem},
has been carried out in double precision. 

\section{Towards a complete deployment of PROP on GPUs} 
\label{gpu-implem}

We now detail how PROP has been progressively deployed on
GPUs in double precision in order to avoid the
expensive memory transfers between the host and the GPU.
Different versions with successive improvements and optimizations are presented.
We use CUDA~\cite{CUDA_ProgGuide} for GPU programming, as well as the
CUBLAS~\cite{CUBLAS} 
and MAGMA \cite{MAGMA} libraries for linear algebra operations.
Since PROP is written in Fortran 90, {\em wrappers} in C are used to
enable calls to CUDA kernels from PROP routines.


\subsection{Computing the output $R$-matrix on GPU}
\label{gpu-RO}

\begin{figure}[h]
  \centering
  \includegraphics[width=0.7\linewidth]{Chapters/chapter15/figures/offdiagonal_nb.pdf}
  \caption{\label{offdiagonal} The six steps of an off-diagonal sector
    evaluation.}
\end{figure}

As mentioned in Algorithm~\ref{prop-algo}, evaluating a sector
mainly consists in constructing local $R$-matrices and in computing
one output $R$-matrix, $\Re^{O}$. In this first step of the porting
process, referred to as GPU V1\label{gpuv1},
we only consider the computation of $\Re^{O}$ on the GPU.
We distinguish the following six steps, related to equations
(\ref{eq_CAPS_1}), (\ref{eq_CAPS_2}) and (\ref{eq_CAPS_3}), and illustrated in
Fig.~\ref{offdiagonal} for an off-diagonal sector.

\begin{description}
\item[Step 1] (``Input copies''):~data are copied from $\Re^{I}$
  to temporary arrays ($A$, $U$, $V$) and to $\Re^{O}$.
  These copies, along with possible scalings or transpositions, are
  implemented as CUDA kernels which can be applied to two
  matrices of any size and starting at any offset. 
  Memory accesses are coalesced \cite{CUDA_ProgGuide} in order to
  provide the best performance for such memory-bound kernels.
\item[Step 2] (``Local copies''):~data are copied from
  local $R$-matrices to temporary arrays ($U$, $V$) and to $\Re^{O}$.
  Moreover data from local $R$-matrix
  $\mathbf{r}_{II}$ 
  is added to matrix $A$ (via a CUDA kernel) and zeroes are written in
   $\Re^{O}$  where required.
\item[Step 3] (``Linear system solving''):~matrix $A$ is factorized
  using the MAGMA DGETRF 
   routine and the result is stored in-place.
\item[Step 4] (``Linear system solving'' cont.):~the matrix system
 of linear equations  $AW$ = $V$ is solved using the MAGMA DGETRS 
routine. The  solution is stored in matrix $V$.
\item[Step 5] (``Output matrix product''):~matrix $U$
  is multiplied by matrix $V$ using the CUBLAS DGEMM 
  routine. The result is stored in a temporary matrix~$t$.
\item[Step 6] (``Output add''):~$t$ is added to $\Re^{O}$ (CUDA
  kernel).
\end{description}

All the involved matrices are stored in the GPU memory. Only the
local $R$-matrices are first constructed on the host and then sent
to the GPU memory, since these matrices vary from sector to sector.
The evaluation of the axis and diagonal sectors is similar.
However, fewer operations and copies are required because of
symmetry considerations \cite{2DRMP}.

\subsection{Constructing the local $R$-matrices on GPU}

\begin{figure}[t]
 \centering
  \includegraphics[width=0.7\linewidth]{Chapters/chapter15/figures/amplitudes_nb.pdf} 
 \caption{\label{amplitudes} Constructing the local $R$-matrix R34
 from the $j$ amplitude array associated with edge 4 and the $i$
 amplitude array associated with edge~3.}
\end{figure}

Local $R$-matrices are constructed using two three dimensional arrays,
$i$ and $j$. Each three dimensional array contains four
matrices corresponding to the surface amplitudes associated with the
four edges of a sector. Those matrices are named {\em amplitude arrays}.
 $j$ amplitude arrays are read from data files and $i$ amplitude arrays
are obtained by scaling each row of the $j$ amplitude arrays. 
The main part of the construction of a local $R$-matrix,
presented in Fig.~\ref{amplitudes},
is a matrix product between
one $i$ amplitude array and one transposed $j$ amplitude array 
which is performed by a single DGEMM 
BLAS call. 
In this version, hereafter referred to as GPU
V2\label{gpuv2}, $i$ and $j$ amplitude arrays are transferred to the
GPU memory and the required matrix multiplications are performed on
the GPU (via CUBLAS routines).


The involved matrices having medium sizes (either $3066 \times 383$ or
$5997 \times 383$),
performing these matrix multiplications
on the GPU is expected to be faster than on the CPU.
However, this implies a greater communication volume
between the CPU and the GPU
since the $i$ and $j$ amplitude arrays are larger than the local
$R$-matrices.
It can be noticed that correction data are also used in the
construction of a local $R$-matrix,
but this is a minor part in the
computation. However, these correction data also have to be
transferred from the CPU to the GPU for each sector.

\subsection{Scaling amplitude arrays on GPU}
It  
should be worthwhile to try to reduce the CPU-GPU data
transfers of the GPU V2, where the $i$ and $j$ amplitude arrays are
constructed on the host and then sent to the GPU memory for each sector. 
In this new version, hereafter referred to as GPU V3\label{gpuv3}, we
transfer only the $j$ amplitude arrays and the
required scaling factors (stored in one 1D array) to the GPU memory,
so that the $i$ amplitude arrays are
then directly computed on the GPU by multiplying the $j$ amplitude
arrays by these scaling factors (via a CUDA kernel).
Therefore, we save the transfer of four $i$ amplitude arrays on
each sector by transferring  only this 1D array of scaling factors. 
Moreover, scaling $j$ amplitude arrays is expected to be faster on the
GPU than on the CPU, thanks to the massively parallel architecture of
the GPU and thanks to its higher internal memory bandwidth.

\subsection{Using double-buffering to overlap I/O and computation}

\begin{figure}[t]
  \centering
  \includegraphics[width=0.8\linewidth]{Chapters/chapter15/figures/C1060_V3_IO_COMP.pdf} 
  \caption{\label{overlapping} Compute and I/O times for the GPU V3 on
    one C1060.}  
\end{figure} 

As described in Algorithm~\ref{prop-algo}, there are two main steps in
the propagation across a sector: reading  amplitude arrays
and correction data from I/O files and
evaluating the current sector.
Fig.~\ref{overlapping} shows the I/O times and the evaluation times
of each sector for the huge case execution (210 sectors, 20 strips) of the GPU V3 
on one C1060. 
Whereas the times required by the off-diagonal sectors are similar
within each of the 20 strips, 
the times for diagonal sectors of each strip
are the shortest ones, the times for the axis sectors being
intermediate.
The I/O times are roughly constant among all strips.
The evaluation time is equivalent to the I/O
time for the first sectors. But this evaluation time grows 
linearly with the strip number, and rapidly exceeds the I/O 
time.

It is thus interesting to use a double-buffering technique to overlap the
I/O time with the evaluation time:
for each sector, the evaluation of sector $n$ is performed
(on GPU) simultaneously with the reading of data for sector
$n+1$ (on CPU). This requires the duplication in the CPU memory of all the
data structures
used for storing data read from I/O files for each sector.
This version, hereafter referred to as GPU
V4\label{gpuv4}, uses POSIX threads. Two threads are
executed concurrently: an I/O thread that reads data from I/O files
for each sector, and a computation thread, dedicated to the propagation
of the global $R$-matrix, that performs successively for each sector
all necessary computations on GPU, 
as well as all required CPU-GPU data transfers.
The evaluation of a sector uses the data read for this sector as well
as the global $R$-matrix computed on the previous sector.
This dependency requires synchronizations between the I/O thread and
the computation thread which are implemented through standard POSIX
thread mechanisms.


\subsection{Matrix padding}
The CUBLAS DGEMM 
performance and the MAGMA DGETRF/DGETRS 
performance is reduced when the sizes (or
the leading dimensions) of the matrix are not multiples of the inner blocking size \cite{NTD10a}.
This inner blocking size can be 32 or 64, depending on the computation
and on the underlying  
GPU architecture \cite{MAGMA,NTD10b}. 
In this version (GPU V5\label{gpuv5}), 
the matrices are therefore padded with $0.0$ (and $1.0$ on the diagonal for the linear systems)
so that their sizes are 
multiples of 64.
This corresponds indeed to the optimal size for the matrix product on the
Fermi architecture \cite{NTD10b}. And as far as linear system solving is
concerned, all the matrices have sizes which are multiples of 383: we
therefore use padding to obtain multiples of 384 (which are 
again  multiples of 64). 
It can be noticed that this padding has to be performed dynamically
as the matrices increase in size during the propagation 
(when possible the
 maximum required storage space is however allocated only once in the
 GPU memory). 



\section{Performance results}
\subsection{PROP deployment on GPU}

\begin{table*}[ht]
\begin{center}
\begin{tabular}{|c||c|c||}
 \hline
  PROP version & \multicolumn{2}{c|}{Execution time} \\
  \hline
  \hline
  CPU version: 1  core & \multicolumn{2}{c|}{201m32s} \\
  \hline
  CPU version: 4  cores &  \multicolumn{2}{c|}{113m28s} \\
  \hline \hline
GPU version  & C1060 & C2050 \\
  \hline\hline
  GPU V1 (\S~\ref{gpuv1}) & 79m25s & 66m22s  \\
  \hline
  GPU V2 (\S~\ref{gpuv2}) & 47m58s & 29m52s \\
  \hline
  GPU V3 (\S~\ref{gpuv3}) & 41m28s & 23m46s \\
  \hline
  GPU V4 (\S~\ref{gpuv4}) & 27m21s & 13m55s\\
  \hline
  GPU V5 (\S~\ref{gpuv5}) & 24m27s & 12m39s  \\
  \hline
\end{tabular}
\caption{\label{table:time} 
Execution time of PROP on CPU and GPU}
\end{center}
\end{table*}

\begin{comment}
\begin{table*}[ht]
\begin{center}
\begin{tabular}{|c||c|c||}
 \hline
  PROP version & \multicolumn{2}{c|}{Execution time} \\
  \hline  \hline
CPU version & 1 core & 4 cores  \\\hline
& {201m32s} & {113m28s} \\ \hline  \hline
GPU version  & C1060 & C2050 \\
  \hline\hline
  GPU V1 (\ref{gpuv1}) & 79m25s & 66m22s  \\
  \hline
  GPU V2 (\ref{gpuv2}) & 47m58s & 29m52s \\
  \hline
  GPU V3 (\ref{gpuv3}) & 41m28s & 23m46s \\
  \hline
  GPU V4 (\ref{gpuv4}) & 27m21s & 13m55s\\
  \hline
  GPU V5 (\ref{gpuv5}) & 24m27s & 12m39s  \\
  \hline
\end{tabular}
\caption{\label{table:time} 
Execution time of the successive GPU versions}
\end{center}
\end{table*}
\end{comment}

\begin{figure}[h]
\centering
\subfigure[Speedup over 1 CPU core]{
\includegraphics*[width=0.76
        \linewidth]{Chapters/chapter15/figures/histo_speedup_1core.pdf}
   \label{fig:speedup_1core}
 }

\subfigure[Speedup over 4 CPU cores]{
\includegraphics*[width=0.76
        \linewidth]{Chapters/chapter15/figures/histo_speedup_4cores.pdf}
 \label{fig:speedup_4cores}
 }
\label{fig:speedup}
\caption{Speedup of the successive GPU versions.}
\end{figure}

Table~\ref{table:time} presents 
the execution times 
of PROP on CPUs and GPUs, 
each version solves the propagation equations in the 
form~(\ref{eq_CAPS_1}-\ref{eq_CAPS_3}) as proposed by CAPS. 
Fig.~\ref{fig:speedup_1core} (respectively \ref{fig:speedup_4cores})
shows the speedup of the successive GPU versions
over one CPU core (respectively four CPU cores). 
We use here Intel Q8200 quad-core CPUs (2.33 GHz), one C1060 GPU and
one C2050 (Fermi) GPU, located at 
 UPMC (Universit\'e Pierre et Marie Curie, Paris, France). 
As a remark, the execution times measured on the C2050 would be the same 
on the C2070 and on  the C2075, the only difference between these GPUs 
being their memory size and their TDP (Thermal Design Power). 
We emphasize that the execution times correspond to the
complete propagation for all six energies of the large case (see
Table~\ref{data-sets}), that is to say to the complete execution of
the PROP program.  
Since energies are independent, execution times for more energies
(e.g. the huge case) should be proportional
to those reported in Table~\ref{table:time}.  

These tests, which have been performed with CUDA 3.2, CUBLAS 3.2 and 
MAGMA 0.2, 
show that the successive GPU versions of PROP offer 
increasing, and at the end interesting, speedups.
More precisely, 
V2 shows that it is worth increasing slightly the
CPU-GPU communication volume in order to perform
large enough matrix products on the GPU. 
 This communication volume can fortunately be
reduced thanks to~V3,  
which also accelerates the computation of
amplitude arrays thanks to the GPU. 
The 
double-buffering technique implemented in V4
 effectively enables the overlapping of 
I/O operations with computation, while the 
padding implemented in V5 also improves the computation times.
It 
is noticed that the padding 
does offer much more performance gain with,  
for example,  CUDA 3.1 and the MAGMA DGEMM~\cite{NTD10b}:  the
speedup with respect to one 
CPU core was increased from 6.3 to 8.1 on C1060, and from 9.5 to 14.3
on C2050. 
Indeed since CUBLAS 3.2 performance has been improved for non block multiple
matrix sizes through MAGMA code~\cite{NTD10a}. 
Although for all versions the C2050 (with its improved
double precision performance) offers up to almost 
double speedup compared to 
the C1060, the performance obtained with both architectures justifies the use of 
the GPU for such an application.

\subsection{PROP execution profile}

\begin{figure}[t]
  \centering
 \includegraphics*[width=0.64\linewidth]{Chapters/chapter15/figures/CPU_1_CORE_TIMES.pdf}
\caption{CPU (1 core)  execution times for the off-diagonal sectors of the large case.}
\label{fig:CPU-timing}
\end{figure}

\begin{figure}[t]
  \centering
  \subfigure[GPU V5 on one C1060]{ 
\includegraphics*[width=0.64\linewidth]{Chapters/chapter15/figures/GPU_V5_C1060_TIMES.pdf}
\label{GPU-timing}} 

  \subfigure[GPU V5 on one C2050]{
\includegraphics*[width=0.64\linewidth]{Chapters/chapter15/figures/GPU_V5_C2050_TIMES.pdf}
\label{fermi-timing}} 
 
\caption{GPU execution times for the off-diagonal sectors of
  the large case.}
\label{fig:profileGPU}
\end{figure}

We detail here the execution profile on 
the CPU and the GPU for the evaluation of all off-diagonal sectors 
(the most representative ones) for a complete energy propagation. 
 Fig.~\ref{fig:CPU-timing} and \ref{fig:profileGPU} show CPU and GPU execution times for the 
171 off-diagonal sectors of  the large case (see Table \ref{data-sets}). 
``Copying, adding, scaling'' corresponds to the amplitude
  array construction (scaling) as well as to steps 1, 2 and 6 in
  Sect.~\ref{gpu-RO}, all implemented via CUDA kernels.
``CPU-GPU transfers'' 
aggregate transfer times for the $j$ amplitude
arrays and the scaling factors, as well as for the correction data.
 ``Amplitude matrix product'' corresponds to the DGEMM call to
 construct the local R-matrices from the $i$ and $j$ amplitude
 arrays. 
``Linear system solving'' and ``Output matrix product'' correspond
respectively to steps 3-4 and to step 5 in Sect.~\ref{gpu-RO}.

On one CPU core (see  Fig.~\ref{fig:CPU-timing}), 
 matrix products for the construction of the local
$R$-matrices require the same 
computation time during the whole propagation. Likewise the CPU time required by
 matrix products for the output $R$-matrix is constant within each
 strip. But as the global $R$-matrix is propagated from strip to
 strip, the sizes of
the matrices $U$ and $V$ increase, so does their multiplication time.
The time required to solve the linear system increases
slightly during the propagation.
These three operations (``Amplitude matrix product'', ``Output matrix
product'' and ``Linear system solving'') are clearly dominant in terms
of computation
time compared to the other remaining operations, which justify our
primary focus on such three linear algebra operations.


On the C1060 (see Fig.~\ref{GPU-timing}), we have
generally managed to obtain a similar speedup for all operations
(around 8, which corresponds to Fig.~\ref{fig:speedup_1core}). Only the linear system solving
presents a lower speedup (around~4). 
The CUDA kernels and the remaining CPU-GPU transfers make a minor contribution
to the overall computation time and do not require
additional improvements.

On the C2050 (see Fig.~\ref{fermi-timing}), additional speedup is
obtained for all operations, except for the 
CPU-GPU transfers and the linear system solving.
The CPU-GPU transfers are mainly due to the $j$ amplitude arrays, and
currently still correspond to minor times. When required, the
double-buffering technique may also be used to overlap such transfers
with computation on the GPU.



\section{Propagation of multiple concurrent energies on GPU}

Finally, we present here an improvement that can  
benefit from the Fermi architecture, as well as from the newest Kepler
architecture, 
both of which enable the concurrent execution of multiple 
CUDA kernels, thus offering additional speedup on  
GPUs for small or medium computation grain kernels.
In our case, the performance gain on the GPU is indeed limited
since most matrices have small or medium sizes. 
By using multiple streams within one CUDA context~\cite{CUDA_ProgGuide},
we can propagate multiple energies
concurrently on the Fermi GPU. 
It can be noticed that all GPU computations for all streams are
launched by the same host thread. We therefore rely here on the {\em legacy
API} of CUBLAS~\cite{CUBLAS} (like MAGMA)
without thread safety problems. 
A {\em breadth first} issue order is used for kernel
launchs \cite{CUDA_stream}: for a given GPU kernel, all kernel launchs
are indeed issued together in the host thread, using one stream for each
concurrent energy, in order to maximize concurrent kernel
execution.  
Of course, the memory available on the GPU must be large enough to
store all data structures required by each energy. 
Moreover, multiple streams are also used within the
propagation of each single energy 
in order to enable concurrent executions among the required kernels.


\begin{table}[t]
\begin{center}
\begin{tabular}{|c|c||c|c|c|c|c|}
  \hline
  \multirow{4}{0.18\linewidth}{Medium case} & Number of & 
  \multirow{2}{0.07\linewidth}{\centering 1} & 
  \multirow{2}{0.07\linewidth}{\centering 2} & 
  \multirow{2}{0.07\linewidth}{\centering 4} & 
  \multirow{2}{0.07\linewidth}{\centering 8} & 
  \multirow{2}{0.07\linewidth}{\centering 16} \\  
  & energies & & & & & \\  
  \cline{2-7}
  & Time (s) &  11.18 & 6.87 & 5.32 & 4.96 & 4.76 \\
  \cline{2-7}
  & Speedup & - & 1.63 & 2.10 & 2.26 & 2.35 \\
  \hline
  \hline
  \multirow{4}{0.18\linewidth}{Large case} & Number of & 
  \multirow{2}{0.07\linewidth}{\centering 1} & 
  \multicolumn{2}{c|}{\multirow{2}{0.07\linewidth}{\centering 2}} & 
  \multicolumn{2}{c|}{\multirow{2}{0.07\linewidth}{\centering 3}} \\
  & energies & & \multicolumn{2}{c|}{~} & \multicolumn{2}{c|}{~}  \\  
  \cline{2-7}
  & Time (s) & 509.51 & \multicolumn{2}{c|}{451.49} & \multicolumn{2}{c|}{436.72}  \\  
  \cline{2-7}
  & Speedup & - & \multicolumn{2}{c|}{1.13} & \multicolumn{2}{c|}{1.17}  \\  
  \hline
\end{tabular}
\caption{\label{t:perfs_V6} Performance results with multiple
  concurrent energies 
  on one C2070 GPU. GPU initialization times are not considered here. }
\end{center}
\end{table}

In order to have enough GPU memory to run two or three concurrent
energies for the large case, we use one C2070 GPU 
(featuring 6GB of memory) 
with one Intel X5650 hex-core CPU, CUDA 4.1 and CUBLAS 4.1, as
well as the latest MAGMA release (version 1.3.0). 
Substantial changes have been required 
in the MAGMA calls with respect to the previous versions of PROP that were using MAGMA 0.2. 
 Table~\ref{t:perfs_V6} presents the speedups 
obtained on the Fermi GPU for multiple concurrent
energies (up to sixteen since this is the maximum number of concurrent
kernel launches currently supported \cite{CUDA_ProgGuide}). 
With the medium case, speedups greater than 2 are obtained with four
concurrent energies or more. 
With the more realistic large case, the performance gain is lower mainly because of
the increase in matrix sizes, which implies a better GPU usage
with only one energy on the GPU. The concurrent execution of multiple
kernels is also limited by other operations on the
GPU \cite{CUDA_ProgGuide,CUDA_stream} and by the current MAGMA code which
prevents concurrent MAGMA calls in different streams. 
Better speedups can be here expected on the latest Kepler GPUs which
offer additional compute power, and whose {\em Hyper-Q} feature may help
improve further the GPU utilization with concurrent energies. 
On the contrary, the same code on the C1060 shows no speedup
 since the concurrent kernel launches are
serialized on this previous GPU architecture. 








\section{Conclusion and future work}
\label{conclusion}


In this chapter, we have presented our methodology and our results in
the deployment on 
a GPU of an application (the PROP program) in atomic physics. 

We have started by studying the numerical stability of PROP using
single precision arithmetic. This has shown that PROP
using single precision, while relatively stable for some small cases,
gives unsatisfactory results for realistic simulation cases above the
ionization threshold where there is a 
significant density of pseudo-states. It is
expected that this may not be the case below the ionization threshold
where the actual target states are less dense. This requires further
investigation. 

We have 
therefore deployed the PROP code in double precision on 
a GPU, with successive improvements. The different GPU versions 
each offer increasing speedups over the CPU version.
Compared to the single (respectively four) core(s) CPU version, the
optimal GPU implementation 
gives a speedup of 8.2 (resp. 4.6) on one C1060 GPU,
and a speedup of 15.9 (resp. 9.0) on one 
C2050 GPU with improved double precision performance.
An additional gain of around 15\% 
can also be obtained on one Fermi GPU
with large memory (C2070) thanks to concurrent kernel execution. 

Such speedups 
cannot be directly compared with the 1.15 speedup 
obtained with the HMPP version, since in our tests the CPUs are
different and the CUBLAS versions are more recent.  
However, the programming effort required 
progressively to deploy PROP on GPUs clearly offers improved and interesting speedups for this 
real-life   application in double precision with varying-size matrices. 


We are currently working on a hybrid CPU-GPU version that spreads the
computations of the independent energies on both the CPU
and the GPU. This will enable 
multiple energy execution on the CPU, with 
one or several core(s) dedicated to each energy (via multi-threaded
BLAS libraries). Multiple 
concurrent energies may also be propagated on each Fermi GPU. 
By merging this work with the current MPI PROP program, we will
obtain a scalable hybrid CPU-GPU version.
This final version will offer an additional level of parallelism
thanks to the MPI
standard in order to exploit multiple
nodes with multiple CPU cores and possibly multiple GPU cards. 

\putbib[Chapters/chapter15/biblio]