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

/* Paralelizovane generovani histogramu normalizovaneho rovnomerneho
   rozdeleni. Normalizovane rovnomerne rozdeleni ma konstantni hustotu
   p-ti f(x) = 1.0 nad intervalem <0.0, 1.0>. Kazdy proces udela svuj
   (dilci) histogram, proces root (ID=0) pak posklada dilci vysledky
   do celkoveho histogramu a vytiskne jednak shromazdena data
   (pole a_all) a vlastni histogram. Histogram musi vyjit tak,
   ze relativni cetnosti jsou priblizne stejne a dohromady davaji
   hodnotu 1.0 */

#define K 10
/* pocet (sub)intervalu histogramu nad intervalem <0,1> */

#define N 10000
/* pocet nahodnych cisel generovanych jednim procesem */

int my_id;     /* MPI cislo procesu v rozmezi 0 az (mpi_N - 1) */
int mpi_N;     /* MPI pocet procesu N */

int main(int argc, char *argv[]) {
  long int *a_loc = NULL;    /* pole citacu histogramu - ma kazdy proces */
  long int *a_all = NULL;    /* pole vsech histogramu - ma jen root */
  long int *sum = NULL;      /* soucet vsech histogramu - ma jen root */
  int i,j, jj, i_where;
  double x;

  MPI_Init (&argc, &argv);
  MPI_Comm_rank (MPI_COMM_WORLD, &my_id);
  MPI_Comm_size (MPI_COMM_WORLD, &mpi_N);

  /* vyrobeni dynamickych poli */
  a_loc = (long int*) malloc(K*sizeof(long int));  /* vsichni */
  if (my_id == 0){    /* jen root */
    a_all = (long int*) malloc(mpi_N * K * sizeof(long int));
    sum = (long int*) malloc(K*sizeof(long int));
  }

  for (j=0; j < K; j++) a_loc[j]=0;  /* nulovani citacu histogramu */

  srand((int)(my_id*371)); /* !!!!rozdilne pocatecni nastaveni generatoru
                          nahodnych cisel pro jednotlive procesy (tj.
                          kazdy proces musi delat jine pokusy), jinak by
                          paralelizace nemela smysl (plati pro metody
                          Monte Carlo obecne) */

  for (i=0; i < N; i++) { /* vypocet histogramu - delaji vsichni */
    x = ((double)rand())/RAND_MAX;
    i_where = ((int)(x*K))%K;
    a_loc[i_where]++;
  }

  /* vyvolani gather() - delaji vsichni */
  MPI_Gather (a_loc, K, MPI_LONG, a_all, K, MPI_LONG,
    0, MPI_COMM_WORLD);

  /* vypocet a tisk celkoveho vysledku - dela jen root */
  if (my_id==0) {
    printf("\n mpi_N = %d\n", mpi_N);
    for (j=0; j<K*mpi_N; j++)
      printf("\n a_all[%d] = %d", j, a_all[j]);
    for (j=0; j<K; j++){
      sum[j] = 0;
      for (jj=0; jj<mpi_N; jj++) sum[j]+= a_all [j + (jj*K)];
      printf("\n rel_cnt[%d] = %2.18f", j, ((double)(sum[j])/(N*mpi_N)));
    }
    printf("\n");
    free (sum);   /* sum a a_all uvolnuje jen root */
    free(a_all);
  }

  free (a_loc);   /* a_loc uvolnuji vsichni */



  /* odhlaseni z MPI */
  MPI_Finalize();
}