/******************************************************************************
* FILE: mpi_prime.c
* DESCRIPTION:
*   Generates prime numbers.  All tasks distribute the work evenly, taking
*   every nth number, where n is the stride computed as:  (rank *2) + 1
*   so that even numbers are automatically skipped.  The method of using
*   stride is preferred over contiguous blocks of numbers, since numbers
*   in the higher range require more work to compute and may result in 
*   load imbalance.  This program demonstrates embarrassing parallelism.
*   Collective communications calls are used to reduce the only two data
*   elements requiring communications: the number of primes found and
*   the largest prime.
* AUTHOR: Blaise Barney 11/25/95 - adapted from version contributed by 
*   Richard Ng &  Wong Sze Cheong during MHPCC Singapore Workshop (8/22/95).
* LAST REVISED: 04/13/05
******************************************************************************/
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define LIMIT     2500000     /* Increase this to find more primes */
#define FIRST     0           /* Rank of first task */

int isprime(int n) {
int i,squareroot;
if (n>10) {
   squareroot = (int) sqrt(n);
      for (i=3; i<=squareroot; i=i+2)
            if ((n%i)==0)
                     return 0;
                        return 1;
                           }
                           /* Assume first four primes are counted elsewhere. Forget everything else */
                           else
                              return 0;
                              }


                              int main (int argc, char *argv[])
                              {
                              int   ntasks,               /* total number of tasks in partitiion */
                                    rank,                 /* task identifier */
                                          n,                    /* loop variable */
                                                pc,                   /* prime counter */
                                                      pcsum,                /* number of primes found by all tasks */
                                                            foundone,             /* most recent prime found */
                                                                  maxprime,             /* largest prime found */
                                                                        mystart,              /* where to start calculating */
                                                                              stride;               /* calculate every nth number */

                                                                              double start_time,end_time;
                                                                                
                                                                                MPI_Init(&argc,&argv);
                                                                                MPI_Comm_rank(MPI_COMM_WORLD,&rank);
                                                                                MPI_Comm_size(MPI_COMM_WORLD,&ntasks);
                                                                                if (((ntasks%2) !=0) || ((LIMIT%ntasks) !=0)) {
                                                                                   printf("Sorry - this exercise requires an even number of tasks.\n");
                                                                                      printf("evenly divisible into %d.  Try 4 or 8.\n",LIMIT);
                                                                                         MPI_Finalize();
                                                                                            exit(0);
                                                                                               }

                                                                                               start_time = MPI_Wtime();   /* Initialize start time */
                                                                                               mystart = (rank*2)+1;       /* Find my starting point - must be odd number */
                                                                                               stride = ntasks*2;          /* Determine stride, skipping even numbers */
                                                                                               pc=0;                       /* Initialize prime counter */
                                                                                               foundone = 0;               /* Initialize */

                                                                                               /******************** task with rank 0 does this part ********************/
                                                                                               if (rank == FIRST) {
                                                                                                  printf("Using %d tasks to scan %d numbers\n",ntasks,LIMIT);
                                                                                                     pc = 4;                  /* Assume first four primes are counted here */
                                                                                                        for (n=mystart; n<=LIMIT; n=n+stride) {
                                                                                                              if (isprime(n)) {
                                                                                                                       pc++;
                                                                                                                                foundone = n;
                                                                                                                                         /***** Optional: print each prime as it is found
                                                                                                                                                  printf("%d\n",foundone);
                                                                                                                                                           *****/
                                                                                                                                                                    }
                                                                                                                                                                          }
                                                                                                                                                                             MPI_Reduce(&pc,&pcsum,1,MPI_INT,MPI_SUM,FIRST,MPI_COMM_WORLD);
                                                                                                                                                                                MPI_Reduce(&foundone,&maxprime,1,MPI_INT,MPI_MAX,FIRST,MPI_COMM_WORLD);
                                                                                                                                                                                   end_time=MPI_Wtime();
                                                                                                                                                                                      printf("Done. Largest prime is %d Total primes %d\n",maxprime,pcsum);
                                                                                                                                                                                         printf("Wallclock time elapsed: %.2lf seconds\n",end_time-start_time);
                                                                                                                                                                                            }


                                                                                                                                                                                            /******************** all other tasks do this part ***********************/
                                                                                                                                                                                            if (rank > FIRST) {
                                                                                                                                                                                               for (n=mystart; n<=LIMIT; n=n+stride) {
                                                                                                                                                                                                     if (isprime(n)) {
                                                                                                                                                                                                              pc++;
                                                                                                                                                                                                                       foundone = n;
                                                                                                                                                                                                                                /***** Optional: print each prime as it is found
                                                                                                                                                                                                                                         printf("%d\n",foundone);
                                                                                                                                                                                                                                                  *****/
                                                                                                                                                                                                                                                           }
                                                                                                                                                                                                                                                                 }
                                                                                                                                                                                                                                                                    MPI_Reduce(&pc,&pcsum,1,MPI_INT,MPI_SUM,FIRST,MPI_COMM_WORLD);
                                                                                                                                                                                                                                                                       MPI_Reduce(&foundone,&maxprime,1,MPI_INT,MPI_MAX,FIRST,MPI_COMM_WORLD);
                                                                                                                                                                                                                                                                          }

                                                                                                                                                                                                                                                                          MPI_Finalize();
                                                                                                                                                                                                                                                                          }