UPC Version of the Sieve of Erathosthenes

Phil Merkey

Introduction

The Seive of Erathosthenes is not a good application for UPC, but it does give us an excuse to talk about different ways to control multiple threads.

The seive is simple program. Write down as many integers as you want starting at two, then cross out multiples of the primes. The simplicity of the algorithm is fact that the next smallest thing not crossed out is a prime. One crosses out the evens, then the next thing not crossed out is 3. Cross out all multiples of 3 and the next thing available is 5 and so on.

To find the primes less than a hundred with a C program we could write:

<serial-sieve.c>=
#include <stdio.h>

#define SqrtN 10
#define MaxN (SqrtN * SqrtN)

int z[MaxN];

main ()
{
    int i, sieve;

    for (i = 2; i < MaxN; i++)
        z[i] = i;

    sieve = 1;
    while( 1 )
    {
        for (sieve++; sieve <= SqrtN; sieve++) {
            if (z[sieve] > 0)
                break;
        }
        if( sieve >= SqrtN )
            break;
        for (i = 2 * sieve; i < MaxN; i += sieve)
            z[i] = 0;
    };

    for (i = 2; i < MaxN; i++) {
        if (z[i] > 0)
            printf("%d\n", z[i]);
    }
}

Sieve of Erathosthenes in UPC

The problem of ``parallelizing'' the sieve is to decide where the elements in the array are going to live and how the work is distributed among the threads. We will make the array a share array and distribute the work by having each thread sieve (ie, zero out the appropriate elements to which it has affinity). We have thread 0 find the next sieve and then since the variable sieve is shared all the other can use it.

<upc-sieve-barrier.c>=
#include <stdio.h>
#include <upc.h>

#define SqrtN 10
#define MaxN (SqrtN*SqrtN)

shared int z[MaxN];
shared int sieve;

main ()
{
  int i;

  upc_forall( i=2; i< MaxN; i++; &z[i] )
      z[i] =  i;

  upc_barrier(1);

  sieve = 1;
  while(1){
    if( MYTHREAD == 0 ) {    // find the sieve prime
       for( sieve++; sieve<=SqrtN; sieve++){
          if( z[sieve] > 0 )    // found it
             break;
       }
    }
    upc_barrier(1);  // everybody waits for the new sieve
    if( sieve >= SqrtN )  // we are done
       break;

    // here is the parallel part
    upc_forall( i = 2*sieve; i<MaxN; i += sieve ; &z[i])
      z[i] = 0;
    upc_barrier(2);
  }

  if( MYTHREAD == 0 ){  // have thread 0 print the result
    for(i=2; i<MaxN; i++) {
      if( z[i] != 0 )
         printf("%d\n", z[i]);
    }
  }
}

In the next example we will have each thread find it own prime to sieve with. The basic idea is to have each thread working independently. A thread get a lock on the critical section of code where finds the next prime and then it marks that element by making in negative. It then releases the lock and then zeros out all multiples of the prime it just found.

<upc-sieve-lock.c>=
#include <stdio.h>
#include <upc.h>

#define MaxN 100
#define SqrtN 10
shared int z[MaxN];

upc_lock_t *instr_lock;

main ()
{
  int i, j;
  int sieve;

  instr_lock = upc_all_lock_alloc();
  upc_lock_init( instr_lock );

  upc_forall( i=2; i< MaxN; i++; &z[i] )
      z[i] =  i;
  upc_barrier(1);

  do {
    // find my sieve prime
    upc_lock (instr_lock );
      sieve = 0;
      for( i=2; i<SqrtN; i++){
         sieve = z[i];
         if( sieve > 0 ) {  // found next candidate
            for ( j=i-1; j>1 ; j-- )
                if( z[j] < 0 ) {  // check against this prime
                   if( sieve % z[j] == 0 ){ // somebody hasn't got here yet
                                  // this works even with (z[j] < 0) 
                      sieve = 0;   // indicator to keep looking
                      break;
                   } 
                } 
                if( sieve > 0 ) { // found one
                    printf("%d: sieve is %d\n", MYTHREAD, sieve);
                    z[i] = -1 * z[i];
                    break;
                }
         }
      }
      upc_unlock (instr_lock );

      // now sieve across all the elements in the array
      for( i = 2*sieve; sieve && i<MaxN; i += sieve )
          z[i] = 0;

  } while( sieve );

  upc_barrier(2);

  if( MYTHREAD == 0 )
    for(i=2; i<MaxN; i++) {
      if( z[i] != 0 )
         printf("%d\n", z[i]);
    }
}

Of course the Sieve of Erathosthenes is really only interesting as a theoretical algorithm and doing it in parallel is even more silly then doing it serially. We choose it as a example just to show how one could have a critical section of code protected by a barrier or by a lock.