#include "stdio.h"
#include "math.h"
#include "mpi.h"

float func(float x[4]) 
{   return pow(x[0],4)+pow(x[1],2)+pow(x[2],2)+pow(x[3],4);
}

float vnorm(float x[4]) 
{   return  pow((x[0]*x[0]+x[1]*x[1]+x[2]*x[2]+x[3]*x[3]),0.5);
}

int main(int argc, char *argv[]) 
{/* mpi control declarations */
   int numprocs, numslaves, master, slave, myrank, tag;
   MPI_Status status;  
   float done, donetest, donetop;
/* other declarations */   
   int i, j, count, maxcount;
   float xa[5][4],x[5],xs[5],ya[5],y0,y,dy[4],ndy,udy[4]; 
   float delta,eps,h,m;   
/* mpi data initalizations */
   master = 0;
   tag = 1; 
   done = 0.0;
   donetest = 1.0;
   donetop = 2.0;
/* master node parameter and data initalizations */
   delta = 0.01;
   h = 1;
   eps = pow(10,-6);          
   maxcount = 1000;
   for (j=0; j <= 3; j=j+1)
      {xa[0][j] = 1;
       x[j] = xa[0][j];
      }         
   y0 = 10000;
   y  = 0.0;
   count = 0;   
/* mpi driver */
   MPI_Init(&argc,&argv);
   MPI_Comm_size( MPI_COMM_WORLD, &numprocs );
   MPI_Comm_rank( MPI_COMM_WORLD, &myrank );
/* check if number of processors is right */
   if (!(numprocs == 6))
   {  if (myrank == master)
        printf("Error: the number of processors should be 6, but was allocated %d \n",numprocs);
      MPI_Finalize();  
      return 0;
   }   
   if (myrank == master)      
     /* master node code */  
   { while (done < donetest)
     {/* set seek directions */
       for (i=1; i<=4;i=i+1)
       {  for (j=0; j<=3;j=j+1) 
             xa[i][j]= xa[0][j];
          xa[i][i-1] = xa[i][i-1]+delta;
       }
      /* compute values in seek directions and then gradient */        
       for (i=0; i<=4;i=i+1)
       {  for (j=0; j<=3;j=j+1)
             xs[j] = xa[i][j];
          xs[4] = done;
          MPI_Send(&xs,5,MPI_FLOAT,i+1,tag,MPI_COMM_WORLD);  
       }       
         
       for (i=0; i<=4;i=i+1)
           MPI_Recv(&ya[i],1,MPI_FLOAT,i+1,tag,MPI_COMM_WORLD,&status);       
       for (i=0; i<=3;i=i+1)
           dy[i] = (ya[i+1]-ya[0])/delta;
       ndy = vnorm(dy);
       printf("count = %d, ndy = %8.3f, \n",count,ndy);
       for (j=0; j<=3;j=j+1)
           udy[j] = dy[j]/ndy;   
       /* update seek point */
       if (ndy > 1)
          m = 1;
       else
          m = ndy;
       for (j=0; j<=3;j=j+1)
       {   x[j] = xa[0][j];
           xa[0][j] = xa[0][j]-m*h*udy[j];
       }    
       y0 = y;
       y = ya[0];
       count = count+1;
      /* check if done */ 
       if ((eps > fabs(y-y0))||(count == maxcount)) 
       { done = donetop;
         for (i=0; i<=4;i=i+1)   
         {  xs[4] = done;    
            MPI_Send(&xs,5,MPI_FLOAT,i+1,tag,MPI_COMM_WORLD);
         }  
       }  
     }
      printf("x-values %11.8f, %11.8f, %11.8f, %11.8f, \n",x[0],x[1],x[2],x[3]);
      printf("y-values %11.8f, \n",y);   
   }  
   else
    /* slave node code */  
    while (done < donetest)
    {  MPI_Recv(&xs,5,MPI_FLOAT,master,tag,MPI_COMM_WORLD,&status);
       if (xs[4] < donetest) 
       {   y = func(xs);
           MPI_Send(&y,1,MPI_FLOAT,master,tag,MPI_COMM_WORLD);
       }         
       else 
           done = donetop;               
    }  
/* mpi cleanup */
   MPI_Finalize();  
   return 0;
}