/* L-25 MCS 572 Monday 13 March 2006 : Newton's method */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "dcmplx.h"
#include "mpi.h"

#define tag 100 /* tag for sending a number */

dcmplx eval ( int n, dcmplx f[n], dcmplx x );
/*
 * DESCRIPTION :
 *   Applies Horner's method to evaluate the polynomial f at x.
 *
 * ON ENTRY :
 *   n        #coefficients of the polynomial equals one plus its degree;
 *   f        complex coefficients of the polynomial,
 *            f[i] is the coefficient with the monomial x^i;
 *   x        where to evaluate the polynomial. 
 *
 * ON RETURN :
 *   y        y = eval() the function value of f at x. */

dcmplx* diff ( int n, dcmplx f[n] );
/*
 * DESCRIPTION :
 *   Allocates memory and computes the coefficients of the derivative of f. */

int newton ( int n, dcmplx f[n], dcmplx *df, dcmplx *x, float eps, int m,
             int verbose );
/*
 * DESCRIPTION :
 *   Performs at most m steps of Newton's method to find a root of f.
 *
 * ON ENTRY :
 *   n       number of coefficients in f, equals one plus the degree of f;
 *   f       complex coefficients of a polynomial in one variable;
 *   df      derivative, given as array of n-1 complex coefficients;
 *   x       initial approximation for the root;
 *   eps     accuracy requirement;
 *   m       maximal number of iterations;
 *   verbose should be 1 if output written to screen, 0 for no output.
 *
 * ON RETURN :
 *   x       improved approximation for the root;
 *   v       v = newton() return value of newton is number of iterations,
 *           if v = m+1, then the desired accuracy was not reached. */

int update_roots ( int n, dcmplx r[n], int *cnt, dcmplx x, float tol,
                   int verbose );
/*
 * DESCRIPTION :
 *   Updates r with a new approximate root x.
 *
 * ON ENTRY :
 *   n       degree of the polynomial, total size of r;
 *   r       roots currently been found;
 *   cnt     number of roots currently in r;
 *   x       new approximate root;
 *   tol     to decide whether two numbers are equal;
 *   verbose should be 1 if intermediate output is wanted, 0 otherwise.
 *
 * ON RETURN :
 *   r       updated sequence of roots;
 *   cnt     updated counter;
 *   v       v = update_roots() is index of x in r. */

double* convert_to_doubles ( int n, dcmplx c[n] );
/*
 * DESCRIPTION :
 *   Returns an array of 2*n doubles with consecutive real and
 *   imaginary parts of the n complex numbers in c. */

dcmplx* convert_to_dcmplx ( int n, double d[n] );
/*
 * DESCRIPTION :
 *   Returns an array of n/2 complex numbers, defined by the
 *   consecutive real and imaginary parts in the array d. */

int main ( int argc, char *argv[] )
{
   int n,i,j,k,max,nit,cnt,p,myid,stages;
   dcmplx *r,*f,*df,x,y;
   double *fd,*dfd;
   float tol;
   char ans;
   MPI_Status status;

   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD,&p);
   MPI_Comm_rank(MPI_COMM_WORLD,&myid);

   if(myid==0)
   {
      srand(time(NULL));
      printf("Give the degree of the polynomial : "); scanf("%d",&n);
      r = (dcmplx*)calloc(n,sizeof(dcmplx));    /* space for roots */
      f = (dcmplx*)calloc(++n,sizeof(dcmplx));  /* coefficient vector */
      for(i=0; i<n; i++) f[i] = random_dcmplx1();
      printf("The random coefficients : \n");
      for(i=0; i<n; i++,printf("\n")) write_dcmplx(f[i]);
      df = diff(n,f);
      printf("Give the maximal number of steps : "); scanf("%d",&max);
      printf("Give the required accuracy : "); scanf("%f",&tol);
      fd = convert_to_doubles(n,f);
      dfd = convert_to_doubles(n-1,df);
      printf("Give the number of stages : "); scanf("%d",&stages);
   }
   MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
   if(myid!=0)
   {
      fd = (double*)calloc(2*n,sizeof(double));
      dfd = (double*)calloc(2*(n-1),sizeof(double));
   }
   MPI_Bcast(fd,2*n,MPI_DOUBLE,0,MPI_COMM_WORLD);
   MPI_Bcast(dfd,2*(n-1),MPI_DOUBLE,0,MPI_COMM_WORLD);
   if(myid!=0)
   {
      f = convert_to_dcmplx(2*n,fd);
      df = convert_to_dcmplx(2*(n-1),dfd);
   }
   MPI_Bcast(&max,1,MPI_INT,0,MPI_COMM_WORLD);
   MPI_Bcast(&tol,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
   MPI_Bcast(&stages,1,MPI_INT,0,MPI_COMM_WORLD);
   for(k=0; k<stages; k++)
   {
      if(myid==0)
      {
         cnt = 0;
         for(i=1; i<p; i++)
         {
            x = random_dcmplx1();
            printf("%d: Sending %.3le %.3le to node %d.\n",k,x.re,x.im,i);
            fflush(stdout);
            MPI_Send(&x.re,1,MPI_DOUBLE,i,tag,MPI_COMM_WORLD);
            MPI_Send(&x.im,1,MPI_DOUBLE,i,tag,MPI_COMM_WORLD);
            MPI_Recv(&nit,1,MPI_INT,i,tag,MPI_COMM_WORLD,&status);
            MPI_Recv(&x.re,1,MPI_DOUBLE,i,tag,MPI_COMM_WORLD,&status);
            MPI_Recv(&x.im,1,MPI_DOUBLE,i,tag,MPI_COMM_WORLD,&status);
            printf("%d: Received %.3le %.3le and %d steps from %d.\n",
                   k,x.re,x.im,nit,i); fflush(stdout);
            if(nit<max) j = update_roots(n-1,r,&cnt,x,tol,1);
            else printf("%d: Rejected root coming from %d.\n",k,i);
            fflush(stdout);
         }
      }
      else
      {
         MPI_Recv(&x.re,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
         MPI_Recv(&x.im,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD,&status);
         nit = newton(n,f,df,&x,tol,max,0);
         printf("%d: Node %d finds %.3le %.3le after %d steps.\n",
                k,myid,x.re,x.im,nit); fflush(stdout);
         MPI_Send(&nit,1,MPI_INT,0,tag,MPI_COMM_WORLD);
         MPI_Send(&x.re,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
         MPI_Send(&x.im,1,MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
      }
   }
   MPI_Finalize();
   return 0;
}

dcmplx eval ( int n, dcmplx f[n], dcmplx x )
{
   dcmplx y = f[n-1];
   int i;

   for(i=n-2; i>=0; i--)
   {
      y = mul_dcmplx(y,x);         /* y = y*x */
      y = add_dcmplx(y,f[i]);      /* y = y + f[i] */
   }
   return y;
}

dcmplx* diff ( int n, dcmplx f[n] )
{
   dcmplx *y = (dcmplx*)calloc(n-1,sizeof(dcmplx));
   int i;

   for(i=1; i<n; i++)
      y[i-1] = mul_double(f[i],((double) i));

   return y;
}

int newton ( int n, dcmplx f[n], dcmplx *df, dcmplx *x, float eps, int m,
             int verbose )
{
   int i;
   dcmplx y,dy,dx;

   y = eval(n,f,*x);            /* y = f(x) */

   for(i=1; i<=m; i++)
   {
      dy = eval(n-1,df,*x);     /* dy = f'(x) */
      dx = div_dcmplx(y,dy);    /* dx = f(x)/f'(x) */
      *x = sub_dcmplx(*x,dx);   /* x = x - f(x)/f'(x) */
      y = eval(n,f,*x);         /* y = f(x), for stop criterium */
      if(verbose > 0)
      {
         printf("************ results at step %d ******************\n",i);
         printf("    Newton update : "); write_dcmplx(dx); printf("\n");
         printf("new approximation : "); write_dcmplx(*x); printf("\n");
         printf("   function value : "); write_dcmplx(y);  printf("\n");
      }
      if((modulus(y)<=eps) || (modulus(dx) <= eps)) return i;
   }

   return i;
}

int update_roots ( int n, dcmplx r[n], int *cnt, dcmplx x, float tol,
                   int verbose )
{
   int i;

   for(i=0; i<*cnt; i++)
      if(equal_dcmplx(r[i],x,tol) == 1)
      {
         if(verbose > 0)
         {
            printf("Found "); write_dcmplx(r[i]); printf(" at %d.\n",i);
         }
         return i;
      }
   r[(*cnt)++] = x;
   if(verbose > 0)
   {
      printf("Added "); write_dcmplx(x); printf(" as root %d.\n",*cnt);
   }
   return (*cnt-1);
}

double* convert_to_doubles ( int n, dcmplx c[n] )
{
   double *d = (double*)calloc(2*n,sizeof(double));
   int i,j;

   for(i=0,j=0; i<n; i++)
   {
      d[j++] = c[i].re;
      d[j++] = c[i].im;
   }
   return d;
}

dcmplx* convert_to_dcmplx ( int n, double d[n] )
{
   int n2 = n/2;
   dcmplx *c = (dcmplx*)calloc(n2,sizeof(dcmplx));
   int i,j;

   for(i=0,j=0; i<n2; i++)
   {
      c[i].re = d[j++];
      c[i].im = d[j++];
   }
   return c;
}