]> AND Private Git Repository - kahina_paper2.git/blob - code.c
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
d681fa0a63a4adbf805f4cc39b48121992074f84
[kahina_paper2.git] / code.c
1
2
3
4
5 //Normal version of the Ehrlich-Aberth method
6 cuDoubleComplex FirstH_EA(int i, cuDoubleComplex *Z) {
7   
8   cuDoubleComplex  result;
9   cuDoubleComplex  C,F,Fp;
10   int              j;
11   cuDoubleComplex  sum;
12   cuDoubleComplex  un;
13
14   //evaluate the polynomial
15   F = Fonction(Z[i]);
16   //evaluate the derivative of the poly.
17   Fp=FonctionD(Z[i]);  
18   
19   double mod=Cmodule(F);
20   sum.x=0;sum.y=0;
21   un.x=1;un.y=0;
22   C=Cdiv(F,Fp);       //P(z)/P'(z)
23
24   //for all roots, compute the sum
25   //for the Ehrlich-Aberth iteration 
26   for ( j=0 ; j<P.PolyDegre ; j++ )
27   {
28     if ( i != j )
29     {
30       sum=Cadd(sum,Cdiv(un,Csub(Z[i],Z[j])));
31     }
32   }
33   sum=Cdiv(C,Csub(un,Cmul(C,sum)));   //C/(1-Csum) 
34   result=Csub(Z[i], sum);
35   return (result);
36 }
37
38
39 //Log Exp version of the Ehrlich-Aberth method
40 cuDoubleComplex NewH_EA(int i, cuDoubleComplex *Z) {
41
42   cuDoubleComplex result;
43   cuDoubleComplex F,Fp;
44   cuDoubleComplex one,denominator,sum;
45   int             j;
46   one.x=1;one.y=0;
47   sum.x=0;
48   sum.y=0;
49
50   //evaluate the polynomial with
51   //the LogExp version
52   Fp = LogFonctionD(Z[i]);
53   //evaluate the derivative of the polynomial
54   //with the LogExp version
55   F = LogFonction(Z[i]); 
56   
57   cuDoubleComplex FdivFp=Csub(F,Fp);
58
59   //for all roots, compute the sum
60   //for the Ehrlich-Aberth iteration 
61   for ( j=0 ; j<P.degrePolynome ; j++ )
62   {
63     if ( i != j )
64     {
65       sum=Cadd(sum,Cdiv(un,Csub(Z[i],Z[j])));
66     }
67   }
68
69   //then terminate the computation
70   //of the Ehrlich-Aberth method
71   denominator=Cln(Csub(un,Cexp(Cadd(FdivFp,Cln(sum)))));
72   result=Csub(FdivFp,denominator);
73   result=Csub(Z[i],Cexp(res));
74   
75   return   result;
76   
77 }
78
79
80 //kernels to update a root i
81 cuDoubleComplex H_EA(int i, cuDoubleComplex *Z) {
82   cuDoubleComplex c;
83   //if the root needs to be updated
84   if(!finished[i]) {
85     //according to the module of the root
86     if (Cmodule(Z[i])<=maxRadius)
87       //selects the normal version
88       c = FirstH_EA(i,Z);  
89     else
90       //of the Log Exp version
91       c = NewH_EA(i,Z);
92     return c;
93   }
94   else return Z[i];
95 }
96