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