]> AND Private Git Repository - kahina_paper1.git/blob - Root.tex
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
this is the commit of the parallel Aberth algorithm
[kahina_paper1.git] / Root.tex
1 \documentclass[11pt,a4paper]{article}\r
2 \usepackage[latin1]{inputenc}\r
3 \usepackage{amsmath}\r
4 \usepackage{amsfonts}\r
5 \usepackage{amssymb}\r
6 %%\usepackage{algorithm2e}\r
7 \usepackage[ruled,vlined]{algorithm2e}\r
8 %%\usepackage{algo}\r
9 \author{ghidouche}\r
10 \title{Paper1_kahina}\r
11 \begin{document}\r
12 \section{Root finding problem}\r
13 We consider a polynomial of degree \textit{n} having coefficients\r
14 in the complex \textit{C} and zeros $\alpha_{i},\textit{i=1,...,n}$. \\\r
15 \begin{center}\r
16 \begin{equation}\r
17      {\Large p(x)=\sum{a_{i}x^{i}}=a_{n}\prod(x-\alpha_{i}),a_{0} a_{n}\neq 0}\r
18 \end{equation}\r
19 \end{center}\r
20 \r
21  the root finding problem consist to find\r
22 all n root of \textit{p(x)}. the problem of finding a root is\r
23 equivalent to the problem of finding a fixed-point. To see this\r
24 consider the fixed-point problem of finding the n-dimensional\r
25 vector x such that\r
26 \begin{center}\r
27 $x=g(x).  $\r
28 \end{center}\r
29 where $g : C^{n}\longrightarrow C^{n}$. Note that we can easily\r
30 rewrite this fixed-point problem as a root-finding problem by\r
31 setting $f (x) = x-g(x)$ and likewise we can recast the\r
32 root-finding problem into a fixed-point problem by setting\r
33 \begin{center}\r
34 $g(x)= f(x)-x$\r
35 \end{center}\r
36 Often it will not be possible to solve such nonlinear equation\r
37 root-finding problems analytically. When this occurs we turn to\r
38 numerical methods to approximate the solution. Generally speaking,\r
39 algorithms for solving problems numerically can be divided into\r
40 two main groups: direct methods and iterative methods.\r
41 \\\r
42  Direct methods exist only for $n \leqslant4$,solved in closed form by G. Cardano\r
43 in the mid-16th century. However, N.H. Abel in the early 19th\r
44 century showed that polynomials of degree five or more could not\r
45 be solved by  directs methods. Since then researchers have\r
46 concentrated on numerical (iterative) methods such as the famous\r
47 Newton s method, Bernoulli s method of the 18th, and Graeffe s.\r
48 With the advent of electronic computers, different methods has\r
49 been developed such as the Jenkins-Traub method, Larkin s method,\r
50 Muller s method, and several methods for simultaneous\r
51 approximation of all the roots, starting with the Durand-Kerner\r
52 method:\r
53 \begin{center}\r
54 \begin{equation} Z_{i}=Z_{i}-\frac{P(Z_{i})}{\prod_{i\neq j}(z_{i}-z_{j})}\r
55 \end{equation}\r
56 \end{center}\r
57 \r
58 This formula is mentioned for the first time from\r
59 Weiestrass~\cite{Weierstrass03} as part of the fundamental theorem\r
60 of Algebra and is rediscovered from Ilieff~\cite{Ilie50},\r
61 Docev~\cite{Docev62}, Durand~\cite{Durand60},\r
62 Kerner~\cite{Kerner66}. Another method discovered from\r
63 Borsch-Supan~\cite{ Borch-Supan63} and also described and brought\r
64 in the following form from Ehrlich~\cite{Ehrlich67} and\r
65 Aberth~\cite{Aberth73}.\r
66 \begin{center}\r
67 \begin{equation}\r
68  Z_{i}=Z_{i}-\frac{1}{{\frac {P'(Z_{i})} {P(Z_{i})}}-{\sum_{i\neq j}(z_{i}-z_{j})}}\r
69 \end{equation}\r
70 \end{center}\r
71 \r
72 Aberth, Ehrlich and Farmer-Loizou~\cite{Loizon83} have proved that\r
73 the above method has cubic order of convergence for simple roots.\r
74 \r
75 \r
76 Iterative methods raise several problem when implemented e.g.\r
77 specific sizes of numbers must be used to deal with this\r
78 difficulty.Moreover,the convergence time of iterative methods\r
79 drastically increase like the degrees of high polynomials. The\r
80 parallelization of these algorithms will improve the convergence\r
81 time.\r
82 \r
83 Many authors have treated the problem of parallelization of\r
84 simultaneous methods. Freeman~\cite{Freeman89} has tested the DK\r
85 method, EA method and another method of the fourth order proposed\r
86 from Farmer and Loizou~\cite{Loizon83},on a 8- processor linear\r
87 chain, for polynomial of degree up to 8. The third method often\r
88 diverges, but the first two methods have speed-up 5.5\r
89 (speed-up=(Time on one processor)/(Time on p processors)). Later\r
90 Freeman and Bane~\cite{Freemanall90}  consider asynchronous\r
91 algorithms, in which each processor continues to update its\r
92 approximations even although the latest values of other $z_i((k))$\r
93 have not received from the other processors, in difference with\r
94 the synchronous version where it would wait.\r
95 in~\cite{Raphaelall01}proposed two methods of parallelization for\r
96 architecture with shared memory and distributed memory,it able to\r
97 compute the root of polynomial degree  10000 on 430 s with only 8\r
98 pc and 2 communications per iteration. Compare to the sequential\r
99 it take 3300 s to obtain the same results.\r
100 \r
101 After this few works discuses this problem until the apparition of\r
102 the Compute Unified Device Architecture (CUDA)~\cite{CUDA10},a\r
103 parallel computing platform and a programming model invented by\r
104 NVIDIA. the computing ability of GPU has exceeded the counterpart\r
105 of CPU. It is a waste of resource to be just a graphics card for\r
106 GPU.  CUDA adopts a totally new computing architecture to use the\r
107 hardware resources provided by GPU in order to offer a stronger\r
108 computing ability to the massive data computing.\r
109 \r
110 \r
111 Indeed,~\cite{Kahinall14}proposed the implementation of the\r
112 Durand-Kerner method on GPU (Graphics Processing Unit). The main\r
113 result prove that a parallel implementation is 10 times as fast as\r
114 the sequential implementation on a single CPU for high degree\r
115 polynomials that is greater than about 48000.\r
116 \paragraph{}\r
117 The mean part of our work is to implement the Aberth method on GPU\r
118 and compare it with the Durand Kerner\r
119 implementation.................To be continued..................\r
120 \r
121 \r
122 \section{Aberth method and difficulties}\r
123 A cubically convergent iteration method for finding zeros of\r
124 polynomials was proposed by O.Aberth~\cite{Aberth73}.The Aberth\r
125 method is a purely algebraic derivation.To illustrate the\r
126 derivation, we let $w_{i}(z)$ be the product of linear factor $\r
127 w_{i}(z)=\prod_{j=1,j \neq i}^{n} (z-x_{j})$\r
128 \r
129 and rational function $R_{i}(z)$ be the correction term of\r
130 Weistrass method~\cite{Weierstrass03}:\r
131 \r
132 \begin{equation}\r
133 R_{i}(z)=\dfrac{p(z)}{w_{i}(Z)} , i=1,2,...,n.\r
134 \end{equation}\r
135 \r
136 Differentiating the rational function $R_{i}(z)$ and applying the\r
137 Newton method, we have\r
138 \r
139 \begin{equation}\r
140 \dfrac{R_{i}(z)}{R_{i}^{'}(z)}=\r
141 \dfrac{p(z)}{p^{'}(z)-p(z)\dfrac{w_{i}(z)}{w_{i}^{'}(z)}}=\r
142 \dfrac{p(z)}{p^{'}(z)-p(z) \sum _{j=1,j \neq\r
143 i}^{n}\dfrac{1}{z-x_{i}}}, i=1,2,...,n\r
144 \end{equation}\r
145 \r
146 Substituting $x_{j}$ for z we obtain the Aberth iteration method\r
147 \r
148 Let present the means stages of Aberth's method.\r
149 \r
150 \subsection{Polynomials Initialization}\r
151  The initialization of polynomial P(z) with complex coefficients\r
152  are given by:\r
153 \r
154 \begin{equation}\r
155   p(z)=\sum{a_{i}z^{n-i}}. where a_{n} \neq 0,a_{0}=1, a_{i}\subset C\r
156 \end{equation}\r
157 \r
158 \r
159 \subsection{Vector $Z^{0)}$ Initialization}\r
160 \r
161 The choice of the initial points $z^{(0)}_{i} , i = 1, . . . , n,$\r
162 from which starting the iteration  (2) or (3), is rather delicate\r
163 since the number of steps needed by the iterative method to reach\r
164 a given approximation strongly depends on it.\r
165 In~\cite{Aberth73}the Aberth iteration is started by selecting n\r
166 equispaced points on a circle of center 0 and radius r, where r is\r
167 an upper bound to the moduli of the zeros. After,~\cite{Bini96}\r
168 performs this choice by selecting complex numbers along different\r
169 circles and relies on the result of~\cite{Ostrowski41}.\r
170 \r
171 \begin{equation}\r
172 \sigma_{0}=\frac{u+v}{2};u=\frac{\sum_{i=1}^{n}u_{i}}{n.max_{i=1}^{n}u_{i}};\r
173 v=\frac{\sum_{i=0}^{n-1}v_{i}}{n.min_{i=0}^{n-1}v_{i}};u_{i}=2.|a_{i}|^{\frac{1}{i}};\r
174 v_{i}=\frac{|\frac{a_{n}}{a_{i}}|^{\frac{1}{n-i}}}{2}\r
175 \end{equation}\r
176 \r
177 \subsection{Iterative Function Hi}\r
178 The operator used with Aberth's method is corresponding to the\r
179 following equation which will enable the convergence towards\r
180 polynomial solutions, provided all the roots are distinct.\r
181 \r
182 \begin{equation}\r
183 H_{i}(z)=z_{i}-\frac{1}{\frac{P^{'}(z_{i})}{P(z_{i})}-\sum_{j\neq\r
184 i}{\frac{1}{z_{i}-z_{j}}}}\r
185 \end{equation}\r
186 \r
187 \subsection{Convergence condition}\r
188 determines the success of the termination. It consists in stopping\r
189 the iterative function $H_{i}(z)$ when the are stable,the method\r
190 converge sufficiently:\r
191 \r
192 \begin{equation}\r
193 \forall i \in\r
194 [1,n];\frac{z_{i}^{(k)}-z_{i}^{(k-1)}}{z_{i}^{(k)}}<\xi\r
195 \end{equation}\r
196 \r
197 \r
198 \section{Difficulties and amelioration}\r
199 the Aberth method implementation suffer of overflow problems. This\r
200 situation occurs, for instance, in the case where a polynomial\r
201 having positive coefficients and large degree is computed at a\r
202 point $\xi$ where $|\xi| > 1$.Indeed the limited number in the\r
203 mantissa of floating takings the computation of P(z) wrong when z\r
204 is large. for example $(10^{50}) +1+ (- 10_{50})$ will give result\r
205 0 instead of 1 in reality.consequently we can't compute the roots\r
206 for large polynomial's degree. This problem was discuss in\r
207 ~\cite{Karimall98} for the Durand-Kerner method, the authors\r
208 propose to use the logratihm and the exponential of a complex:\r
209 \r
210 \begin{equation}\r
211  \forall(x,y)\in R^{*2}; \ln (x+i.y)=\ln(x^{2}+y^{2})\r
212 2+i.\arcsin(y\sqrt{x^{2}+y^{2}})_{\left] -\pi, \pi\right[ }\r
213 \end{equation}\r
214 %%\begin{equation}\r
215 \begin{align}\r
216  \forall(x,y)\in R^{*2}; \exp(x+i.y)&= \exp(x).\exp(i.y)\\\r
217                                     &=\exp(x).\cos(y)+i.\exp(x).\sin(y)\r
218 \end{align}\r
219 %%\end{equation}\r
220 \r
221 The application of logarithm can replace any multiplications and\r
222 divisions with additions and subtractions; consequently it\r
223 manipulates lower absolute values and can be compute the roots for\r
224 large polynomial's degree exceed~\cite{Karimall98}.\r
225 \r
226 Applying this solution for the Aberth method we obtain the\r
227 iteration function with logarithm:\r
228 %%$$ \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}})$$\r
229 \begin{equation}\r
230 H_{i}(z)=z_{i}^{k}-\exp \left(\ln \left(\r
231 p(z_{k})\right)-\ln\left(p(z_{k}^{'})\right)- \ln\r
232 \left(1-Q(z_{k})\right)\right)\r
233 \end{equation}\r
234 where:\r
235 \r
236 \begin{equation}\r
237 Q(z_{k})=\exp\left( \ln (p(z_{k}))-\ln(p(z_{k}^{'}))+\ln \left(\r
238 \sum_{k\neq j}^{n}\frac{1}{z_{k}-z_{j}}\right)\right)\r
239 \end{equation}\r
240 \r
241 \r
242 this solution is applying when it is necessary\r
243 \r
244 \section{The implementation of simultaneous methods in a parallel computer}\r
245     The main problem of the simultaneous methods is that the necessary\r
246 time needed for the convergence is increased with the increasing\r
247 of the degree of the polynomial. The parallelization of these\r
248 algorithms will improve the convergence time. Researchers usually\r
249 adopt one of the two following approaches to parallelize root\r
250 finding algorithms. One approach is to reduce the total number of\r
251 iterations as implemented by Miranker\r
252 ~\cite{Mirankar68,Mirankar71}, Schedler~\cite{Schedler72} and\r
253 Winogard~\cite{Winogard72}. Another approach is to reduce the\r
254 computation time per iteration, as reported\r
255 in~\cite{Benall68,Jana06,Janall99,Riceall06}. There are many\r
256 schemes for simultaneous approximations of all roots of a given\r
257 polynomial. Several works on different methods and issues of root\r
258 finding have been reported in~\cite{Azad07,Gemignani07,Kalantari08\r
259 ,Skachek08,Zhancall08,Zhuall08}. However, Durand-Kerner and\r
260 Ehrlich methods are the most practical choices among\r
261 them~\cite{Bini04}. These two methods have been extensively\r
262 studied for parallelization due to their following advantages. The\r
263 computation involved in these methods has some inherent\r
264 parallelism that can be suitably exploited by SIMD machines.\r
265 Moreover, they have fast rate of convergence (quadratic for the\r
266 Durand-Kerner method and cubic for the Ehrlich). Various parallel\r
267 algorithms reported for these methods can be found\r
268 in~\cite{Cosnard90, Freeman89,Freemanall90,,Jana99,Janall99}.\r
269 Freeman and Bane~\cite{Freemanall90} presented two parallel\r
270 algorithms on a local memory MIMD computer with the compute-to\r
271 communication time ratio O(n). However, their algorithms require\r
272 each processor to communicate its current approximation to all\r
273 other processors at the end of each iteration. Therefore they\r
274 cause a high degree of memory conflict. Recently the author\r
275 in~\cite{Mirankar71} proposed two versions of parallel algorithm\r
276 for the Durand-Kerner method, and Aberth method on an on model of\r
277 Optoelectronic Transpose Interconnection System (OTIS).The\r
278 algorithms are mapped on an OTIS-2D torus using N processors. This\r
279 solution need N processors to compute N roots, that it is not\r
280 practical (is not suitable to compute large polynomial's degrees).\r
281 Until then, the related works are not able to compute the root of\r
282 the large polynomial's degrees (higher then 1000) and with small\r
283 time.\r
284 \r
285     Finding polynomial roots rapidly and accurately it is our\r
286 objective, with the apparition of the CUDA(Compute Unified Device\r
287 Architecture), finding the roots of polynomials becomes rewarding\r
288 and very interesting, CUDA adopts a totally new computing\r
289 architecture to use the hardware resources provided by GPU in\r
290 order to offer a stronger computing ability to the massive data\r
291 computing.in~\cite{Kahinall14} we proposed the first implantation\r
292 of the root finding polynomials method on GPU (Graphics Processing\r
293 Unit),which is the Durand-Kerner method. The main result prove\r
294 that a parallel implementation is 10 times as fast as the\r
295 sequential implementation on a single CPU for high degree\r
296 polynomials that is greater than about 48000. Indeed, in this\r
297 paper we present a parallel implementation of Aberth's method on\r
298 GPU, more details are discussed in the following of this paper.\r
299 \r
300 \section {A parallel implementation of Aberth's method}\r
301 \subsection{Background on the GPU architecture}\r
302 A GPU is viewed as an accelerator for the data-parallel and\r
303 intensive arithmetic computations. It draws its computing power\r
304 from the parallel nature of its hardware and software\r
305 architectures. A GPU is composed of hundreds of Streaming\r
306 Processors (SPs) organized in several blocks called Streaming\r
307 Multiprocessors (SMs). It also has a memory hierarchy. It has a\r
308 private read-write local memory per SP, fast shared memory and\r
309 read-only constant and texture caches per SM and a read-write\r
310 global memory shared by all its SPs~\cite{NVIDIA10}\r
311 \r
312     On a CPU equipped with a GPU, all the data-parallel and intensive\r
313 functions of an application running on the CPU are off-loaded onto\r
314 the GPU in order to accelerate their computations. A similar\r
315 data-parallel function is executed on a GPU as a kernel by\r
316 thousands or even millions of parallel threads, grouped together\r
317 as a grid of thread blocks. Therefore, each SM of the GPU executes\r
318 one or more thread blocks in SIMD fashion (Single  Instruction,\r
319 Multiple Data) and in turn each SP of a GPU SM runs one or more\r
320 threads within a block in SIMT fashion (Single Instruction,\r
321 Multiple threads). Indeed at any given clock cycle, the threads\r
322 execute the same instruction of a kernel, but each of them\r
323 operates on different data.\r
324  GPUs only work on data filled in their\r
325 global memories and the final results of their kernel executions\r
326 must be communicated to their CPUs. Hence, the data must be\r
327 transferred in and out of the GPU. However, the speed of memory\r
328 copy between the GPU and the CPU is slower than the memory\r
329 bandwidths of the GPU memories and, thus, it dramatically affects\r
330 the performances of GPU computations. Accordingly, it is necessary\r
331 to limit data transfers between the GPU and its CPU during the\r
332 computations.\r
333 \subsection{Background on the CUDA Programming Model}\r
334 \r
335 The CUDA programming model is similar in style to a single program\r
336 multiple-data (SPMD) softwaremodel. The GPU is treated as a\r
337 coprocessor that executes data-parallel kernel functions. CUDA\r
338 provides three key abstractions, a hierarchy of thread groups,\r
339 shared memories, and barrier synchronization. Threads have a three\r
340 level hierarchy. A grid is a set of thread blocks that execute a\r
341 kernel function. Each grid consists of blocks of threads. Each\r
342 block is composed of hundreds of threads. Threads within one block\r
343 can share data using shared memory and can be synchronized at a\r
344 barrier. All threads within a block are executed concurrently on a\r
345 multithreaded architecture.The programmer specifies the number of\r
346 threads per block, and the number of blocks per grid. A thread in\r
347 the CUDA programming language is much lighter weight than a thread\r
348 in traditional operating systems. A thread in CUDA typically\r
349 processes one data element at a time. The CUDA programming model\r
350 has two shared read-write memory spaces, the shared memory space\r
351 and the global memory space. The shared memory is local to a block\r
352 and the global memory space is accessible by all blocks. CUDA also\r
353 provides two read-only memory spaces, the constant space and the\r
354 texture space, which reside in external DRAM, and are accessed via\r
355 read-only caches\r
356 \r
357 \subsection{A parallel implementation of the Aberth's method }\r
358 %%\subsection{A CUDA implementation of the Aberth's method }\r
359 %%\subsection{A GPU implementation of the Aberth's method }\r
360 \r
361 \r
362 \r
363 \subsubsection{A sequential Aberth algorithm}\r
364 The means steps of Aberth's method can expressed as an algorithm\r
365 like:\r
366   \r
367 \begin{algorithm}[H]\r
368 \LinesNumbered\r
369 \caption{Algorithm to find root polynomial with Aberth method}\r
370 \r
371 \KwIn{$Z^{0}$(Initial root's vector),$\varepsilon$ (error\r
372 tolerance threshold),P(Polynomial to solve)}\r
373 \r
374 \KwOut {Z(The solution root's vector)}\r
375 \r
376 \BlankLine\r
377 \r
378 Initialization of the parameter of the polynomial to solve\;\r
379 Initialization of the solution vector $Z^{0}$\;\r
380 \r
381 \While {$\Delta z_{max}\succ \epsilon$}{\r
382  Let $\Delta z_{max}=0$\;\r
383 \For{$j \gets 0 $ \KwTo $n$}{\r
384 $ZPrec\left[j\right]=Z\left[j\right]$\;\r
385 $Z\left[j\right]=H\left(j,Z\right)$\;\r
386 }\r
387 \r
388 \For{$i \gets 0 $ \KwTo $n-1$}{\r
389 $c=\frac{\left|Z\left[i\right]-ZPrec\left[i\right]\right|}{Z\left[i\right]}$\;\r
390 \If{$c\succ\Delta z_{max}$ }{\r
391 $\Delta z_{max}$=c\;}\r
392 }\r
393 }\r
394 \end{algorithm}\r
395 ~\\ \r
396 ~\\ \r
397 In this sequential algorithm one thread CPU execute all steps, let see the step 3 the execution of the iterative function , 2 instructions are needed, the first instruction \textit{save} the solution vector for the previous iteration, the second instruction \textit{update} or compute a new values of the roots.\r
398 We have two manner to execute the iterative function, taking a Jacobi iteration who need all the previous value $z^{(k)}_{i}$ to compute the new value $z^{(k+1)}_{i}$we have:\r
399 \r
400 \begin{equation}\r
401 H(i,z^{k+1})=\frac{p(z^{(k)}_{i})}{p'(z^{(k)}_{i})-p(z^{(k)}_{i})\sum^{n}_{j=1 j\neq i}\frac{1}{z^{(k)}_{i}-z^{(k)}_{j}}}, i=1,...,n.\r
402 \end{equation}\r
403 \r
404 Or with the Gauss-seidel iteration, we have:\r
405 \begin{equation}\r
406 H(i,z^{k+1})=\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)}_{j}}+\sum^{n}_{j=i+1}\frac{1}{z^{(k)}_{i}-z^{(k)}_{j}}}, i=1,...,n.\r
407 \end{equation}\r
408 \r
409 In formula(16) the iteration function use the $z^{k+1}_{i}$ computed in the current iteration to compute the rest of the roots, which take him to converge more quickly compare to the jacobi iteration (it's well now that the Gauss-seidel iteration converge more quickly because they used the most fresh computed root, so we used Gauss-seidel iteration.)\r
410 \r
411 The steps 4 of the Aberth's method compute the convergence of the roots, using(9) formula.\r
412 Both steps 3 and 4 use 1 thread to compute N roots on CPU, which is faster and hard, it make the algorithm faster and hard for the large polynomial's roots finding.\r
413 \r
414 \paragraph{The execution time}\r
415 Let $T_{i}(N)$: the time to compute one new root's value of the step 3,$T_{i}$ depend on the polynomial's degrees N, when N increase $T_{i}$ increase to.We need $N.T_{i}(N)$ to compute all the new root's value in one iteration on the step 3.\r
416 \r
417 Let $T_{j}$: the time to compute one root's convergence value of the step 4, we need $N.T_{j}$ to compute all the root's convergence value in one iteration on the step 4.\r
418 \r
419 The execution time for both steps 3 and 4 can see like:\r
420 \begin{equation}\r
421 T_{exe}=N(T_{i}(N)+T_{j})+O(n).\r
422 \end{equation}\r
423 Let Nbr\_iter the number of iteration necessary to compute all the roots,so the total execution time $Total\_time_{exe}$ can give like:\r
424 \r
425 \begin{equation}\r
426 Total\_time_{exe}=\left[N\left(T_{i}(N)+T_{j}\right)+O(n)\right].Nbr\_iter.\r
427 \end{equation}\r
428 The execution time increase with the increasing of the polynomial's root, which take necessary to parallilize this step to reduce the execution time. In the following paper you explain how we parrallelize  this step using GPU architecture with CUDA platform.\r
429 \r
430 \subsubsection{Parralelize the steps on GPU }\r
431 On the CPU Aberth algorithm both steps 3 and 4 contain the loop \verb=for= , it use one thread to execute all the instruction in the loop N times.Here we explain how the GPU architecture can compute this loop and reduce the execution time.\r
432 The GPU architecture affect the execution of this loop to a groups of parallel threads organized as a grid of blocks each block contain a number of threads. All threads within a block are executed concurrently in parallel. the instruction are executed as a kernel.\r
433 \r
434 Let nbr\_thread be the number of threads executed in parallel, so you can easily transform the (18)formula like this: \r
435 \r
436 \begin{equation}\r
437 Total\_time_{exe}=\left[\frac{N}{nbr\_thread}\left(T_{i}(N)+T_{j}\right)+O(n)\right].Nbr\_iter.\r
438 \end{equation}\r
439 \r
440 In theory, the $Total\_time_{exe}$ on GPU is speed up nbr\_thread times as a $Total\_time_{exe}$ on CPU. We show more details in the experiment part. \r
441 ~\\\r
442 ~\\\r
443 In CUDA platform, All the instruction of the loop \verb=for= are executed by the GPU as a kernel form. A kernel is a procedure written in CUDA and defined by a heading \verb=__global__=, which means that it is to be executed by the GPU.the following algorithm see the kernels created for the step 3 and 4 of the Aberth algorithm:\r
444 \r
445 \begin{algorithm}[H]\r
446 \LinesNumbered\r
447 \caption{Algorithm to find root polynomial with Aberth method}\r
448 \r
449 \KwIn{$Z^{0}$(Initial root's vector),$\varepsilon$ (error\r
450 tolerance threshold),P(Polynomial to solve)}\r
451 \r
452 \KwOut {Z(The solution root's vector)}\r
453 \r
454 \BlankLine\r
455 \r
456 Initialization of the parameter of the polynomial to solve\;\r
457 Initialization of the solution vector $Z^{0}$\;\r
458 \r
459 \While {$\Delta z_{max}\succ \epsilon$}{\r
460  Let $\Delta z_{max}=0$\;\r
461 $\prec DimGrid,DimBloc\succ kernel\_save(d\_Z^{k-1});$\r
462 $\prec DimGrid,DimBloc \succ kernel\_update(d\_z^{k});$\r
463 \r
464 $\prec DimGrid,DimBloc\succ kernel\_testConverge (d_?z_{max},d_Z^{k},d_Z^{k-1});$\r
465 }\r
466 \end{algorithm}\r
467 ~\\ \r
468 \r
469 \r
470 \r
471 here we need to create two kernel for the step 3 \textit{Kernel\_save} is used to save vector $Z^{K-1}$ and \textit{kernel\_update} is used to update the $Z^{k}$ vector. In phase 4 a kernel is created to test the convergence of the method\r
472 \r
473 \r
474 \subsubsection{the kernel corresponding }\r
475 \subsubsection{Comparison between sequential algorithm and GPU algorithm }\r
476 \bibliographystyle{plain}\r
477 \bibliography{biblio}\r
478 %% \begin{thebibliography}{2}\r
479 \r
480 %% \bibitem [1] {1} O. Aberth, Iteration Methods for Finding\r
481 \r
482 %% all Zeros of a Polynomial Simultaneously, Math. Comput. 27, 122\r
483 %% (1973) 339\96344.\r
484 \r
485 %% \bibitem [2] {2} Ilieff, L. (1948-50), On the approximations of Newton, Annual\r
486 %% Sofia Univ. 46, 167-171.\r
487 \r
488 %% \bibitem [3] {3} Docev, K. (1962), An alternative method of Newton for\r
489 %% simultaneous calculation of all the roots of a given algebraic\r
490 %% equation, Phys. Math. J., Bulg. Acad. Sci. 5, 136-139.\r
491 \r
492 %% \bibitem [4]{4} Durand, E. (1960), Solution Numerique des Equations\r
493 %% Algebriques, Vol. 1, Equations du Type F(x)=0, Racines d'une\r
494 %% Polynome. Masson, Paris.\r
495 \r
496 %% \bibitem [4]  {4} Aberth, O. (1973), Iterative methods for finding all zeros of\r
497 %% a polynomial simultaneously, Math. Comp. 27, 339-344.\r
498 \r
499 %% \bibitem [5] {5} Kerner, I.O. (1966), Ein Gesamtschritteverfahren zur\r
500 %% Berechnung der Nullstellen von Polynomen, Numer. Math. 8, 290-294.\r
501 \r
502 %% \bibitem [6]{6} Borch-Supan, W. (1963), A posteriori error for the zeros of\r
503 %% polynomials, Numer. Math. 5, 380-398.\r
504 \r
505 %% \bibitem [7] {7} Ehrlich, L. W. (1967), A modified Newton method for\r
506 %% polynomials, Comm. Ass. Comput. Mach. 10, 107-108.\r
507 \r
508 \r
509 \r
510 %% \bibitem [10] {10}Loizon, G. (1983), Higher-order iteration functions for\r
511 %% simultaneously approximating polynomial zeros, Intern. J. Computer\r
512 %% Math. 14, 45-58.\r
513 \r
514 %% \bibitem [11]{11} E. Durand, Solutions numŽeriques des Žequations algŽebriques,\r
515 %% Tome 1: Equations du type F(X) = 0; Racines d\92un polyn\88ome,\r
516 %% Masson, Paris 1960.\r
517 \r
518 %% \bibitem [12] {12} Weierstrass, K. (1903), Neuer Beweis des Satzes, dass\r
519 %% jede ganze rationale function einer veranderlichen dagestellt\r
520 %% werden kann als ein product aus linearen functionen derselben\r
521 %% veranderlichen, Ges. Werke 3, 251-269.\r
522 \r
523 \r
524 %%\r
525 \r
526 %% \bibitem [16]{16} Kahina, G. Raphaël, C. Abderrahmane, S. A\r
527 %% parallel implementation of the Durand-Kerner algorithm for\r
528 %% polynomial root-finding on GPU. In INDS 2014, Int. Conf. on\r
529 %% advanced Networking, Distributed Systems and Applications, Bejaia,\r
530 %% Algeria, pages 53--57, June 2014. IEEE\r
531 \r
532 %% \bibitem [17] {17} Karim Rhofir, François Spies, and Jean-Claude Miellou.\r
533 %%Perfectionnements de la méthode asynchrone de Durand-Kerner pour\r
534 %%les polynômes complexes. Calculateurs Parallèles, 10(4):449-- 458,\r
535 %%1998.\r
536 %% \bibitem [18] {18} Bini, D. A. Numerical computation of polynomial zeros by\r
537 %%means of Aberth\92s method. Numerical Algorithms 13 (1996), 179\96\r
538 %%200.\r
539 %% \bibitem [19] {19} A. Ostrowski, On a Theorem by J.L. Walsh Concerning the Moduli of Roots of Algebraic Equations,\r
540 %%Bull. A.M.S., 47 (1941) 742\96746.\r
541 \r
542 %%\bibitem [20] {20} Mirankar WL (1968) Parallel methods for\r
543 %%approximating the roots of a function. IBM Res Dev 297\96 301 30.\r
544 %%\bibitem [21] {21} Mirankar WL (1971) A survey of parallelism in\r
545 %%numerical analysis. SIAM Rev 524\96547\r
546 \r
547 %ù\bibitem [22] {22}Bini DA, Gemignani L (2004) Inverse power and\r
548 %%Durand\96Kerner iterations for univariate polynomial root-finding.\r
549 %%Comput Math Appl 47:447\96459\r
550 \r
551 %%\bibitem [23] {23}Ben-Or M, Feig E, Kozzen D, Tiwary P (1968) A fast parallel\r
552 %algorithm for determining all roots of a polynomial with real\r
553 %%roots. In: Proc of ACM, pp 340\96349\r
554 \r
555 %%\bibitem [24] {24}Zhanc X, Wan M, Yi Z (2008) A constrained learning algorithm for\r
556 %%finding multiple real roots of polynomial. In: Proc of the 2008\r
557 %%intl symposium on computational intelligence and design, pp 38\9641\r
558 \r
559 %%\bibitem [25] {25}Kalantari B (2008) Polynomial root finding and polynomiography.\r
560 %%World Scientific, New Jersey\r
561 \r
562 %%\bibitem [27] {27} Gemignani L (2007) Structured matrix methods for polynomial root\r
563 %%finding. In: Proc of the 2007 Intl symposium on symbolic and\r
564 %%algebraic computation, pp 175\96180 Skachek V, Roth RM (2008)\r
565 \r
566 %%\bibitem [28] {28}Probabilistic algorithm for finding roots of linearized\r
567 %%polynomials. Design, codes and cryptography. Kluwer, Norwell\r
568 \r
569 %%\bibitem [29] {29}Schedler GS (1967) Parallel numerical methods for the solution of\r
570 %%equations. Commun ACM 286\96 290 Ben-Or M, Feig E, Kozzen D, Tiwary\r
571 \r
572 %%\bibitem [30] {30}P (1968) A fast parallel algorithm for determining all roots of a\r
573 %%polynomial with real roots. In: Proc of ACM, pp 340\96349\r
574 \r
575 %%\bibitem [31] {31}Rice TA, Jamieson LH (1989) A highly parallel algorithm for root\r
576 %%extraction. IEEE Trans Comp 38(3):443\96449 20. Jana PK (2006)\r
577 \r
578 %%\bibitem [32] {32}Winogard S (1972) Parallel iteration methods in complexity of\r
579 %%computer communications. Plenum, New York\r
580 \r
581 %ù\bibitem [33] {33} Cosnard M, Fraigniaud P (1990) Finding the roots of a polynomial\r
582 %%on an MIMD multicomputer. Parallel Comput 15:75\9685\r
583 \r
584 %%\bibitem [41] {41} Jana PK (1999) Finding polynomial zeroes on a Multi-mesh of trees\r
585 %%(MMT). In: Proc of the 2nd int conference on information\r
586 %%technology, Bhubaneswar, December 20\9622, pp 202\96206\r
587 \r
588 %%\bibitem [42] {42}Zhu W, Zeng Z, Lin D (2008) An adaptive algorithm finding\r
589 %%multiple roots of polynomials. Lect Notes Comput Sci 5262:674\96681\r
590 \r
591 \r
592 \r
593 %%\bibitem [43] {43}Polynomial interpolation and polynomial root finding on OTIS-Mesh.\r
594 %%Parallel Comput 32:301\96312\r
595 \r
596 %%\bibitem [44] {44}Jana PK, Sinha BP, Datta Gupta R (1999) Efficient parallel\r
597 %%algorithms for finding polynomial zeroes. In: Proc of the 6th int\r
598 %%conference on advance computing, CDAC, Pune University Campus,\r
599 %%India, December 14\9616, pp 189\96196\r
600 \r
601 \r
602 \r
603 %% \end{thebibliography}\r
604 \end{document}\r