--- trunk/OOPSE/libmdtools/mpiSimulation.cpp 2003/03/26 19:34:49 406 +++ trunk/OOPSE/libmdtools/mpiSimulation.cpp 2003/03/27 23:33:40 432 @@ -1,7 +1,8 @@ #ifdef IS_MPI - +#include #include #include +#include #include #include @@ -10,6 +11,7 @@ #include "fortranWrappers.hpp" #include "randomSPRNG.hpp" +#define BASE_SEED 123456789 mpiSimulation* mpiSim; @@ -23,7 +25,6 @@ mpiSimulation::mpiSimulation(SimInfo* the_entryPlug) MolToProcMap = new int[entryPlug->n_mol]; MolComponentType = new int[entryPlug->n_mol]; - AtomToProcMap = new int[entryPlug->n_atoms]; mpiSim = this; @@ -33,6 +34,10 @@ mpiSimulation::~mpiSimulation(){ mpiSimulation::~mpiSimulation(){ + delete[] MolToProcMap; + delete[] MolComponentType; + delete[] AtomToProcMap; + delete mpiPlug; // perhaps we should let fortran know the party is over. @@ -44,7 +49,7 @@ int* mpiSimulation::divideLabor( void ){ int nComponents; MoleculeStamp** compStamps; - randomSPRNG myRandom; + randomSPRNG *myRandom; int* componentsNmol; int* AtomsPerProc; @@ -58,8 +63,11 @@ int* mpiSimulation::divideLabor( void ){ int molIndex, atomIndex, compIndex, compStart; int done; int nLocal, molLocal; - int i, index; + int i, j, loops, which_proc, nmol_local, natoms_local; + int nmol_global, natoms_global; + int local_index, index; int smallDiff, bigDiff; + int baseSeed = BASE_SEED; int testSum; @@ -75,7 +83,7 @@ int* mpiSimulation::divideLabor( void ){ mpiPlug->nSRIGlobal = entryPlug->n_SRI; mpiPlug->nMolGlobal = entryPlug->n_mol; - myRandom = new randomSPRNG(); + myRandom = new randomSPRNG( baseSeed ); a = (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal; @@ -120,34 +128,17 @@ int* mpiSimulation::divideLabor( void ){ // Pick a processor at random - which_proc = (int) (myRandom.getRandom() * mpiPlug->numberProcessors); + 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 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. @@ -164,13 +155,30 @@ int* mpiSimulation::divideLabor( void ){ MolToProcMap[i] = which_proc; AtomsPerProc[which_proc] += add_atoms; for (j = 0 ; j < add_atoms; j++ ) { - atomIndex++; - AtomToProcMap[atomIndex] = which_proc; + AtomToProcMap[atomIndex] = which_proc; + atomIndex++; } done = 1; continue; } + if (old_atoms >= nTarget) 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 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 @@ -180,15 +188,15 @@ int* mpiSimulation::divideLabor( void ){ // where a = 1 / (average atoms per molecule) x = (double) (new_atoms - nTarget); - y = myRandom.getRandom(); + 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; - } + AtomToProcMap[atomIndex] = which_proc; + atomIndex++; + } done = 1; continue; } else { @@ -200,32 +208,34 @@ int* mpiSimulation::divideLabor( void ){ // Spray out this nonsense to all other processors: - MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, + MPI::COMM_WORLD.Bcast(MolToProcMap, mpiPlug->nMolGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, + MPI::COMM_WORLD.Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, + MPI::COMM_WORLD.Bcast(MolComponentType, mpiPlug->nMolGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors, + 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::COMM_WORLD.Bcast(MolToProcMap, mpiPlug->nMolGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, + MPI::COMM_WORLD.Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, + MPI::COMM_WORLD.Bcast(MolComponentType, mpiPlug->nMolGlobal, MPI_INT, 0); - MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors, + MPI::COMM_WORLD.Bcast(AtomsPerProc, mpiPlug->numberProcessors, MPI_INT, 0); + + } @@ -245,6 +255,8 @@ int* mpiSimulation::divideLabor( void ){ } } + std::cerr << "proc = " << mpiPlug->myNode << " atoms = " << natoms_local << "\n"; + 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); @@ -277,21 +289,12 @@ int* mpiSimulation::divideLabor( void ){ local_index = 0; for (i = 0; i < mpiPlug->nAtomsGlobal; i++) { if (AtomToProcMap[i] == mpiPlug->myNode) { - local_index++; globalIndex[local_index] = i; + local_index++; } } - - - - index = mpiPlug->myAtomStart; -// for( i=0; imyNlocal; i++){ -// globalIndex[i] = index; -// index++; -// } - -// return globalIndex; + return globalIndex; }