Add support for Charm++'s change to use 64 bit ID for load balancing
[namd.git] / src / NamdCentLB.C
1 /*****************************************************************************
2  * $Source: /home/cvs/namd/cvsroot/namd2/src/NamdCentLB.C,v $
3  * $Author: jim $
4  * $Date: 2017/03/30 20:06:17 $
5  * $Revision: 1.125 $
6  *****************************************************************************/
7
8 #if !defined(WIN32) || defined(__CYGWIN__)
9 #include <unistd.h>
10 #endif
11 #include <fcntl.h>
12
13 #include "InfoStream.h"
14 #include "NamdCentLB.h"
15 #include "NamdCentLB.def.h"
16 #include "Node.h"
17 #include "PatchMap.h"
18 #include "ComputeMap.h"
19 #include "LdbCoordinator.h"
20
21 // #define DUMP_LDBDATA 1
22 // #define LOAD_LDBDATA 1
23
24 double *cpuloads = NULL;
25
26 void CreateNamdCentLB() {
27   // CkPrintf("[%d] creating NamdCentLB %d\n",CkMyPe(),loadbalancer);
28   loadbalancer = CProxy_NamdCentLB::ckNew();
29   // CkPrintf("[%d] created NamdCentLB %d\n",CkMyPe(),loadbalancer);
30   if (CkMyRank() == 0 && cpuloads == NULL) {    
31     cpuloads = new double[CkNumPes()];
32     CmiAssert(cpuloads != NULL);
33     for (int i=0; i<CkNumPes(); i++) cpuloads[i] = 0.0;
34   }
35 }
36
37 NamdCentLB *AllocateNamdCentLB() {
38   return new NamdCentLB((CkMigrateMessage*)NULL);
39 }
40
41 /**
42  * Migratable Object Constructor.
43  */
44 NamdCentLB::NamdCentLB(CkMigrateMessage *msg): CentralLB(msg) {
45   processorArray = 0;
46   patchArray = 0;
47   computeArray = 0;
48
49
50 NamdCentLB::NamdCentLB(): CentralLB(CkLBOptions(-1))
51 {
52   //  if (CkMyPe()==0)
53   //   CkPrintf("[%d] NamdCentLB created\n",CkMyPe());
54   processorArray = 0;
55   patchArray = 0;
56   computeArray = 0;
57 }
58
59 /*
60 NamdCentLB::~NamdCentLB()
61 {
62   delete [] processorArray;
63   delete [] patchArray;
64   delete [] computeArray;
65 }
66 */
67
68 bool NamdCentLB::QueryBalanceNow(int _step)
69 {
70   //  CkPrintf("[%d] Balancing on step %d\n",CkMyPe(),_step);
71   if ( LdbCoordinator::Object()->takingLdbData ) {
72     return true;
73   } else {
74     return false;
75   }
76 }
77
78 bool NamdCentLB::QueryDumpData()
79 {
80 #if 0
81   if (LdbCoordinator::Object()->ldbCycleNum == 1)  return true;
82   if (LdbCoordinator::Object()->ldbCycleNum == 2)  return true;
83 #endif
84   return false;
85 }
86
87 CLBMigrateMsg* NamdCentLB::Strategy(LDStats* stats)
88 {
89   //  CkPrintf("LDB: All statistics received at %f, %f\n",
90   //  CmiTimer(),CmiWallTimer());
91
92   int numProcessors = stats->nprocs();
93   int numPatches = PatchMap::Object()->numPatches();
94   ComputeMap *computeMap = ComputeMap::Object();
95   const int numComputes = computeMap->numComputes();
96   const SimParameters* simParams = Node::Object()->simParameters;
97
98   // these sizes should never change
99   if ( ! processorArray ) processorArray = new processorInfo[numProcessors];
100   if ( ! patchArray ) patchArray = new patchInfo[numPatches];
101   if ( ! computeArray ) computeArray = new computeInfo[numComputes];
102
103   int nMoveableComputes = buildData(stats);
104
105 #if LDB_DEBUG
106 #define DUMP_LDBDATA 1
107 #define LOAD_LDBDATA 1
108 #endif
109
110 #if DUMP_LDBDATA 
111   dumpDataASCII("ldbd_before", numProcessors, numPatches, nMoveableComputes);
112 #elif LOAD_LDBDATA
113   loadDataASCII("ldbd_before.5", numProcessors, numPatches, nMoveableComputes);
114   // CkExit();
115 #endif
116
117   double averageLoad = 0.;
118   double avgCompute = 0.;
119   if ( nMoveableComputes ) {
120    int i;
121    double total = 0.;
122    double maxCompute = 0.;
123    int maxi = 0;
124    for (i=0; i<nMoveableComputes; i++) {
125       double load = computeArray[i].load;
126       total += load;
127       if ( load > maxCompute ) { maxCompute = load;  maxi = i; }
128    }
129    avgCompute = total / nMoveableComputes;
130
131     int P = stats->nprocs();
132    int numPesAvailable = 0;
133    for (i=0; i<P; i++) {
134       if (processorArray[i].available) {
135         ++numPesAvailable;
136         total += processorArray[i].backgroundLoad;
137       }
138    }
139    if (numPesAvailable == 0)
140      NAMD_die("No processors available for load balancing!\n");
141
142    averageLoad = total/numPesAvailable;
143    CkPrintf("LDB: Largest compute %d load %f is %.1f%% of average load %f\n",
144             LdbIdField(computeArray[maxi].handle.id, 0),
145             maxCompute, 100. * maxCompute / averageLoad, averageLoad);
146    CkPrintf("LDB: Average compute %f is %.1f%% of average load %f\n",
147             avgCompute, 100. * avgCompute / averageLoad, averageLoad);
148   }
149
150   if ( step() == 1 ) {
151     // compute splitting only
152     // partitions are stored as char but mostly limited by
153     // high load noise at low outer-loop iteration counts
154     int maxParts = 10;
155 #ifdef NAMD_CUDA
156 //split LCPO compute very small, else CUDA compute is delayed
157     if (simParams->LCPOOn) {
158       maxParts = 20;
159     }
160 #endif
161     int totalAddedParts = 0;
162     double maxCompute = averageLoad / 10.;
163     if ( maxCompute < 2. * avgCompute ) maxCompute = 2. * avgCompute;
164     if ( simParams->ldbRelativeGrainsize > 0. ) {
165       maxCompute = averageLoad * simParams->ldbRelativeGrainsize;
166     }
167     CkPrintf("LDB: Partitioning computes with target load %f\n", maxCompute);
168     double maxUnsplit = 0.;
169     for (int i=0; i<nMoveableComputes; i++) {
170       computeArray[i].processor = computeArray[i].oldProcessor;
171       const int cid = LdbIdField(computeArray[i].handle.id, 0);
172       const double load = computeArray[i].load;
173       if ( computeMap->numPartitions(cid) == 0 ) {
174         if ( load > maxUnsplit ) maxUnsplit = load;
175         continue;
176       }
177       int nparts = (int) ceil(load / maxCompute);
178       if ( nparts > maxParts ) nparts = maxParts;
179       if ( nparts < 1 ) nparts = 1;
180       if ( 0 && nparts > 1 ) {
181         CkPrintf("LDB: Partitioning compute %d with load %f by %d\n",
182                   cid, load, nparts);
183       }
184       computeMap->setNewNumPartitions(cid,nparts);
185       totalAddedParts += nparts - 1;
186     }
187     CkPrintf("LDB: Increased migratable compute count from %d to %d\n",
188               nMoveableComputes,nMoveableComputes+totalAddedParts);
189     CkPrintf("LDB: Largest unpartitionable compute is %f\n", maxUnsplit);
190   } else if (simParams->ldbStrategy == LDBSTRAT_DEFAULT) { // default
191     if (step() < 4)
192       TorusLB(computeArray, patchArray, processorArray,
193                   nMoveableComputes, numPatches, numProcessors);
194     else
195       RefineTorusLB(computeArray, patchArray, processorArray,
196                   nMoveableComputes, numPatches, numProcessors, 1);
197   } else if (simParams->ldbStrategy == LDBSTRAT_COMPREHENSIVE) {
198     TorusLB(computeArray, patchArray, processorArray,
199                   nMoveableComputes, numPatches, numProcessors);
200   } else if (simParams->ldbStrategy == LDBSTRAT_REFINEONLY) {
201     RefineTorusLB(computeArray, patchArray, processorArray,
202                   nMoveableComputes, numPatches, numProcessors, 1);
203   } else if (simParams->ldbStrategy == LDBSTRAT_OLD) {
204     if (step() < 4)
205       Alg7(computeArray, patchArray, processorArray,
206                   nMoveableComputes, numPatches, numProcessors);
207     else
208       RefineOnly(computeArray, patchArray, processorArray, 
209                   nMoveableComputes, numPatches, numProcessors);
210   }
211
212 #if LDB_DEBUG && USE_TOPOMAP
213   TopoManager tmgr;
214   int pe1, pe2, pe3, hops=0;
215   /* This is double counting the hops
216   for(int i=0; i<nMoveableComputes; i++)
217   {
218     pe1 = computeArray[i].processor;
219     pe2 = patchArray[computeArray[i].patch1].processor;
220     pe3 = patchArray[computeArray[i].patch2].processor;
221     hops += tmgr.getHopsBetweenRanks(pe1, pe2);
222     if(computeArray[i].patch1 != computeArray[i].patch2)
223       hops += tmgr.getHopsBetweenRanks(pe1, pe3);  
224   }*/
225   for (int i=0; i<numPatches; i++)  {
226     //int num = patchArray[i].proxiesOn.numElements();
227     pe1 = patchArray[i].processor;
228     Iterator nextProc;
229     processorInfo *p = (processorInfo *)patchArray[i].proxiesOn.iterator((Iterator *)&nextProc);
230     while (p) {
231       pe2 = p->Id;
232       hops += tmgr.getHopsBetweenRanks(pe1, pe2);
233       p = (processorInfo *)patchArray[i].proxiesOn.next((Iterator*)&nextProc);
234     }
235   }
236   CkPrintf("Load Balancing: Number of Hops: %d\n", hops);
237 #endif
238
239 #if DUMP_LDBDATA
240   dumpDataASCII("ldbd_after", numProcessors, numPatches, nMoveableComputes);
241 #elif LOAD_LDBDATA
242   dumpDataASCII("ldbd_after.5", numProcessors, numPatches, nMoveableComputes);
243   // loadDataASCII("ldbd_after", numProcessors, numPatches, nMoveableComputes);
244   // CkExit();
245 #endif
246
247   // For error checking:
248   // Count up computes, to see if somebody doesn't have any computes
249   int i;
250 #if 0
251   int* computeCount = new int[numProcessors];
252   for(i=0; i<numProcessors; i++)
253     computeCount[i]=0;
254   for(i=0; i<nMoveableComputes; i++)
255     computeCount[computeArray[i].processor]++;
256   for(i=0; i<numProcessors; i++) {
257     if (computeCount[i]==0)
258       iout << iINFO <<"Warning: Processor " << i 
259            << " has NO moveable computes.\n" << endi;
260   }
261   delete [] computeCount;
262 #endif
263   
264   CkVec<MigrateInfo *> migrateInfo;
265   for(i=0;i<nMoveableComputes;i++) {
266     if (computeArray[i].processor != computeArray[i].oldProcessor) {
267       //      CkPrintf("[%d] Obj %d migrating from %d to %d\n",
268       //               CkMyPe(),computeArray[i].handle.id.id[0],
269       //               computeArray[i].processor,computeArray[i].oldProcessor);
270       MigrateInfo *migrateMe = new MigrateInfo;
271       migrateMe->obj = computeArray[i].handle;
272       migrateMe->from_pe = computeArray[i].oldProcessor;
273       migrateMe->to_pe = computeArray[i].processor;
274       migrateInfo.insertAtEnd(migrateMe);
275
276       // sneak in updates to ComputeMap
277       computeMap->setNewNode(LdbIdField(computeArray[i].handle.id, 0),
278                                 computeArray[i].processor);
279     }
280   }
281   
282   int migrate_count=migrateInfo.length();
283   // CkPrintf("NamdCentLB migrating %d elements\n",migrate_count);
284   CLBMigrateMsg* msg = new(migrate_count,CkNumPes(),CkNumPes(),0) CLBMigrateMsg;
285
286   msg->n_moves = migrate_count;
287   for(i=0; i < migrate_count; i++) {
288     MigrateInfo* item = migrateInfo[i];
289     msg->moves[i] = *item;
290     delete item;
291     migrateInfo[i] = 0;
292   }
293
294   for (i=0; i<numProcessors; i++) {
295     cpuloads[i] = processorArray[i].load;
296   }
297
298   delete [] processorArray;
299   delete [] patchArray;
300   delete [] computeArray;
301
302   processorArray = NULL;
303   patchArray = NULL;
304   computeArray = NULL;
305   
306   return msg;
307 };
308
309 #ifndef WIN32
310
311 void NamdCentLB::dumpDataASCII(char *file, int numProcessors,
312                                int numPatches, int numComputes)
313 {
314   char filename[128];
315   sprintf(filename, "%s.%d", file, step());
316   FILE* fp = fopen(filename,"w");
317   if (fp == NULL){
318      perror("dumpLDStatsASCII");
319      return;
320   }
321   CkPrintf("***** DUMP data to file: %s ***** \n", filename);
322   fprintf(fp,"%d %d %d\n",numProcessors,numPatches,numComputes);
323
324   int i;
325   for(i=0;i<numProcessors;i++) {
326     processorInfo* p = processorArray + i;
327     fprintf(fp,"%d %e %e %e %e\n",p->Id,p->load,p->backgroundLoad,p->computeLoad,p->idleTime);
328   }
329
330   for(i=0;i < numPatches; i++) {
331     patchInfo* p = patchArray + i;
332     fprintf(fp,"%d %e %d %d\n",p->Id,p->load,p->processor,p->numAtoms);
333   }
334     
335   for(i=0; i < numComputes; i++) {
336     computeInfo* c = computeArray + i;
337     fprintf(fp,"%d %e %d %d %d %d",c->Id,c->load,c->patch1,c->patch2,
338             c->processor,c->oldProcessor);
339     fprintf(fp, "\n");
340   }
341
342   // dump patchSet
343   for (i=0; i< numProcessors; i++) {
344       int num = processorArray[i].proxies.numElements();
345       fprintf(fp, "%d %d: ", i, num);
346       Iterator nextProxy;
347       patchInfo *p = (patchInfo *)processorArray[i].proxies.
348         iterator((Iterator *)&nextProxy);
349       while (p) {
350           fprintf(fp, "%d ", p->Id);
351           p = (patchInfo *)processorArray[i].proxies.
352             next((Iterator*)&nextProxy);
353       }
354       fprintf(fp, "\n");
355   }
356   // dump proxiesOn
357   for (i=0; i<numPatches; i++)  {
358     int num = patchArray[i].proxiesOn.numElements();
359     fprintf(fp, "%d %d: ", i, num);
360       Iterator nextProc;
361       processorInfo *p = (processorInfo *)patchArray[i].proxiesOn.
362         iterator((Iterator *)&nextProc);
363       while (p) {
364         fprintf(fp, "%d ", p->Id);
365         p = (processorInfo *)patchArray[i].proxiesOn.
366           next((Iterator*)&nextProc);
367       }
368       fprintf(fp, "\n");
369   }
370
371   fclose(fp);
372   //CkExit();
373 }
374
375 void NamdCentLB::loadDataASCII(char *file, int &numProcessors,
376                                int &numPatches, int &numComputes)
377 {
378   char filename[128];
379   //sprintf(filename, "%s.%d", file, step());
380   sprintf(filename, "%s", file);
381
382   CkPrintf("***** Load ascii data from file: %s ***** \n", filename);
383
384   FILE* fp = fopen(filename, "r");
385   if (fp == NULL){
386      perror("loadDataASCII");
387      return;
388   }
389
390   fscanf(fp,"%d %d %d",&numProcessors,&numPatches,&numComputes);
391
392   printf("numProcs: %d numPatches: %d numComputes: %d\n", numProcessors,numPatches, numComputes);
393
394   delete [] processorArray;
395   delete [] patchArray;
396   delete [] computeArray;
397   processorArray = new processorInfo[numProcessors];
398   patchArray = new patchInfo[numPatches];
399   computeArray = new computeInfo[numComputes];
400
401   int i;
402   for(i=0;i<numProcessors;i++) {
403     processorInfo* p = processorArray + i;
404     fscanf(fp,"%d %le %le %le", &p->Id, &p->load, &p->backgroundLoad, &p->computeLoad);
405     fscanf(fp,"%le\n", &p->idleTime);
406     if (p->Id != i) CmiAbort("Reading processorArray error!");
407 //    p->backgroundLoad = 0.0;
408   }
409
410   for(i=0;i < numPatches; i++) {
411     patchInfo* p = patchArray + i;
412     fscanf(fp,"%d %le %d %d\n",&p->Id,&p->load,&p->processor,&p->numAtoms);
413     if (p->Id != i || p->processor > numProcessors || p->processor < 0) 
414       CmiAbort("Reading patchArray error!");
415   }
416     
417   for(i=0; i < numComputes; i++) {
418     computeInfo* c = computeArray + i;
419     fscanf(fp,"%d %le %d %d %d %d",&c->Id,&c->load,&c->patch1,&c->patch2,
420             &c->processor,&c->oldProcessor);
421
422     if (c->patch1 < 0 || c->patch1 > numPatches || c->patch2 < 0 || c->patch2 > numPatches)
423       CmiAbort("Reading computeArray error!");
424   // printf("%d %e %d %d %d %d\n", c->Id,c->load,c->patch1,c->patch2,c->processor,c->oldProcessor);
425   }
426
427   // dump patchSet
428   for (i=0; i< numProcessors; i++) {
429       int num, curp;
430       fscanf(fp,"%d %d: ",&curp, &num);
431       if(curp != i)
432         CmiAbort("Reading patchsSet error!");
433       for (int j=0; j<num; j++) {
434           int id;
435           fscanf(fp,"%d",&id);
436           processorArray[i].proxies.unchecked_insert(&patchArray[id]);
437       }
438   }
439   // dump proxiesOn
440   for (i=0; i<numPatches; i++)  {
441       int num, curp;
442       fscanf(fp,"%d %d: ",&curp, &num);
443       if(curp != i)
444         CmiAbort("Reading proxiesOn error!");
445       for (int j=0; j<num; j++) {
446           int id;
447           fscanf(fp,"%d",&id);
448           patchArray[i].proxiesOn.insert(&processorArray[id]);
449       }
450   }
451
452   fclose(fp);
453 }
454 #endif
455
456 extern int isPmeProcessor(int); 
457 #ifdef MEM_OPT_VERSION
458 extern int isOutputProcessor(int); 
459 #endif
460 #if defined(NAMD_MIC)
461 extern int isMICProcessor(int);
462 #endif
463
464 int NamdCentLB::buildData(LDStats* stats)
465 {
466   int n_pes = stats->nprocs();
467
468   PatchMap* patchMap = PatchMap::Object();
469   ComputeMap* computeMap = ComputeMap::Object();
470   const SimParameters* simParams = Node::Object()->simParameters;
471
472   BigReal bgfactor = simParams->ldbBackgroundScaling;
473   BigReal pmebgfactor = simParams->ldbPMEBackgroundScaling;
474   BigReal homebgfactor = simParams->ldbHomeBackgroundScaling;
475   int pmeOn = simParams->PMEOn;
476   int unLoadPme = simParams->ldbUnloadPME;
477   int pmeBarrier = simParams->PMEBarrier;
478   int unLoadZero = simParams->ldbUnloadZero;
479   int unLoadOne = simParams->ldbUnloadOne;
480   int unLoadIO= simParams->ldbUnloadOutputPEs;
481   int i;
482   for (i=0; i<n_pes; ++i) {
483     processorArray[i].Id = i;
484     processorArray[i].available = true;
485     if ( pmeOn && isPmeProcessor(i) ) {
486       processorArray[i].backgroundLoad = pmebgfactor * stats->procs[i].bg_walltime;
487     } else if (patchMap->numPatchesOnNode(i) > 0) {
488       processorArray[i].backgroundLoad = homebgfactor * stats->procs[i].bg_walltime;
489     } else {
490       processorArray[i].backgroundLoad = bgfactor * stats->procs[i].bg_walltime;
491     }
492     processorArray[i].idleTime = stats->procs[i].idletime;
493     processorArray[i].load = processorArray[i].computeLoad = 0.0;
494   }
495
496 /* *********** this code is defunct *****************
497 #if 0
498   double bgfactor = 1.0 + 1.0 * CkNumPes()/1000.0;
499   if ( bgfactor > 2.0 ) bgfactor = 2.0;
500   iout << iINFO << "Scaling background load by " << bgfactor << ".\n" << endi;
501   int i;
502   for (i=0; i<n_pes; i++) {
503     processorArray[i].Id = i;
504     processorArray[i].backgroundLoad = bgfactor * stats[i].bg_walltime;
505   }
506
507   double bg_weight = 0.7;
508
509   int i;
510   for (i=0; i<n_pes; i++) {
511     processorArray[i].Id = i;
512     if (patchMap->numPatchesOnNode(i) > 0)
513       processorArray[i].backgroundLoad = bg_weight * stats->procs[i].bg_walltime;
514     else 
515       processorArray[i].backgroundLoad = stats[i].bg_walltime;
516   }
517   
518   //Modification to reduce the coputeload on PME processors
519   const SimParameters* simParams = Node::Object()->simParameters;  
520   
521   // CkPrintf("BACKGROUND LOAD\n");
522   if(simParams->PMEOn) {
523     double bgfactor = 1.0 + 1.0 * CkNumPes()/1000.0;
524     if ( bgfactor > 2.0 ) bgfactor = 2.0;
525     for (i=0; i<n_pes; i++) {
526       // CkPrintf("BG[%d] =  %5.5lf,", i, processorArray[i].backgroundLoad);
527       if(isPmeProcessor(i)) {
528         processorArray[i].backgroundLoad *= bgfactor;
529       }
530       // CkPrintf("%5.5lf;  ", processorArray[i].backgroundLoad);
531     }
532   }
533   // CkPrintf("\n");
534 #endif  
535 *********** end of defunct code *********** */
536
537   if (unLoadZero) processorArray[0].available = false;
538   if (unLoadOne) processorArray[1].available = false;
539
540   // if all pes are Pme, disable this flag
541   if (pmeOn && unLoadPme) {
542     for (i=0; i<n_pes; i++) {
543       if (!isPmeProcessor(i))  break;
544     }
545     if (i == n_pes) {
546       iout << iINFO << "Turned off unLoadPme flag!\n"  << endi;
547       unLoadPme = 0;
548     }
549   }
550   
551   if (pmeOn && unLoadPme) {
552     for (i=0; i<n_pes; i++) {
553       if ((pmeBarrier && i==0) || isPmeProcessor(i)) 
554         processorArray[i].available = false;
555     }
556   }
557   // if all pes are output, disable this flag
558 #ifdef MEM_OPT_VERSION
559
560   if (unLoadIO) {
561       if (simParams->numoutputprocs == n_pes) {
562           iout << iINFO << "Turned off unLoadIO flag!\n"  << endi;
563           unLoadIO = 0;
564       }
565   }
566   if (unLoadIO){
567     iout << iINFO << "Testing for output processors!\n"  << endi;
568       for (i=0; i<n_pes; i++) {
569           if (isOutputProcessor(stats->procs[i].pe)) 
570             {
571               //              iout << iINFO << "Removed output PE "<< stats->procs[i].pe <<" from available list!\n"  << endi;
572               processorArray[i].available = false;
573             }
574           else
575             {
576               //              iout << iINFO << "Nonoutput PE "<< stats->procs[i].pe <<" is in available list!\n"  << endi;
577             }
578       }
579   }
580 #endif
581
582   // Unload PEs driving MIC devices, if need be
583   #if defined(NAMD_MIC)
584     if (simParams->mic_unloadMICPEs != 0) {
585       for (i = 0; i < n_pes; i++) {
586         if (isMICProcessor(i) != 0) { processorArray[i].available = false; }
587       }
588     }
589   #endif
590
591   int nMoveableComputes=0;
592   int nProxies = 0;             // total number of estimated proxies
593   int nIdleComputes = 0;
594
595   int j;
596   for (j=0; j < stats->n_objs; j++) {
597       const LDObjData &this_obj = stats->objData[j];
598       int frompe = stats->from_proc[j];
599
600       // filter out non-NAMD managed objects (like PME array)
601       if (this_obj.omID().id.idx != 1) {
602         // CkPrintf("non-NAMD object %d on pe %d with walltime %lf\n",
603         // this_obj.id().id[0], stats->from_proc[j], this_obj.wallTime);
604         processorArray[stats->from_proc[j]].backgroundLoad += this_obj.wallTime;
605         continue;
606       }
607
608       if (LdbIdField(this_obj.id(), 1) == PATCH_TYPE) { // Its a patch
609         const int pid = LdbIdField(this_obj.id(), 0);
610         int neighborNodes[PatchMap::MaxOneAway + PatchMap::MaxTwoAway];
611
612         patchArray[pid].Id = pid;
613         patchArray[pid].numAtoms = 0;
614         patchArray[pid].processor = stats->from_proc[j];
615         const int numProxies = 
616 #if USE_TOPOMAP
617         requiredProxiesOnProcGrid(pid,neighborNodes);
618 #else
619         requiredProxies(pid, neighborNodes);
620 #endif
621
622         nProxies += numProxies;
623
624         for (int k=0; k<numProxies; k++) {
625           processorArray[neighborNodes[k]].proxies.unchecked_insert(&patchArray[pid]);
626           patchArray[pid].proxiesOn.unchecked_insert(&processorArray[neighborNodes[k]]);
627         }
628         processorArray[stats->from_proc[j]].backgroundLoad += this_obj.wallTime;
629       } else if (LdbIdField(this_obj.id(), 1) == BONDED_TYPE) { // Its a bonded compute
630         processorArray[stats->from_proc[j]].backgroundLoad += this_obj.wallTime;
631       } else if (this_obj.migratable) { // Its a compute
632        if ( this_obj.wallTime == 0. ) { // don't migrate idle computes
633          ++nIdleComputes;
634        } else {
635         const int cid = LdbIdField(this_obj.id(), 0);
636         const int p0 = computeMap->pid(cid,0);
637
638         // For self-interactions, just return the same pid twice
639         int p1;
640         if (computeMap->numPids(cid) > 1)
641           p1 = computeMap->pid(cid,1);
642         else p1 = p0;
643         computeArray[nMoveableComputes].Id = cid;
644         computeArray[nMoveableComputes].oldProcessor = stats->from_proc[j];
645         processorArray[stats->from_proc[j]].computeLoad += this_obj.wallTime;
646         computeArray[nMoveableComputes].processor = -1;
647         computeArray[nMoveableComputes].patch1 = p0;
648         computeArray[nMoveableComputes].patch2 = p1;
649         computeArray[nMoveableComputes].handle = this_obj.handle;
650         computeArray[nMoveableComputes].load = this_obj.wallTime;
651         nMoveableComputes++;
652        }
653       } else {
654         processorArray[stats->from_proc[j]].backgroundLoad += this_obj.wallTime;
655       }
656     }
657
658    if ( nIdleComputes )
659      CkPrintf("LDB: %d computes have load of zero\n", nIdleComputes);
660
661 /* *********** this code is defunct *****************
662 #if 0
663   int averageProxy = nProxies / n_pes;
664   CkPrintf("total proxies: %d, avervage: %d\n", nProxies, averageProxy);
665   for (i=0; i<n_pes; i++) {
666     // too many proxies on this node, weight the background load
667     int proxies = processorArray[i].proxies.numElements();
668     if (proxies > averageProxy) {
669       double factor = 1.0*(proxies-averageProxy)/nProxies;
670       processorArray[i].backgroundLoad *= (1.0 + factor);
671       CkPrintf("On [%d]: too many proxies: %d, increased bg load by %f\n", i, nProxies, factor);
672     }
673   }
674 #endif
675 *********** end of defunct code *********** */
676
677   for (i=0; i<n_pes; i++) {
678     processorArray[i].load = processorArray[i].backgroundLoad + processorArray[i].computeLoad;
679   }
680   stats->clear();
681   return nMoveableComputes;
682 }
683
684 // Figure out which proxies we will definitely create on other
685 // nodes, without regard for non-bonded computes.  This code is swiped
686 // from ProxyMgr, and changes there probable need to be propagated here.
687
688 int NamdCentLB::requiredProxies(PatchID id, int neighborNodes[])
689 {
690   PatchMap* patchMap = PatchMap::Object();
691   int myNode = patchMap->node(id);
692   int nProxyNodes = 0;
693
694 #define IF_NEW_NODE \
695     int j; \
696     for ( j=0; j<nProxyNodes && neighborNodes[j] != proxyNode; ++j ); \
697     if ( j == nProxyNodes )
698
699   PatchID neighbors[1 + PatchMap::MaxOneAway + PatchMap::MaxTwoAway];
700   neighbors[0] = id;
701   int numNeighbors = 1 + patchMap->downstreamNeighbors(id,neighbors+1);
702   for ( int i = 0; i < numNeighbors; ++i ) {
703     const int proxyNode = patchMap->basenode(neighbors[i]);
704     if ( proxyNode != myNode ) {
705       IF_NEW_NODE {
706         neighborNodes[nProxyNodes] = proxyNode;
707         nProxyNodes++;
708       }
709     }
710   }
711
712   // Distribute initial default proxies across empty processors.
713   // This shouldn't be necessary, but may constrain the load balancer
714   // and avoid placing too many proxies on a single processor.  -JCP
715
716   // This code needs to be turned off when the creation of ST is
717   // shifted to the load balancers -ASB
718
719 #if 1
720   int numPes = CkNumPes();
721   int numPatches = patchMap->numPatches();
722   int emptyNodes = numPes - numPatches;
723   if ( emptyNodes > numPatches ) {
724     int nodesPerPatch = nProxyNodes + 1 + (emptyNodes-1) / numPatches;
725     int maxNodesPerPatch = PatchMap::MaxOneAway + PatchMap::MaxTwoAway;
726     if ( nodesPerPatch > maxNodesPerPatch ) nodesPerPatch = maxNodesPerPatch;
727     int proxyNode = (myNode + 1) % numPes;
728     while ( nProxyNodes < nodesPerPatch &&
729                         ! patchMap->numPatchesOnNode(proxyNode) ) {
730       if ( proxyNode != myNode ) {
731         IF_NEW_NODE {
732           neighborNodes[nProxyNodes] = proxyNode;
733           nProxyNodes++;
734         }
735       }
736       proxyNode = (proxyNode + 1) % numPes;
737     }
738     proxyNode = (myNode - 1 + numPes) % numPes;
739     while ( nProxyNodes < nodesPerPatch &&
740                         ! patchMap->numPatchesOnNode(proxyNode) ) {
741       if ( proxyNode != myNode ) {
742         IF_NEW_NODE {
743           neighborNodes[nProxyNodes] = proxyNode;
744           nProxyNodes++;
745         }
746       }
747       proxyNode = (proxyNode - 1 + numPes) % numPes;
748     }
749     proxyNode = (myNode + 1) % numPes;
750     int count = 0;
751     while ( nProxyNodes < nodesPerPatch ) {
752       if ( ! patchMap->numPatchesOnNode(proxyNode) && proxyNode != myNode ) {
753         IF_NEW_NODE {
754           neighborNodes[nProxyNodes] = proxyNode;
755           nProxyNodes++;
756         }
757       }
758       proxyNode = (proxyNode + 1) % numPes;
759       count ++; if (count == numPes) break;   // we looped all
760     }
761   } else {
762     int proxyNode = myNode - 1;
763     if ( proxyNode >= 0 && ! patchMap->numPatchesOnNode(proxyNode) ) {
764       if ( proxyNode != myNode ) {
765         IF_NEW_NODE {
766           neighborNodes[nProxyNodes] = proxyNode;
767           nProxyNodes++;
768         }
769       }
770     }
771     proxyNode = myNode + 1;
772     if ( proxyNode < numPes && ! patchMap->numPatchesOnNode(proxyNode) ) {
773       if ( proxyNode != myNode ) {
774         IF_NEW_NODE {
775           neighborNodes[nProxyNodes] = proxyNode;
776           nProxyNodes++;
777         }
778       }
779     }
780   }
781 #endif
782
783   return nProxyNodes;
784 }
785
786 #if USE_TOPOMAP 
787 // Figure out which proxies we will definitely create on other nodes,
788 // without regard for non-bonded computes.  This code is swiped from
789 // ProxyMgr, and changes there probable need to be propagated here.
790 // The proxies are placed on nearby processors on the 3d-grid along
791 // the X, Y, Z and T dimensions
792
793 int NamdCentLB::requiredProxiesOnProcGrid(PatchID id, int neighborNodes[])
794 {
795   enum proxyHere { No, Yes };
796   int numPes = CkNumPes();
797   proxyHere *proxyNodes = new proxyHere[numPes];
798   int nProxyNodes;
799   int i, j, k, l;
800
801   int xsize = 0, ysize = 0, zsize = 0, tsize = 0;
802   int my_x = 0, my_y = 0, my_z = 0, my_t = 0;
803
804   PatchMap* patchMap = PatchMap::Object();
805   int myNode = patchMap->node(id);
806     
807   TopoManager tmgr;
808   xsize = tmgr.getDimNX();
809   ysize = tmgr.getDimNY();
810   zsize = tmgr.getDimNZ();
811   tsize = tmgr.getDimNT();
812   
813   tmgr.rankToCoordinates(myNode, my_x, my_y, my_z, my_t);
814   
815   if(xsize * ysize * zsize * tsize != CkNumPes()) {
816     delete [] proxyNodes;
817     return requiredProxies(id, neighborNodes);
818   }  
819
820   // Note all home patches.
821   for ( i = 0; i < numPes; ++i )
822   {
823     proxyNodes[i] = No;
824   }
825   nProxyNodes = 0;
826
827   // Check all two-away neighbors.
828   // This is really just one-away neighbors, since 
829   // two-away always returns zero: RKB
830   PatchID neighbors[1 + PatchMap::MaxOneAway + PatchMap::MaxTwoAway];
831
832   // Assign a proxy to all your neighbors. But dont increment counter
833   // because these have to be there anyway.
834   neighbors[0] = id;  
835   int numNeighbors = 1 + patchMap->downstreamNeighbors(id,neighbors+1);
836   
837   // Small Flag chooses between different loadbalancing schemes.
838   // Small Flag == true, patches are close to each other
839   // false, patches are far from each other
840   bool smallFlag = false;
841   double pnodes = CkNumPes();
842   pnodes *= 0.25;    
843   smallFlag = (patchMap->numPatches() > pnodes )?1:0;
844
845   //If there are lot of patches its likely they will all be neighbors, 
846   //so all we need to do is to place proxies on downstream patches.
847   //if (smallFlag) {
848   for ( i = 1; i < numNeighbors; ++i )
849     {
850       int proxyNode = patchMap->basenode(neighbors[i]);
851       
852       if (proxyNode != myNode)
853         if (proxyNodes[proxyNode] == No)
854           {
855             proxyNodes[proxyNode] = Yes;
856             neighborNodes[nProxyNodes] = proxyNode;
857             nProxyNodes++;
858           }
859     }
860   //}
861  
862   if (step() > 2) {
863     delete [] proxyNodes;
864     return nProxyNodes;
865   }
866  
867   // Place numPesPerPatch proxies on the 3d torus neighbors of a processor
868
869   int numPatches = patchMap->numPatches();
870   int emptyNodes = numPes - numPatches;
871   //if ( emptyNodes > numPatches ) {
872   
873   int nodesPerPatch = nProxyNodes + 4 * (emptyNodes-1) / numPatches + 1;
874   int proxyNode = 0 ;
875   int proxy_x=0, proxy_y=0, proxy_z=0;
876   
877   //Choose from the 26 neighbors of mynode.
878   //CkAssert(nodesPerPatch - nProxyNodes <= 26);  
879   //Too few patches otherwise, try twoaway?
880   
881   for(k=-1; k<= 1; k++) {
882     proxy_z = (my_z + k + zsize) % zsize;
883     for(j=-1; j <= 1; j++) {
884       proxy_y = (my_y + j + ysize) % ysize;
885       for(i = -1; i <= 1; i++) {
886         proxy_x = (my_x + i + xsize) % xsize;
887         for(l = 0; l < tsize; l++) {
888           if(i == 0 && j == 0 && k == 0 && l == 0)
889             continue;
890
891           proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z, l);
892
893           if((! patchMap->numPatchesOnNode(proxyNode) || !smallFlag) &&
894              proxyNodes[proxyNode] == No) {
895             proxyNodes[proxyNode] = Yes;
896             neighborNodes[nProxyNodes] = proxyNode;
897             nProxyNodes++;
898           }
899           
900           if(nProxyNodes >= nodesPerPatch || 
901              nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
902             break;
903         } // end for
904
905         if(nProxyNodes >= nodesPerPatch || 
906            nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
907           break;
908       } // end for
909       
910       if(nProxyNodes >= nodesPerPatch || 
911          nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
912         break;    
913     } // end for
914
915     if(nProxyNodes >= nodesPerPatch || 
916        nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
917       break;      
918   } // end for
919
920 #if 1
921   if(!smallFlag) {
922     for(k=-2; k<= 2; k+=2) {
923       proxy_z = (my_z + k + zsize) % zsize;
924       for(j=-2; j <= 2; j+=2) {
925         proxy_y = (my_y + j + ysize) % ysize;
926         for(i = -2; i <= 2; i+=2) {
927           proxy_x = (my_x + i + xsize) % xsize;
928           for(l = 0; l < tsize; l++) {
929             if(i == 0 && j == 0 && k == 0 && l == 0)
930               continue;
931           
932             proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z, l);
933           
934             if((! patchMap->numPatchesOnNode(proxyNode) || !smallFlag) &&
935                proxyNodes[proxyNode] == No) {
936               proxyNodes[proxyNode] = Yes;
937               neighborNodes[nProxyNodes] = proxyNode;
938               nProxyNodes++;
939             }
940             
941             if(nProxyNodes >= nodesPerPatch || 
942                nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
943               break;
944           } // end for
945
946           if(nProxyNodes >= nodesPerPatch || 
947              nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
948             break;
949         } // end for
950         
951         if(nProxyNodes >= nodesPerPatch || 
952            nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
953           break;          
954       } // end for
955
956       if(nProxyNodes >= nodesPerPatch || 
957          nProxyNodes >= PatchMap::MaxOneAway + PatchMap::MaxTwoAway)
958         break;    
959     } // end for
960   }
961
962 #else
963   #if 0
964   const SimParameters* params = Node::Object()->simParameters;
965
966   if(!smallFlag) {
967     //Add two-away proxies
968     if(patchMap->numaway_a() == 2) {
969       proxy_y = (my_y + 2) % ysize;
970       proxy_x = my_x  % xsize;
971       proxy_z = my_z  % zsize;
972       
973       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
974       if(proxyNodes[proxyNode] == No) {
975         proxyNodes[proxyNode] = Yes;
976         neighborNodes[nProxyNodes] = proxyNode;
977       nProxyNodes++;
978       }
979       
980       proxy_y = (my_y - 2 + ysize) % ysize;
981       proxy_x = my_x  % xsize;
982       proxy_z = my_z % zsize;
983       
984       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
985       if(proxyNodes[proxyNode] == No) {
986         proxyNodes[proxyNode] = Yes;
987         neighborNodes[nProxyNodes] = proxyNode;
988         nProxyNodes++;
989       }
990     }
991     
992     //Add two away proxies
993     if(patchMap->numaway_b() == 2) {
994       proxy_y = my_y  % ysize;
995       proxy_x = my_x  % xsize;
996       proxy_z = (my_z + 2) % zsize;
997       
998       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
999       if(proxyNodes[proxyNode] == No) {
1000         proxyNodes[proxyNode] = Yes;
1001         neighborNodes[nProxyNodes] = proxyNode;
1002         nProxyNodes++;
1003       }
1004       
1005       proxy_y = my_y  % ysize;
1006       proxy_x = my_x  % xsize;
1007       proxy_z = (my_z - 2 + zsize) % zsize;
1008       
1009       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
1010       if(proxyNodes[proxyNode] == No) {
1011         proxyNodes[proxyNode] = Yes;
1012         neighborNodes[nProxyNodes] = proxyNode;
1013         nProxyNodes++;
1014       }
1015     }
1016     
1017     //Add two away proxies
1018     if(patchMap->numaway_c() == 2) {
1019       proxy_y = my_y  % ysize;
1020       proxy_x = (my_x + 2) % xsize;
1021       proxy_z = my_z  % zsize;
1022       
1023       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
1024       if(proxyNodes[proxyNode] == No) {
1025         proxyNodes[proxyNode] = Yes;
1026         neighborNodes[nProxyNodes] = proxyNode;
1027       nProxyNodes++;
1028       }
1029       
1030       proxy_y = my_y  % ysize;
1031       proxy_x = (my_x  - 2 + xsize) % xsize;
1032       proxy_z = my_z % zsize;
1033       
1034       proxyNode = tmgr.coordinatesToRank(proxy_x, proxy_y, proxy_z);
1035       if(proxyNodes[proxyNode] == No) {
1036         proxyNodes[proxyNode] = Yes;
1037         neighborNodes[nProxyNodes] = proxyNode;
1038         nProxyNodes++;
1039       }
1040     }
1041   }
1042   #endif
1043 #endif
1044   
1045   // CkPrintf("Returning %d proxies\n", nProxyNodes);
1046
1047   delete [] proxyNodes;
1048   return nProxyNodes;
1049 }
1050
1051 #endif