--- trunk/OOPSE/libmdtools/mpiSimulation.cpp 2003/03/26 16:12:53 404 +++ trunk/OOPSE/libmdtools/mpiSimulation.cpp 2003/03/26 18:02:18 405 @@ -8,10 +8,9 @@ #include "mpiSimulation.hpp" #include "simError.h" #include "fortranWrappers.hpp" +#include "randomSPRNG.hpp" - - mpiSimulation* mpiSim; mpiSimulation::mpiSimulation(SimInfo* the_entryPlug) @@ -21,7 +20,12 @@ mpiSimulation::mpiSimulation(SimInfo* the_entryPlug) 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 ); } @@ -34,19 +38,21 @@ mpiSimulation::~mpiSimulation(){ } - - 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; @@ -60,7 +66,8 @@ int* mpiSimulation::divideLabor( void ){ 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; @@ -68,135 +75,222 @@ int* mpiSimulation::divideLabor( void ){ 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. - numerator = (double) entryPlug->n_atoms; - denominator = (double) mpiPlug->numberProcessors; - precast = numerator / denominator; - nTarget = (int)( precast + 0.5 ); - - molIndex = 0; - atomIndex = 0; - compIndex = 0; - compStart = 0; - for( i=0; i<(mpiPlug->numberProcessors-1); i++){ - - done = 0; - nLocal = 0; - molLocal = 0; + if (old_atoms >= nTarget) continue; - if( i == mpiPlug->myNode ){ - mpiPlug->myMolStart = molIndex; - mpiPlug->myAtomStart = atomIndex; - } + add_atoms = compStamps[MolComponentType[i]]->getNatoms(); + new_atoms = old_atoms + add_atoms; - while( !done ){ - - if( (molIndex-compStart) >= componentsNmol[compIndex] ){ - compStart = molIndex; - compIndex++; - 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++ ) { + atomIndex++; + AtomToProcMap[atomIndex] = which_proc; + } + done = 1; + continue; + } - nLocal += compStamps[compIndex]->getNAtoms(); - atomIndex += compStamps[compIndex]->getNAtoms(); - molIndex++; - molLocal++; - - if ( nLocal == nTarget ) done = 1; - - else if( nLocal < nTarget ){ - smallDiff = nTarget - nLocal; - } - else if( nLocal > nTarget ){ - bigDiff = nLocal - nTarget; - - if( bigDiff < smallDiff ) done = 1; - else{ - molIndex--; - molLocal--; - atomIndex -= compStamps[compIndex]->getNAtoms(); - nLocal -= compStamps[compIndex]->getNAtoms(); - done = 1; - } + // 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: - if( i == mpiPlug->myNode ){ - mpiPlug->myMolEnd = (molIndex - 1); - mpiPlug->myAtomEnd = (atomIndex - 1); - mpiPlug->myNlocal = nLocal; - mpiPlug->myMol = molLocal; - } + MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, + MPI_INT, 0); - numerator = (double)( entryPlug->n_atoms - atomIndex ); - denominator = (double)( mpiPlug->numberProcessors - (i+1) ); - precast = numerator / denominator; - nTarget = (int)( precast + 0.5 ); + 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); } - - if( mpiPlug->myNode == mpiPlug->numberProcessors-1 ){ - mpiPlug->myMolStart = molIndex; - mpiPlug->myAtomStart = atomIndex; - nLocal = 0; - molLocal = 0; - while( compIndex < nComponents ){ - - if( (molIndex-compStart) >= componentsNmol[compIndex] ){ - compStart = molIndex; - compIndex++; - continue; - } - nLocal += compStamps[compIndex]->getNAtoms(); - atomIndex += compStamps[compIndex]->getNAtoms(); - molIndex++; - molLocal++; - } - - mpiPlug->myMolEnd = (molIndex - 1); - mpiPlug->myAtomEnd = (atomIndex - 1); - mpiPlug->myNlocal = nLocal; - mpiPlug->myMol = molLocal; + // 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( &nLocal, &testSum, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); + 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( mpiPlug->myNode == 0 ){ - if( testSum != entryPlug->n_atoms ){ - sprintf( painCave.errMsg, - "The summ of all nLocals, %d, did not equal the total number of atoms, %d.\n", - testSum, entryPlug->n_atoms ); - painCave.isFatal = 1; - simError(); - } + 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(); - // lets create the identity array + mpiPlug->myNMol = nmol_local; + mpiPlug->myNlocal = natoms_local; globalIndex = new int[mpiPlug->myNlocal]; - index = mpiPlug->myAtomStart; - for( i=0; imyNlocal; i++){ - globalIndex[i] = index; - index++; + local_index = 0; + for (i = 0; i < mpiPlug->nAtomsGlobal; i++) { + if (AtomToProcMap[i] == mpiPlug->myNode) { + globalIndex[local_index] = + } } + - return globalIndex; + + + index = mpiPlug->myAtomStart; +// for( i=0; imyNlocal; i++){ +// globalIndex[i] = index; +// index++; +// } + +// return globalIndex; }