OOPSE/libmdtools/mpiSimulation.cpp

#ifdef IS_MPI

#include <cstdlib>
#include <cstring>
#include <mpi.h>
#include <mpi++.h>

#include "mpiSimulation.hpp"
#include "simError.h"
#include "fortranWrappers.hpp"
#include "randomSPRNG.hpp"


mpiSimulation* mpiSim;

mpiSimulation::mpiSimulation(SimInfo* the_entryPlug)
{
  entryPlug = the_entryPlug;
  mpiPlug = new mpiSimData;
  
  mpiPlug->numberProcessors = MPI::COMM_WORLD.Get_size();
  mpiPlug->myNode = worldRank;

  MolToProcMap = new int[entryPlug->n_mol];
  MolComponentType = new int[entryPlug->n_mol];

  AtomToProcMap = new int[entryPlug->n_atoms];

  mpiSim = this;
  wrapMeSimParallel( this );
}


mpiSimulation::~mpiSimulation(){
  
  delete mpiPlug;
  // perhaps we should let fortran know the party is over.
  
}

int* mpiSimulation::divideLabor( void ){

  int* globalIndex;

  int nComponents;
  MoleculeStamp** compStamps;
  randomSPRNG myRandom;
  int* componentsNmol;
  int* AtomsPerProc;

  double numerator;
  double denominator;
  double precast;
  double x, y, a;
  int old_atoms, add_atoms, new_atoms;

  int nTarget;
  int molIndex, atomIndex, compIndex, compStart;
  int done;
  int nLocal, molLocal;
  int i, index;
  int smallDiff, bigDiff;

  int testSum;

  nComponents = entryPlug->nComponents;
  compStamps = entryPlug->compStamps;
  componentsNmol = entryPlug->componentsNmol;
  AtomsPerProc = new int[mpiPlug->numberProcessors];
  
  mpiPlug->nAtomsGlobal = entryPlug->n_atoms;
  mpiPlug->nBondsGlobal = entryPlug->n_bonds;
  mpiPlug->nBendsGlobal = entryPlug->n_bends;
  mpiPlug->nTorsionsGlobal = entryPlug->n_torsions;
  mpiPlug->nSRIGlobal = entryPlug->n_SRI;
  mpiPlug->nMolGlobal = entryPlug->n_mol;

  myRandom = new randomSPRNG();

  a = (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal;

  // Initialize things that we'll send out later:
  for (i = 0; i < mpiPlug->numberProcessors; i++ ) {
    AtomsPerProc[i] = 0;
  }
  for (i = 0; i < mpiPlug->nMolGlobal; i++ ) {
    // default to an error condition:
    MolToProcMap[i] = -1;
    MolComponentType[i] = -1;
  }
  for (i = 0; i < mpiPlug->nAtomsGlobal; i++ ) {
    // default to an error condition:
    AtomToProcMap[i] = -1;
  }
    
  if (mpiPlug->myNode == 0) {
    numerator = (double) entryPlug->n_atoms;
    denominator = (double) mpiPlug->numberProcessors;
    precast = numerator / denominator;
    nTarget = (int)( precast + 0.5 );

    // Build the array of molecule component types first
    molIndex = 0;
    for (i=0; i < nComponents; i++) {
      for (j=0; j < componentsNmol[i]; j++) {        
        MolComponentType[molIndex] = i;
        molIndex++;
      }
    }

    atomIndex = 0;

    for (i = 0; i < molIndex; i++ ) {

      done = 0;
      loops = 0;

      while( !done ){
        loops++;
        
        // Pick a processor at random

        which_proc = (int) (myRandom.getRandom() * mpiPlug->numberProcessors);

        // How many atoms does this processor have?
        
        old_atoms = AtomsPerProc[which_proc];

        // If the processor already had too many atoms, just skip this
        // processor and try again.

        if (old_atoms >= nTarget) continue;

        add_atoms = compStamps[MolComponentType[i]]->getNatoms();
        new_atoms = old_atoms + add_atoms;
    
        // If we can add this molecule to this processor without sending
        // it above nTarget, then go ahead and do it:
    
        if (new_atoms <= nTarget) {
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            atomIndex++;
            AtomToProcMap[atomIndex] = which_proc;
          }
          done = 1;
          continue;
        }

        // If we've been through this loop too many times, we need
        // to just give up and assign the molecule to this processor
        // and be done with it. 
        
        if (loops > 100) {          
          sprintf( painCave.errMsg,
                   "I've tried 100 times to assign molecule %d to a "
                   " processor, but can't find a good spot.\n"  
                   "I'm assigning it at random to processor %d.\n",
                   i, which_proc);
          painCave.isFatal = 0;
          simError();
          
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            atomIndex++;
            AtomToProcMap[atomIndex] = which_proc;
          }
          done = 1;
          continue;
        }

        // The only situation left is where old_atoms < nTarget, but
        // new_atoms > nTarget.   We want to accept this with some
        // probability that dies off the farther we are from nTarget

        // roughly:  x = new_atoms - nTarget
        //           Pacc(x) = exp(- a * x)
        // where a = 1 / (average atoms per molecule)

        x = (double) (new_atoms - nTarget);
        y = myRandom.getRandom();
        
        if (exp(- a * x) > y) {
          MolToProcMap[i] = which_proc;
          AtomsPerProc[which_proc] += add_atoms;
          for (j = 0 ; j < add_atoms; j++ ) {
            atomIndex++;
            AtomToProcMap[atomIndex] = which_proc;
          }
          done = 1;
          continue;
        } else {
          continue;
        }       
        
      }
    }

    // Spray out this nonsense to all other processors:

    MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, 
                          MPI_INT, 0);

    MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, 
                          MPI_INT, 0);

    MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, 
                          MPI_INT, 0);

    MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors,
                          MPI_INT, 0);    
  } else {

    // Listen to your marching orders from processor 0:
    
    MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, 
                          MPI_INT, 0);
    
    MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, 
                          MPI_INT, 0);

    MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, 
                          MPI_INT, 0);
    
    MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors,
                          MPI_INT, 0);
  }


  // Let's all check for sanity:

  nmol_local = 0;
  for (i = 0 ; i < mpiPlug->nMolGlobal; i++ ) {
    if (MolToProcMap[i] == mpiPlug->myNode) {
      nmol_local++;
    }
  }

  natoms_local = 0;
  for (i = 0; i < mpiPlug->nAtomsGlobal; i++) {
    if (AtomToProcMap[i] == mpiPlug->myNode) {
      natoms_local++;      
    }
  }

  MPI::COMM_WORLD.Allreduce(&nmol_local,&nmol_global,1,MPI_INT,MPI_SUM);
  MPI::COMM_WORLD.Allreduce(&natoms_local,&natoms_global,1,MPI_INT,MPI_SUM);
  
  if( nmol_global != entryPlug->n_mol ){
    sprintf( painCave.errMsg,
             "The sum of all nmol_local, %d, did not equal the "
             "total number of molecules, %d.\n",
             nmol_global, entryPlug->n_mol );
    painCave.isFatal = 1;
    simError();
  }
  
  if( natoms_global != entryPlug->n_atoms ){
    sprintf( painCave.errMsg,
             "The sum of all natoms_local, %d, did not equal the "
             "total number of atoms, %d.\n",
             natoms_global, entryPlug->n_atoms );
    painCave.isFatal = 1;
    simError();
  }

  sprintf( checkPointMsg,
           "Successfully divided the molecules among the processors.\n" );
  MPIcheckPoint();

  mpiPlug->myNMol = nmol_local;
  mpiPlug->myNlocal = natoms_local;

  globalIndex = new int[mpiPlug->myNlocal];
  local_index = 0;
  for (i = 0; i < mpiPlug->nAtomsGlobal; i++) {
    if (AtomToProcMap[i] == mpiPlug->myNode) {
      local_index++;
      globalIndex[local_index] = i;
    }
  }
  

   index = mpiPlug->myAtomStart;
//   for( i=0; i<mpiPlug->myNlocal; i++){
//     globalIndex[i] = index;
//     index++;
//   }

//   return globalIndex;
}


void mpiSimulation::mpiRefresh( void ){

  int isError, i;
  int *globalIndex = new int[mpiPlug->myNlocal];

  for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex();

  
  isError = 0;
  setFsimParallel( mpiPlug, &(entryPlug->n_atoms), globalIndex, &isError );
  if( isError ){

    sprintf( painCave.errMsg,
             "mpiRefresh errror: fortran didn't like something we gave it.\n" );
    painCave.isFatal = 1;
    simError();
  }

  delete[] globalIndex;

  sprintf( checkPointMsg,
           " mpiRefresh successful.\n" );
  MPIcheckPoint();
}
 

#endif // is_mpi
Revision:	406
Committed:	Wed Mar 26 19:34:49 2003 UTC (21 years, 3 months ago) by chuckv
File size:	8526 byte(s)
Log Message:	Finished globalIndex.
#	User	Rev	Content
1	mmeineke	377	#ifdef IS_MPI
2
3			#include <cstdlib>
4			#include <cstring>
5			#include <mpi.h>
6			#include <mpi++.h>
7
8			#include "mpiSimulation.hpp"
9			#include "simError.h"
10			#include "fortranWrappers.hpp"
11	gezelter	405	#include "randomSPRNG.hpp"
12	mmeineke	377
13
14			mpiSimulation* mpiSim;
15
16			mpiSimulation::mpiSimulation(SimInfo* the_entryPlug)
17			{
18			entryPlug = the_entryPlug;
19			mpiPlug = new mpiSimData;
20
21			mpiPlug->numberProcessors = MPI::COMM_WORLD.Get_size();
22			mpiPlug->myNode = worldRank;
23	gezelter	405
24			MolToProcMap = new int[entryPlug->n_mol];
25			MolComponentType = new int[entryPlug->n_mol];
26
27			AtomToProcMap = new int[entryPlug->n_atoms];
28
29	mmeineke	377	mpiSim = this;
30			wrapMeSimParallel( this );
31			}
32
33
34			mpiSimulation::~mpiSimulation(){
35
36			delete mpiPlug;
37			// perhaps we should let fortran know the party is over.
38
39			}
40
41			int* mpiSimulation::divideLabor( void ){
42
43			int* globalIndex;
44
45			int nComponents;
46			MoleculeStamp** compStamps;
47	gezelter	405	randomSPRNG myRandom;
48	mmeineke	377	int* componentsNmol;
49	gezelter	405	int* AtomsPerProc;
50	mmeineke	377
51			double numerator;
52			double denominator;
53			double precast;
54	gezelter	405	double x, y, a;
55			int old_atoms, add_atoms, new_atoms;
56	mmeineke	377
57			int nTarget;
58			int molIndex, atomIndex, compIndex, compStart;
59			int done;
60			int nLocal, molLocal;
61			int i, index;
62			int smallDiff, bigDiff;
63
64			int testSum;
65
66			nComponents = entryPlug->nComponents;
67			compStamps = entryPlug->compStamps;
68			componentsNmol = entryPlug->componentsNmol;
69	gezelter	405	AtomsPerProc = new int[mpiPlug->numberProcessors];
70
71	mmeineke	377	mpiPlug->nAtomsGlobal = entryPlug->n_atoms;
72			mpiPlug->nBondsGlobal = entryPlug->n_bonds;
73			mpiPlug->nBendsGlobal = entryPlug->n_bends;
74			mpiPlug->nTorsionsGlobal = entryPlug->n_torsions;
75			mpiPlug->nSRIGlobal = entryPlug->n_SRI;
76			mpiPlug->nMolGlobal = entryPlug->n_mol;
77
78	gezelter	405	myRandom = new randomSPRNG();
79	chuckv	404
80	gezelter	405	a = (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal;
81	chuckv	404
82	gezelter	405	// Initialize things that we'll send out later:
83			for (i = 0; i < mpiPlug->numberProcessors; i++ ) {
84			AtomsPerProc[i] = 0;
85			}
86			for (i = 0; i < mpiPlug->nMolGlobal; i++ ) {
87			// default to an error condition:
88			MolToProcMap[i] = -1;
89			MolComponentType[i] = -1;
90			}
91			for (i = 0; i < mpiPlug->nAtomsGlobal; i++ ) {
92			// default to an error condition:
93			AtomToProcMap[i] = -1;
94			}
95
96			if (mpiPlug->myNode == 0) {
97			numerator = (double) entryPlug->n_atoms;
98			denominator = (double) mpiPlug->numberProcessors;
99			precast = numerator / denominator;
100			nTarget = (int)( precast + 0.5 );
101	chuckv	404
102	gezelter	405	// Build the array of molecule component types first
103			molIndex = 0;
104			for (i=0; i < nComponents; i++) {
105			for (j=0; j < componentsNmol[i]; j++) {
106			MolComponentType[molIndex] = i;
107			molIndex++;
108			}
109			}
110	chuckv	404
111	gezelter	405	atomIndex = 0;
112	chuckv	404
113	gezelter	405	for (i = 0; i < molIndex; i++ ) {
114	chuckv	404
115	gezelter	405	done = 0;
116			loops = 0;
117	chuckv	404
118	gezelter	405	while( !done ){
119			loops++;
120
121			// Pick a processor at random
122	chuckv	404
123	gezelter	405	which_proc = (int) (myRandom.getRandom() * mpiPlug->numberProcessors);
124	chuckv	404
125	gezelter	405	// How many atoms does this processor have?
126
127			old_atoms = AtomsPerProc[which_proc];
128	chuckv	404
129	gezelter	405	// If the processor already had too many atoms, just skip this
130			// processor and try again.
131	chuckv	404
132	gezelter	405	if (old_atoms >= nTarget) continue;
133	mmeineke	377
134	gezelter	405	add_atoms = compStamps[MolComponentType[i]]->getNatoms();
135			new_atoms = old_atoms + add_atoms;
136	mmeineke	377
137	gezelter	405	// If we can add this molecule to this processor without sending
138			// it above nTarget, then go ahead and do it:
139
140			if (new_atoms <= nTarget) {
141			MolToProcMap[i] = which_proc;
142			AtomsPerProc[which_proc] += add_atoms;
143			for (j = 0 ; j < add_atoms; j++ ) {
144			atomIndex++;
145			AtomToProcMap[atomIndex] = which_proc;
146			}
147			done = 1;
148			continue;
149			}
150	mmeineke	377
151	gezelter	405	// If we've been through this loop too many times, we need
152			// to just give up and assign the molecule to this processor
153			// and be done with it.
154
155			if (loops > 100) {
156			sprintf( painCave.errMsg,
157			"I've tried 100 times to assign molecule %d to a "
158			" processor, but can't find a good spot.\n"
159			"I'm assigning it at random to processor %d.\n",
160			i, which_proc);
161			painCave.isFatal = 0;
162			simError();
163
164			MolToProcMap[i] = which_proc;
165			AtomsPerProc[which_proc] += add_atoms;
166			for (j = 0 ; j < add_atoms; j++ ) {
167			atomIndex++;
168			AtomToProcMap[atomIndex] = which_proc;
169			}
170			done = 1;
171			continue;
172			}
173
174			// The only situation left is where old_atoms < nTarget, but
175			// new_atoms > nTarget. We want to accept this with some
176			// probability that dies off the farther we are from nTarget
177
178			// roughly: x = new_atoms - nTarget
179			// Pacc(x) = exp(- a * x)
180			// where a = 1 / (average atoms per molecule)
181
182			x = (double) (new_atoms - nTarget);
183			y = myRandom.getRandom();
184
185			if (exp(- a * x) > y) {
186			MolToProcMap[i] = which_proc;
187			AtomsPerProc[which_proc] += add_atoms;
188			for (j = 0 ; j < add_atoms; j++ ) {
189			atomIndex++;
190			AtomToProcMap[atomIndex] = which_proc;
191			}
192			done = 1;
193			continue;
194			} else {
195			continue;
196			}
197
198	mmeineke	377	}
199			}
200	gezelter	405
201			// Spray out this nonsense to all other processors:
202
203			MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal,
204			MPI_INT, 0);
205
206			MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal,
207			MPI_INT, 0);
208
209			MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal,
210			MPI_INT, 0);
211
212			MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors,
213			MPI_INT, 0);
214			} else {
215
216			// Listen to your marching orders from processor 0:
217	mmeineke	377
218	gezelter	405	MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal,
219			MPI_INT, 0);
220	mmeineke	377
221	gezelter	405	MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal,
222			MPI_INT, 0);
223
224			MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal,
225			MPI_INT, 0);
226
227			MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors,
228			MPI_INT, 0);
229	mmeineke	377	}
230
231
232	gezelter	405	// Let's all check for sanity:
233
234			nmol_local = 0;
235			for (i = 0 ; i < mpiPlug->nMolGlobal; i++ ) {
236			if (MolToProcMap[i] == mpiPlug->myNode) {
237			nmol_local++;
238			}
239	mmeineke	377	}
240
241	gezelter	405	natoms_local = 0;
242			for (i = 0; i < mpiPlug->nAtomsGlobal; i++) {
243			if (AtomToProcMap[i] == mpiPlug->myNode) {
244			natoms_local++;
245			}
246			}
247	mmeineke	377
248	gezelter	405	MPI::COMM_WORLD.Allreduce(&nmol_local,&nmol_global,1,MPI_INT,MPI_SUM);
249			MPI::COMM_WORLD.Allreduce(&natoms_local,&natoms_global,1,MPI_INT,MPI_SUM);
250	mmeineke	377
251	gezelter	405	if( nmol_global != entryPlug->n_mol ){
252			sprintf( painCave.errMsg,
253			"The sum of all nmol_local, %d, did not equal the "
254			"total number of molecules, %d.\n",
255			nmol_global, entryPlug->n_mol );
256			painCave.isFatal = 1;
257			simError();
258	mmeineke	377	}
259	gezelter	405
260			if( natoms_global != entryPlug->n_atoms ){
261			sprintf( painCave.errMsg,
262			"The sum of all natoms_local, %d, did not equal the "
263			"total number of atoms, %d.\n",
264			natoms_global, entryPlug->n_atoms );
265			painCave.isFatal = 1;
266			simError();
267			}
268	mmeineke	377
269			sprintf( checkPointMsg,
270			"Successfully divided the molecules among the processors.\n" );
271			MPIcheckPoint();
272
273	gezelter	405	mpiPlug->myNMol = nmol_local;
274			mpiPlug->myNlocal = natoms_local;
275	mmeineke	377
276			globalIndex = new int[mpiPlug->myNlocal];
277	gezelter	405	local_index = 0;
278			for (i = 0; i < mpiPlug->nAtomsGlobal; i++) {
279			if (AtomToProcMap[i] == mpiPlug->myNode) {
280	chuckv	406	local_index++;
281			globalIndex[local_index] = i;
282	gezelter	405	}
283	mmeineke	377	}
284	gezelter	405
285	mmeineke	377
286	gezelter	405
287
288			index = mpiPlug->myAtomStart;
289			// for( i=0; i<mpiPlug->myNlocal; i++){
290			// globalIndex[i] = index;
291			// index++;
292			// }
293
294			// return globalIndex;
295	mmeineke	377	}
296
297
298			void mpiSimulation::mpiRefresh( void ){
299
300			int isError, i;
301			int *globalIndex = new int[mpiPlug->myNlocal];
302
303			for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex();
304
305
306			isError = 0;
307			setFsimParallel( mpiPlug, &(entryPlug->n_atoms), globalIndex, &isError );
308			if( isError ){
309
310			sprintf( painCave.errMsg,
311			"mpiRefresh errror: fortran didn't like something we gave it.\n" );
312			painCave.isFatal = 1;
313			simError();
314			}
315
316			delete[] globalIndex;
317
318			sprintf( checkPointMsg,
319			" mpiRefresh successful.\n" );
320			MPIcheckPoint();
321			}
322
323
324			#endif // is_mpi