<< " [-d step_size_control_parameter]\n"                           << "         [-e diagnostics_interval]"                           << " [-o output_interval]\n"                           << "         [-t total_duration]"                           << " [-i (start output at t = 0)]\n"                           << "         [-x (extra debugging diagnostics)]"                           << endl;                      return false;         // execution should stop after help            case 'd': dt_param = atof(optarg);                      break;            case 'e': dt_dia = atof(optarg);                      break;            case 'i': i_flag = true;                      break;            case 'o': dt_out = atof(optarg);                      break;            case 't': dt_tot = atof(optarg);                      break;            case 'x': x_flag = true;                      break;            case '?': cerr << "usage: " << argv[0]                           << " [-h (for help)]"                           << " [-d step_size_control_parameter]\n"                           << "         [-e diagnostics_interval]"                           << " [-o output_interval]\n"                           << "         [-t total_duration]"                           << " [-i (start output at t = 0)]\n"                           << "         [-x (extra debugging diagnostics)]"                           << endl;                      return false;        // execution should stop after error            }    return true;                         // ready to continue program execution}/*----------------------------------------------------------------------------- *  get_snapshot  --  reads a single snapshot with the root processor  *                    from the input stream cin. * *  note: in this implementation, particle number, time and the data for *        individual particles are read in. *        All communication goes via the root (0) process. *  note: The routine also declared the MPI_Datatype particletype which is *        the contianer for the particle data. *        This structure is used for further MPI comminication. *----------------------------------------------------------------------------- */Particle *get_snapshot(int &n, real &t,		       MPI_Datatype &particletype){  int rank, size;  MPI_Comm_rank( MPI_COMM_WORLD, &rank );  MPI_Comm_size( MPI_COMM_WORLD, &size );  Particle *p_tmp;  if(rank==root) {    cerr << "Reading snapshot" << endl;    cin >> n;    cin >> t;    p_tmp = new Particle[n];    for(int i=0; i<n; i++) {      p_tmp[i].id = i;      cin >> p_tmp[i].mass;                         // mass of particle i      for (int k = 0; k < NDIM; k++)	cin >> p_tmp[i].pos[k];                 // position of particle i      for (int k = 0; k < NDIM; k++)	cin >> p_tmp[i].vel[k];                 // velocity of particle i    }  }  MPI_Bcast(&n,1,MPI_INT,root,MPI_COMM_WORLD); // broadcasts particle number  int n_local = (int)(floor(1.0*n/size));  if(n != n_local*size && rank==root) {    cerr << "WARNING: Paticle number in input is not a mulitple of the number of processors." 	 << endl;    cerr << "         Action: Reduce particle number to n = " 	 << n_local*size << "." << endl;  }  Particle *p = new Particle[n_local];  // defining the particletype  int inputblockcounts[2] = {1, 13};  MPI_Datatype ntypes[] = {MPI::INT, MPI::DOUBLE};  MPI::Aint displs[2];  MPI_Address(&p[0].id, &displs[0]);  MPI_Address(&p[0].mass, &displs[1]);  displs[1] -= displs[0];        //make them relative  displs[0] = 0;  MPI_Type_struct(2, inputblockcounts, displs, ntypes, &particletype);  MPI_Type_commit(&particletype);  // Distribute the particles over the processors  MPI_Scatter(p_tmp,n_local,particletype,p,n_local,particletype,root,MPI_COMM_WORLD);  n = n_local;  MPI_Bcast(&t,1,MPI_DOUBLE,root,MPI_COMM_WORLD);  delete []p_tmp;  return p;}/*----------------------------------------------------------------------------- *  put_snapshot  --  writes a single snapshot on the output stream cout. * *  note: Communication goes via the root (0) proceesor *----------------------------------------------------------------------------- */void put_snapshot(Particle p[], int n, real t,		  MPI_Datatype particletype){    int rank;    MPI_Comm_rank( MPI_COMM_WORLD, &rank );    int ntot;    MPI_Allreduce(&n, &ntot, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );    Particle *p_all;    if(rank==root) {      p_all = new Particle[ntot];    }    MPI_Gather(p,n,particletype,p_all,n,particletype,root,MPI_COMM_WORLD);    cout.precision(16);                       // full double precision    if(rank==root) {      cout << ntot << endl;                     // N, total particle number      cout << t << endl;                        // current time    for (int i = 0; i < ntot ; i++){      cout << p_all[i].mass;                      // mass of particle i        for (int k = 0; k < NDIM; k++)            cout << ' ' << p_all[i].pos[k];         // position of particle i        for (int k = 0; k < NDIM; k++)            cout << ' ' << p_all[i].vel[k];         // velocity of particle i        cout << endl;    }    delete []p_all;    }}    /*----------------------------------------------------------------------------- *  write_diagnostics  --  writes diagnostics on the error stream cerr: *                         current time; number of integration steps so far; *                         kinetic, potential, and total energy; absolute and *                         relative energy errors since the start of the run. *                         If x_flag (x for eXtra data) is true, all internal *                         data are dumped for each particle (mass, position, *                         velocity, acceleration, and jerk). * *  note: the kinetic energy is calculated here, while the potential energy is *        calculated in the function get_acc_jerk_pot_coll() and broadcasted *        by the other processors. * *  note: Communication goes via the root (0) proceesor *----------------------------------------------------------------------------- */void write_diagnostics(Particle p[], int n, real t, real epot_local,                       int nsteps, real & einit, bool init_flag,                       bool x_flag, real &tcpu){    int rank;    MPI_Comm_rank( MPI_COMM_WORLD, &rank );    real ekin_local = 0;       // kinetic energy of the n-body system    for (int i = 0; i < n ; i++)        for (int k = 0; k < NDIM ; k++)            ekin_local += 0.5 * p[i].mass * p[i].vel[k] * p[i].vel[k];    real ekin;    MPI_Allreduce(&ekin_local, &ekin, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );    real epot;    MPI_Allreduce(&epot_local, &epot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );    epot *= 0.5;  // against double counting    real etot = ekin + epot;             // total energy of the n-body system    if (init_flag)                       // at first pass, pass the initial        einit = etot;                    // energy back to the calling function    tcpu = MPI_Wtime() - tcpu;    if(rank==0) {    cerr << "at time t = " << t << ", after " << nsteps         << " steps (CPU = " << tcpu << "): \n  E_kin = " << ekin         << " , E_pot = " << epot         << " , E_tot = " << etot << endl;    cerr << "                "         << "absolute energy error: E_tot - E_init = "         << etot - einit << endl;    cerr << "                "         << "relative energy error: (E_tot - E_init) / E_init = "         << (etot - einit) / einit << endl;    }    if (x_flag){        cerr << "  for debugging purposes, here is the internal data "             << "representation:\n";        for (int i = 0; i < n ; i++){            cerr << "    internal data for particle " << i+1 << " : " << endl;            cerr << "      ";            cerr << p[i].id << " " << p[i].mass;            for (int k = 0; k < NDIM; k++)                cerr << ' ' << p[i].pos[k];            for (int k = 0; k < NDIM; k++)                cerr << ' ' << p[i].vel[k];            for (int k = 0; k < NDIM; k++)                cerr << ' ' << p[i].acc[k];            for (int k = 0; k < NDIM; k++)                cerr << ' ' << p[i].jerk[k];            cerr << endl;        }    }} /*----------------------------------------------------------------------------- *  evolve  --  integrates an N-body system, for a total duration dt_tot. *              Snapshots are sent to the standard output stream once every *              time interval dt_out.  Diagnostics are sent to the standard *              error stream once every time interval dt_dia. * *  note: the integration time step, shared by all particles at any given time, *        is variable.  Before each integration step we use coll_time (short *        for collision time, an estimate of the time scale for any significant *        change in configuration to happen), multiplying it by dt_param (the *        accuracy parameter governing the size of dt in units of coll_time), *        to obtain the new time step size. * *  Before moving any particles, we start with an initial diagnostics output *  and snapshot output if desired.  In order to write the diagnostics, we *  first have to calculate the potential energy, with get_acc_jerk_pot_coll(). *  That function also calculates accelerations, jerks, and an estimate for the *  collision time scale, all of which are needed before we can enter the main *  integration loop below. *       In the main loop, we take as many integration time steps as needed to *  reach the next output time, do the output required, and continue taking *  integration steps and invoking output this way until the final time is *  reached, which triggers a `break' to jump out of the infinite loop set up *  with `while(true)'. *----------------------------------------------------------------------------- */void evolve(Particle p[],             int n, real & t, real dt_param, real dt_dia, real dt_out,            real dt_tot, bool init_out, bool x_flag,	    void *pipe, MPI_Datatype particletype){    int rank, size;    MPI_Comm_rank( MPI_COMM_WORLD, &rank );    MPI_Comm_size( MPI_COMM_WORLD, &size );    if(rank==root)       cerr << "Starting a Hermite integration for a " << n*size	   << "-body system,\n  from time t = " << t 	   << " with time step control parameter dt_param = " << dt_param	   << "  until time " << t + dt_tot 	   << " ,\n  with diagnostics output interval dt_dia = "