\subsection{Parallel implementation with CUDA }
On the CPU, both steps 3 and 4 contain the loop \verb=for= and a single thread executes all the instructions in the loop $n$ times. In this subsection, we explain how the GPU architecture can compute this loop and reduce the execution time.
-In the GPU, the schduler assigns the execution of this loop to a group of threads organised as a grid of blocks with block containing a number of threads. All threads within a block are executed concurrently in parallel. The instructions run on the GPU are grouped in special function called kernels. It's up to the programmer, to describe the execution context, that is the size of the Grid, the number of blocks and the number of threads per block upon the call of a given kernel, according to a special syntax defined by CUDA.
+In the GPU, the scheduler assigns the execution of this loop to a
+group of threads organised as a grid of blocks with block containing a
+number of threads. All threads within a block are executed
+concurrently in parallel. The instructions run on the GPU are grouped
+in special function called kernels. With CUDA, a programmer must
+describe the kernel execution context: the size of the Grid, the number of blocks and the number of threads per block.
-In CUDA programming, all the instructions of the \verb=for= loop are executed by the GPU as a kernel. A kernel is a function written in CUDA and defined by the \verb=__global__= qualifier added before a usual \verb=C= function, which instructs the compiler to generate appropriate code to pass it to the CUDA runtime in order to be executed on the GPU.
+%In CUDA programming, all the instructions of the \verb=for= loop are executed by the GPU as a kernel. A kernel is a function written in CUDA and defined by the \verb=__global__= qualifier added before a usual \verb=C= function, which instructs the compiler to generate appropriate code to pass it to the CUDA runtime in order to be executed on the GPU.
Algorithm~\ref{alg2-cuda} shows a sketch of the Ehrlich-Aberth algorithm using CUDA.
%\LinesNumbered
\caption{CUDA Algorithm to find roots with the Ehrlich-Aberth method}
-\KwIn{$Z^{0}$(Initial root's vector),$\varepsilon$ (error tolerance threshold), P(Polynomial to solve), $\Delta z_{max}$ (maximum value of stop condition)}
+\KwIn{$Z^{0}$ (Initial root's vector), $\varepsilon$ (error tolerance threshold), P(Polynomial to solve), $\Delta z_{max}$ (maximum value of stop condition)}
-\KwOut {Z(The solution root's vector)}
+\KwOut {Z (The solution root's vector)}
\BlankLine
-Initialization of the coeffcients of the polynomial to solve\;
+Initialization of the coefficients of the polynomial to solve\;
Initialization of the solution vector $Z^{0}$\;
Allocate and copy initial data to the GPU global memory\;
k=0\;
After the initialisation step, all data of the root finding problem to be solved must be copied from the CPU memory to the GPU global memory, because the GPUs only access data already present in their memories. Next, all the data-parallel arithmetic operations inside the main loop \verb=(do ... while(...))= are executed as kernels by the GPU. The first kernel named \textit{save} in line 6 of Algorithm~\ref{alg2-cuda} consists in saving the vector of polynomial's root found at the previous time-step in GPU memory, in order to check the convergence of the roots after each iteration (line 8, Algorithm~\ref{alg2-cuda}).
-The second kernel executes the iterative function $H$ and updates $z^{k}$, according to Algorithm~\ref{alg3-update}. We notice that the update kernel is called in two forms, separated with the value of \emph{R} which determines the radius beyond which we apply the logarithm computation of the power of a complex.
+The second kernel executes the iterative function $H$ and updates
+$z^{k}$, according to Algorithm~\ref{alg3-update}. We notice that the
+update kernel is called in two forms, separated with the value of
+\emph{R} which determines the radius beyond which we apply the
+exponential logarithm algorithm.
\begin{algorithm}[H]
\label{alg3-update}
\eIf{$(\left|Z^{(k)}\right|<= R)$}{
$kernel\_update(d\_z^{k})$\;}
{
-$kernel\_update\_Log(d\_z^{k})$\;
+$kernel\_update\_ExpoLog(d\_z^{k})$\;
}
\end{algorithm}
-The first form executes formula \ref{eq:SimplePolynome} if the modulus of the current complex is less than the a certain value called the radius i.e. ($ |z^{k}_{i}|<= R$), else the kernel executes formulas (Eq.~\ref{deflncomplex},Eq.~\ref{defexpcomplex}). The radius $R$ is evaluated as :
+The first form executes formula \ref{eq:SimplePolynome} if the modulus
+of the current complex is less than the a certain value called the
+radius i.e. ($ |z^{k}_{i}|<= R$), else the kernel executes the EA.EL
+function Eq.~\ref{Log_H2}
+( with Eq.~\ref{deflncomplex}, Eq.~\ref{defexpcomplex}). The radius $R$ is evaluated as :
$$R = \exp( \log(DBL\_MAX) / (2*n) )$$ where $DBL\_MAX$ stands for the maximum representable double value.
