+
+\subsubsection{Polynomials Initialization}
+The initialization of a polynomial $p(z)$ is done by setting each of the $n$ complex coefficients $a_{i}$:
+
+\begin{equation}
+\label{eq:SimplePolynome}
+ p(z)=\sum{a_{i}z^{n-i}} , a_{n} \neq 0,a_{0}=1, a_{i}\subset C
+\end{equation}
+
+
+\subsubsection{Vector $Z^{(0)}$ Initialization}
+\label{sec:vec_initialization}
+As for any iterative method, we need to choose $n$ initial guess points $z^{0}_{i}, i = 1, . . . , n.$
+The initial guess is very important since the number of steps needed by the iterative method to reach
+a given approximation strongly depends on it.
+In~\cite{Aberth73} the Ehrlich-Aberth iteration is started by selecting $n$
+equi-spaced points on a circle of center 0 and radius r, where r is
+an upper bound to the moduli of the zeros. Later, Bini and al.~\cite{Bini96}
+performed this choice by selecting complex numbers along different
+circles which relies on the result of~\cite{Ostrowski41}.
+
+\begin{equation}
+\label{eq:radiusR}
+%%\begin{align}
+\sigma_{0}=\frac{u+v}{2};u=\frac{\sum_{i=1}^{n}u_{i}}{n.max_{i=1}^{n}u_{i}};
+v=\frac{\sum_{i=0}^{n-1}v_{i}}{n.min_{i=0}^{n-1}v_{i}};\\
+%%\end{align}
+\end{equation}
+Where:
+\begin{equation}
+u_{i}=2.|a_{i}|^{\frac{1}{i}};
+v_{i}=\frac{|\frac{a_{n}}{a_{i}}|^{\frac{1}{n-i}}}{2}.
+\end{equation}
+
+\subsubsection{Iterative Function}
+The operator used by the Aberth method is corresponding to the
+following equation~\ref{Eq:EA1} which will enable the convergence towards
+polynomial solutions, provided all the roots are distinct.
+
+%Here we give a second form of the iterative function used by the Ehrlich-Aberth method:
+
+\begin{equation}
+\label{Eq:EA1}
+EA: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
+{1-\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}\sum_{j=1,j\neq i}^{j=n}{\frac{1}{(z_{i}^{k}-z_{j}^{k})}}}, i=1,. . . .,n
+\end{equation}
+
+\subsubsection{Convergence Condition}
+The convergence condition determines the termination of the algorithm. It consists in stopping the iterative function when the roots are sufficiently stable. We consider that the method converges sufficiently when:
+
+\begin{equation}
+\label{eq:Aberth-Conv-Cond}
+\forall i \in [1,n];\vert\frac{z_{i}^{k}-z_{i}^{k-1}}{z_{i}^{k}}\vert<\xi
+\end{equation}
+
+
+%\begin{figure}[htbp]
+%\centering
+ % \includegraphics[angle=-90,width=0.5\textwidth]{EA-Algorithm}
+%\caption{The Ehrlich-Aberth algorithm on single GPU}
+%\label{fig:03}
+%\end{figure}
+
+%the Ehrlich-Aberth method is an iterative method, contain 4 steps, start from the initial approximations of all the
+%roots of the polynomial,the second step initialize the solution vector $Z$ using the Guggenheimer method to assure the distinction of the initial vector roots, than in step 3 we apply the the iterative function based on the Newton's method and Weiestrass operator[...,...], wich will make it possible to converge to the roots solution, provided that all the root are different. At the end of each application of the iterative function, a stop condition is verified consists in stopping the iterative process when the whole of the modules of the roots
+%are lower than a fixed value $ε$
+
+
+\subsection{EA parallel implementation on CUDA}
+Like any parallel code, a GPU parallel implementation first
+requires to determine the sequential tasks and the
+parallelizable parts of the sequential version of the
+program/algorithm. In our case, all the operations that are easy
+to execute in parallel must be made by the GPU to accelerate
+the execution of the application, like the step 3 and step 4. On the other hand, all the
+sequential operations and the operations that have data
+dependencies between threads or recursive computations must
+be executed by only one CUDA or CPU thread (step 1 and step 2). Initially we specifies the organization of threads in parallel, need to specify the dimension of the grid Dimgrid: the number of block per grid and block by DimBlock: the number of threads per block required to process a certain task.
+
+we create the kernel, for step 3 we have two kernels, the
+first named \textit{save} is used to save vector $Z^{K-1}$ and the kernel
+\textit{update} is used to update the $Z^{K}$ vector. In step 4 a kernel is
+created to test the convergence of the method. In order to
+compute function H, we have two possibilities: either to use
+the Jacobi method, or the Gauss-Seidel method which uses the
+most recent computed roots. It is well known that the Gauss-
+Seidel mode converges more quickly. So, we used the Gauss-Seidel mode of iteration. To
+parallelize the code, we created kernels and many functions to
+be executed on the GPU for all the operations dealing with the
+computation on complex numbers and the evaluation of the
+polynomials. As said previously, we managed both functions
+of evaluation of a polynomial: the normal method, based on
+the method of Horner and the method based on the logarithm
+of the polynomial. All these methods were rather long to
+implement, as the development of corresponding kernels with
+CUDA is longer than on a CPU host. This comes in particular
+from the fact that it is very difficult to debug CUDA running
+threads like threads on a CPU host. In the following paragraph
+Algorithm 1 shows the GPU parallel implementation of Ehrlich-Aberth method.
+
+Algorithm~\ref{alg2-cuda} shows a sketch of the Ehrlich-Aberth method using CUDA.
+
+\begin{enumerate}
+\begin{algorithm}[htpb]
+\label{alg1-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), Pu (Derivative of P), $n$ (Polynomial degrees), $\Delta z_{max}$ (Maximum value of stop condition)}
+
+\KwOut {$Z$ (Solution root's vector), $ZPrec$ (Previous solution root's vector)}
+
+%\BlankLine
+
+\item Initialization of the of P\;
+\item Initialization of the of Pu\;
+\item Initialization of the solution vector $Z^{0}$\;
+\item Allocate and copy initial data to the GPU global memory\;
+\item k=0\;
+\While {$\Delta z_{max} > \epsilon$}{
+\item Let $\Delta z_{max}=0$\;
+\item $ kernel\_save(ZPrec,Z)$\;
+\item k=k+1\;
+\item $ kernel\_update(Z,P,Pu)$\;
+\item $kernel\_testConverge(\Delta z_{max},Z,ZPrec)$\;
+
+}
+\item Copy results from GPU memory to CPU memory\;
+\end{algorithm}
+\end{enumerate}
+~\\
+
+
+