code.c

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