#include<stdio.h>
#include<mpi.h>

#define nmax 100
#define pmax 10

void vec_gen(a, b, dim)
  int dim; 
  double a[nmax],b[nmax];
/*
Generira vectori a i b s razmernost dim 
v procesora, koito izvikva funkciata
*/
{ int i;

  for(i=0; i<dim; i++)
  {  a[i]=2. ;
     b[i]=i+1. ;
  }
}

double dot_pr(a,b,dim,np,myid)
  int dim, np, myid;
  double a[nmax],b[nmax];
/*
Presmyata skalarnoto proizvedenie na a i b s np procesora.
Vseki ot procesorite poluchava chast ot a i ot b (generirani 
sa samo v proces 0) i presmyata saotvetnata suma, sled koeto 
ya izprashta otnovo na proces 0.  
Krainiat rezultat se izprashta vav vsichki procesori.
*/
{ int i, m, mod, rem, myid1, in1, atype, btype;
  
  double res=0.;
  double prod=0.;
  double a1[nmax], b1[nmax], vprod[pmax];
  MPI_Comm com=MPI_COMM_WORLD;
  MPI_Status status;

  atype=1;
  btype=2;

  mod=dim/np;
  rem=dim%np;

  if (myid<rem) 
  {  m=mod+1;  }
  else 
  {  m=mod;  } 

  if (rem==0)
  {  MPI_Scatter(a, m, MPI_DOUBLE,a1,m,MPI_DOUBLE ,0, com);
     MPI_Scatter(b, m, MPI_DOUBLE,b1,m,MPI_DOUBLE ,0, com);
  }
  else  
  {  if (myid==0)
     {  for(i=0;i<m; i++)
        {  a1[i]=a[i];
           b1[i]=b[i];
        }
        
        for(i=1; i<=rem-1; i++)
        {  in1=i*(mod+1);

           MPI_Send(&(a[in1]), mod+1, MPI_DOUBLE, i,atype,com);
           MPI_Send(&(b[in1]), mod+1, MPI_DOUBLE, i,btype,com);

        }
        for(i=rem; i<np; i++)
        {  in1=rem*(mod+1)+(i-rem)*mod;

           MPI_Send(&(a[in1]), mod, MPI_DOUBLE, i,atype,com);

           MPI_Send(&(b[in1]), mod, MPI_DOUBLE, i,btype,com);
        }
     }
     else
     {
        MPI_Recv(a1, m, MPI_DOUBLE, 0, atype, com, &status);
        MPI_Recv(b1, m, MPI_DOUBLE, 0, btype, com, &status);

     }
  }
 

  for(i=0; i<m; i++)
  {  prod=prod+a1[i]*b1[i];
  } 
  MPI_Gather(&prod,1,MPI_DOUBLE,vprod,1,MPI_DOUBLE,0, com);

  if (myid==0) 
  {  for(i=0; i<np; i++)
     {  res=res+vprod[i];
     }
  }  
  MPI_Bcast(&res, 1, MPI_DOUBLE, 0, com);

  return  res;
}  

int main(argc, argv)
  int argc;
  char **argv;
{ int myid, np, dim, i;
  double a[nmax], b[nmax], c;
  MPI_Comm com=MPI_COMM_WORLD;
  MPI_Status status;

  MPI_Init(&argc, &argv);

  MPI_Comm_size(com, &np);
  MPI_Comm_rank(com, &myid);

  if(myid==0)
  {  printf("Enter the dimension of the vectors: \n");
     scanf("%d",&dim);
     vec_gen(a,b,dim);
  }
  MPI_Bcast(&dim, 1, MPI_INT, 0, com);

  printf("I am %d of %d \n", myid, np);


  c=dot_pr(a,b,dim,np,myid);

  if (myid==0) 
  {  printf("Dot product of a and b is: %f \n", c );
  }

  MPI_Finalize();
  return 0;
}