Python modul mpi4py: kolektivna komunikacija

Ova značajka je dio standarda MPI-1.

• operacije broadcast, scatter, gather, allgather, alltoall -- jednostavno objašnjenje dano je na slajdovima 30 i 31 prezentacije MPI for Python od Lisandra Dalcina

• drastično pojednostavljuje način rada u situaciji više procesa izvodi vrlo sličan ili identičan niz naredbi

  • komunikacija točka-do-točke: u pravilu se slanja i primanja kodiraju posebno za svaki proces
  • kolektivna komunikacija: komunikacijske operacije pišu se jednom za sve procese koji sudjeluju u komunikaciji

Operacija barrier

• comm.Barrier() -- sinkronizacija korištenjem barijere; čeka se da svi procesi dođu do barijere

import numpy as np
# import random
import time
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

sleep_time = 10 * np.random.random()
# sleep_time = 10 * random.random()
print("Proces ranga", rank, "spavat će", sleep_time, "sekundi")
time.sleep(sleep_time)

comm.Barrier()

print("Proces ranga", rank, "prošao je barijeru")

Operacija broadcast

• comm.bcast(obj, root) vrijednost varijable sa korijenskog procesa šalje svima

  • obj je Python objekt koje si šalje
  • root je rang korijenskog procesa; podatak s njega bit će poslan svim ostalima i prepisati preko podataka na ostalim procesima

# varijanta s Python objektima
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

if rank == 0:
    sendmsg = (7, "abc", [1.0, 2+3j], {3: 4})
else:
    sendmsg = None

recvmsg = comm.bcast(sendmsg, root=0)

• comm.Bcast(array, root) -- vrijednost varijable sa korijenskog procesa šalje svima

  • array je numpy polje koje se šalje; mora biti istog oblika na svim procesima
  • root je rang korijenskog procesa; podatak s njega bit će poslan svim ostalima i prepisati preko podataka na ostalim procesima

# varijanta s NumPy poljima
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

if rank == 0:
    msg = np.array([1.5, 2.5, 3.5, 4.5])
else:
    msg = np.empty(4)

print("Prije broadcasta proces ranga", rank, "ima vrijednost varijable", msg)
comm.Bcast(msg, root=0)
print("Nakon broadcasta proces ranga", rank, "ima vrijednost varijable", msg)

Operacija scatter

• comm.scatter(sendmsg, root) listu od n elemenata raspršuje na n procesa, tako da svaki proces dobiva po jedan element i to točno onaj s indeksom koliki je njegov rang

  • sendmsg je Python lista koja se raspršuje; mora biti iste veličine kao broj procesa na korijenskom procesu, i None na ostalima
  • root je rang korijenskog procesa

# varijanta s Python objektima
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

if rank == 0:
    sendmsg = [i**2 for i in range(size)]
else:
    sendmsg = None

recvmsg = comm.scatter(sendmsg, root=0)

• comm.Scatter(sendmsg, recvmsg, root) -- listu od n elemenata raspršuje na n procesa, tako da svaki proces dobiva po jedan element i to točno onaj s indeksom koliki je njegov rang

  • sendmsg je numpy polje koje se raspršuje; mora biti iste veličine kao broj procesa na korijenskom procesu, i None na ostalima
  • recvmsg je numpy polje veličine jednog elementa u koje se sprema rezultat raspršenja
  • root je rang korijenskog procesa

# varijanta s NumPy poljima
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

if rank == 0:
    sendmsg = np.arange(size, dtype=np.int32)
else:
    sendmsg = None

recvmsg = np.empty(1, dtype=np.int32)
comm.Scatter(sendmsg, recvmsg, root=0)
print("Proces ranga", rank, "ima vrijednost poruke za slanje", sendmsg, "primljena poruka", recvmsg)

Operacije gather i allgather

• comm.gather(sendmsg, root) radi obrnuto od scatter(), po jednu varijablu sa svakog procesa skuplja u listu na korijenskom procesu i to tako da vrijednost postaje element s indeksom koliki je rang pripadnog procesa

  • sendmsg je Python objekt koji se šalje korijenskom procesu
  • root je rang korijenskog procesa

• comm.allgather(sendmsg) vrši istu operaciju kao gather(), osim što sada rezultirajuću listu dobivaju svi procesi, a ne samo korijenski

  • sendmsg je Python objekt koji se šalje

# varijanta s Python objektima
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

sendmsg = rank**2

recvmsg1 = comm.gather(sendmsg, root=0)

recvmsg2 = comm.allgather(sendmsg)

• comm.Gather(sendmsg, recvmsg, root) -- vrši obrnut proces od Scatter(), po jednu varijablu sa svakog procesa skuplja u polje na korijenskom procesu i to tako da se vrijednost postavlja na indeks koliki je rang pripadnog procesa

  • sendmsg je numpy polje koje se šalje korijenskom procesu
  • recvmsg je numpy polje u koje se sakupljaju primljene vrijednosti; ono mora biti veličine koliki je broj procesa na korijenskom procesu, i None na ostalima
  • root je rang korijenskog procesa

• comm.Allgather(sendmsg, recvmsg) -- vrši operaciju sakupljanja kao i Gather(), osim što sada rezultirajuću listu dobivaju svi procesi, a ne samo korijenski

  • sendmsg je numpy polje koje se šalje
  • recvmsg je numpy polje u koje se sakupljaju primljene vrijednosti

# varijanta s NumPy poljima
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

sendmsg = np.array([rank**2], dtype=np.int32)

print("Proces ranga", rank, "ima vrijednost poruke za slanje", sendmsg)

if rank == 0:
    recvmsg1 = np.empty(size, dtype=np.int32)
else:
    recvmsg1 = None

comm.Gather(sendmsg, recvmsg1, root=0)

recvmsg2 = np.empty(size, dtype=np.int32)

comm.Allgather(sendmsg, recvmsg2)

print("Proces ranga", rank, "ima vrijednost prve primljene poruke", recvmsg1, "i druge primljene poruke", recvmsg2)

Operacije reduce i allreduce

• comm.reduce(sendmsg, op, root) vrši proces redukcije (suma, produkt, maksimum, minimum) na danom skupu elemenata koji čini po jedna vrijednost sa svakog od procesa i rezultat redukcije sprema u varijablu na korijenskom procesu

  • sendmsg je Python objekt na kojem se vrši redukcija
  • op je operator redukcije: MAX, MIN, SUM, PROD, ... (čitav popis predefiniranih operatora moguće je pronaći u standardu MPI 1.1)

• comm.allreduce(sendmsg, op) vrši istu operaciju kao reduce(), osim što sada rezultat dobivaju svi procesi, a ne samo korijenski

  • sendmsg je Python objekt na temelju kojeg se vrši redukcija
  • op je operator redukcije

# varijanta s Python objektima
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

sendmsg = rank

recvmsg1 = comm.reduce(sendmsg, op=MPI.SUM, root=0)

recvmsg2 = comm.allreduce(sendmsg)

• comm.Reduce(sendmsg, recvmsg, op, root) -- vrši proces redukcije (suma, produkt, maksimum, minimum) na danom skupu elemenata koji čini po jedna vrijednost sa svakog od procesa i rezultat redukcije sprema u varijablu na korijenskom procesu

  • sendmsg je numpy polje na kojem se vrši redukcija
  • recvmsg je numpy polje u koje se sprema rezultat; ono mora biti odgovarajuće veličine na korijenskom procesu, i None na ostalima
  • op je operator redukcije: MAX, MIN, SUM, PROD, ... (čitav popis predefiniranih operatora moguće je pronaći u standardu MPI 1.1)
  • root je rang korijenskog procesa

• comm.Allreduce(sendmsg, recvmsg, op) -- vrši istu operaciju kao Reduce(), osim što sada rezultat dobivaju svi procesi, a ne samo korijenski

  • sendmsg je numpy polje na temelju kojeg se vrši redukcija
  • recvmsg je numpy polje u koje se sprema rezultat
  • op je operator redukcije

# varijanta s NumPy poljima
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

sendmsg = np.array([rank**2], dtype=np.int32)

if rank == 0:
    recvmsg1 = np.empty(1, dtype=np.int32)
else:
    recvmsg1 = None

comm.Reduce(sendmsg, recvmsg1, op=MPI.SUM, root=0)

recvmsg2 = np.empty(1, dtype=np.int32)

comm.Allreduce(sendmsg, recvmsg2, op=MPI.SUM)

print("Proces ranga", rank, "ima vrijednost prve primljene poruke", recvmsg1, "i druge primljene poruke", recvmsg2)

Operacija alltoall

• comm.alltoall(sendmsg) -- element vektora sa indeksom j iz sendmsg na procesu i postaje element vektora sa indeksom i u recvmsg na procesu j

  • sendmsg je lista koja se šalje, veličine jednake ukupnom broju procesa

# varijanta s Python objektima
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

sendmsg = [rank * size + i for i in range(size)]

recvmsg = comm.alltoall(sendmsg)

print("Proces ranga", rank, "ima vrijednost poslane poruke", sendmsg, "vrijednost primljene poruke", recvmsg)

• comm.Alltoall(sendmsg, recvmsg) -- element vektora sa indeksom j iz sendmsg na procesu i postaje element vektora sa indeksom i u recvmsg na procesu j

  • sendmsg je numpy polje koje se šalje, veličine jednake ukupnom broju procesa
  • recvmsg je numpy polje koje se prima, veličine jednake ukupnom broju procesa

# varijanta s NumPy poljima
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

sendmsg = np.array([rank * size + i for i in range(size)], dtype=np.int32)

recvmsg = np.empty(size, dtype=np.int32)

comm.Alltoall(sendmsg, recvmsg)

print("Proces ranga", rank, "ima vrijednost poslane poruke", sendmsg, "vrijednost primljene poruke", recvmsg)

Primjeri korištenja kolektivne komunikacije

Zbroj vektora

• sekvencijalni kod je vrlo jednostavan smo, kao što smo već vidjeli

• korištenje kolektivne komunikacije tipa scatter-gather značajno pojednostavljuje kod u odnosu na komunikaciju tipa točka-do-točke; pretpostavka je da se kod izvodi u 4 procesa, inače scatter()/Scatter() i gather()/Gather() javljaju grešku

  # varijanta s Python objektima
  from mpi4py import MPI
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = [1, 2, 3, 9]
      b = [4, 5, 6, 7]
  else:
      a = None
      b = None
  
  komponenta_od_a = comm.scatter(a, root=0)
  komponenta_od_b = comm.scatter(b, root=0)
  zbroj = komponenta_od_a + komponenta_od_b
  print("Proces ranga", rank, "izračunao je zbroj", zbroj)
  
  zbroj_vektor = comm.gather(zbroj, root=0)
  if rank==0: # ostali procesi će ispisati None, slično kao kod vektora za scatter()
      print("Zbroj vektora a, b i c iznosi", zbroj_vektor)
  

  # varijanta s NumPy poljima
  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = np.array([1, 2, 3, 9], dtype=np.int32)
      b = np.array([4, 5, 6, 7], dtype=np.int32)
      zbroj_vektor = np.empty(4, dtype=np.int32)
  else:
      a = None
      b = None
      zbroj_vektor = None
  
  element_a = np.empty(1, dtype=np.int32)
  comm.Scatter(a, element_a, root=0)
  element_b = np.empty(1, dtype=np.int32)
  comm.Scatter(b, element_b, root=0)
  
  element_zbroj = element_a + element_b
  print("Proces ranga", rank, "izračunao je zbroj", element_zbroj)
  
  comm.Gather(element_zbroj, zbroj_vektor, root=0)
  if rank == 0:
      print("Zbroj vektora a i b iznosi", zbroj_vektor)
  

# varijanta s Python objektima
a = [[1, 2, 3],
     [4, 5, 6],
     [7, 8, 9]]
b = [[10, 20, 30],
     [40, 50, 60],
     [70, 80, 90]]
zbroj = [[0, 0, 0],
         [0, 0, 0],
         [0, 0, 0]]

for i in range(3):
    for j in range(3):
        zbroj[i][j] = a[i][j] + b[i][j]

# varijanta s NumPy poljima
a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.int32)
b = np.array([[10, 20, 30], [40, 50, 60], [70, 80, 90]], dtype=int32)
zbroj = np.zeros((3, 3), dtype=np.int32)

for i in range(3):
    for j in range(3):
        zbroj[i][j] = a[i][j] + b[i][j]

# zbroj = a + b

Računanje skalarnog produkta vektora

Skalarni produkt vektora a i b veličine n elemenata je izraz oblika

skalarni_produkt = a[0] * b[0] + a[1] * b[1] + ... a[n] * b[n]

• serijski kod za skalarni produkt je oblika

  # varijanta s Python objektima
  a = [1, 2, 3, 9, 10, 11]
  b = [4, 5, 6, 7, 8, 24]
  
  skalarni_produkt = 0
  for i in range(6):
      skalarni_produkt += a[i] + b[i]
  
  # varijanta s NumPy poljima
  a = np.array([1, 2, 3, 9, 10, 11], dtype=np.float32)
  b = np.array([4, 5, 6, 7, 8, 24], dtype=np.float32)
  
  skalarni_produkt = 0.0
  for i in range(6):
      skalarni_produkt += a[i] * b[i]
  
  # skalarni_produkt = np.dot(a, b)
  

• paralelni kod za koji koristi kolektivnu komunikaciju tipa scatter-gather je oblika

  # varijanta s Python objektima
  from mpi4py import MPI
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = [1, 2, 3, 9, 10, 11]
      b = [4, 5, 6, 7, 8, 24]
  else:
      a = None
      b = None
  
  komponenta_od_a = comm.scatter(a, root=0)
  komponenta_od_b = comm.scatter(b, root=0)
  produkt = komponenta_od_a * komponenta_od_b
  
  produkt_vektora = comm.gather(produkt, root=0)
  if rank==0:
      suma = 0
      for x in produkt_vektora:
          suma += x
      print("Skalarni produkt iznosi", suma)
  

  # varijanta s NumPy poljima
  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = np.array([1, 2, 3, 9, 10, 11], dtype=np.float32)
      b = np.array([4, 5, 6, 7, 8, 24], dtype=np.float32)
  else:
      a = None
      b = None
  
  element_a = np.empty(1, dtype=np.float32)
  comm.Scatter(a, element_a, root=0)
  element_b = np.empty(1, dtype=np.float32)
  comm.Scatter(b, element_b, root=0)
  produkt_elemenata = element_a * element_b
  
  vektor_produkata_elemenata = np.empty(6, dtype=np.float32)
  comm.Gather(produkt_elemenata, vektor_produkata_elemenata, root=0)
  
  if rank == 0:
      print("Skalarni produkt je", vektor_produkata_elemenata.sum())
  

• paralelni kod za koji koristi kolektivnu komunikaciju tipa scatter-reduce je oblika

  # varijanta s Python objektima
  from mpi4py import MPI
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = [1, 2, 3, 9, 10, 11]
      b = [4, 5, 6, 7, 8, 24]
  else:
      a = None
      b = None
  
  komponenta_od_a = comm.scatter(a, root=0)
  komponenta_od_b = comm.scatter(b, root=0)
  produkt = komponenta_od_a * komponenta_od_b
  
  skalarni_produkt_vektora = comm.reduce(produkt, op=MPI.SUM, root=0)
  if rank == 0:
      print("Skalarni produkt je", skalarni_produkt_vektora)
  

  # varijanta s NumPy poljima
  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      a = np.array([1, 2, 3, 9, 10, 11], dtype=np.float32)
      b = np.array([4, 5, 6, 7, 8, 24], dtype=np.float32)
  else:
      a = None
      b = None
  
  element_a = np.empty(1, dtype=np.float32)
  comm.Scatter(a, element_a, root=0)
  element_b = np.empty(1, dtype=np.float32)
  comm.Scatter(b, element_b, root=0)
  produkt_elemenata = element_a * element_b
  
  skalarni_produkt = np.empty(1, dtype=np.float32)
  comm.Reduce(produkt_elemenata, skalarni_produkt, op=MPI.SUM, root=0)
  if rank == 0:
      print("Skalarni produkt je", skalarni_produkt)
  

Produkt matrice i vektora

• sekvencijalni produkt matrice i vektora

  # varijanta s NumPy poljima
  import numpy as np
  
  A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.float32)
  x = np.array([10, 11, 12], dtype=np.float32)
  
  y = np.empty(3, dtype=np.float32)
  
  for i in range(3):
      y[i] = np.dot(A[i], x)
  
  print("Produkt matrice\n", A)
  print("i vektora\n", x)
  print("iznosi\n", y)
  
  # y = np.dot(A, x)
  

• paralelni kod korištenjem Scatter() i Reduce()

  # varijanta s NumPy poljima
  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  if rank == 0:
      A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [11, 22, 33]], dtype=np.float32)
      x = np.array([10, 11, 12], dtype=np.float32)
      y = np.empty(4, dtype=np.float32)
  else:
      A = None
      x = np.empty(3, dtype=np.float32)
      y = None
  
  redak_A = np.empty(3, dtype=np.float32)
  comm.Scatter(A, redak_A, root=0)
  comm.Bcast(x, root=0)
  
  element_y = np.dot(redak_A, x)
  comm.Gather(element_y, y, root=0)
  
  if rank == 0:
      print("Produkt matrice\n", A)
      print("i vektora\n", x)
      print("iznosi\n", y)
  

Aproksimativno računanje broja \(\pi\) kao sume niza

Broj \(\pi\) može se aproksimirati formulom

• sekvencijalni kod je oblika

  # varijanta s Python objektima
  import math
  
  n_iterations = 10
  
  def compute_pi(n):
      h = 1.0 / n
      s = 0.0
      for i in range(n):
          x = h * (i + 0.5)
          s += 4.0 / (1.0 + x**2)
      return s * h
  
  pi = compute_pi(n_iterations)
  error = abs(pi - math.pi)
  print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

  # varijanta s NumPy poljima
  import numpy as np
  
  n_iterations = 10
  
  def compute_pi(n):
      h = 1.0 / n
      s = 0.0
      for i in range(n):
          x = h * (i + 0.5)
          s += 4.0 / (1.0 + x**2)
      return s * h
  
  pi = compute_pi(n_iterations)
  error = abs(pi - np.math.pi)
  print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

• paralelni kod je oblika

  # varijanta s Python objektima
  from mpi4py import MPI
  import math
  
  n_iterations = 10
  
  def compute_pi(n, start=0, step=1):
      h = 1.0 / n
      s = 0.0
      for i in range(start, n, step):
          x = h * (i + 0.5)
          s += 4.0 / (1.0 + x**2)
      return s * h
  
  comm = MPI.COMM_WORLD
  size = comm.Get_size()
  rank = comm.Get_rank()
  
  if rank == 0:
      n = 10
  else:
      n = None
  
  n = comm.bcast(n, root=0)
  
  pi_part = compute_pi(n, rank, size)
  
  pi = comm.reduce(pi_part, op=MPI.SUM, root=0)
  
  if rank == 0:
     error = abs(pi - math.pi)
     print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

  # varijanta s NumPy poljima
  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  size = comm.Get_size()
  rank = comm.Get_rank()
  
  n_iterations = 1000
  
  def compute_pi(n, start=0, step=1):
      h = 1.0 / n
      s = 0.0
      for i in range(start, n, step):
          x = h * (i + 0.5)
          s += 4.0 / (1.0 + x**2)
      return s * h
  
  pi_part = np.empty(1)
  pi_part[0] = compute_pi(n_iterations, start=rank, step=size)
  
  if rank == 0:
      pi = np.empty(1)
  else:
      pi = None
  
  comm.Reduce(pi_part, pi, op=MPI.SUM, root=0)
  
  if rank == 0:
      error = abs(pi - np.math.pi)
      print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

Aproksimativno računanje broja \(\pi\) korištenjem Monte Carlo metode

Metode Monte Carlo su skupina metoda koja na temelju slučajnih brojeva i velikog broja pokusa daju aproksimaciju određene vrijednosti.

Broj \(\pi\) može se aproksimativno računati i korištenjem Monte Carlo metode. Da bi to vidjeli, uzmimo da je u koordinatnom sustavu ucrtan jedinični krug unutar jediničnog kvadrata. Promotrimo njegovu četvrtinu u prvom kvadrantu:

• četvrtina kvadrata je površine 1 (obzirom da su mu obje stranice duljine 1),
• čevrtina kruga je površine \(\frac{\pi}{4}\) (obzirom da mu je radijus 1).

Slučajno odabrana točka unutar četvrtine kvadrata ima vjerojatnost \(\frac{\pi}{4}\) da upadne unutar četvrtine kruga. Dakle, ako s \(n\) označimo broj slučajno odabranih točaka, a s \(h\) broj točaka koje se od \(n\) slučajno odabranih nalaze unutar četvrtine kruga, aproksimativno možemo odrediti broj \(\pi\) kao

Povećavajući \(n\) dobivamo točniju aproksimaciju, a navedena metoda naziva se Monte Carlo metoda. Programski kod koji aproksimira broj \(\pi\) korištenjem Monte Carlo metode je oblika:

• sekvencijalni kod

  import numpy as np
  # import random
  # import math
  
  n_random_choices = 10000
  hits = 0
  throws = 0
  
  for i in range (0, n_random_choices):
      throws += 1
      x = np.random.random() # x = random.random()
      y = np.random.random() # y = random.random()
      dist = np.math.sqrt(x * x + y * y) # dist = math.sqrt(x * x + y * y)
      if dist <= 1.0:
          hits += 1
  
  pi = 4 * (hits / throws)
  
  error = abs(pi - np.math.pi) # error = abs(pi - math.pi)
  print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

• paralelni kod korištenjem Reduce()

  import numpy as np
  from mpi4py import MPI
  
  comm = MPI.COMM_WORLD
  rank = comm.Get_rank()
  size = comm.Get_size()
  
  n_random_choices = 100000
  hits = 0
  throws = 0
  
  for i in range (0, n_random_choices):
      throws += 1
      x = np.random.random()
      y = np.random.random()
      dist = np.math.sqrt(x * x + y * y)
      if dist <= 1.0:
          hits += 1
  
  pi_part = np.empty(1)
  pi_part[0] = 4 * (hits / throws)
  
  if rank == 0:
      pi_reduced = np.empty(1)
  else:
      pi_reduced = None
  
  comm.Reduce(pi_part, pi_reduced, op=MPI.SUM, root=0)
  
  if rank == 0:
      pi = pi_reduced[0] / size
      error = abs(pi - np.math.pi)
      print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))
  

brojevi = [[1, 2, 3, ..., 125000],
           [125001, 125002, ..., 250000],
           [250001, 250002, ..., 375000],
           [375001, 375002, ..., 500000]]

Author: Vedran Miletić