1 \documentclass[review]{elsarticle}
3 \usepackage{lineno,hyperref}
4 \usepackage[utf8]{inputenc}
5 %%\usepackage[T1]{fontenc}
6 %%\usepackage[french]{babel}
8 \usepackage{amsmath,amsfonts,amssymb}
9 \usepackage[ruled,vlined]{algorithm2e}
10 %\usepackage[french,boxed,linesnumbered]{algorithm2e}
11 \usepackage{array,multirow,makecell}
12 \usepackage[textsize=footnotesize]{todonotes}
14 \newcommand{\RC}[2][inline]{%
15 \todo[color=blue!10,#1]{\sffamily\textbf{RC:} #2}\xspace}
16 \newcommand{\KG}[2][inline]{%
17 \todo[color=green!10,#1]{\sffamily\textbf{KG:} #2}\xspace}
19 \newcommand{\AS}[2][inline]{%
20 \todo[color=orange!10,#1]{\sffamily\textbf{AS:} #2}\xspace}
25 \newcolumntype{R}[1]{>{\raggedleft\arraybackslash }b{#1}}
26 \newcolumntype{L}[1]{>{\raggedright\arraybackslash }b{#1}}
27 \newcolumntype{C}[1]{>{\centering\arraybackslash }b{#1}}
32 \journal{Journal of \LaTeX\ Templates}
34 %%%%%%%%%%%%%%%%%%%%%%%
35 %% Elsevier bibliography styles
36 %%%%%%%%%%%%%%%%%%%%%%%
37 %% To change the style, put a % in front of the second line of the current style and
38 %% remove the % from the second line of the style you would like to use.
39 %%%%%%%%%%%%%%%%%%%%%%%
42 %\bibliographystyle{model1-num-names}
44 %% Numbered without titles
45 %\bibliographystyle{model1a-num-names}
48 %\bibliographystyle{model2-names.bst}\biboptions{authoryear}
51 %\usepackage{numcompress}\bibliographystyle{model3-num-names}
53 %% Vancouver name/year
54 %\usepackage{numcompress}\bibliographystyle{model4-names}\biboptions{authoryear}
57 %\bibliographystyle{model5-names}\biboptions{authoryear}
60 %\usepackage{numcompress}\bibliographystyle{model6-num-names}
62 %% `Elsevier LaTeX' style
63 \bibliographystyle{elsarticle-num}
64 %%%%%%%%%%%%%%%%%%%%%%%
70 \title{Efficient high degree polynomial root finding using GPU}
72 %% Group authors per affiliation:
73 %\author{Elsevier\fnref{myfootnote}}
74 %\address{Radarweg 29, Amsterdam}
75 %\fntext[myfootnote]{Since 1880.}
77 %% or include affiliations in footnotes:
78 \author[mymainaddress]{Kahina Ghidouche}
79 %%\ead[url]{kahina.ghidouche@univ-bejaia.dz}
80 \cortext[mycorrespondingauthor]{Corresponding author}
81 \ead{kahina.ghidouche@univ-bejaia.dz}
83 \author[mysecondaryaddress]{Raphaël Couturier\corref{mycorrespondingauthor}}
84 %%\cortext[mycorrespondingauthor]{Corresponding author}
85 \ead{raphael.couturier@univ-fcomte.fr}
87 \author[mymainaddress]{Abderrahmane Sider}
88 %%\cortext[mycorrespondingauthor]{Corresponding author}
89 \ead{ar.sider@univ-bejaia.dz}
91 \address[mymainaddress]{Laboratoire LIMED, Faculté des sciences
92 exactes, Université de Bejaia, 06000, Algeria}
93 \address[mysecondaryaddress]{FEMTO-ST Institute, University of
94 Bourgogne Franche-Comte, France }
97 Polynomials are mathematical algebraic structures that play a great
98 role in science and engineering. Finding the roots of high degree
99 polynomials is computationally demanding. In this paper, we present
100 the results of a parallel implementation of the Ehrlich-Aberth
101 algorithm for the root finding problem for high degree polynomials on
102 GPU architectures. The main result of this
103 work is to be able to solve high degree polynomials (up
104 to 1,000,000) efficiently. We also compare the results with a
105 sequential implementation and the Durand-Kerner method on full and
110 Polynomial root finding, Iterative methods, Ehrlich-Aberth, Durand-Kerner, GPU
117 \section{The problem of finding the roots of a polynomial}
118 Polynomials are mathematical algebraic structures used in science and engineering to capture physical phenomena and to express any outcome in the form of a function of some unknown variables. Formally speaking, a polynomial $p(x)$ of degree \textit{n} having $n$ coefficients in the complex plane \textit{C} is :
121 {\Large p(x)=\sum_{i=0}^{n}{a_{i}x^{i}}}.
125 The root finding problem consists in finding the values of all the $n$ values of the variable $x$ for which \textit{p(x)} is nullified. Such values are called zeros of $p$. If zeros are $\alpha_{i},\textit{i=1,...,n}$ the $p(x)$ can be written as :
127 {\Large p(x)=a_{n}\prod_{i=1}^{n}(x-\alpha_{i}), a_{0} a_{n}\neq 0}.
130 The problem of finding a root is equivalent to that of solving a fixed-point problem. To see this, consider the fixed-point problem of finding the $n$-dimensional
131 vector $x$ such that :
135 where $g : C^{n}\longrightarrow C^{n}$. Usually, we can easily
136 rewrite this fixed-point problem as a root-finding problem by
137 setting $f(x) = x-g(x)$ and likewise we can recast the
138 root-finding problem into a fixed-point problem by setting :
143 It is often impossible to solve such nonlinear equation
144 root-finding problems analytically. When this occurs, we turn to
145 numerical methods to approximate the solution.
146 Generally speaking, algorithms for solving problems can be divided into
147 two main groups: direct methods and iterative methods.
149 Direct methods only exist for $n \leq 4$, solved in closed form
150 by G. Cardano in the mid-16th century. However, N. H. Abel in the early 19th
151 century proved that polynomials of degree five or more could not
152 be solved by direct methods. Since then, mathematicians have
153 focussed on numerical (iterative) methods such as the famous
154 Newton method, the Bernoulli method of the 18th century, and the Graeffe method.
156 Later on, with the advent of electronic computers, other methods have
157 been developed such as the Jenkins-Traub method, the Larkin method,
158 the Muller method, and several other methods for the simultaneous
159 approximation of all the roots, starting with the Durand-Kerner (DK)
164 DK: z_i^{k+1}=z_{i}^{k}-\frac{P(z_i^{k})}{\prod_{i\neq j}(z_i^{k}-z_j^{k})}, i = 1, . . . , n,
167 where $z_i^k$ is the $i^{th}$ root of the polynomial $p$ at the
171 This formula is mentioned for the first time by
172 Weiestrass~\cite{Weierstrass03} as part of the fundamental theorem
173 of Algebra and is rediscovered by Ilieff~\cite{Ilie50},
174 Docev~\cite{Docev62}, Durand~\cite{Durand60},
175 Kerner~\cite{Kerner66}. Another method discovered by
176 Borsch-Supan~\cite{ Borch-Supan63} and also described and brought
177 in the following form by Ehrlich~\cite{Ehrlich67} and
178 Aberth~\cite{Aberth73} uses a different iteration formula given as:
182 EA: z_i^{k+1}=z_i^{k}-\frac{1}{{\frac {P'(z_i^{k})} {P(z_i^{k})}}-{\sum_{i\neq j}\frac{1}{(z_i^{k}-z_j^{k})}}}, i = 1, . . . , n,
185 where $p'(z)$ is the polynomial derivative of $p$ evaluated in the
188 Aberth, Ehrlich and Farmer-Loizou~\cite{Loizou83} have proved that
189 the Ehrlich-Aberth method (EA) has a cubic order of convergence for simple roots whereas the Durand-Kerner has a quadratic order of convergence.
192 Moreover, the convergence times of iterative methods
193 drastically increases like the degrees of high polynomials. It is expected that the
194 parallelization of these algorithms will reduce the execution times.
196 Many authors have dealt with the parallelization of
197 simultaneous methods, i.e. that find all the zeros simultaneously.
198 Freeman~\cite{Freeman89} implemented and compared DK, EA and another method of the fourth order proposed
199 by Farmer and Loizou~\cite{Loizou83}, on an 8-processor linear
200 chain, for polynomials of degree 8. The third method often
201 diverges, but the first two methods have speed-ups equal to 5.5. Later,
202 Freeman and Bane~\cite{Freemanall90} considered asynchronous
203 algorithms, in which each processor continues to update its
204 approximations even though the latest values of other roots
205 have not yet been received from the other processors. In contrast,
206 synchronous algorithms wait the computation of all roots at a given
207 iterations before making a new one.
208 Couturier and al.~\cite{Raphaelall01} proposed two methods of parallelization for
209 a shared memory architecture and for distributed memory one. They were able to
210 compute the roots of sparse polynomials of degree 10,000 in 430 seconds with only 8
211 personal computers and 2 communications per iteration. Compared to sequential implementations
212 where it takes up to 3,300 seconds to obtain the same results, the
213 authors' work experiment show an interesting speedup.
215 Few works have been conducted after those works until the appearance of
216 the Compute Unified Device Architecture (CUDA)~\cite{CUDA10}, a
217 parallel computing platform and a programming model invented by
218 NVIDIA. The computing power of GPUs (Graphics Processing Unit) has exceeded that of CPUs. However, CUDA adopts a totally new computing architecture to use the
219 hardware resources provided by GPU in order to offer a stronger
220 computing ability to the massive data computing.
223 Ghidouche and al~\cite{Kahinall14} proposed an implementation of the
224 Durand-Kerner method on GPU. Their main
225 result showed that a parallel CUDA implementation is about 10 times faster than
226 the sequential implementation on a single CPU for sparse
227 polynomials of degree 48,000.
230 In this paper, we focus on the implementation of the Ehrlich-Aberth
231 method for high degree polynomials on GPU. We propose an adaptation of
232 the exponential logarithm in order to be able to solve sparse and full
233 polynomial of degree up to $1,000,000$. The paper is organized as
234 follows. Initially, we recall the Ehrlich-Aberth method in
235 Section~\ref{sec1}. Improvements for the Ehrlich-Aberth method are
236 proposed in Section \ref{sec2}. Related work to the implementation of
237 simultaneous methods using a parallel approach is presented in Section
238 \ref{secStateofArt}. In Section~\ref{sec5} we propose a parallel
239 implementation of the Ehrlich-Aberth method on GPU and discuss
240 it. Section~\ref{sec6} presents and investigates our implementation
241 and experimental study results. Finally, Section~\ref{sec7} concludes
242 this paper and gives some hints for future research directions in this
245 \section{Ehrlich-Aberth method}
247 A cubically convergent iteration method to find zeros of
248 polynomials was proposed by O. Aberth~\cite{Aberth73}. The
249 Ehrlich-Aberth method contains 4 main steps, presented in what
252 %The Aberth method is a purely algebraic derivation.
253 %To illustrate the derivation, we let $w_{i}(z)$ be the product of linear factors
256 %w_{i}(z)=\prod_{j=1,j \neq i}^{n} (z-x_{j})
259 %And let a rational function $R_{i}(z)$ be the correction term of the
260 %Weistrass method~\cite{Weierstrass03}
263 %R_{i}(z)=\frac{p(z)}{w_{i}(z)} , i=1,2,...,n.
266 %Differentiating the rational function $R_{i}(z)$ and applying the
267 %Newton method, we have:
270 %\frac{R_{i}(z)}{R_{i}^{'}(z)}= \frac{p(z)}{p^{'}(z)-p(z)\frac{w_{i}(z)}{w_{i}^{'}(z)}}= \frac{p(z)}{p^{'}(z)-p(z) \sum _{j=1,j \neq i}^{n}\frac{1}{z-x_{j}}}, i=1,2,...,n
272 %where R_{i}^{'}(z)is the rational function derivative of F evaluated in the point z
273 %Substituting $x_{j}$ for $z_{j}$ we obtain the Aberth iteration method.%
276 \subsection{Polynomials Initialization}
277 The initialization of a polynomial $p(z)$ is done by setting each of the $n$ complex coefficients $a_{i}$:
280 \label{eq:SimplePolynome}
281 p(z)=\sum{a_{i}z^{n-i}} , a_{n} \neq 0,a_{0}=1, a_{i}\subset C
285 \subsection{Vector $Z^{(0)}$ Initialization}
286 \label{sec:vec_initialization}
287 As for any iterative method, we need to choose $n$ initial guess points $z^{0}_{i}, i = 1, . . . , n.$
288 The initial guess is very important since the number of steps needed by the iterative method to reach
289 a given approximation strongly depends on it.
290 In~\cite{Aberth73} the Ehrlich-Aberth iteration is started by selecting $n$
291 equi-spaced points on a circle of center 0 and radius r, where r is
292 an upper bound to the moduli of the zeros. Later, Bini and al.~\cite{Bini96}
293 performed this choice by selecting complex numbers along different
294 circles which relies on the result of~\cite{Ostrowski41}.
299 \sigma_{0}=\frac{u+v}{2};u=\frac{\sum_{i=1}^{n}u_{i}}{n.max_{i=1}^{n}u_{i}};
300 v=\frac{\sum_{i=0}^{n-1}v_{i}}{n.min_{i=0}^{n-1}v_{i}};\\
305 u_{i}=2.|a_{i}|^{\frac{1}{i}};
306 v_{i}=\frac{|\frac{a_{n}}{a_{i}}|^{\frac{1}{n-i}}}{2}.
309 \subsection{Iterative Function}
310 %The operator used by the Aberth method is corresponding to the
311 %following equation~\ref{Eq:EA} which will enable the convergence towards
312 %polynomial solutions, provided all the roots are distinct.
314 Here we give a second form of the iterative function used by the Ehrlich-Aberth method:
318 EA2: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
319 {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
321 It can be noticed that this equation is equivalent to Eq.~\ref{Eq:EA},
322 but we prefer the latter one because we can use it to improve the
323 Ehrlich-Aberth method and find the roots of high degree polynomials. More
324 details are given in Section~\ref{sec2}.
325 \subsection{Convergence Condition}
326 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:
329 \label{eq:Aberth-Conv-Cond}
330 \forall i \in [1,n];\vert\frac{z_{i}^{k}-z_{i}^{k-1}}{z_{i}^{k}}\vert<\xi
334 \section{Improving the Ehrlich-Aberth Method for high degree polynomials with exp-log formulation}
336 With high degree polynomial, the Ehrlich-Aberth method implementation,
337 as well as the Durand-Kerner implementation, suffers from overflow problems. This
338 situation occurs, for instance, in the case where a polynomial,
339 having positive coefficients and a large degree, is computed at a
340 point $\xi$ where $|\xi| > 1$, where $|z|$ stands for the modulus of a complex $z$. Indeed, the limited number in the
341 mantissa of floating points representations makes the computation of $p(z)$ wrong when z
342 is large. For example $(10^{50}) +1+ (- 10^{50})$ will give the wrong result
343 of $0$ instead of $1$. Consequently, we can not compute the roots
344 for large degrees. This problem was discussed earlier in
345 ~\cite{Karimall98} for the Durand-Kerner method. The authors
346 propose to use the logarithm and the exponential of a complex in order to compute the power at a high exponent.
350 \forall(x,y)\in R^{*2}; \ln (x+i.y)=\ln(x^{2}+y^{2})
351 2+i.\arcsin(y\sqrt{x^{2}+y^{2}})_{\left] -\pi, \pi\right[ }
355 \label{defexpcomplex}
356 \forall(x,y)\in R^{*2}; \exp(x+i.y) & = \exp(x).\exp(i.y)\\
357 & =\exp(x).\cos(y)+i.\exp(x).\sin(y)\label{defexpcomplex1}
361 Using the logarithm (eq.~\ref{deflncomplex}) and the exponential (eq.~\ref{defexpcomplex}) operators, we can replace any multiplications and divisions with additions and subtractions. Consequently, computations
362 manipulate lower absolute values and the roots for large polynomial degrees can be looked for successfully~\cite{Karimall98}.
364 Applying this solution for the iteration function Eq.~\ref{Eq:Hi} of
365 Ehrlich-Aberth method, we obtain the following iteration function with exponential and logarithm:
366 %%$$ \exp \bigl( \ln(p(z)_{k})-ln(\ln(p(z)_{k}^{'}))- \ln(1- \exp(\ln(p(z)_{k})-ln(\ln(p(z)_{k}^{'})+\ln\sum_{i\neq j}^{n}\frac{1}{z_{k}-z_{j}})$$
369 EA.EL: z^{k+1}_{i}=z_{i}^{k}-\exp \left(\ln \left(
370 p(z_{i}^{k})\right)-\ln\left(p'(z^{k}_{i})\right)- \ln\left(1-Q(z^{k}_{i})\right)\right),
377 Q(z^{k}_{i})=\exp\left( \ln (p(z^{k}_{i}))-\ln(p'(z^{k}_{i}))+\ln \left(
378 \sum_{i\neq j}^{n}\frac{1}{z^{k}_{i}-z^{k}_{j}}\right)\right)i=1,...,n,
381 This solution is applied when the root except the circle unit, represented by the radius $R$ evaluated in C language as :
384 R = exp(log(DBL\_MAX)/(2*n) );
390 %R = \exp( \log(DBL\_MAX) / (2*n) )
392 where \verb=DBL_MAX= stands for the maximum representable \verb=double= value.
394 \section{Implementation of simultaneous methods in a parallel computer}
395 \label{secStateofArt}
396 The main problem of simultaneous methods is that the
397 time needed for convergence is increased when we increase
398 the degree of the polynomial. The parallelization of these
399 algorithms is expected to improve the convergence time.
400 Authors usually adopt one of the two following approaches to parallelize root
401 finding algorithms. The first approach aims at reducing the total number of
402 iterations as in Miranker
403 ~\cite{Mirankar68,Mirankar71}, Schedler~\cite{Schedler72} and
404 Winograd~\cite{Winogard72}. The second approach aims at reducing the
405 computation time per iteration, as reported
406 in~\cite{Benall68,Jana06,Janall99,Riceall06}.
408 There are many schemes for the simultaneous approximation of all roots of a given
409 polynomial. Several works on different methods and issues of root
410 finding have been reported in~\cite{Azad07, Gemignani07, Kalantari08,
411 Zhancall08, Zhuall08}. However, the Durand-Kerner and the Ehrlich-Aberth methods are the most practical choices among
412 them~\cite{Bini04}. These two methods have been extensively
413 studied for parallelization due to their intrinsic parallelism, i.e. the
414 computations involved in both methods have some inherent
415 parallelism that can be suitably exploited by SIMD machines.
416 Moreover, they have fast a rate of convergence (quadratic for the
417 Durand-Kerner and cubic for the Ehrlich-Aberth). Various parallel
418 algorithms reported for these methods can be found
419 in~\cite{Cosnard90, Freeman89,Freemanall90,Jana99,Janall99}.
420 Freeman and Bane~\cite{Freemanall90} presented two parallel
421 algorithms on a local memory MIMD computer with the compute-to
422 communication time ratio O(n). However, their algorithms require
423 each processor to communicate its current approximation to all
424 other processors at the end of each iteration (synchronous). Therefore they
425 cause a high degree of memory conflict. Recently the author
426 in~\cite{Mirankar71} proposed two versions of parallel algorithm
427 for the Durand-Kerner method, and the Ehrlich-Aberth method on a model of
428 Optoelectronic Transpose Interconnection System (OTIS). The
429 algorithms are mapped on an OTIS-2D torus using $N$ processors. This
430 solution needs $N$ processors to compute $N$ roots, which is not
431 practical for solving large degree polynomials.
433 %Until very recently, the literature did not mention implementations
434 %able to compute the roots of large degree polynomials (higher then
435 %1000) and within small or at least tractable times.
437 Finding polynomial roots rapidly and accurately is the main objective of our work.
438 With the advent of CUDA (Compute Unified Device
439 Architecture), finding the roots of polynomials receives a new attention because of the new possibilities to solve higher degree polynomials in less time.
440 In~\cite{Kahinall14} we already proposed the first implementation
441 of a root finding method on GPUs, that of the Durand-Kerner method. The main result showed
442 that a parallel CUDA implementation is 10 times as fast as the
443 sequential implementation on a single CPU for high degree
444 polynomials of 48,000.
445 %In this paper we present a parallel implementation of Ehrlich-Aberth
446 %method on GPUs for sparse and full polynomials with high degree (up
450 %% \section {A CUDA parallel Ehrlich-Aberth method}
451 %% In the following, we describe the parallel implementation of Ehrlich-Aberth method on GPU
452 %% for solving high degree polynomials (up to $1,000,000$). First, the hardware and software of the GPUs are presented. Then, the CUDA parallel Ehrlich-Aberth method is presented.
455 %% \subsection{Background on the CUDA Programming Model}
457 %% The CUDA programming model is similar in style to a single program
458 %% multiple-data (SPMD) software model. The GPU is viewed as a
459 %% coprocessor that executes data-parallel kernel functions. CUDA
460 %% provides three key abstractions, a hierarchy of thread groups,
461 %% shared memories, and barrier synchronization. Threads have a three
462 %% level hierarchy. A grid is a set of thread blocks that execute a
463 %% kernel function. Each grid consists of blocks of threads. Each
464 %% block is composed of hundreds of threads. Threads within one block
465 %% can share data using shared memory and can be synchronized at a
466 %% barrier. All threads within a block are executed concurrently on a
467 %% multithreaded architecture.The programmer specifies the number of
468 %% threads per block, and the number of blocks per grid. A thread in
469 %% the CUDA programming language is much lighter weight than a thread
470 %% in traditional operating systems. A thread in CUDA typically
471 %% processes one data element at a time. The CUDA programming model
472 %% has two shared read-write memory spaces, the shared memory space
473 %% and the global memory space. The shared memory is local to a block
474 %% and the global memory space is accessible by all blocks. CUDA also
475 %% provides two read-only memory spaces, the constant space and the
476 %% texture space, which reside in external DRAM, and are accessed via
479 \section{ GPU Implementation of the Ehrlich-Aberth method}
481 \KG{In the following, we describe the parallel implementation on GPU of the Ehrlich-Aberth method, for solving high degree polynomials. First, the hardware and software of the GPUs are presented. Then, the principal steps of the implementation are described.}
483 \subsection{Background on the GPU architecture}
484 \KG{A GPU is viewed as an accelerator for the data-parallel and
485 intensive arithmetic computations. It draws its computing power
486 from the parallel nature of its hardware and software
487 architectures. A GPU is composed of hundreds of Streaming
488 Processors (SPs) organized in several blocks called Streaming
489 Multiprocessors (SMs). It also has a memory hierarchy. It has a
490 private read-write local memory per SP, fast shared memory and
491 read-only constant and texture caches per SM and a read-write
492 global memory shared by all its SPs~\cite{NVIDIA10}.
494 On a CPU equipped with a GPU, all the data-parallel and intensive
495 functions of an application running on the CPU are off-loaded onto
496 the GPU in order to accelerate their computations. A similar
497 data-parallel function is executed on a GPU as a kernel by
498 thousands or even millions of parallel threads, grouped together
499 as a grid of thread blocks. Therefore, each SM of the GPU executes
500 one or more thread blocks in SIMD fashion (Single Instruction,
501 Multiple Data) and in turn each SP of a GPU SM runs one or more
502 threads within a block in SIMT fashion (Single Instruction,
503 Multiple threads). Indeed at any given clock cycle, the threads
504 execute the same instruction of a kernel, but each of them
505 operates on different data.
506 GPUs only work on data filled in their
507 global memories and the final results of their kernel executions
508 must be communicated to their CPUs. Hence, the data must be
509 transferred in and out of the GPU. However, the speed of memory
510 copy between the GPU and the CPU is slower than the memory
511 bandwidths of the GPU memories and, thus, it dramatically affects
512 the performances of GPU computations. Accordingly, it is necessary
513 to limit as much as possible, data transfers between the GPU and its CPU during the
515 \subsection{Background on the CUDA Programming Model}
516 \KG{The CUDA programming model is similar in style to a single program
517 multiple-data (SPMD) software model. The GPU is viewed as a
518 coprocessor that executes data-parallel kernel functions. CUDA
519 provides three key abstractions, a hierarchy of thread groups,
520 shared memories, and barrier synchronization. Threads have a three
521 level hierarchy. A grid is a set of thread blocks that execute a
522 kernel function. Each grid consists of blocks of threads. Each
523 block is composed of hundreds of threads. Threads within one block
524 can share data using shared memory and can be synchronized at a
525 barrier. All threads within a block are executed concurrently on a
526 multi-threaded architecture.The programmer specifies the number of
527 threads per block, and the number of blocks per grid. A thread in
528 the CUDA programming language is much lighter weight than a thread
529 in traditional operating systems. A thread in CUDA typically
530 processes one data element at a time. The CUDA programming model
531 has two shared read-write memory spaces, the shared memory space
532 and the global memory space. The shared memory is local to a block
533 and the global memory space is accessible by all blocks. CUDA also
534 provides two read-only memory spaces, the constant space and the
535 texture space, which reside in external DRAM, and are accessed via
540 %%\subsection{A CUDA implementation of the Aberth's method }
541 %%\subsection{A GPU implementation of the Aberth's method }
545 %% \subsection{Sequential Ehrlich-Aberth algorithm}
546 %% The main steps of Ehrlich-Aberth method are shown in Algorithm.~\ref{alg1-seq} :
548 %% \begin{algorithm}[H]
551 %% \caption{A sequential algorithm to find roots with the Ehrlich-Aberth method}
553 %% \KwIn{$Z^{0}$ (Initial root's vector), $\varepsilon$ (error tolerance
554 %% threshold), $P$ (Polynomial to solve),$Pu$ (the derivative of P) $\Delta z_{max}$ (maximum value
555 %% of stop condition), k (number of iteration), n (Polynomial's degrees)}
556 %% \KwOut {$Z$ (The solution root's vector), $ZPrec$ (the previous solution root's vector)}
560 %% Initialization of $P$\;
561 %% Initialization of $Pu$\;
562 %% Initialization of the solution vector $Z^{0}$\;
563 %% $\Delta z_{max}=0$\;
566 %% \While {$\Delta z_{max} > \varepsilon$}{
567 %% Let $\Delta z_{max}=0$\;
568 %% \For{$j \gets 0 $ \KwTo $n$}{
569 %% $ZPrec\left[j\right]=Z\left[j\right]$;// save Z at the iteration k.\
571 %% $Z\left[j\right]=H\left(j, Z, P, Pu\right)$;//update Z with the iterative function.\
575 %% \For{$i \gets 0 $ \KwTo $n-1$}{
576 %% $c= testConverge(\Delta z_{max},ZPrec\left[j\right],Z\left[j\right])$\;
577 %% \If{$c > \Delta z_{max}$ }{
578 %% $\Delta z_{max}$=c\;}
585 %% In this sequential algorithm, one CPU thread executes all the steps. Let us look to the $3^{rd}$ step i.e. the execution of the iterative function, 2 sub-steps are needed. The first sub-step \textit{save}s the solution vector of the previous iteration, the second sub-step \textit{update}s or computes the new values of the roots vector.
587 \subsection{Parallel implementation with CUDA }
589 %In order to implement the Ehrlich-Aberth method in CUDA, it is
590 %possible to use the Jacobi scheme or the Gauss-Seidel one. With the
591 %Jacobi iteration, at iteration $k+1$ we need all the previous values
592 %$z^{k}_{i}$ to compute the new values $z^{k+1}_{i}$, that is :
595 %EAJ: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
596 %{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.
599 %With the Gauss-Seidel iteration, we have:
601 %\label{eq:Aberth-H-GS}
602 %EAGS: z^{k+1}_{i}=\frac{p(z^{k}_{i})}{p'(z^{k}_{i})-p(z^{k}_{i})(\sum^{i-1}_{j=1}\frac{1}{z^{k}_{i}-z^{k+1}_{j}}+\sum^{n}_{j=i+1}\frac{1}{z^{k}_{i}-z^{k}_{j}})}, i=1,...,n.
606 %\label{eq:Aberth-H-GS}
607 %%EAGS: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
608 %%{1-\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}(\sum^{i-1}_{j=1}\frac{1}{z^{k}_{i}-z^{k+1}_{j}}+\sum_{j=i+1}^{j=n}{\frac{1}{(z_{i}^{k}-z_{j}^{k})}})}, i=1,. . . .,n
611 %Using Eq.~\ref{eq:Aberth-H-GS} to update the vector solution
612 %\textit{Z}, we expect the Gauss-Seidel iteration to converge more
613 %quickly because, just as any Jacobi algorithm (for solving linear
614 %systems of equations), it uses the freshest computed roots $z^{k+1}_{i}$.
616 %The $4^{th}$ step of the algorithm checks the convergence condition using Eq.~\ref{eq:Aberth-Conv-Cond}.
617 %Both steps 3 and 4 use 1 thread to compute all the $n$ roots on CPU, which is very harmful for performance in case of the large degree polynomials.
621 %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.
622 %In the GPU, the scheduler assigns the execution of this loop to a
623 %group of threads organised as a grid of blocks with block containing a
624 %number of threads. All threads within a block are executed
625 %concurrently in parallel. The instructions run on the GPU are grouped
626 %in special function called kernels. With CUDA, a programmer must
627 %describe the kernel execution context: the size of the Grid, the number of blocks and the number of threads per block.
629 %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.
631 Algorithm~\ref{alg2-cuda} defines the main key points for finding roots polynomials with Ehrlich-Aberth method on GPU. Where $P$, $Pu$ and $Z$ are, respectively, the polynomial to solve, Derivative of P and the solution root vector.
633 %%Algorithm~\ref{alg2-cuda} shows a sketch of the Ehrlich-Aberth method using CUDA.
639 \caption{CUDA Algorithm to find roots with the Ehrlich-Aberth method}
641 \KwIn{$Z^{0}$ (Initial root's vector), $\varepsilon$ (Error tolerance
642 threshold), P (Polynomial to solve), Pu (Derivative of P), $n$ (Polynomial degrees), $\Delta z_{max}$ (Maximum value of stop condition)}
644 \KwOut {$Z$ (Solution root's vector), $ZPrec$ (Previous solution root's vector)}
648 %\item Initialization of P\;
649 %\item Initialization of Pu\;
650 %\item Initialization of the solution vector $Z^{0}$\;
651 \item Initialization of the parameters of roots finding problem (P, Pu, $Z^{0}$);
652 \item Allocate and copy initial data to the GPU global memory\;
654 \While {$\Delta z_{max} > \varepsilon$}{
655 \item Let $\Delta z_{max}=0$\;
656 \item $ save(ZPrec,Z)$\;
658 \item $ Find(Z,P,Pu)$\;
659 \item $testConverge(\Delta z_{max},Z,ZPrec)$\;
662 \item Copy results from GPU memory to CPU memory\;
667 After the initialization step, all data of the root finding problem
668 must be copied from the CPU memory to the GPU global memory, because
669 the GPUs only work on the data filled in their memories. Next, the algorithm uses an iterative method for finding polynomial root,defined in the function \textit{Find()} in line 7. The iterative method used in this algorithm is the Ehrlich-Aberth method corresponding to Eq.~\ref{Eq:EA}, At every time step, the initial guess for the iterative method is set to the solution found
670 at the previous time step ($ZPrec$ )defined in the function \textit{Save()} in line 7.
671 The iterative function terminates its computations when the error tolerance
672 threshold, $\varepsilon$ have been achieved, and/or all the roots have converged. Finally, the solution of the root finding problem is copied back from the GPU global memory to the CPU memory. All the data-parallel arithmetic operations inside the main loop \verb=(while(...))= are executed as kernels by the GPU.
674 The Ehrlich-Aberth is an iterative method for finding root polynomial. it is based on arithmetic vector operations that are easy to implement on parallel computers and, thus, on GPUs. Indeed, the GPU executes the vector operations as kernels and the CPU executes the sequential operations, launches the kernels and supplies the GPU with data.
676 % Algorithm 2 shows the main key points of the iterative function. %All the data-parallel arithmetic operations inside the main loop (repeat ... until(...)) are executed as kernels by the GPU.
678 In order to implement the Ehrlich-Aberth method in CUDA, it is
679 possible to use the Jacobi scheme or the Gauss-Seidel one. With the Jacobi iteration, at iteration $k+1$ we need all the previous values $z^{k}_{i}$ to compute the new values $z^{k+1}_{i}$, that is :
682 EAJ: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
683 {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.
686 With the Gauss-Seidel iteration, we have:
689 \label{eq:Aberth-H-GS}
690 EAGS: z^{k+1}_{i}=z_{i}^{k}-\frac{\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}}
691 {1-\frac{p(z_{i}^{k})}{p'(z_{i}^{k})}(\sum^{i-1}_{j=1}\frac{1}{z^{k}_{i}-z^{k+1}_{j}}+\sum_{j=i+1}^{j=n}{\frac{1}{(z_{i}^{k}-z_{j}^{k})}})}, i=1,. . . .,n
694 Using Eq.~\ref{eq:Aberth-H-GS} to update the vector solution
695 \textit{Z}, we expect the Gauss-Seidel iteration to converge more
696 quickly because, just as any Jacobi algorithm (for solving linear
697 systems of equations), it uses the freshest computed roots $z^{k+1}_{i}$.
699 Algorithm~\ref{alg3-update} shows the main key points of the iterative function implemented as a kernel.
703 %Next, all the data-parallel arithmetic operations inside the main loop
704 %\verb=(while(...))= are executed as kernels by the GPU. The
705 %first kernel named \textit{save} in line 7 of
706 %Algorithm~\ref{alg2-cuda} consists in saving the vector of
707 %polynomial roots found at the previous time-step in GPU memory, in
708 %order to check the convergence of the roots after each iteration (line
709 %10, Algorithm~\ref{alg2-cuda}).
711 %The second kernel executes the iterative function and updates
712 %$Z$, according to Algorithm~\ref{alg3-update}. We notice that the
713 %update kernel is called in two forms, according to the value
714 %\emph{R} which determines the radius beyond which we apply the
715 %exponential logarithm algorithm.
720 \caption{Kernel update for the iterative function}
722 \eIf{$(\left|Z\right|<= R)$}{
723 $kernel\_update(Z,P,Pu)$\;}
725 $kernel\_update\_ExpoLog(Z,P,Pu)$\;
729 We notice that the update kernel is called in two forms, according to the value
730 \emph{R} which determines the radius beyond which we apply the exponential logarithm algorithm.
731 If the modulus of the current complex is less than a given value called the
732 radius i.e. ($ |z^{k}_{i}|<= R$), then the classical form of the EA
733 function Eq.~\ref{Eq:Hi} is executed, else the EA.EL function Eq.~\ref{Log_H2} is executed.
734 (with Eq.~\ref{deflncomplex}, Eq.~\ref{defexpcomplex}). The radius $R$ is evaluated as in Eq.~\ref{R.EL}. \KG{experimentally it is difficult to solve the high degree polynomial with the classical Ehrlich-Aberth method, so if the root are under to the unit circle ($R$)the kernel \textit{update} is called in the EA.EL function Eq.~\ref{Log_H2} (with Eq.~\ref{deflncomplex}, Eq.~\ref{defexpcomplex}) line 3, to put into account the limited of grander floating manipulated by processors and compute more roots}.
735 \KG{} We notice that we used \verb= cuDoubleComplex= to exploit the complex number in CUDA, and the functions of the CUBLAS library to implement some vector operations on the GPU. We use the following functions:
738 \item \verb= cublasIdamax()= for the
739 \item \verb= cublasGetVector()= for the
743 \KG{Listing ~\ref{lst:01} shows a simplified version of second kernel
744 code (some parameters in the kernels have been simplified in order to
745 increase the readability). As can be seen this
746 kernel calls multiple kernels, all the kernels for complex numbers and
747 kernels for the evaluation of a polynomial are not detailed.}
750 \lstinputlisting[label=lst:01,caption=Kernels to update the roots]{code.c}
752 \KG{}The Kernel \verb= __global__EA_update()= is a code called from the host and executed on the device, need to call \verb=__device__FirstH_EA()= on the device to compute the root who are under the unit circle ($R$), the \verb= __device__NewH_EA()= is called to compute root with the exp-log version. The Horner schema are used to evaluate the polynomial and his derivative in $Z[i]$ \verb=__device__Fonction()= and \verb=__deviceFonctionD()= respectively. Its exp-log version need to be implemented \verb= __device__LogFonction()=, \verb= __device__LogFonctionD()=
753 and used in the exp-log version of the Ehrlich-Aberth method
755 \item \verb= __device__LogFonction()= to
756 \item \verb= __device__LogFonctionD()= for the
758 %need to call some function on device
760 %The Kernel \verb= EA_update= is a code implemented to executed on GPU and lanced in CPU,
761 This kernel is executed by a large number of GPU threads such that each thread is in charge of the computation of one component of the iterate vector $Z$. We set the size of a thread block, \textit{Threads}, to 512 threads and the number of thread blocks of a kernel, \textit{Blocks}, is computed so as each GPU thread is in charge of one vector element:
762 \[ Blocks=\frac{N+Threads-1}{Threads}, N: Polynomial size \]
763 %$ Blocks=\frac{N+Threads-1}{Threads}, N: Polynomial size$
765 Each GPU threads in grid compute one root en parallel, if the polynomial size exceed the capacity of the grid the G.S schema are finely executed, like the grid can only compute << Blocks,Threads>> roots at the same time, if we need to compute more roots, the grid can used the roots previously executed to compute other root ih the same iteration, like the following schema:
767 %\begin{figure}[htbp]
769 % \includegraphics[width=0.8\textwidth]{figures/G.S}
770 %\caption{Gauss Seidel iteration}
774 The last kernel checks the convergence of the roots after each update
775 of $Z^{k}$, according to formula Eq.~\ref{eq:Aberth-Conv-Cond}. We used the functions of the CUBLAS Library (CUDA Basic Linear Algebra Subroutines) to implement this kernel.
776 This listing shows the kernel \textit{TestConvergence} code
778 \lstinputlisting[label=lst:01,caption=Kernel to test the convergence of the roots]{code1.c}
781 The kernel terminates its computations when all the roots have
782 converged. It should be noticed that, as blocks of threads are
783 scheduled automatically by the GPU, we have absolutely no control on
784 the order of the blocks. Consequently, our algorithm is executed more
785 or less with the asynchronous iteration model, where blocks of roots
786 are updated in a non deterministic way. As the Durand-Kerner method
787 has been proved to converge with asynchronous iterations \KG{Ajouter la reference qui montre que DK converge with asynchronous iteration}, we think it
788 is similar with the Ehrlich-Aberth method, but we did not try to prove
789 this in that paper. Another consequence of that, is that several
790 executions of our algorithm with the same polynomial do not
791 necessarily give the same result (but roots have the same accuracy)
792 and the same number of iterations (even if the variation is not very
799 %%HIER END MY REVISIONS (SIDER)
800 \section{Experimental study}
802 %\subsection{Definition of the used polynomials }
803 We study two categories of polynomials: sparse polynomials and the full polynomials.\\
804 {\it A sparse polynomial} is a polynomial for which only some
805 coefficients are not null. In this paper, we consider sparse polynomials for which the roots are distributed on 2 distinct circles:
807 \forall \alpha_{1} \alpha_{2} \in C,\forall n_{1},n_{2} \in N^{*}; P(z)= (z^{n_{1}}-\alpha_{1})(z^{n_{2}}-\alpha_{2})
808 \end{equation}\noindent
809 {\it A full polynomial} is, in contrast, a polynomial for which
810 all the coefficients are not null. A full polynomial is defined by:
812 %%\forall \alpha_{i} \in C,\forall n_{i}\in N^{*}; P(z)= \sum^{n}_{i=1}(z^{n^{i}}.a_{i})
816 {\Large \forall a_{i} \in C, i\in N; p(x)=\sum^{n}_{i=0} a_{i}.x^{i}}
818 %With this form, we can have until \textit{n} non zero terms whereas the sparse ones have just two non zero terms.
820 %\subsection{The study condition}
821 %Two parameters are studied are
822 %the polynomial degree and the execution time of our program
823 %to converge on the solution. The polynomial degree allows us
824 %to validate that our algorithm is powerful with high degree
825 %polynomials. The execution time remains the
826 %element-key which justifies our work of parallelization.
827 For our tests, a CPU Intel(R) Xeon(R) CPU
828 E5620@2.40GHz and a GPU K40 (with 6 Go of ram) are used.
831 %\subsection{Comparative study}
832 %First, performances of the Ehrlich-Aberth method of root finding polynomials
833 %implemented on CPUs and on GPUs are studied.
835 We performed a set of experiments on the sequential and the parallel algorithms, for both sparse and full polynomials of different sizes. We took into account the execution times, the polynomial size and the number of threads per block performed by sum or each experiment on CPU and on GPU.
837 All experimental results obtained from the simulations are made in
838 double precision data, the convergence threshold of the methods is set
840 %Since we were more interested in the comparison of the
841 %performance behaviors of Ehrlich-Aberth and Durand-Kerner methods on
842 %CPUs versus on GPUs.
843 The initialization values of the vector solution
844 of the methods are given in Section~\ref{sec:vec_initialization}.
846 \subsection{Comparison of execution times of the Ehrlich-Aberth method
847 on a CPU with OpenMP (1 core and 4 cores) vs. on a Tesla GPU}
851 \includegraphics[width=0.8\textwidth]{figures/openMP-GPU}
852 \caption{Comparison of execution times of the Ehrlich-Aberth method
853 on a CPU with OpenMP (1 core, 4 cores) and on a Tesla GPU}
856 %%Figure 1 %%show a comparison of execution time between the parallel
857 %%and sequential version of the Ehrlich-Aberth algorithm with sparse
858 %%polynomial exceed 100000,
860 In Figure~\ref{fig:01}, we report the execution times of the
861 Ehrlich-Aberth method on one core of a Quad-Core Xeon E5620 CPU, on
862 four cores on the same machine with \textit{OpenMP} and on a Nvidia
863 Tesla K40 GPU. We chose different sparse polynomials with degrees
864 ranging from 100,000 to 1,000,000. We can see that the implementation
865 on the GPU is faster than those implemented on the CPU.
866 However, the execution time for the
867 CPU (4 cores) implementation exceed 5,000s for 250,000 degrees
868 polynomials. On the other hand, the GPU implementation for the same
869 polynomials do not take more 100s. With the GPU
870 we can solve high degree polynomials very quickly up to degree 1,000,000. We can also notice that the GPU implementation are
871 almost 40 times faster then the implementation on the CPU (4 cores).
876 %This reduction of time allows us to compute roots of polynomial of more important degree at the same time than with a CPU.
878 %We notice that the convergence precision is a round $10^{-7}$ for the both implementation on CPU and GPU. Consequently, we can conclude that Ehrlich-Aberth on GPU are faster and accurately then CPU implementation.
880 \subsection{Influence of the number of threads on the execution times of different polynomials (sparse and full)}
881 To optimize the performances of an algorithm on a GPU, it is necessary to maximize the use of cores GPU (maximize the number of threads executed in parallel). In fact, it is interesting to see the influence of the number of threads per block on the execution time of Ehrlich-Aberth algorithm.
882 For that, we noticed that the maximum number of threads per block for
883 the Nvidia Tesla K40 GPU is 1,024, so we varied the number of threads
884 per block from 8 to 1,024. We took into account the execution time for
885 10 different sparse and full polynomials of degree 50,000 and of degree 500,000.
889 \includegraphics[width=0.8\textwidth]{figures/influence_nb_threads}
890 \caption{Influence of the number of threads on the execution times of different polynomials (sparse and full)}
894 Figure~\ref{fig:02} shows that, the best execution time for both
895 sparse and full polynomial are given when the threads number varies
896 between 64 and 256 threads per block. We notice that with small
897 polynomials the best number of threads per block is 64, whereas for large polynomials the best number of threads per block is
898 256. However, in the following experiments we specify that the number
899 of threads per block is 256.
902 \subsection{Influence of exp-log solution to compute high degree polynomials}
904 In this experiment we report the performance of the exp-log solution described in Section~\ref{sec2} to compute high degree polynomials.
907 \includegraphics[width=0.8\textwidth]{figures/sparse_full_explog}
908 \caption{The impact of exp-log solution to compute high degree polynomials}
913 Figure~\ref{fig:03} shows a comparison between the execution time of
914 the Ehrlich-Aberth method using the exp-log solution and the
915 execution time of the Ehrlich-Aberth method without this solution,
916 with full and sparse polynomials degrees. We can see that the
917 execution times for both algorithms are the same with full polynomials
918 degree inferior to 4,000 and sparse polynomials inferior to 150,000. We
919 also clearly show that the classical version (without exp-log) of
920 Ehrlich-Aberth algorithm does not converge after these degrees with
921 sparse and full polynomials. On the contrary, the new version of
922 the Ehrlich-Aberth algorithm with the exp-log solution can solve
923 high degree polynomials.
925 %in fact, when the modulus of the roots are up than \textit{R} given in ~\ref{R},this exceed the limited number in the mantissa of floating points representations and can not compute the iterative function given in ~\ref{eq:Aberth-H-GS} to obtain the root solution, who justify the divergence of the classical Ehrlich-Aberth algorithm. However, applying exp-log solution given in ~\ref{sec2} took into account the limit of floating using the iterative function in(Eq.~\ref{Log_H1},Eq.~\ref{Log_H2} and allows to solve a very large polynomials degrees .
930 \subsection{Comparison of the Durand-Kerner and the Ehrlich-Aberth methods}
932 In this part, we compare the Durand-Kerner and the Ehrlich-Aberth
933 methods on GPU. We took into account the execution times, the number of iterations and the polynomials size for both sparse and full polynomials.
937 \includegraphics[width=0.8\textwidth]{figures/EA_DK}
938 \caption{Execution times of the Durand-Kerner and the Ehrlich-Aberth methods on GPU}
942 Figure~\ref{fig:04} shows the execution times of both methods with
943 sparse polynomial degrees ranging from 1,000 to 1,000,000. We can see
944 that the Ehrlich-Aberth algorithm is faster than Durand-Kerner
945 algorithm, being on average 25 times faster. Then, when degrees of
946 polynomials exceed 500,000 the execution times with DK are very long.
948 %with double precision not exceed $10^{-5}$.
952 \includegraphics[width=0.8\textwidth]{figures/EA_DK_nbr}
953 \caption{The number of iterations to converge for the Ehrlich-Aberth
954 and the Durand-Kerner methods}
958 Figure~\ref{fig:05} shows the evaluation of the number of iterations according
959 to the degree of polynomials for both EA and DK algorithms. We can see
960 that the number of iterations of DK is of order 100 while EA is of order
961 10. Indeed the computation of the derivative of P in the iterative function (Eq.~\ref{Eq:Hi}) executed by EA
962 allows the algorithm to converge faster. On the contrary, the
963 DK operator (Eq.~\ref{DK}) needs low operations, consequently low
964 execution times per iteration, but it needs more iterations to converge.
969 \section{Conclusion and perspectives}
971 In this paper we have presented the parallel implementation
972 Ehrlich-Aberth method on GPU for the problem of finding roots
973 polynomial. Moreover, we have improved the classical Ehrlich-Aberth
974 method which suffers from overflow problems, the exp-log solution
975 applied to the iterative function allows to solve high degree
978 We have performed many experiments with the Ehrlich-Aberth method in
979 GPU. These experiments highlight that this method is more efficient in
980 GPU than all the other implementations. The improvement with
981 the exponential logarithm solution allows us to solve sparse and full
982 high degree polynomials up to 1,000,000 degree. Hence, it may be
983 possible to consider using polynomial root finding methods in other
984 numerical applications on GPU.
987 In future works, we plan to investigate the possibility of using
988 several multiple GPUs simultaneously, either with a multi-GPU machine or
989 with a cluster of GPUs. It may also be interesting to study the
990 implementation of other root finding polynomial methods on GPU.
994 \bibliography{mybibfile}