/******************************************************************************
* TotalView MPI Example 
* The number of MPI tasks must be equal to 4.
****************************************************************************/
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#define  ARRAYSIZE	160000
#define  MASTER		0

double *data;

int main (int argc, char *argv[])

{
int   numtasks, taskid, rc=0, dest, myoffset, offset, i, j, 
      tag1=1, tag2=2, source, chunksize; 
double mysum, sum, isum;
double update(int myoffset, int chunk, int myid);
double *Mdata;
MPI_Status status;

/***** Initializations *****/
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
if (numtasks != 4) {
   printf("Quitting. Number of MPI tasks must equal 4.\n");
   MPI_Abort(MPI_COMM_WORLD, rc);
   exit(0);
   }
MPI_Comm_rank(MPI_COMM_WORLD,&taskid);
printf ("MPI task %d has started...\n", taskid);
chunksize = (ARRAYSIZE / numtasks);


/***** Master task only ******/
if (taskid == MASTER){
  Mdata=(double*)malloc(sizeof(double)*ARRAYSIZE);

  /* Initialize the array */
  isum = 0.0;
  for(i=0; i<ARRAYSIZE; i++) {
    Mdata[i] =  i * 1.0;
    isum = isum + Mdata[i];
    }
  printf("*** Initial sum = %e\n",isum);
}

  /* Alloc memory for each Task */
  data=(double*)malloc(sizeof(double)*chunksize);

if (taskid == MASTER){
  /* Send each task its portion of the array - master keeps 1st part */
  offset = chunksize;
  for (dest=1; dest<numtasks; dest++) {
    MPI_Send(&offset, 2, MPI_INT, dest, tag1, MPI_COMM_WORLD);
    printf("Send %d elements to task %d offset= %d\n",chunksize,dest,offset);
    offset = offset + chunksize;
    }
  offset = 0;
} else {
  /* Receive my offset from the master task */
  source = MASTER;
  MPI_recv(&offset, 1, MPI_INT, source, tag1, MPI_COMM_WORLD, &status);
}

/* Distribute data */
    MPI_scatter(Mdata, ARRAYSIZE, MPI_DOUBLE, data,chunksize, MPI_DOUBLE, 
                MASTER,  MPI_COMM_WORLD);

  if (taskid == MASTER) free(Mdata);

  mysum = update( offset, chunksize, taskid);
  sum=0.0;

  /* Get final sum and print sample results */  
  MPI_Reduce(&mysum, &sum, 1, MPI_DOUBLE, MPI_SUM, MASTER, MPI_COMM_WORLD);

if (taskid == MASTER){
  printf("*** Final sum= %e *** (Twice Initial %e*%.1f) \n",sum,isum,sum/isum);
}

MPI_Finalize();

}   /* end of main */


double update(int myoffset, int chunk, int myid) {
  int i; 
  double mysum;
  /* Perform addition to each of my array elements and keep my sum */
  mysum = 0.0;
  for(i=0; i < chunk; i++) {
    data[i] = data[i] + (i+myoffset) * 1.0;
    mysum = mysum + data[i];
    }
  printf("Task %d mysum = %e\n",myid,mysum);
  return(mysum);
  }