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

/* common globals for parallel processing
   nproc = number of real or virtual processors obtained
   lgrp = current communicator
   mreal = default datatype for reals
   mint = default datatype for integers               */

int nproc;
MPI_Comm lgrp;
MPI_Datatype mreal, mint;

FILE *unit2;

void ppinit(int argc, char *argv[], int *idproc, int *nvp);
void ppexit(void);
void timera(int icntrl, char *chr, float *time);

float fltmax(float val1, float val2) {
   return (val1 > val2 ? val1 : val2);
}

/* this is a network performance ping-pong test
   where one processor sends data and another echoes it back
   on each iteration, the length of data sent doubles         */

/* iter = number of different message lengths */
#define iter   21
/* nter = number of interations within timing loop */
#define nter    1
/* mter = number of times test is repeated to accumulate statistics */
#define mter   10

#define maxl   (1L<<iter-1)

int isbuf[maxl], irbuf[maxl];

int main(int argc, char *argv[])
{
   float rate[iter], var[iter], rmin[iter], rmax[iter];
   int nerr[iter], mist[iter];
   char label[17], name[37];
   MPI_Status istat;
   MPI_Request msid;
   int i, j, k;
   int lint, idproc, nvp, id0, id1, lmsg, ierr;
   float abytes, altime, time, r, anorm;
   FILE *unit7, *unit8;

/* SET MACHINE SPECIFIC DATA
   lint = number of bytes in an integer
   label = machine name and type     */

   lint = sizeof(int);
   strcpy(label,"Mac G3 Cluster");
   printf("Ping Pong test initializing...\n"); 

/* initialize data */

   abytes = (float)(2*lint*nter)*1.0e-6;
   for (i = 0; i < maxl; i++) {
      isbuf[i] = i+1;
      irbuf[i] = 0;
   }

/* initialize for parallel processing */

   nvp = 0;
   ppinit(argc, argv, &idproc, &nvp);

/* set ids for participating processors */

   id0 = 2*(idproc/2);
   id1 = id0 + 1;

/* suppress last node if number of nodes is odd */

   if (id1==nvp) 
      id0 = id1;

/* perform ping-pong */

   altime = 0.;

/* clear out counters */

   for (i = 0; i < iter; i++) {
      rate[i] = 0.;
      var[i] = 0.;
      rmin[i] = -1.;
      rmax[i] = 0.;
      nerr[i] = 0;
      mist[i] = 0;
   }

   for (k = 0; k < mter; k++) {

   sprintf(name,"Ping-Pong Test, Pass %2d\n",k+1);
   Logname(name);

/* loop over different message lengths */

      for (i = 0; i < iter; i++) {
         lmsg = 1L<<i;
/*       if (idproc==0) 
            printf("data length = %d\n",lmsg); */
         timera(-1,"pingpong",&time);
         for (j = 0; j < nter; j++) {
            if (idproc==id0) {
               ierr = MPI_Send(isbuf,lmsg,mint,id1,i+1,lgrp);
               ierr = MPI_Recv(irbuf,lmsg,mint,id1,i+2,lgrp,&istat);
            }
            else if (idproc==id1) {
               ierr = MPI_Recv(irbuf,lmsg,mint,id0,i+1,lgrp,&istat);
               ierr = MPI_Send(irbuf,lmsg,mint,id0,i+2,lgrp);
            }
         }
         timera(1,"pingpong",&time);
         if (idproc==0) {

/* compare result with expected sum */

            ierr = 0;
            for (j = 0; j < lmsg; j++) {
               if (isbuf[j] != irbuf[j])
                  ierr++;
            }
            nerr[i] += ierr;

/* calculate result */

            if (time > 0.) {
               r = (abytes*(float) lmsg)/time;
               rate[i] += r;
               var[i] += r*r;
               if (rmin[i] < 0)
                  rmin[i] = r;
               else
                  rmin[i] = r < rmin[i] ? r : rmin[i];
               rmax[i] = r > rmax[i] ? r : rmax[i];
            }
            else
               mist[i] += 1;
/*          printf("rate(MBS) = %f, errors = %d\n",r,ierr); */
            if (i==0)
               altime += time;
         }
      }
   }

/* write out results to file */

   if (idproc==0) {
      unit7 = fopen("PING.DAT","w");
      if (mter > mist[0])
         anorm = 1./(float) (mter - mist[0]);
      else
         anorm = 0;
      altime = 500.*altime*anorm/(float) nter;
      fprintf(unit7,"# Ping-Pong Test on %s\n",label);
      fprintf(unit7,"# Average latency(msec) = %f\n",altime);
      fprintf(unit7,"Rate(MBS) variance ln2(data) data_length min max errors\n");
      for (i = 0; i < iter; i++) {
         if (mter > mist[i])
            anorm = 1./(float) (mter - mist[i]);
         else
            anorm = 0;
         rate[i] = rate[i]*anorm;
         var[i] = sqrt(fltmax((var[i]*anorm - rate[i]*rate[i]),0.));
         lmsg = 1L<<i;
         fprintf(unit7,"%f %f %d %d %f %f %d\n",
                 rate[i],var[i],i,lmsg,rmin[i],rmax[i],nerr[i]);
      }
   }

/* perform swap */

   altime = 0.;

/* clear out counters */

   for(i = 0; i < iter; i++) {
      rate[i] = 0.;
      var[i] = 0.;
      rmin[i] = -1.;
      rmax[i] = 0.;
      nerr[i] = 0;
      mist[i] = 0;
   }

   for (k = 0; k < mter; k++) {

   sprintf(name,"Swap Test, Pass %2d\n",k+1);
   Logname(name);

/* loop over different message lengths */

      for (i = 0; i < iter; i++) {
         lmsg = 1L<<i;
/*       if (idproc==0)
            printf("data length = %d\n",lmsg); */
         timera(-1,"swap    ",&time);
         for (j = 0; j < nter; j++) {
            if (idproc==id0) {
               ierr = MPI_Irecv(irbuf,lmsg,mint,id1,i+2,lgrp,&msid);
               ierr = MPI_Send(isbuf,lmsg,mint,id1,i+1,lgrp);
               ierr = MPI_Wait(&msid,&istat);
            }
            else if (idproc==id1) {
               ierr = MPI_Irecv(irbuf,lmsg,mint,id0,i+1,lgrp,&msid);
               ierr = MPI_Send(isbuf,lmsg,mint,id0,i+2,lgrp);
               ierr = MPI_Wait(&msid,&istat);
            }
         }
         timera(1,"swap    ",&time);
         if (idproc==0) {

/* compare result with expected sum */
            ierr = 0;
            for (j = 0; j < lmsg; j++) {
               if (isbuf[j] != irbuf[j]) 
                  ierr++;
            }
            nerr[i] += ierr;

/* calculate result */

            if (time > 0.) {
               r = (abytes*(float) lmsg)/time;
               rate[i] += r;
               var[i] += r*r;
               if (rmin[i] < 0)
                  rmin[i] = r;
               else
                  rmin[i] = r < rmin[i] ? r : rmin[i];
               rmax[i] = r > rmax[i] ? r : rmax[i];
            }
            else
               mist[i] += 1;
/*         printf("rate(MBS) = %f, errors = %d\n", r,ierr); */
           if (i==0) 
              altime += time;
         }
      }
   }

/* write out results to file */

   if (idproc==0) {
      unit8 = fopen("SWAP.DAT","w");
      if (mter > mist[0])
         anorm = 1./(float) (mter - mist[0]);
      else
         anorm = 0;
      altime = 500.*altime*anorm/(float) nter;
      fprintf(unit8,"# Swap Test on %s\n",label);
      fprintf(unit8,"# Average latency(msec) = %f\n",altime);
      fprintf(unit8,"Rate(MBS) variance ln2(data) data_length min max errors\n");
      for (i = 0; i < iter; i++) {
         if (mter > mist[i])
            anorm = 1./(float) (mter - mist[i]);
         else
            anorm = 0;
         rate[i] = rate[i]*anorm;
         var[i] = sqrt(fltmax((var[i]*anorm - rate[i]*rate[i]),0.));
         lmsg = 1L<<i;
         fprintf(unit8,"%f %f %d %d %f %f %d\n",
                 rate[i],var[i],i,lmsg,rmin[i],rmax[i],nerr[i]);
      }
   }

/* terminate parallel processing */

   ppexit();
   return 0;

#undef iter
#undef nter
#undef mter
#undef maxl
}

/*--------------------------------------------------------------------*/

void ppinit(int argc, char *argv[], int *idproc, int *nvp)
{
/* this subroutine initializes parallel processing
   input: nvp, output: idproc
   idproc = processor id
   nvp = number of real or virtual processors requested

common globals for parallel processing
   nproc = number of real or virtual processors obtained
   lgrp = current communicator
   mreal = default datatype for reals
   mint = default datatype for integers

local data
   ndprec = (0,1) = (no,yes) use (normal,autodouble) precision */

   int ierr, ndprec = 0;

   unit2 = fopen("FOR002.DAT","w");

/* initialize the MPI execution environment */

   ierr = MPI_Init(&argc, &argv);
   if (ierr)
      exit(1);
   lgrp = MPI_COMM_WORLD;

/* determine the rank of the calling process in the communicator */

   ierr = MPI_Comm_rank(lgrp,idproc);

/* determine the size of the group associated with a communicator */

   ierr = MPI_Comm_size(lgrp,&nproc);

/* set default datatypes */

   mint = MPI_INT;

/* single precision */

   if (ndprec==0)
      mreal = MPI_FLOAT;

/* double precision */

   else
      mreal = MPI_DOUBLE;

/* requested number of processors not obtained */

   if (*nvp==0)
      *nvp = nproc;
   else if (nproc != *nvp) {
      fprintf(unit2,"processor number error: nvp, nproc=%d %d\n",
              *nvp,nproc);
      ppexit();
      exit(1);
   }
}

/*--------------------------------------------------------------------*/
void ppexit(void)
{
/* this subroutine terminates parallel processing

common globals for parallel processing
   lgrp = current communicator                 */

   int ierr;

/* synchronize processes */

   ierr = MPI_Barrier(lgrp);

/* terminate MPI execution environment */

   ierr = MPI_Finalize();

/* delete file if empty */

   if (!fseek(unit2,0,SEEK_END)) {
      ierr = ftell(unit2);
      fclose(unit2);
      if (!ierr)
         remove("FOR002.DAT");
   }
}
/*--------------------------------------------------------------------*/
void timera(int icntrl, char *chr, float *time)
{
/* this subroutine performs timing
   input: icntrl, chr
   icntrl = (-1,0,1) = (initialize,ignore,read) clock
   clock should be initialized before it is read!
   chr = string variable for labeling timings
   time = elapsed time in seconds
   written for mpi

common globals for parallel processing
   lgrp = current communicator
   mreal = default datatype for reals

local data                                          */

   int idproc, ierr;
   float nclock, mclock;
   static double jclock = 0.;
   if (icntrl==0)
      return;

/* read clock and write time difference from last clock initialization */

   if (icntrl==1) {
      nclock = (float) (MPI_Wtime() - jclock);
      ierr = MPI_Allreduce(&nclock,time,1,mreal,MPI_MIN,lgrp);
      mclock = *time;
      ierr = MPI_Allreduce(&nclock,time,1,mreal,MPI_MAX,lgrp);
      ierr = MPI_Comm_rank(lgrp,&idproc);
/*    if (idproc==0)
         printf("%s max/min real time = %14.7e %14.7e sec\n",
                chr,*time,mclock); */
   }

/* initialize clock */

   else {
      ierr = MPI_Barrier(lgrp);
      jclock = MPI_Wtime();
   }
}