OOPSE/libmdtools/mpiSimulation.cpp

#ifdef IS_MPI
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cmath>
#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;
  
  MPI_Comm_size(MPI_COMM_WORLD, &(mpiPlug->numberProcessors) );
  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[] MolToProcMap;
  delete[] MolComponentType;
  delete[] AtomToProcMap;

  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;
  int done;
  int i, j, loops, which_proc, nmol_local, natoms_local;
  int nmol_global, natoms_global;
  int local_index;
  int baseSeed = entryPlug->getSeed();

  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( baseSeed );

  a = 3.0 * (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];
        add_atoms = compStamps[MolComponentType[i]]->getNAtoms();
        new_atoms = old_atoms + add_atoms;

        // 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++ ) {
            AtomToProcMap[atomIndex] = which_proc;
            atomIndex++;
          }
          done = 1;
          continue;
        }
    
        // 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++ ) {
            AtomToProcMap[atomIndex] = which_proc;
            atomIndex++;
          }
          done = 1;
          continue;
        }


        // The only situation left is when 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 = penalty / (average atoms per molecule)

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

    // Spray out this nonsense to all other processors:

    MPI_Bcast(MolToProcMap, mpiPlug->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(MolComponentType, mpiPlug->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(AtomsPerProc, mpiPlug->numberProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);    
  } else {

    // Listen to your marching orders from processor 0:
    
    MPI_Bcast(MolToProcMap, mpiPlug->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);
    
    MPI_Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Bcast(MolComponentType, mpiPlug->nMolGlobal, 
              MPI_INT, 0, MPI_COMM_WORLD);
    
    MPI_Bcast(AtomsPerProc, mpiPlug->numberProcessors,
              MPI_INT, 0, MPI_COMM_WORLD);


  }


  // 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_Allreduce(&nmol_local,&nmol_global,1,MPI_INT,MPI_SUM, 
                MPI_COMM_WORLD);
  MPI_Allreduce(&natoms_local,&natoms_global,1,MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);
  
  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) {
      globalIndex[local_index] = i;
      local_index++;
    }
  }
  
  return globalIndex;
}


void mpiSimulation::mpiRefresh( void ){

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

  // Fortran indexing needs to be increased by 1 in order to get the 2 languages to
  // not barf

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

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