Buffonovo vprašanje

Klasičen primer Monte Carlo (MC) metode je metanje kemičnega svinčnika na tla s ploščicami. Želimo ugotoviti kolikšna je verjetnost, da se kemik ustavi na fugi ko ga vržemo na tla. Problem je nekoliko zahtevnejši od računanja razmerja med površino kroga in kvadrata, saj robni pogoji niso tako očitni. Za začetek poenostavimo problem in predpostavimo, da nas zanima samo prekrivanje vzdolžnih črt, ki so z razmaknjene za razdaljo t. Zanima nas konkreten primer ko je dolžina kemika polovico širine med fugami (l = t/2). Verjetnost je, da pade na fugo je očitno manjša kot 1. Koliko torej je? Naredimo numerični eksperiment z mentanjem tako da nas na koncu zanima razmerje med vsemi vrženimi (N) in številom, ko je bila fuga prekrita. Kolikšno je to razmerje.

Napotki:

1. Uporabite uniformni generator nakljčnih števil kot je npr drand48() in njegove variante za izdelavo serijskega programa. Pri izdelavi serijskega programa se osredotočite na pravilnost delovanja in ne na hitrost in s tem optimiranje kode do nerazumljivosti. Če želite si lahko za l/t izberete enote, ki so vam najbolj razumljive.
2. Ker v praksi ne moremo preveriti kakšno je končno razmerje uporabite dolga števila za štetje. Predlagamo N=300000000
3. Izmerite čas serijskega programa na enem procesorju poljubnega vozlišča z ukazom single ./buffon
4. Program razširite z OpenMP ukazi ter uporabo nitno varne funkcije drand48_r()
5. Na vozlišču prevedite in poženite z različnim številom niti (od 1 do 24) z ukazi

   gcc -O2 -fopenmp -o buffonmp buffonmp.c -lm
   
   node time OMP_NUM_THREADS=24 ./buffonmp
   

6. Za sestavljanje posameznih seštevkov zadetkov uporabite redukcijo

7. Predelajte taisti program še za MPI, ki ga poganjate na več procesorjih s klasičnim načinom

   bsub -n 4,24 -oo buffonmpi.log mpirun buffonmpi # ali z
   
   node time mpirun -np 4 buffonmpi # za ročno merjenje časa
   

8. Naredite si tabelo s stenčasi.
9. Poskusite s hibridnim programom OpenMPI/MP, kjer prevajamo z

   mpicc -fopenmp -o buffonhy buffonhy.c -lm
   
   bsub -n 24 -R "span[ptile=12]" -oo buffonhy.log mpirun  buffonhy 
   

10. Koliko metov/ našega časa rabimo da dobimo še malo boljši rezultat?
11. Ko ste vse to naredili, Vam je šef rekel, da ga ne zanima poenostavljen primer ampak tisti, ki ima še horizontalne fuge.

Serijski program

#include 

#include 

#include 

#define N 300000000

#define t 1.0

#define l 0.5



int main(int argc, char *argv[])

{

  long i, count = 0;

  double x, phi;

  for(i = 0; i < N; ++i)

    {

      x = t*drand48() - t/2;

      phi = 2*M_PI*drand48();

      double x1 = x-l/2*cos(phi);

      double x2 = x+l/2*cos(phi);

      if (x1 < 0 && x2 > 0 || x1 > 0 && x2 < 0)

        ++count;

    }

  printf("Razmerje = %lf\n", (double)N/count);

}

OpenMP

Zaradi deljenega spomina je pri OpenMP potrebno paziti, da se uporabljajo rutine iz knjižnic, ki so thread safe saj ne smejo vsebovati statični spremeljivk stanja.

#include 

#include 

#include 

#include 



#define N 300000000

#define t 1.0

#define l 0.5



int main(int argc, char *argv[])

{

  long i, count = 0;

  struct drand48_data buffer;



#pragma omp parallel private(buffer)

  {

    double x, phi;

    srand48_r(omp_get_thread_num(), &buffer);

#pragma omp for private (x, phi) reduction(+:count)

    for(i = 0; i < N; ++i)

      {

        drand48_r(&buffer, &x);

        x -= 0.5;

        drand48_r(&buffer, &phi);

        phi *= 2*M_PI;

        double x1 = x-l/2*cos(phi);

        double x2 = x+l/2*cos(phi);

        if (x1 < 0 && x2 > 0 || x1 > 0 && x2 < 0)

          ++count;

      }

  }

  printf("Razmerje = %f\n", (float)N/count);

}

MPI

#include 

#include 

#include 

#include 



#define N 300000000

#define t 1.0

#define l 0.5



int main(int argc, char *argv[])

{

  long i, count = 0;

  struct drand48_data buffer;

  int id, procesov;



  MPI_Init(&argc, &argv);

  MPI_Comm_rank(MPI_COMM_WORLD, &id);

  MPI_Comm_size(MPI_COMM_WORLD, &procesov);



  double x, phi;

  srand48_r(id, &buffer);

  for(i = id; i < N; i += procesov)

    {

      drand48_r(&buffer, &x);

      x -= 0.5;

      drand48_r(&buffer, &phi);

      phi *= 2*M_PI;

      double x1 = x-l/2*cos(phi);

      double x2 = x+l/2*cos(phi);

      if (x1 < 0 && x2 > 0 || x1 > 0 && x2 < 0)

        ++count;

    }

  long total=0;

  MPI_Reduce(&count, &total, 1, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD);

  if (id == 0)

    printf("Razmerje = %f\n", (float)N/total);

  MPI_Finalize();

}