C++ PatchMap类代码示例

void Cartographer::fixupResourceMaps(const ResourceType *rt, const Vec2i &pos) {
    PF_TRACE();
    const Map &map = *world->getMap();
    Vec2i junk;
    for (set<ResourceMapKey>::iterator it = resourceMapKeys.begin(); it != resourceMapKeys.end(); ++it) {
        if (it->resourceType == rt) {
            PatchMap<1> *pMap = resourceMaps[*it];
            const int &size = it->workerSize;
            const Field &field = it->workerField;

            Vec2i tl = pos + OrdinalOffsets[odNorthWest] * size;
            Vec2i br(tl.x + size + 2, tl.y + size + 2);

            Util::RectIterator iter(tl, br);
            while (iter.more()) {
                Vec2i cur = iter.next();
                if (map.isInside(cur) && masterMap->canOccupy(cur, size, field)
                        && map.isResourceNear(cur, size, rt, junk)) {
                    pMap->setInfluence(cur, 1);
                } else {
                    pMap->setInfluence(cur, 0);
                }
            }
        }
    }
}

	foreach (KeyList, it, storeMaps) {
		PatchMap<1> *pMap = g_world.getCartographer()->getStoreMap(*it, false);
		if (pMap && pMap->getInfluence(cell)) {
			colour = Vec4f(0.f, 1.f, 0.3f, 0.7f);
			return true;
		}
	}

void Sync::releaseComputes()
{
  PatchMap *patchMap = PatchMap::Object();

  traceUserEvent(eventReleaseComputes);

  for (int i= 0; i<cnum; i++) {
    int &pid = clist[i].pid;
    if (pid == -1) continue;
    if (clist[i].step != step) {
      continue;
    }
    //         CkPrintf(" %d-%d-%d ",
    //	 clist[i].pid, clist[i].step,
    //      patchMap->patch(pid)->flags.sequence);

    ComputePtrListIter cid = clist[i].cid;

    int compute_count = 0;
    for(cid = cid.begin(); cid != cid.end(); cid++) {
      compute_count++;
      (*cid)->patchReady(pid,clist[i].doneMigration,step);
    }
    if (compute_count == 0 && patchMap->node(pid) != CkMyPe()) {
	   iout << iINFO << "PATCH_COUNT-Sync step " << step
		<< "]: Patch " << pid << " on PE " 
		<< CkMyPe() <<" home patch " 
		<< patchMap->node(pid) << " does not have any computes\n" 
		<< endi;
    }
    pid = -1;
  }
//  CkPrintf("\n");
}

bool ResourceMapOverlay::operator()(const Vec2i &cell, Vec4f &colour) {
	ResourceMapKey mapKey(rt, Field::LAND, 1);
	PatchMap<1> *pMap = g_world.getCartographer()->getResourceMap(mapKey);
	if (pMap && pMap->getInfluence(cell)) {
		colour = Vec4f(1.f, 1.f, 0.f, 0.7f);
		return true;
	}
	return false;
}

//! \brief Free up the memory allocated in the PatchMap, making sure to avoid double free's
//! \param map The map to free up
void freePatchMap(PatchMap& map)
{
  std::vector<RealArray*> free;
  for (PatchMap::iterator it = map.begin(); it != map.end(); ++it) {
    for (size_t i = 0; i < it->second.Patch.size(); ++i)
      delete it->second.Patch[i];
    free.insert(free.end(), it->second.FakeModel.begin(), it->second.FakeModel.end());
  }

  std::vector<RealArray*>::iterator it = std::unique(free.begin(), free.end());
  for (std::vector<RealArray*>::iterator it2 = free.begin(); it2 != it; ++it2)
    delete *it2;

  map.clear();
}

 PatchMapIt _add(const int eTag)
 {
   PatchMapIns insPatch = patch.insert(PatchMapVal(eTag, PatchDataListIt()));
   if(insPatch.second) {
     insPatch.first->second =
       patchData.insert(patchData.end(), PatchData(eTag));
   }
   return insPatch.first;
 }

int NamdHybridLB::requiredProxies(PatchID id, int neighborNodes[])
{
  PatchMap* patchMap = PatchMap::Object();
  int myNode = patchMap->node(id);
  int nProxyNodes = 0;

#define IF_NEW_NODE \
    int j; \
    for ( j=0; j<nProxyNodes && neighborNodes[j] != proxyNode; ++j ); \
    if ( j == nProxyNodes )

  PatchID neighbors[1 + PatchMap::MaxOneAway + PatchMap::MaxTwoAway];
  neighbors[0] = id;
  int numNeighbors = 1 + patchMap->downstreamNeighbors(id,neighbors+1);
  for ( int i = 0; i < numNeighbors; ++i ) {
    const int proxyNode = patchMap->basenode(neighbors[i]);
    if ( proxyNode != myNode ) {
      IF_NEW_NODE {
        neighborNodes[nProxyNodes] = proxyNode;
        nProxyNodes++;
      }
    }
  }

void Sync::triggerCompute()
{
  PatchMap *patchMap = PatchMap::Object();
  const int numHomePatches = patchMap->numHomePatches();

  if (numPatches == -1) 
    numPatches = ProxyMgr::Object()->numProxies() + numHomePatches;

// if (CkMyPe()<=8) CkPrintf("SYNC[%d]: PATCHREADY:%d %d patches:%d %d\n", CkMyPe(), counter, numHomePatches, nPatcheReady, numPatches);
//  CkPrintf("SYNC[%d]: PATCHREADY:%d %d patches:%d %d\n", CkMyPe(), counter, PatchMap::Object()->numHomePatches(), nPatcheReady, numPatches);

  if (homeReady == 0 && counter >= numHomePatches) {
    homeReady = 1;
 // if (CkMyPe()<=8) CkPrintf("HOMEREADY[%d]\n", CkMyPe());
    if (!useProxySync)  releaseComputes();
  }

  if (homeReady && nPatcheReady == numPatches)
  {
// if (CkMyPe()<=8) CkPrintf("TRIGGERED[%d]\n", CkMyPe());
//     CkPrintf("TRIGGERED[%d]\n", CkMyPe());
    if (useProxySync) releaseComputes();

    // reset counter
    numPatches = -1;
    step++;
    nPatcheReady = 0;
    for (int i= 0; i<cnum; i++) {
      if (clist[i].pid != -1 && clist[i].step == step) ++nPatcheReady;
    }
    homeReady = 0;
    if ( numHomePatches ) {
      counter -= numHomePatches;
      if (counter >= numHomePatches) triggerCompute();
    }
  }
}

/**
 * @brief Builds the data structures required for the load balancing strategies in NAMD.
 */ 
int NamdHybridLB::buildData(LDStats* stats) {
  int n_pes = stats->nprocs();

  PatchMap* patchMap = PatchMap::Object();
  ComputeMap* computeMap = ComputeMap::Object();
  const SimParameters* simParams = Node::Object()->simParameters;

  BigReal bgfactor = simParams->ldbBackgroundScaling;
  BigReal pmebgfactor = simParams->ldbPMEBackgroundScaling;
  BigReal homebgfactor = simParams->ldbHomeBackgroundScaling;
  int pmeOn = simParams->PMEOn;
  int unLoadPme = simParams->ldbUnloadPME;
  int pmeBarrier = simParams->PMEBarrier;
  int unLoadZero = simParams->ldbUnloadZero;
  int unLoadOne = simParams->ldbUnloadOne;
  int unLoadIO= simParams->ldbUnloadOutputPEs;
  // traversing the list of processors and getting their load information
  int i, pe_no;
  for (i=0; i<n_pes; ++i) {
    pe_no = stats->procs[i].pe;

    // BACKUP processorArray[i].Id = i; 
    processorArray[i].Id = pe_no;               // absolute pe number
    processorArray[i].available = true;
    // BACKUP if ( pmeOn && isPmeProcessor(i) )
    if ( pmeOn && isPmeProcessor(pe_no) ) {
      processorArray[i].backgroundLoad = pmebgfactor * stats->procs[i].bg_walltime;
    // BACKUP } else if (patchMap->numPatchesOnNode(i) > 0) {
    } else if (patchMap->numPatchesOnNode(pe_no) > 0) {
      processorArray[i].backgroundLoad = homebgfactor * stats->procs[i].bg_walltime;
    } else {
      processorArray[i].backgroundLoad = bgfactor * stats->procs[i].bg_walltime;
    }
    processorArray[i].idleTime = stats->procs[i].idletime;
    processorArray[i].load = processorArray[i].computeLoad = 0.0;
  }

  // If I am group zero, then offload processor 0 and 1 in my group
  if(stats->procs[0].pe == 0) {
    if(unLoadZero) processorArray[0].available = false;
    if(unLoadOne) processorArray[1].available = false;
  }

  // if all pes are Pme, disable this flag
  if (pmeOn && unLoadPme) {
    for (i=0; i<n_pes; i++) {
      if(!isPmeProcessor(stats->procs[i].pe))  break;
    }
    if (i == n_pes) {
      iout << iINFO << "Turned off unLoadPme flag!\n"  << endi;
      unLoadPme = 0;
    }
  }

  if (pmeOn && unLoadPme) {
    for (i=0; i<n_pes; i++) {
      if ((pmeBarrier && i==0) || isPmeProcessor(stats->procs[i].pe)) 
	processorArray[i].available = false;
    }
  }

  // if all pes are output, disable this flag
#ifdef MEM_OPT_VERSION
  if (unLoadIO) {
      if (simParams->numoutputprocs == n_pes) {
	  iout << iINFO << "Turned off unLoadIO flag!\n"  << endi;
	  unLoadIO = 0;
      }
  }
  if (unLoadIO){
      for (i=0; i<n_pes; i++) {
	  if (isOutputProcessor(stats->procs[i].pe)) 
	      processorArray[i].available = false;
      }
  }
#endif

  // need to go over all patches to get all required proxies
  int numPatches = patchMap->numPatches();
  int totalLocalProxies = 0;
  int totalProxies = 0;
  for ( int pid=0; pid<numPatches; ++pid ) {
	int neighborNodes[PatchMap::MaxOneAway + PatchMap::MaxTwoAway];

	patchArray[pid].Id = pid;
	patchArray[pid].numAtoms = 0;
	patchArray[pid].processor = patchMap->node(pid);

	const int numProxies = 
#if 0 // USE_TOPOMAP - this function needs to be there for the hybrid case
	requiredProxiesOnProcGrid(pid,neighborNodes);
#else
	requiredProxies(pid, neighborNodes);
#endif

        int numLocalProxies = 0;
	for (int k=0; k<numProxies; k++) {
//.........这里部分代码省略.........

void OrbLB::mapPartitionsToNodes()
{
  int i,j;
#if 1
  if (!_lb_args.ignoreBgLoad()) {
      // processor mapping has already been determined by the background load pe
    for (i=0; i<npartition; i++) partitions[i].node = partitions[i].bkpes[0];
  }
  else {
    int n = 0;
    for (i=0; i<P; i++) { 
      if (!statsData->procs[i].available) continue;
      partitions[n++].node = i;
    }
  }
#else
  PatchMap *patchMap = PatchMap::Object();

  int **pool = new int *[P];
  for (i=0; i<P; i++) pool[i] = new int[P];
  for (i=0; i<P; i++) for (j=0; j<P; j++) pool[i][j] = 0;

  // sum up the number of nodes that patches of computes are on
  for (i=0; i<numComputes; i++)
  {
    for (j=0; j<P; j++)
      if (computeLoad[i].refno == partitions[j].refno) 
      {
	int node1 = patchMap->node(computes[i].patch1);
	int node2 = patchMap->node(computes[i].patch2);
	pool[j][node1]++;
	pool[j][node2]++;
      }
  }
#ifdef DEBUG
  for (i=0; i<P; i++) {
    for (j=0; j<P; j++) CmiPrintf("%d ", pool[i][j]);
    CmiPrintf("\n");
  }
#endif
  while (1)
  {
    int index=-1, node=0, eager=-1;
    for (j=0; j<npartition; j++) {
      if (partitions[j].node != -1) continue;
      int wantmost=-1, maxnodes=-1;
      for (k=0; k<P; k++) if (pool[j][k] > maxnodes && !partitions[k].mapped) {wantmost=k; maxnodes = pool[j][k];}
      if (maxnodes > eager) {
	index = j; node = wantmost; eager = maxnodes;
      }
    }
    if (eager == -1) break;
    partitions[index].node = node;
    partitions[node].mapped = 1;
  }

  for (i=0; i<P; i++) delete [] pool[i];
  delete [] pool;
#endif

/*
  if (_lb_args.debug()) {
    CmiPrintf("partition load: ");
    for (i=0; i<npartition; i++) CmiPrintf("%f ", partitions[i].load);
    CmiPrintf("\n");
    CmiPrintf("partitions to nodes mapping: ");
    for (i=0; i<npartition; i++) CmiPrintf("%d ", partitions[i].node);
    CmiPrintf("\n");
  }
*/
  if (_lb_args.debug()) {
    CmiPrintf("After partitioning: \n");
    for (i=0; i<npartition; i++) {
      double bgload = 0.0;
      if (!_lb_args.ignoreBgLoad())
        bgload = statsData->procs[partitions[i].bkpes[0]].bg_walltime;
      CmiPrintf("[%d=>%d] (%d,%d,%d) (%d,%d,%d) load:%f count:%d objload:%f\n", i, partitions[i].node, partitions[i].origin[0], partitions[i].origin[1], partitions[i].origin[2], partitions[i].corner[0], partitions[i].corner[1], partitions[i].corner[2], partitions[i].load, partitions[i].count, partitions[i].load-bgload);
    }
    for (i=npartition; i<P; i++) CmiPrintf("[%d] --------- \n", i);
  }

}

int main (int argc, char** argv)
{
  int format = 1;
  int n[3] = { 2, 2, 2 };
  int dims = 3;
  int skip=1;
  int start=0;
  int end=-1;
  bool last=false;
  char* infile = 0;
  char* vtffile = 0;
  float starttime = -1, endtime = -1;

  for (int i = 1; i < argc; i++)
    if (!strcmp(argv[i],"-format") && i < argc-1) {
      if (!strcasecmp(argv[++i],"ascii"))
        format = 0;
      else if (!strcasecmp(argv[i],"binary"))
        format = 1;
      else
        format = atoi(argv[i]);
    }
    else if (!strcmp(argv[i],"-nviz") && i < argc-1)
      n[0] = n[1] = n[2] = atoi(argv[++i]);
    else if (!strcmp(argv[i],"-1D"))
      dims = 1;
    else if (!strcmp(argv[i],"-2D"))
      dims = 2;
    else if (!strcmp(argv[i],"-last"))
      last = true;
    else if (!strcmp(argv[i],"-start") && i < argc-1)
      start = atoi(argv[++i]);
    else if (!strcmp(argv[i],"-starttime") && i < argc-1)
      starttime = atof(argv[++i]);
    else if (!strcmp(argv[i],"-end") && i < argc-1)
      end = atoi(argv[++i]);
    else if (!strcmp(argv[i],"-endtime") && i < argc-1)
      endtime = atof(argv[++i]);
    else if (!strcmp(argv[i],"-ndump") && i < argc-1)
      skip = atoi(argv[++i]);
    else if (!infile)
      infile = argv[i];
    else if (!vtffile)
      vtffile = argv[i];
    else
      std::cerr <<"  ** Unknown option ignored: "<< argv[i] << std::endl;

  if (!infile) {
    std::cout <<"usage: "<< argv[0]
              <<" <inputfile> [<vtffile>|<vtufile>] [-nviz <nviz>] \n"
              << "[-ndump <ndump>] [-last] [-start <level>] [-end <level>]\n"
              << "[-starttime <time>] [-endtime <time>] [-1D|-2D]\n"
              << "[-format <0|1|ASCII|BINARY>]\n";
    return 0;
  }
  else if (!vtffile)
    vtffile = infile;

  std::cout <<"\n >>> IFEM HDF5 to VT[F|U] converter <<<"
            <<"\n ==================================\n"
            <<"\nInput file: " << infile;

  std::cout <<"\nOutput file: "<< vtffile
            <<"\nNumber of visualization points: "
            << n[0] <<" "<< n[1] << " " << n[2] << std::endl;

  VTF* myVtf;
  if (strstr(vtffile,".vtf"))
    myVtf = new VTF(vtffile,format);
  else
    myVtf = new VTU(vtffile,last?1:0);

  // Process XML - establish fields and collapse bases
  PatchMap patches;

  HDF5Writer hdf(strtok(infile,"."),ProcessAdm(),true,true);
  XMLWriter xml(infile,ProcessAdm());
  xml.readInfo();

  int levels = xml.getLastTimeLevel();
  std::cout <<"Reading "<< infile <<": Time levels = "<< levels << std::endl;

  const std::vector<XMLWriter::Entry>& entry = xml.getEntries();
  std::vector<XMLWriter::Entry>::const_iterator it;

  ProcessList processlist;
  for (it = entry.begin(); it != entry.end(); ++it) {
    if (!it->basis.empty() && it->type != "restart") {
      processlist[it->basis].push_back(*it);
      std::cout << it->name <<"\t"<< it->description <<"\tnc="<< it->components
                <<"\t"<< it->basis << std::endl;
    }
    if (it->type == "eigenmodes") {
      levels = it->components-1;
      processlist[it->basis].back().components = 1;
    }
    if (it->type == "nodalforces")
      processlist["nodalforces"].push_back(*it);
  }

//.........这里部分代码省略.........

void ComputeGlobal::recvResults(ComputeGlobalResultsMsg *msg) {
  DebugM(3,"Receiving results (" << msg->aid.size() << " forces, "
	 << msg->newgdef.size() << " new group atoms) on client\n");

  // set the forces only if we aren't going to resend the data
  int setForces = !msg->resendCoordinates;

  if(setForces) { // we are requested to 
    // Store forces to patches
    PatchMap *patchMap = PatchMap::Object();
    int numPatches = patchMap->numPatches();
    AtomMap *atomMap = AtomMap::Object();
    const Lattice & lattice = patchList[0].p->lattice;
    ResizeArrayIter<PatchElem> ap(patchList);
    Force **f = new Force*[numPatches];
    FullAtom **t = new FullAtom*[numPatches];
    for ( int i = 0; i < numPatches; ++i ) { f[i] = 0; t[i] = 0; }
    Force extForce = 0.;
    Tensor extVirial;

    for (ap = ap.begin(); ap != ap.end(); ap++) {
      (*ap).r = (*ap).forceBox->open();
      f[(*ap).patchID] = (*ap).r->f[Results::normal];
      t[(*ap).patchID] = (*ap).p->getAtomList().begin();
    }

    AtomIDList::iterator a = msg->aid.begin();
    AtomIDList::iterator a_e = msg->aid.end();
    ForceList::iterator f2 = msg->f.begin();
    for ( ; a != a_e; ++a, ++f2 ) {
      DebugM(1,"processing atom "<<(*a)<<", F="<<(*f2)<<"...\n");
      /* XXX if (*a) is out of bounds here we get a segfault */
      LocalID localID = atomMap->localID(*a);
      if ( localID.pid == notUsed || ! f[localID.pid] ) continue;
      Force f_atom = (*f2);
      f[localID.pid][localID.index] += f_atom;
      Position x_orig = t[localID.pid][localID.index].position;
      Transform trans = t[localID.pid][localID.index].transform;
      Position x_atom = lattice.reverse_transform(x_orig,trans);
      extForce += f_atom;
      extVirial += outer(f_atom,x_atom);
    }
    DebugM(1,"done with the loop\n");

  // calculate forces for atoms in groups
    Molecule *mol = Node::Object()->molecule;
    AtomIDList::iterator g_i, g_e;
    g_i = gdef.begin(); g_e = gdef.end();
    ResizeArray<BigReal>::iterator gm_i = gmass.begin();
    ForceList::iterator gf_i = msg->gforce.begin();
    //iout << iDEBUG << "recvResults\n" << endi;
    for ( ; g_i != g_e; ++g_i, ++gm_i, ++gf_i ) {
      //iout << iDEBUG << *gf_i << '\n' << endi;
      Vector accel = (*gf_i) / (*gm_i);
      for ( ; *g_i != -1; ++g_i ) {
	//iout << iDEBUG << *g_i << '\n' << endi;
	LocalID localID = atomMap->localID(*g_i);
	if ( localID.pid == notUsed || ! f[localID.pid] ) continue;
	Force f_atom = accel * mol->atommass(*g_i);
	f[localID.pid][localID.index] += f_atom;
        Position x_orig = t[localID.pid][localID.index].position;
        Transform trans = t[localID.pid][localID.index].transform;
        Position x_atom = lattice.reverse_transform(x_orig,trans);
        extForce += f_atom;
        extVirial += outer(f_atom,x_atom);
      }
    }
    DebugM(1,"done with the groups\n");

    for (ap = ap.begin(); ap != ap.end(); ap++) {
      (*ap).forceBox->close(&((*ap).r));
    }

    delete [] f;
    delete [] t;

    ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
    ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
    reduction->submit();
  }
  // done setting the forces

  // Get reconfiguration if present
  if ( msg->reconfig ) configure(msg->newaid, msg->newgdef);

  // send another round of data if requested

  if(msg->resendCoordinates) {
    DebugM(3,"Sending requested data right away\n");
    sendData();
  }

  delete msg;
  DebugM(3,"Done processing results\n");
}

void registerUserEventsForAllComputeObjs()
{
#ifdef TRACE_COMPUTE_OBJECTS
    ComputeMap *map = ComputeMap::Object();
    PatchMap *pmap = PatchMap::Object();     
    char user_des[50];
    int p1, p2;
    int adim, bdim, cdim;
    int t1, t2;
    int x1, y1, z1, x2, y2, z2;
    int dx, dy, dz;
    for (int i=0; i<map->numComputes(); i++)
    {
        memset(user_des, 0, 50);
        switch ( map->type(i) )
        {
        case computeNonbondedSelfType:
            sprintf(user_des, "computeNonBondedSelfType_%d_pid_%d", i, map->pid(i,0));
            break;
        case computeLCPOType:
            sprintf(user_des, "computeLCPOType_%d_pid_%d", i, map->pid(i,0));
            break;
        case computeNonbondedPairType:
            adim = pmap->gridsize_a();
            bdim = pmap->gridsize_b();
            cdim = pmap->gridsize_c();
            p1 = map->pid(i, 0);
            t1 = map->trans(i, 0);
            x1 = pmap->index_a(p1) + adim * Lattice::offset_a(t1);
            y1 = pmap->index_b(p1) + bdim * Lattice::offset_b(t1);
            z1 = pmap->index_c(p1) + cdim * Lattice::offset_c(t1);
            p2 = map->pid(i, 1);
            t2 = map->trans(i, 1);
            x2 = pmap->index_a(p2) + adim * Lattice::offset_a(t2);
            y2 = pmap->index_b(p2) + bdim * Lattice::offset_b(t2);
            z2 = pmap->index_c(p2) + cdim * Lattice::offset_c(t2);
            dx = abs(x1-x2);
            dy = abs(y1-y2);
            dz = abs(z1-z2);
            sprintf(user_des, "computeNonBondedPairType_%d(%d,%d,%d)", i, dx,dy,dz);
            break;
        case computeExclsType:
            sprintf(user_des, "computeExclsType_%d", i);
            break;
        case computeBondsType:
            sprintf(user_des, "computeBondsType_%d", i);
            break;
        case computeAnglesType:
            sprintf(user_des, "computeAnglesType_%d", i);
            break;
        case computeDihedralsType:
            sprintf(user_des, "computeDihedralsType_%d", i);
            break;
        case computeImpropersType:
            sprintf(user_des, "computeImpropersType_%d", i);
            break;
        case computeTholeType:
            sprintf(user_des, "computeTholeType_%d", i);
            break;
        case computeAnisoType:
            sprintf(user_des, "computeAnisoType_%d", i);
            break;
        case computeCrosstermsType:
            sprintf(user_des, "computeCrosstermsType_%d", i);
            break;
        case computeSelfExclsType:
            sprintf(user_des, "computeSelfExclsType_%d", i);
            break;
        case computeSelfBondsType:
            sprintf(user_des, "computeSelfBondsType_%d", i);
            break;
        case computeSelfAnglesType:
            sprintf(user_des, "computeSelfAnglesType_%d", i);
            break;
        case computeSelfDihedralsType:
            sprintf(user_des, "computeSelfDihedralsType_%d", i);
            break;
        case computeSelfImpropersType:
            sprintf(user_des, "computeSelfImpropersType_%d", i);
            break;
        case computeSelfTholeType:
            sprintf(user_des, "computeSelfTholeType_%d", i);
            break;
        case computeSelfAnisoType:
            sprintf(user_des, "computeSelfAnisoType_%d", i);
            break;
        case computeSelfCrosstermsType:
            sprintf(user_des, "computeSelfCrosstermsType_%d", i);
            break;
#ifdef DPMTA
        case computeDPMTAType:
            sprintf(user_des, "computeDPMTAType_%d", i);
            break;
#endif
#ifdef DPME
        case computeDPMEType:
            sprintf(user_des, "computeDPMEType_%d", i);
            break;
#endif
        case computePmeType:
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    } else if ( patch[0]->flags.usePairlists && patch[1]->flags.usePairlists ) {
      if ( ! pairlistsValid ||
           ( patch[0]->flags.maxAtomMovement +
             patch[1]->flags.maxAtomMovement > pairlistTolerance ) ) {
        reductionData[pairlistWarningIndex] += 1;
      } else {
        params.usePairlists = 1;
      }
    }
    if ( ! params.usePairlists ) {
      pairlistsValid = 0;
    }
    params.plcutoff = cutoff;
    params.groupplcutoff = cutoff +
	patch[0]->flags.maxGroupRadius + patch[1]->flags.maxGroupRadius;
    if ( params.savePairlists ) {
      pairlistsValid = 1;
      pairlistTolerance = patch[0]->flags.pairlistTolerance +
                          patch[1]->flags.pairlistTolerance;
      params.plcutoff += pairlistTolerance;
      params.groupplcutoff += pairlistTolerance;
    }


    const Lattice &lattice = patch[0]->lattice;
    params.offset = lattice.offset(trans[a]) - lattice.offset(trans[b]);

    // Atom Sorting : If we are sorting the atoms along the line connecting
    //   the patch centers, then calculate a normalized vector pointing from
    //   patch a to patch b (i.e. outer loop patch to inner loop patch).
    #if NAMD_ComputeNonbonded_SortAtoms != 0

      // Center of patch a (outer-loop; i-loop) and patch b (inner-loop; j/k-loop)
      PatchMap* patchMap = PatchMap::Object();
      ScaledPosition p_a_center, p_b_center;
      p_a_center.x = patchMap->min_a(patchID[a]);
      p_a_center.y = patchMap->min_b(patchID[a]);
      p_a_center.z = patchMap->min_c(patchID[a]);
      p_b_center.x = patchMap->min_a(patchID[b]);
      p_b_center.y = patchMap->min_b(patchID[b]);
      p_b_center.z = patchMap->min_c(patchID[b]);
      p_a_center.x += patchMap->max_a(patchID[a]);
      p_a_center.y += patchMap->max_b(patchID[a]);
      p_a_center.z += patchMap->max_c(patchID[a]);
      p_b_center.x += patchMap->max_a(patchID[b]);
      p_b_center.y += patchMap->max_b(patchID[b]);
      p_b_center.z += patchMap->max_c(patchID[b]);
      p_a_center *= (BigReal)0.5;
      p_b_center *= (BigReal)0.5;
      p_a_center = lattice.unscale(p_a_center);
      p_b_center = lattice.unscale(p_b_center);

      // Adjust patch a's center by the offset
      p_a_center.x += params.offset.x;
      p_a_center.y += params.offset.y;
      p_a_center.z += params.offset.z;

      // Calculate and fill in the projected line vector
      params.projLineVec = p_b_center - p_a_center;
      params.projLineVec /= params.projLineVec.length(); // Normalize the vector

    #endif

      params.p[0] = p[a];
      params.p[1] = p[b];
      params.pExt[0] = pExt[a];