C++ Comm类代码示例

void Teuchos::updateParametersFromXmlFileAndBroadcast(
  const std::string &xmlFileName,
  const Ptr<ParameterList> &paramList,
  const Comm<int> &comm
  )
{
  if (comm.getSize()==1)
    updateParametersFromXmlFile(xmlFileName, paramList);
  else {
    if (comm.getRank()==0) {
      XMLParameterListReader xmlPLReader;
      xmlPLReader.setAllowsDuplicateSublists( false );
      FileInputSource xmlFile(xmlFileName);
      XMLObject xmlParams = xmlFile.getObject();
      std::string xmlString = toString(xmlParams);
      int strsize = xmlString.size();
      broadcast<int, int>(comm, 0, &strsize);
      broadcast<int, char>(comm, 0, strsize, &xmlString[0]);
      updateParametersFromXmlString(xmlString, paramList);
    }
    else {
      int strsize;
      broadcast<int, int>(comm, 0, &strsize);
      std::string xmlString;
      xmlString.resize(strsize);
      broadcast<int, char>(comm, 0, strsize, &xmlString[0]);
      updateParametersFromXmlString(xmlString, paramList);
    }
  }
}

void PotJMEAMSpline::fast_compute_densities(const Comm& comm, Vector<double> *density_ptr)
{
  // Untrap procs if necessary
  if (is_trapped_) {
    int flag = 4;
    comm.bcast(&flag, 1, MPI_INT, comm.get_root());
  }

  // Initialize potential on all procs
  initialize_pot(comm);

  // Setup local variables for potential function
  std::vector<T *> rho_fns;
  for (Basis*& fn : pot_fns_[1].fns)
    rho_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> f_fns;
  for (Basis*& fn : pot_fns_[3].fns)
    f_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> g_fns;
  for (Basis*& fn : pot_fns_[4].fns)
    g_fns.push_back(static_cast<T *>(fn));

  // Make list of densities for each atom
  int natoms = mmz->config->total_natoms;
  Vector<double> densities(natoms, 0.0);

  // Loop over all atoms in atomvec
  for (Atom*& atom_i_ptr : mmz->atomvec->atoms) {
    AtomJMEAMSpline &atom_i = *(static_cast<AtomJMEAMSpline *>(atom_i_ptr));  // tmp atom

    for (Pair*& pair_ij_ptr : atom_i.pairs) {
      PairJMEAMSpline &pair_ij = *(static_cast<PairJMEAMSpline *>(pair_ij_ptr));  // tmp pair

      if (pair_ij.rho_knot != -1)  // pair distance is inside density function
        densities[atom_i.global_idx] += rho_fns[pair_ij.rho_idx]->T::splint(pair_ij.rho_knot, pair_ij.rho_shift);

      if (pair_ij.f_knot != -1)  // radial distance inside f-potential
        pair_ij.f = f_fns[pair_ij.f_idx]->T::splint(pair_ij.f_knot, pair_ij.f_shift);
      else
        pair_ij.f = 0.0;
    } // END LOOP OVER PAIRS

    for (Triplet*& triplet_ijk_ptr : atom_i.triplets) {
      TripletJMEAMSpline &triplet_ijk = *(static_cast<TripletJMEAMSpline *>(triplet_ijk_ptr));  // tmp triplet
      PairJMEAMSpline &pair_ij = *(static_cast<PairJMEAMSpline *>(triplet_ijk.pair_ij));  // tmp pairs
      PairJMEAMSpline &pair_ik = *(static_cast<PairJMEAMSpline *>(triplet_ijk.pair_ik));

      double g_val = g_fns[triplet_ijk.g_idx]->T::splint(triplet_ijk.g_knot, triplet_ijk.g_shift);

      densities[atom_i.global_idx] += pair_ij.f * pair_ik.f * g_val;
    } // END LOOP OVER TRIPLETS
  } // END LOOP OVER ATOMS


  // Gather up densities from all procs
  std::vector<double> densities_final(natoms, 0.0);
  comm.reduce(&densities[0], &densities_final[0], natoms, MPI_DOUBLE, MPI_SUM, comm.get_root());

  if (comm.is_root()) density_ptr->swap(densities_final);
}

void s4u::Actor::send(Mailbox &chan, void *payload, size_t simulatedSize) {
  Comm c = Comm::send_init(this,chan);
  c.setRemains(simulatedSize);
  c.setSrcData(payload);
  // c.start() is optional.
  c.wait();
}

        void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
        {
            Dune::InverseOperatorResult result;
            // Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
            if (parallelInformation_.type() == typeid(ParallelISTLInformation))
            {
                const size_t size = opA.getmat().N();
                const ParallelISTLInformation& info =
                    boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);

                // As we use a dune-istl with block size np the number of components
                // per parallel is only one.
                info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
                                  size, 1);
                // Construct operator, scalar product and vectors needed.
                constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
            }
            else
#endif
            {
                OPM_THROW(std::logic_error,"this method if for parallel solve only");
            }

            checkConvergence( result );
        }

bool CommConfigurable::run(Messages & messages, Feedback & feedback) {
  // Run all the subclasses
  bool result = false;
  std::vector<Comm *>::iterator iter;
  for (iter = comm.begin(); iter != comm.end(); iter++) {
    Comm *subclass = (Comm *)(*iter);
    if (subclass->run(messages, feedback)) {
      result = true;
    }
  }

  return result;
}

void
mergeCounterNames (const Comm<int>& comm,
                   const Array<std::string>& localNames,
                   Array<std::string>& globalNames,
                   const ECounterSetOp setOp)
{
    const int myRank = comm.getRank();
    const int left = 0;
    const int right = comm.getSize() - 1;
    Array<std::string> theGlobalNames;
    mergeCounterNamesHelper (comm, myRank, left, right,
                             localNames, theGlobalNames, setOp);

    // Proc 0 has the list of counter names.  Now broadcast it back to
    // all the procs.
    broadcastStrings (comm, theGlobalNames);

    // "Transactional" semantics ensure strong exception safety for
    // output.
    globalNames.swap (theGlobalNames);
}

  MachineRepresentation(const Comm<int> &comm):
    networkDim(0), numProcs(comm.getSize()), myRank(comm.getRank()),
    procCoords(NULL)
  {
    // WIll need this constructor to be specific to RAAMP (MD).
    // Will need a default constructor using, e.g., GeometricGenerator
    // or nothing at all, for when RAAMP is not available as TPL.
    //
    // (AG) In addition, need to be able to run without special
    // privileges in system (e.g., on hopper).
    // Notes:  For now, all cores connected to same NIC will get the
    // same coordinates; later, we could add extra coordinate dimensions
    // to represent nodes or dies (using hwloc info through RAAMP
    // data object).

    // (MD) will modify mapping test to use machine representation
    // #ifdef HAVE_ZOLTAN2_OVIS

    // Call initializer for RAAMP data object (AG)

    //get network dimension.
    //TODO change.
    // Call RAAMP Data Object to get the network dimension (AG)
    networkDim = 3;

    //allocate memory for processor coordinates.
    procCoords = new nCoord_t *[networkDim];
    for (int i = 0; i < networkDim; ++i){
      procCoords[i] = new nCoord_t [numProcs];
      memset (procCoords[i], 0, sizeof(nCoord_t) * numProcs);
    }
    //obtain the coordinate of the processor.
    this->getMyCoordinate(/*nCoord_t &xyz[networkDim]*/);
    // copy xyz into appropriate spot in procCoords. (MD)  // KDD I agree with this

    //reduceAll the coordinates of each processor.
    this->gatherMachineCoordinates();
  }

int PotGMEAM::rescale_3body(const Comm& comm, std::ostream *out, int flag)
{
  if (!flag) return 0; // Don't run rescaling if flag is 0

  int ntypes = mmz->potlist->get_ntypes();

  int f_idx = 0;
  for (Basis*& f_fn : pot_fns_[3].fns) {
    double max_f_mag = f_fn->get_max_y_mag();
    double b = 1.0/max_f_mag;

    // Scale f-pot
    *f_fn *= b;

    // Scale g-pot
    std::vector<Basis *> g_fns = pot_fns_[4].fns;
    for (int i=0; i<ntypes; ++i) {
      for (int j=0; j<ntypes; ++j) {
        for (int k=j; k<ntypes; ++k) {
          int ij_idx = pot_fns_[3].get_2body_alloy_idx(i, j);
          int ik_idx = pot_fns_[3].get_2body_alloy_idx(i, k);
          int ijk_idx = pot_fns_[4].get_3body_alloy_idx(i, j, k);

          if (f_idx == ij_idx)
            *g_fns[ijk_idx] /= b;
          if (f_idx == ik_idx)
            *g_fns[ijk_idx] /= b;
        }
      }
    }

    // Output scaling factor to screen
    if (comm.is_root() && out)
      *out << "MEAM potential scaling factor (b_" << f_idx << ") " << std::fixed << b << std::endl;

    ++f_idx;
  }

  return 1;
}

        void prepareSolver(Operator& wellOpA, Comm& comm)
        {

            Vector& istlb = *(this->rhs_);
            comm.copyOwnerToAll(istlb, istlb);

            const double relax = this->parameters_.ilu_relaxation_;
            const MILU_VARIANT ilu_milu  = this->parameters_.ilu_milu_;

            // TODO: revise choice of parameters
            // int coarsenTarget = 4000;
            int coarsenTarget = 1200;
            Criterion criterion(15, coarsenTarget);
            criterion.setDebugLevel( this->parameters_.cpr_solver_verbose_ ); // no debug information, 1 for printing hierarchy information
            criterion.setDefaultValuesIsotropic(2);
            criterion.setNoPostSmoothSteps( 1 );
            criterion.setNoPreSmoothSteps( 1 );
            //new guesses by hmbn
            //criterion.setAlpha(0.01); // criterion for connection strong 1/3 is default
            //criterion.setMaxLevel(2); //
            //criterion.setGamma(1); //  //1 V cycle 2 WW

            // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
            using AmgType           = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
                                                                BlackoilAmgType, ParallelBlackoilAmgType>::type;
            using SpType            = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
                                                                Dune::SeqScalarProduct<Vector>,
                                                                ParallelScalarProduct >::type;
            using OperatorType      = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
                                                                MatrixAdapter, ParallelMatrixAdapter>::type;
            typedef typename AmgType::Smoother Smoother;
            typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments  SmootherArgs;
            SmootherArgs  smootherArgs;
            smootherArgs.iterations = 1;
            smootherArgs.relaxationFactor = relax;
            const Opm::CPRParameter& params(this->parameters_); // strange conversion
            ISTLUtility::setILUParameters(smootherArgs, ilu_milu);

            auto& opARef = reinterpret_cast<OperatorType&>(*opA_);
            int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
            bool update_preconditioner = false;

            if (this->parameters_.cpr_reuse_setup_ < 1) {
                update_preconditioner = true;
            }
            if (this->parameters_.cpr_reuse_setup_ < 2) {
                if (newton_iteration < 1) {
                    update_preconditioner = true;
                }
            }
            if (this->parameters_.cpr_reuse_setup_ < 3) {
                if (this->iterations() > 10) {
                    update_preconditioner = true;
                }
            }

            if ( update_preconditioner or (amg_== 0) ) {
                amg_.reset( new AmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
            } else {
                if (this->parameters_.cpr_solver_verbose_) {
                    std::cout << " Only update amg solver " << std::endl;
                }
                reinterpret_cast<AmgType*>(amg_.get())->updatePreconditioner(opARef, smootherArgs, comm);
            }
            // Solve.
            //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result);
            //references seems to do something els than refering

            int verbosity_linsolve = 0;
            if (comm.communicator().rank() == 0) {
                verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
            }

            linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, reinterpret_cast<SpType&>(*sp_), reinterpret_cast<AmgType&>(*amg_),
                                                             this->parameters_.linear_solver_reduction_,
                                                             this->parameters_.linear_solver_maxiter_,
                                                             verbosity_linsolve));
        }

int main(int argc, char** argv)
{
  In in;
  in.datafile = NULL;
  int me = 0;                   //local MPI rank
  int nprocs = 1;               //number of MPI ranks
  int num_threads = 1;		//number of OpenMP threads
  int num_steps = -1;           //number of timesteps (if -1 use value from lj.in)
  int system_size = -1;         //size of the system (if -1 use value from lj.in)
  int nx = -1;
  int ny = -1;
  int nz = -1;
  int check_safeexchange = 0;   //if 1 complain if atom moves further than 1 subdomain length between exchanges
  int do_safeexchange = 0;      //if 1 use safe exchange mode [allows exchange over multiple subdomains]
  int use_sse = 0;              //setting for SSE variant of miniMD only
  int screen_yaml = 0;          //print yaml output to screen also
  int yaml_output = 0;          //print yaml output
  int halfneigh = 1;            //1: use half neighborlist; 0: use full neighborlist; -1: use original miniMD version half neighborlist force
  int teams = 1;
  int device = 0;
  int neighbor_size = -1;
  char* input_file = NULL;
  int ghost_newton = 1;
  int sort = -1;

  for(int i = 0; i < argc; i++) {
    if((strcmp(argv[i], "-i") == 0) || (strcmp(argv[i], "--input_file") == 0)) {
      input_file = argv[++i];
      continue;
    }
  }

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &me);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

  int error = 0;

  if(input_file == NULL)
    error = input(in, "in.lj.miniMD");
  else
    error = input(in, input_file);

  if(error) {
    MPI_Finalize();
    exit(0);
  }

  for(int i = 0; i < argc; i++) {
    if((strcmp(argv[i], "-t") == 0) || (strcmp(argv[i], "--num_threads") == 0)) {
      num_threads = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "--teams") == 0)) {
      teams = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0))  {
      num_steps = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-s") == 0) || (strcmp(argv[i], "--size") == 0)) {
      system_size = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-nx") == 0)) {
      nx = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-ny") == 0)) {
      ny = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-nz") == 0)) {
      nz = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-b") == 0) || (strcmp(argv[i], "--neigh_bins") == 0))  {
      neighbor_size = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "--half_neigh") == 0))  {
      halfneigh = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "-sse") == 0))  {
      use_sse = atoi(argv[++i]);
      continue;
    }

    if((strcmp(argv[i], "--check_exchange") == 0))  {
//.........这里部分代码省略.........

double PotJMEAMSpline::fast_compute(const Comm& comm, ErrorVec *error_vec)
{
  // Untrap procs if necessary
  if (is_trapped_) {
    if (error_vec) {
      int flag = 3;
      comm.bcast(&flag, 1, MPI_INT, comm.get_root());
    } else {
      int flag = 2;
      comm.bcast(&flag, 1, MPI_INT, comm.get_root());
    }
  }

  // Initialize potential on all procs
  initialize_pot(comm, error_vec);

  // Initialize potential by resetting forces
  initialize_compute(comm);

  // Setup local variables for potential functions
  std::vector<T *> phi_fns;
  for (Basis*& fn : pot_fns_[0].fns)
    phi_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> rho_fns;
  for (Basis*& fn : pot_fns_[1].fns)
    rho_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> F_fns;
  for (Basis*& fn : pot_fns_[2].fns)
    F_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> f_fns;
  for (Basis*& fn : pot_fns_[3].fns)
    f_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> g_fns;
  for (Basis*& fn : pot_fns_[4].fns)
    g_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> p_fns;
  for (Basis*& fn : pot_fns_[5].fns)
    p_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> q_fns;
  for (Basis*& fn : pot_fns_[6].fns)
    q_fns.push_back(static_cast<T *>(fn));

  // Set up constraint error (error from density going out of bounds of embedding function)
  Vector<double> constraint_err(mmz->config->ncells,0.0);

  // Loop over all atoms in atomvec
  for (Atom*& atom_i_ptr : mmz->atomvec->atoms) {
    // Make temporary atom and cell
    AtomJMEAMSpline &atom_i = *(static_cast<AtomJMEAMSpline *>(atom_i_ptr));
    Cell &cell = *mmz->config->cells[atom_i.cell_idx];

    double rho_val = 0.0; // initialize density for this atom
    double dF = 0.0;      // initialize gradient of embedding fn for this atom

    // Loop over pairs for this atom
    for (Pair*& pair_ij_ptr : atom_i.pairs) {
      PairJMEAMSpline &pair_ij = *(static_cast<PairJMEAMSpline *>(pair_ij_ptr));  // tmp pair

      // Check that neighbor length lies in pair potential radius
      if (pair_ij.phi_knot != -1) {
        AtomJMEAMSpline &atom_j = *(static_cast<AtomJMEAMSpline *>(pair_ij.neigh));  // tmp atom

        // Compute phi(r_ij) and its gradient in one step
        double phigrad;
        double phival = 0.5 * phi_fns[pair_ij.phi_idx]->T::splint_comb(pair_ij.phi_knot, pair_ij.phi_shift, &phigrad);

        phigrad *= 0.5; // only half of the gradient/energy contributes to the force/energy since we are double counting

        cell.energy += phival;  // add in piece contributed by neighbor to energy

        Vect tmp_force = pair_ij.dist * phigrad;  // compute tmp force values
        atom_i.force += tmp_force;  // add in force on atom i from atom j
        atom_j.force -= tmp_force;  // subtract off force on atom j from atom i (Newton's law: action = -reaction)

        // Compute stress on cell
        tmp_force *= pair_ij.r;
        cell.stress -= pair_ij.dist & tmp_force;
      } // END IF STMNT: PAIR LIES INSIDE CUTOFF FOR PAIR POTENTIAL

      // Check that neighbor length lies in rho potential (density function) radius
      if (pair_ij.rho_knot != -1) {
        // Compute density and its gradient in one step
        rho_val += rho_fns[pair_ij.rho_idx]->T::splint_comb(pair_ij.rho_knot, pair_ij.rho_shift, &pair_ij.drho);
      } else {
        pair_ij.drho = 0.0;
      } // END IF STMNT: PAIR LIES INSIDE CUTOFF FOR RHO POTENTIAL

      // Check that neighbor length lies in f- potential radius
      if (pair_ij.f_knot != -1) {
        pair_ij.f = f_fns[pair_ij.f_idx]->T::splint_comb(pair_ij.f_knot, pair_ij.f_shift, &pair_ij.df);
      } else {
        pair_ij.f = 0.0;
        pair_ij.df = 0.0;
      } // END IF STMNT: PAIR LIES INSIDE CUTOFF FOR f- POTENTIAL

      // Check that neighbor length lies in p- potential radius
      if (pair_ij.p_knot != -1) {
        pair_ij.p = p_fns[pair_ij.p_idx]->T::splint_comb(pair_ij.p_knot, pair_ij.p_shift, &pair_ij.dp);
      } else {
        pair_ij.p = 0.0;
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
  //Common miniMD settings
  In in;
  in.datafile = NULL;
  int me=0;                     //local MPI rank
  int nprocs=1;                 //number of MPI ranks
  int num_threads=32;		    //number of Threads per Block threads
  int num_steps=-1;             //number of timesteps (if -1 use value from lj.in)
  int system_size=-1;           //size of the system (if -1 use value from lj.in)
  int check_safeexchange=0;     //if 1 complain if atom moves further than 1 subdomain length between exchanges
  int do_safeexchange=0;        //if 1 use safe exchange mode [allows exchange over multiple subdomains]
  int use_sse=0;                //setting for SSE variant of miniMD only
  int screen_yaml=0;            //print yaml output to screen also
  int yaml_output=0;            //print yaml output
  int halfneigh=0;              //1: use half neighborlist; 0: use full neighborlist; -1: use original miniMD version half neighborlist force
  char* input_file = NULL;
  int ghost_newton = 0;
  int skip_gpu = 999;
  int neighbor_size = -1;

  //OpenCL specific
  int platform = 0;
  int device = 0;
  int subdevice = -1;
  int ppn = 2;
  int use_tex = 0;
  int threads_per_atom = 1;
  int map_device=0;

  for(int i = 0; i < argc; i++) {
    if((strcmp(argv[i], "-i") == 0) || (strcmp(argv[i], "--input_file") == 0)) {
      input_file = argv[++i];
      continue;
    }
    if((strcmp(argv[i],"-p")==0)||(strcmp(argv[i],"--platform")==0)) {platform=atoi(argv[++i]); continue;}
    if((strcmp(argv[i],"-d")==0)||(strcmp(argv[i],"--device")==0)) {device=atoi(argv[++i]); continue;}
	if((strcmp(argv[i],"-sd")==0)||(strcmp(argv[i],"--subdevice")==0)) {subdevice=atoi(argv[++i]); continue;}
	if((strcmp(argv[i],"-sd_map")==0)||(strcmp(argv[i],"--subdevice_mapping")==0)) {subdevice=1-me%ppn; continue;}
	if((strcmp(argv[i],"-ng")==0)||(strcmp(argv[i],"--num_gpus")==0)) {ppn=atoi(argv[++i]); continue;}
	if((strcmp(argv[i],"-dm")==0)||(strcmp(argv[i],"--device_map")==0)) {map_device=1; continue;}
  }

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &me);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

  if(map_device) {device = me%ppn; if(device>=skip_gpu) device++;}

  OpenCLWrapper* opencl = new OpenCLWrapper;
  if( me == 0)
  printf("# Platforms: %i\n",opencl->num_platforms);
  printf("# Proc: %i using device %i\n",me,device);
  opencl->Init(argc,argv,device,device+1,NULL,platform,subdevice);

  int error = 0;
  if(input_file == NULL)
    error = input(in, "in.lj.miniMD");
  else
	error = input(in, input_file);

  if (error)
  {
	  MPI_Finalize();
	  exit(0);
  }

  for(int i=0;i<argc;i++)
  {
     if((strcmp(argv[i],"-t")==0)||(strcmp(argv[i],"--num_threads")==0)) {num_threads=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"-n")==0)||(strcmp(argv[i],"--nsteps")==0))  {num_steps=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"-s")==0)||(strcmp(argv[i],"--size")==0))  {system_size=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"--half_neigh")==0))  {halfneigh=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"-sse")==0))  {use_sse=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"--check_exchange")==0))  {check_safeexchange=1; continue;}
     if((strcmp(argv[i],"-o")==0)||(strcmp(argv[i],"--yaml_output")==0))  {yaml_output=atoi(argv[++i]); continue;}
     if((strcmp(argv[i],"--yaml_screen")==0))  {screen_yaml=1; continue;}
     if((strcmp(argv[i], "-f") == 0) || (strcmp(argv[i], "--data_file") == 0)) {
       if(in.datafile == NULL) in.datafile = new char[1000];

       strcpy(in.datafile, argv[++i]);
       continue;
     }
     if((strcmp(argv[i], "-u") == 0) || (strcmp(argv[i], "--units") == 0)) {
       in.units = strcmp(argv[++i], "metal") == 0 ? 1 : 0;
       continue;
     }

     if((strcmp(argv[i], "-p") == 0) || (strcmp(argv[i], "--force") == 0)) {
       in.forcetype = strcmp(argv[++i], "eam") == 0 ? FORCEEAM : FORCELJ;
       continue;
     }
     if((strcmp(argv[i], "-gn") == 0) || (strcmp(argv[i], "--ghost_newton") == 0)) {
       ghost_newton = atoi(argv[++i]);
       continue;
     }
     if((strcmp(argv[i], "--skip_gpu") == 0)) {
       skip_gpu = atoi(argv[++i]);
       continue;
     }
//.........这里部分代码省略.........

int main (int argc, char **argv)
{
  CGMTimers *timers = new CGMTimers ();
  Comm *comm = Comm::getComm (&argc, &argv, timers);
  int tag = 0;
  int p = comm->getNumberOfProcessors ();
  //printf("Number of processors %d\n",p);
  int id = comm->getMyId ();
  //printf("ID%d\n",id);
  int sendTo = (id + 1) % p;
  int receiveFrom = (id - 1 + p) % p;
  int actualSource = -1;
  vector<int> values;
  if (argc > 1)
  {
	char *size=strtok(argv[1],",");
	while (size!=NULL)
	{
		values.push_back(atoi(size));
		size=strtok(NULL,",");
	}
  }
  else
  {
	printf("Falt argumento com os tamanhos das matrices: ex. 1,3,4,5,6\n");
  }
//	printf("size %d\n",values.size());
//	int bloco=(values.size()-1)/p;
	//printf("Bloco %d\n",bloco);
		SimpleCommObject<int> sample(0);
	int matrix_size=values.size()-1;
	int **total_matrix=new int*[matrix_size];
	for (int i=0;i<matrix_size;i++)
		total_matrix[i]=new int[matrix_size];
	for (int row=0;row<matrix_size;row++)
		for (int col=0;col<matrix_size;col++)
			{
				if (row==col)
					total_matrix[row][row]=0;
				else
					total_matrix[row][col]=-1;
			}
	int *blocos=new int[p];
	int q=(values.size()-1)/p;
	if (q<1) q=1;
	int r=(values.size()-1)%p;
	if (r<1) r=0;
	for (int i=0;i<p;i++)
	{
		if (i<r)
			blocos[i]=q+1;
		else
			blocos[i]=q;
	}
	//printf("ID: %d Bloco: %d\n",id,blocos[id]);
	int bloco_offset=0;
	for (int i=0;i<id;i++)
		bloco_offset+=blocos[i];
	//int row_start=blocos[id]*id;
	int row_start=bloco_offset;
	int row_end=row_start+blocos[id]-1;
	//printf("ID %d, row_start %d row_end %d\n",id,row_start, row_end);
	for (int rodada=0;rodada<=p-id-1;rodada++)
	{
		//printf("ID %d RODADA %d\n",id,rodada);
		//int col_start=blocos[id]*(rodada+id);
		int col_start=bloco_offset+(blocos[id]*rodada);
		//int col_end=blocos[id]*(rodada+id+1)-1;
		int col_end=col_start+blocos[id]-1;
		if (col_start>values.size()-2)
			break;
	//	printf("ID %d, col_start %d, col_end %d RODADA: %d\n",id,col_start, col_end,rodada);
		CommObjectList data_to_send(&sample);

		workOnSubMatrix(&total_matrix, &values, row_start,row_end,col_start,col_end, rodada, blocos[id], id);
		//printf("ID: %d SAIU RODADA: %d\n",id,rodada);
		int **submatrix=new int*[row_end-row_start+1];
		for (int i=0;i<row_end-row_start+1;i++)
			submatrix[i]=new int[col_end-col_start+1];
/*if (id==1)
{		
	for (int row=0;row<matrix_size;row++)
	{
		for (int col=0;col<matrix_size;col++)
			{
				printf("R>%d %d,%d=%d ",rodada,row,col,total_matrix[row][col]);
			}
		printf("\n");
	}
}*/

		convertMatrixToList(&total_matrix,row_start, col_end, col_start, col_end, &data_to_send,id,rodada);
		//printf("ID: %d, COPIOU SUBM RODADA: %d\n",id,rodada);

		CommObjectList data_to_receive(&sample);


		if (id!=0)
		{
			//printf("ID: %d IS SENDING RODADA: %d\n",id,rodada);
//.........这里部分代码省略.........

double PotEAMSpline::fast_compute(const Comm& comm, ErrorVec *error_vec)
{
  // Untrap procs if necessary
  if (is_trapped_) {
    if (error_vec) {
      int flag = 3;
      comm.bcast(&flag, 1, MPI_INT, comm.get_root());
    } else {
      int flag = 2;
      comm.bcast(&flag, 1, MPI_INT, comm.get_root());
    }
  }

  // Initialize potential on all procs
  initialize_pot(comm, error_vec);

  // Initialize potential by resetting forces
  initialize_compute(comm);

  // Setup local variables for potential functions
  std::vector<T *> phi_fns;
  for (Basis*& fn : pot_fns_[0].fns)
    phi_fns.push_back(static_cast<T *>(fn));
  std::vector<T *> rho_fns;
  for (Basis*& fn : pot_fns_[1].fns)
    rho_fns.push_back(static_cast<T *>(fn));
  std::vector<Basis *> F_fns;
  for (Basis*& fn : pot_fns_[2].fns) // allow embedding fn to be any basis
    F_fns.push_back(fn);

  // Set up constraint error (error from density going out of bounds of embedding function)
  Vector<double> constraint_err(mmz->config->ncells,0.0);

  // Loop over all atoms in atomvec
  for (Atom*& atom_i_ptr : mmz->atomvec->atoms) {
    // Make temporary atom and cell
    AtomEAMSpline &atom_i = *(static_cast<AtomEAMSpline *>(atom_i_ptr));
    Cell &cell = *mmz->config->cells[atom_i.cell_idx];

    double rho_val = 0.0; // initialize density for this atom
    double dF = 0.0;      // initialize gradient of embedding fn for this atom

    // Loop over pairs for this atom
    for (Pair*& pair_ij_ptr : atom_i.pairs) {
      PairEAMSpline &pair_ij = *(static_cast<PairEAMSpline *>(pair_ij_ptr));  // tmp pair

      // Check that neighbor length lies in pair potential radius
      if (pair_ij.phi_knot != -1) {
        AtomEAMSpline &atom_j = *(static_cast<AtomEAMSpline *>(pair_ij.neigh));  // tmp atom

        // Compute phi(r_ij) and its gradient in one step
        double phigrad;
        double phival = 0.5 * phi_fns[pair_ij.phi_idx]->T::splint_comb(pair_ij.phi_knot, pair_ij.phi_shift, &phigrad);

        phigrad *= 0.5; // only half of the gradient/energy contributes to the force/energy since we are double counting

        cell.energy += phival;  // add in piece contributed by neighbor to energy

        Vect tmp_force = pair_ij.dist * phigrad;  // compute tmp force values
        atom_i.force += tmp_force;  // add in force on atom i from atom j
        atom_j.force -= tmp_force;  // subtract off force on atom j from atom i (Newton's law: action = -reaction)

        // Compute stress on cell
        tmp_force *= pair_ij.r;
        cell.stress -= pair_ij.dist & tmp_force;
      } // END IF STMNT: PAIR LIES INSIDE CUTOFF FOR PAIR POTENTIAL

      // Check that neighbor length lies in rho potential (density function) radius
      if (pair_ij.rho_knot != -1) {
        // Compute density and its gradient in one step
        rho_val += rho_fns[pair_ij.rho_idx]->T::splint_comb(pair_ij.rho_knot, pair_ij.rho_shift, &pair_ij.drho);
      } else {
        pair_ij.drho = 0.0;
      } // END IF STMNT: PAIR LIES INSIDE CUTOFF FOR RHO POTENTIAL
    } // END LOOP OVER PAIRS

    // Compute energy, gradient for embedding function F
    // Punish this potential for having rho lie outside of F
    if ( rho_val < F_fns[atom_i.F_idx]->get_min_rcut() ) {
      double rho_i = F_fns[atom_i.F_idx]->get_min_rcut();
      constraint_err[atom_i.cell_idx] += cell.weight * DUMMY_WEIGHT * 10. * (rho_i - rho_val) * (rho_i - rho_val);
      if (!embed_extrap_)
        rho_val = rho_i;  // set the density to the inner cutoff if we don't extrapolate embedding fn later
    } else if ( rho_val > F_fns[atom_i.F_idx]->get_max_rcut() ) {
      double rho_f = F_fns[atom_i.F_idx]->get_max_rcut();
      constraint_err[atom_i.cell_idx] += cell.weight * DUMMY_WEIGHT * 10. * (rho_val - rho_f) * (rho_val - rho_f);
      if (!embed_extrap_)
        rho_val = rho_f;  // set the density to the outer cutoff if we don't extrapolate embedding fn later
    }

    // Add energy contribution from embedding function and get gradient in one step
    cell.energy += F_fns[atom_i.F_idx]->eval_comb(rho_val, &dF);

    // Loop over pairs for this atom to compute EAM force
    for (Pair*& pair_ij_ptr : atom_i.pairs) {
      PairEAMSpline &pair_ij = *(static_cast<PairEAMSpline *>(pair_ij_ptr));  // tmp pair
      AtomEAMSpline &atom_j = *(static_cast<AtomEAMSpline *>(pair_ij.neigh));  // tmp atom

      Vect tmp_force = pair_ij.dist * pair_ij.drho * dF;  // compute tmp force values
      atom_i.force += tmp_force;  // add in force on atom i from atom j
//.........这里部分代码省略.........

int main() {
    controls.setPC(comm.getPC());
    controls.setup();
//    comm.printPosition();
    comm.printGains();

    controlsInterrupt.attach_us(&controls, &Controls::loop, 1000);

    while(1) {
        controls.updateIMUS();
        comm.check();
        
        
        if (serialCounter++>100) {
//            comm.printPosition();
//            comm.getPC()->printf("%f\n", controls.getTheta1());
//            comm.getPC()->printf("%f", controls.motor.getPWM());
            serialCounter = 0;
           // float z[4] = {1,2,0,0};
//            comm.getPC()->printf("%f\n",controls.target.getTheta2ForTarget(z));
        }
    }
}