Eratosthenes sieve and MPI

Question

I wrote this sieve of Eratosthenes with MPI, but I'm not sure if it's good enough. Should I use MPI_Scatter and MPI_Gather instead of collecting arrays of separate processes in the root process?

Also, I'm not sure if I call

printArray(myArray, (N/size)-1);

right, where myArray is an array of pointers. GCC gives a warning at this point.

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

#define N 60

int main(int argc, char *argv[])
{
long timeCPU;
int rank, size, k, i, v, n, scatterSize;
int **tmpArray, *myArray, firstElement, lastElement;

MPI_Init(&argc,&argv);

MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);

if(N%size != 0)
{
    MPI_Finalize();
    printf("ERROR!");
    exit(0);
}

firstElement = (rank*N)/size;
lastElement = floor(((rank+1)*N)/size);
myArray = (int*) (malloc(((N/size)-1)*sizeof(int)));

for(i=0, v=firstElement+2; i<=(N/size)-1; i++, v++)
{
    myArray[i] = v;
}
    printArray(myArray, (N/size)-1);
/* All processes have to execute above code first */
MPI_Barrier(MPI_COMM_WORLD);

k = 2;

do
{
    markMultiples(k, myArray, lastElement-firstElement);
    k = nextNumber(k, myArray, lastElement-firstElement);
    MPI_Bcast(&k, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
while(k*k <=N);

MPI_Send(myArray, (N/size)-1, MPI_INT, 0, 50, MPI_COMM_WORLD);

MPI_Barrier(MPI_COMM_WORLD);

if (rank == 0)
{
    // Initialize 2D array
    tmpArray = (int**) (malloc(size*sizeof(int)));
    for(i=0;i<size;i++)
    {
        tmpArray[i] = (int*) (malloc(((N/size)-1)*sizeof(int)));
    }

    // Gather from every process his own array
    for (i=0; i<size;i++)
    {
        MPI_Recv(tmpArray[i], (N/size)-1, MPI_INT, i, 50, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        printArray(tmpArray[i], (N/size)-1);
        printf("-----------------------------------------\n");
    }

    // Build 1d array with primes - TODO

}

//printArray(array, lastElement-firstElement);
free(myArray);
MPI_Finalize();
return 0;
}

int nextNumber(int k, int *array, int n)
{
int i;

for(i=0; i<=n; i++)
{
    if(array[i] > k)
        return array[i];
}
return -1;
}

void markMultiples(int k, int *array, int n)
{
int i;

for (i=0; i<=n; i++)
{
    if(array[i] % k == 0 && array[i] != k)
    {
        array[i] = -1;
    }
}
}

void printArray(int *array, int n)
{
int i;

for(i=0; i<=n; i++)
{
    printf("array[%d] = %d\n", i, array[i]);
}
}

Answer 1

Since you tagged this question c++ I assume you would like a C++ approach also.

I would suggest using std::vector<int> rather than an int* and malloc. Using the vector will give you a number of useful advantages over a normal pointer, for example, memory management is greatly simplified, and there is no need to worry about how big an int is. With vector, the order of all elements is guaranteed to be contiguous.

#include <algorithm>
#include <functional>
#include <iostream>
#include <vector>

int nextNumber(int k, const std::vector<int> &array, int n)
{
    std::vector<int>::iterator i = std::find_if(array.begin(), array.begin() + n, std::bind2nd(std::greater<int>(), k));

    return (i == array.end() ? -1 : *i);
}

bool isMultiple(int a, int b)
{
    return (a % k == 0) && (a != k);
}

void markMultiples(int k, std::vector<int> &array, int n)
{
    std::replace_if(array.begin(), array.begin() + n, std::bind2nd(std::ptr_fun(isMultiple), k), -1);
}


void printArray(const std::vector<int> &array)
{
    for(int i = 0; i < array.size(); i++)
        std::cout << "array[" << i << "] = " << array[i] << std::endl;
}

Disclaimer:

I am not very good at C++, so you may want to take this advice with a teaspoon of salt.

Answer 2

• This error message reveals nothing about the actual error:

  if(N%size != 0)
  {
      MPI_Finalize();
      printf("ERROR!");
      exit(0);
  }
  

  Just displaying "ERROR!" won't tell the user what this exact error is and/or why it has happened. Yes, it's certainly important to make sure the input size and number of processes are divisible in MPI to maintain load balance, but the user must still know about this so that either input can be adjusted accordingly.

  It may also help others (and yourself) to display the computed portion size, whether or not this error has occurred. If you get an unreasonable value (usually if the portion size cannot fit inside an int (overflow), which is a commonly used size type by the MPI functions), then the user will know about this and will attempt to fix it for next time.

  This error, and any others, should be printed to stderr via fprintf() instead.

  As this is not a successful exit, you should use a different exit/return value other than 0. This may ideally depend on the language you'll be using (this uses C and C++), but it may be good to use 1 (or EXIT_FAILURE) for this.

  Finally, you can just have process 0 display the error so that it's displayed only once (any normal output should just be printed by this process). Of course, all of the processes should still call MPI_Finalize() and exit()/return().

• Regarding memory allocation:

  You allocate memory for tmpArray, but never deallocate it. Since it's only allocated by process 0, just have that process handle the deallocation, so that it's done only once (very important).

  if (rank == 0) free(tmpArray);
  

  You may also consider terminating the program if any allocation has failed:

  if (myArray == NULL) // nullptr if using C++11
  {
      fprintf(stderr, "could not allocate array\n");
      MPI_Finalize();
      return EXIT_FAILURE; // or use any suitable non-zero error value
  }
  

• Regarding cleanliness:

  Indentation is inconsistent in some places. Also, you don't indent entire function bodies, which could greatly hurt readability if one cannot tell where some code belongs.

  Some lines, such as this one:

  for(i=0, v=firstElement+2; i<=(N/size)-1; i++, v++)
  

  could use a little more whitespace between some operators:

  for (i = 0, v = firstElement+2; i <= (N/size)-1; i++, v++)
  

  It's mostly up to you, but it may be clearer to at least have the equality and assignment operators separated to help identify the operands.