/* this include file is mandatory */
#include <mpi.h>

/* MPI sample program; NOT INTENDED TO BE EFFICIENT as a prime
   finder, either in algorithm or implementation 

   MPI (Message Passing Interface) is a popular package using
   the "message passing" paradigm for communicating between
   processors in parallel applications; as the name implies,
   processors communicate by passing messages using "send" and 
   "receive" functions

   finds and reports the number of primes less than or equal to N

   each node works on a range of Chunk numbers at a time, testing
   them for primeness by straight division; if a potential divisor
   has already been found to be composite (i.e. nonprime), then we
   skip it

*/

#define MAX_N 1000000
#define NEW_MSG 0  /* type of message containing the new composites */
#define COUNT_MSG 0  /* type of message containing the local composite count */

int NNodes,  /* number of nodes in computation*/
    N,  /* find all primes from 2 to N */
    Chunk,  /* chunk size */
    Block,  /* number of integers checked in each iteration */
    Me,  /* my node number */
    Composite[MAX_N],  /* "global" array of composites found so far,
                          code 1 for composite, 0 otherwise */
    NNew,  /* number of new composites found in this iteration*/
    New[MAX_N],  /* new composites found (first component is NNew */
    MyComposites = 0;  /* total number of composites my node has found,
                          in all iterations so far */

double T1,T2;  /* start and finish times */

Init(Argc,Argv)
   int Argc;  char **Argv;

{  int DebugWait;  

   /* mandatory to begin any MPI program */
   MPI_Init(&Argc,&Argv);

   /* do not access command-line arguments until after MPI_Init()
      is called */
   N = atoi(Argv[1]); 
   Chunk = atoi(Argv[2]);
   DebugWait = atoi(Argv[3]);

   /* this loop is there to synchronize all nodes for debugging */
   while (DebugWait) ;

   /* C library function for convenient initialization to 0 */
   bzero(Composite,MAX_N*sizeof(int));

   /* puts the number of nodes in NNodes */
   MPI_Comm_size(MPI_COMM_WORLD,&NNodes);
   /* puts the node number of this node in Me */
   MPI_Comm_rank(MPI_COMM_WORLD,&Me); 
   Block = NNodes * Chunk;
    /* if N is not an even multiple of Block, then leave, and have
       node 0 explain why (we don't want all nodes to print out 
       redundant versions of this message) */
   if (N % Block > 0)  {
      if (Me == 0) printf("N must be divisible by Block\n");
      exit(1);
   }

   /* OK, get started; first record current time in T1 */
   if (Me == 0) T1 = MPI_Wtime();
}

int IsComposite(K)  /* returns 1 if K is found to be composite */
   int K;

{  int J;

   for (J = 2; J*J <= K; J++)  {
      /* if J is known to be composite, don't bother dividing by it */
      if (!Composite[J])  
         if (K % J == 0) return 1;
   }
   return 0;
}

void UpdateComposite()  /* receive the lists of newly-discovered
                           composites and update my copy of the array */

{  int Node,Them,I,Error;
   MPI_Status Status;  
  
   for (Node = 0; Node < NNodes; Node++)  {
      if (Node != Me)  {
         /* MPI_Recv  --  receive a message
               parameters:
	       pointer to place to store message
	       number of items in message (see notes on 
	          this at the end of this file)
	       item type
	       accept message from which node(s)
	       message type ("tag"), programmer-defined; here it 0,
	          so we are telling MPI to wait until a message of
		  type 0 comes in (all others will be queued)
	       node group number (in this case all nodes)
	       status (see notes on this at the end of this file) */
         Error = MPI_Recv(New,Chunk,MPI_INT,MPI_ANY_SOURCE,NEW_MSG,
            MPI_COMM_WORLD,&Status);
         MPI_Get_count(&Status,MPI_INT,&NNew);
         for (I = 0; I < NNew; I++) Composite[New[I]] = 1;
      }
   }
}

PerformIterations()

{  int First,Last,Iter,NIters,I,Node,NNew,Error;

   First = 1 + Me * Chunk;
   Last = First + Chunk - 1;
   NIters = N / Block;
   for (Iter = 0; Iter < NIters; Iter++)  {
      NNew = 0;
      for (I = First; I <= Last; I++)  
         if (IsComposite(I))  {
	    New[NNew++] = I;
	    Composite[I] = 1;
	    MyComposites++;
	 }
      /* now send all our newfound composites to all the other nodes;
         we could use MPI_BROADCAST() for this, but not with any
         appreciable gain in speed, since MPI will still send
	 out separate messages to each node (except in a special
	 implementation), so we will do this explicitly: */
      for (Node = 0; Node < NNodes; Node++)
         if (Node != Me)
	    /* MPI_Send  --  send a message
               parameters:
	       pointer to place where message is to be drawn from
	       number of items in message
	       item type
	       destination node
	       message type ("tag") programmer-defined; here we are
	          sending a message of type 0, which we will signify
		  a list of new composities
	       node group number (in this case all nodes) */
            Error = MPI_Send(New,NNew,MPI_INT,Node,NEW_MSG,MPI_COMM_WORLD);
      UpdateComposite();
      First += Block;  Last += Block;
   }
}

WrapUp()

{  int Count,Node; MPI_Status Status;  

   /* if not node 0, send node 0 our composite counts; if node 0,
      then receive and process them; note we have called this
      message type 1 */
   if (Me != 0)  MPI_Send(&MyComposites,1,MPI_INT,0,COUNT_MSG,MPI_COMM_WORLD);
   else  {
      Count = MyComposites;
      for (Node = 1; Node <NNodes; Node++)  {
         MPI_Recv(New,1,MPI_INT,MPI_ANY_SOURCE,COUNT_MSG,MPI_COMM_WORLD,
            &Status);
         Count += New[0];
      }
   }
   /* print out the results, again with only node 0 doing the printing */
   if (Me == 0)  {
      /* check the time again, and subtract to find run time */
      T2 = MPI_Wtime();
      printf("elapsed time = %f\n",(float)(T2-T1));
      printf("number of primes = %d\n",N-Count-1);
   }
   /* mandatory for all MPI programs */
   MPI_Finalize();
}

main(argc,argv)
   int argc; char **argv;

{  Init(argc,argv);
   PerformIterations();
   WrapUp();
}

/* explanation of "number of items" and "status" arguments at the end 
   of MPI_Recv():

   when receiving a message you must anticipate the longest possible
   message, but the actual received message may be much shorter than
   this; you can call the MPI_Get_count() function on the status
   argument to find out how many items were actually received, as
   we have done here

   the status argument will be a pointer to a struct, containing the
   node number, message type and error status of the received
   message

   say our last parameter is Status; then Status.MPI_SOURCE
   will contain the number of the sending node, and 
   Status.MPI_TAG will contain the message type; these are
   important if used MPI_ANY_SOURCE or MPI_ANY_TAG in our
   node or tag fields but still have to know who sent the
   message or what kind it is */