quick example of using MPI to parallelise trial division primality checking

trial-division.c

#include <mpi.h>

#include <stdio.h>  
#include <string.h>  
#include <stdlib.h>

typedef unsigned long long ull;
// all communication is either input to workers or output from workers
enum {TAG_INPUT = 2001, TAG_RESULT};

ull isDivisible(ull number, ull from, ull to); // work function
ull getUll(void); // get ull from stdin

int numprocs, // number of worker processes
	myid; // id of current process, main is 0

int main(int argc, char **argv)  
{  
	MPI_Status status;
	// done everytime
	MPI_Init(&argc, &argv);  
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);  
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);  

	if(myid == 0)  
	{ // we are now in the main process
		if (numprocs < 2)
		{
			puts("Error: at least one worker process needed.");
			return 1;
		}
		ull input, divisor;

		do {
			printf("Number to primalty-check: ");
			fflush(stdout);
			input = getUll();
		} while (!input);

		for (int id = 1; id < numprocs; id++) // for all worker processes 
		{
			// pass the number and the range that should be checked
			ull n[] = { input, (input / (numprocs-1))*(id-1), (input / (numprocs-1))*id };
			MPI_Send(n, 3, MPI_UNSIGNED_LONG_LONG, id, TAG_INPUT, MPI_COMM_WORLD);
		}

		for(int id = 1; id < numprocs; id++)  // recieve from worker processes
		{
			//read divisor
			MPI_Recv(&divisor, 1, MPI_UNSIGNED_LONG_LONG, MPI_ANY_SOURCE,
					TAG_RESULT, MPI_COMM_WORLD, &status);

			int sender = status.MPI_SOURCE;

			printf("Process %02i reports that ", sender);
			if (divisor == 0)
				printf("no divisors were found in his range\n"); 
			else
			{
				printf("divisor %llu was found in his range\n", divisor);
				break; // one divisor found, no need to check anything else
			}

		}
		if(divisor == 0)
			puts("Number is prime!");
		else
			puts("Number is not a prime.");

		printf("%s\n", "End running!");  
		MPI_Abort(MPI_COMM_WORLD, 0); // very bad way
	}
	else 
	{ // worker processes
		ull n[3];
		// recieve task from main process
		MPI_Recv(&n, 3, MPI_UNSIGNED_LONG_LONG, 0, TAG_INPUT, MPI_COMM_WORLD, &status);

		// check for 0 and 1 beforehand
		ull result = n[0] > 1 ? isDivisible(n[0], n[1], n[2]) : 1;

		// send back divisor found (or 0 if none)
		MPI_Send(&result, 1, MPI_UNSIGNED_LONG_LONG, 0, TAG_RESULT, MPI_COMM_WORLD);
	}
	MPI_Finalize(); // very very bad  

	return 0;  
}  
ull isDivisible(ull number, ull from, ull to)
{
	// sanitize input
	if (from <= 1) from = 2;
	if (number == to) to--;

	// trial divison
	for (ull i = from; i < to; i++)
		if (number % i == 0) return i;

	return 0; // might be a prime
}
ull getUll(void)
{
#if INPUT
	printf("%llu\n", (ull)INPUT);
	return INPUT;
#else
	char str[32];
	fgets(str, sizeof str, stdin);
	return strtoull(str, NULL, 0);
#endif
}

How to build:

• Install (Open-)MPI (apt-get, brew, w/e)
• compile with mpicc -oprime trial-divison.c -O3 -DINPUT=0x7FFFFFFFULL or leave out the #define for stdin usage
• run with mpirun -n 5 ./prime. I guess the best way to fare is (n+1) processes for n processor cores
• you can contrast that with running with 2 processes (dispatcher/main + worker process). My results:

n = 2

real    0m19.257s
user    0m37.786s
sys 0m0.176s

n = 5

real    0m5.528s
user    0m21.193s
sys 0m0.135s

n = 50

real    0m6.181s
user    0m22.347s
sys 0m0.905s