mpi

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <mpi.h>

double f(double x){
  return 4.*sqrt(1.-x*x);
}

double calculeaza_integrala(double a, double b, double c, long N){

  long i;
  const gsl_rng_type *T;
  gsl_rng * r;
  double u;
  double x, y;
  long N1= 0, N2=0;
  double I;

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  for(i = 0; i < N; i++){
    u = gsl_rng_uniform(r);
    x = a+ (b-a)*u;
    u = gsl_rng_uniform(r);
    y = c*u;

    if(y < f(x)){
      N1= N1 +1;
    }
    else {
      N2 = N2 + 1;
    }
  }
  I = (double)N1/(double)N * (b-a)*c;
  return I;
}

int main(int argc, char **argv){
  int Np, rank;

  long i;
  double a, b, c;
  double N;
  double c1, c2;

  long Ni;
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &Np);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  if(rank ==0){
    double s, I;
    a = 0.;
    b = 1.;
    c = 4.;
    N = 1000000;

    for(i = 1; i < Np; i++){
      MPI_Send (&a, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
      MPI_Send (&b, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
      MPI_Send (&c, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);

      Ni = (long) ((double)N/(double)(Np-1));
      MPI_Send(&Ni, 1, MPI_LONG, i, 0, MPI_COMM_WORLD);
    }
    I = 0.;

    for(i = 1; i< Np; i++){
      MPI_Recv(&s, 1, MPI_DOUBLE, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      I = I + s;
    }
    I = I / (double) (Np-1);
    printf("I = %f \n", I);

  }
  else {
    double s;
    MPI_Recv(&a, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(&b, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(&c, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    MPI_Recv(&Ni, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

    printf("%f %f %f %i \n", a, b, c, Ni);
    s = calculeaza_integrala(a, b, c, Ni);
    printf("s=%f \n", s);
    MPI_Send(&s, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);

  }
  MPI_Finalize();
  return 0;
}