Paraprimes

En el artículo anterior se mostraba el cálculo de números primos con la librería GMP y en el anterior a éste, el cálculo con SQL y Perl.

Ahora vamos un paso más allá, distribuyendo el cálculo en procesos paralelos. Para ello, hay múltiples soluciones, como por ejemplo:

Sin embargo, la opción más sencilla probablemente sea OpenMP, un estándar abierto que ya se encuentra integrado en el compilador GCC.

El siguiente código, que se compila con «gcc -fopenmp -o hola hola.c«, es un ejemplo muy sencillo del uso de OpenMP:

hola.c

[c]
#include <omp.h>
#include <stdio.h>

// @author José Cotrino (https://www.cotrino.com)
int main() {
#pragma omp parallel
printf("- Hola caracola desde el hilo %d, hay %d hilos\n",
omp_get_thread_num(), omp_get_num_threads());
const int N = 100000;
const int MAX_THREADS = 20;
int i, a[N], p[MAX_THREADS];
for (i=0;i<MAX_THREADS;i++)
p[i]= 0;
printf("\n- Calculando un bucle 'for' de %d iteraciones en paralelo: \n",N);
printf(" + cada iteración puede ser asignada teóricamente a un procesador\n");
printf(" + el array sobre el que se calcula está en memoria compartida\n");
#pragma omp parallel for
for (i=0;i<N;i++) {
a[i]= 2*i;
p[omp_get_thread_num()] ++;
}
printf("\n- Yatá!\n");
for (i=0;i<MAX_THREADS;i++) {
if (p[i] != 0)
printf(" + Hilo %d se ha ejecutado %d veces\n", i, p[i]);
}
}
[/c]

En un Intel Core 2 Duo obtenemos la siguiente salida:

- Hola caracola desde el hilo 1, hay 2 hilos
- Hola caracola desde el hilo 0, hay 2 hilos
- Calculando un bucle 'for' de 100000 iteraciones en paralelo:
+ cada iteración puede ser asignada teóricamente a un procesador
+ el array sobre el que se calcula está en memoria compartida
- Yatá!
+ Hilo 0 se ha ejecutado 50000 veces
+ Hilo 1 se ha ejecutado 50000 veces

Volviendo al código que teníamos para calcular primos con GMP, lo que haremos es simplemente dividir todo el rango en varios trozos, de forma que cada procesador calcule los números primos dentro de un subrango. Por ejemplo, si tenemos 2 núcleos y queremos calcular todos los primos entre 0 y 1000, en un procesador calculamos entre 0 y 500, y en el otro entre 501 y 1000. Esto es posible porque el test de primalidad que utiliza GMP no necesita conocer los primos anteriores.

El código queda así:

paraprimes.cpp

[cpp]
#include <stdio.h>
#include <time.h>
#include <gmp.h>
#include <omp.h>

// @author José Cotrino (https://www.cotrino.com)
void calculatePrimes ( mpz_t N_max ) {
time_t start,end;
double elapsed;
mpz_t N_primes, N_range;
mpz_init(N_primes);
mpz_init(N_range);
mpz_set_ui(N_primes, 0);
mpz_set(N_range, N_max);

time(&start);

// calculate N/threads
int N_cores = omp_get_num_procs();
mpz_cdiv_q_ui(N_range, N_range, N_cores);

// each thread calculates the primes between N/t*i and N/t*(i+1), where:
// - N is the total range where to search for primes
// - t is the amount of available processors
// - i is the number of the current thread/processor (0, 1, 2, 3...)
#pragma omp parallel shared(N_range, N_primes) {
int thread = omp_get_thread_num();
mpz_t prime, N_max_thread, N_primes_thread;
mpz_init(prime);
mpz_init(N_max_thread);
mpz_init(N_primes_thread);
mpz_set(prime, N_range);
mpz_mul_ui(prime, prime, thread);
mpz_add_ui(prime, prime, 1);
mpz_set(N_max_thread, N_range);
mpz_mul_ui(N_max_thread, N_max_thread, thread+1);
mpz_set_ui(N_primes_thread, 0);

//printf("Thread %d searching in range %s, %s\n", thread, mpz_get_str(NULL,10,prime), mpz_get_str(NULL,10,N_max_thread));
while( mpz_cmp(prime, N_max_thread) < 0 ) {
mpz_nextprime(prime, prime);
//printf("%s, ", mpz_get_str(NULL,10,prime));
mpz_add_ui(N_primes_thread, N_primes_thread, 1);
}
// skip the last prime, which is over the range
mpz_sub_ui(N_primes_thread, N_primes_thread, 1);
// add the primes found by this thread to the whole amount
mpz_add(N_primes, N_primes, N_primes_thread);
}

time(&end);
elapsed = difftime(end, start);

printf("%s primes found below %s in %.0f seconds\n",
mpz_get_str(NULL,10,N_primes),
mpz_get_str(NULL,10,N_max),
elapsed);
//printf("%s %f\n",mpz_get_str(NULL,10,N_max),elapsed);
}

int main() {
int i;
mpz_t N_max;
mpz_init(N_max);
mpz_set_ui(N_max, 100);

printf("%d cores available\n", omp_get_num_procs());

for(i=0; i<8; i++) {
calculatePrimes(N_max);
mpz_mul_ui(N_max, N_max, 10);
}
}
[/cpp]

Y se compila con:
g++ -lgmp -fopenmp paraprimes.cpp -o paraprimes

Los resultados en un Intel Core 2 Duo T5750 son impresionantes:

Máximo número Primos Duración
100 25 0s
1.000 168 0s
10.000 1.229 0s
100.000 9.592 0s
1.000.000 78.498 0s
10.000.000 664.579 7s
100.000.000 5.761.455 66s (~1.1min)
1.000.000.000 50.847.534 703s (~11.7min)

Comparamos de nuevo con las implementaciones anteriores para calcular los primos por debajo de 100 millones:

  • El algoritmo en C utilizando GMP sin paralelizar necesitaba 126 segundos. Con dos procesadores necesitamos 66 segundos, es decir, es 1.91 veces más rápido.
  • El algoritmo en Perl necesitaba 5.929 segundos, que comparados con los 66 segundos actuales nos da que este algoritmo en paralelo es 89.83 veces más rápido.

 

Teóricamente, con este algoritmo y si la ley de Amdahl lo permite, con un quad-core se podrían calcular todos los primos por debajo de 100 millones en apenas 35 segundos, con un octo-core en 17 segundos, y con un Tesla en ná y menos. Me mantengo a la espera de que alguien me regale uno para poder hacer la prueba.