code.c

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