# | Line 1 | Line 1 | |
---|---|---|
1 | #ifdef IS_MPI | |
2 | < | |
2 | > | #include <iostream> |
3 | #include <cstdlib> | |
4 | #include <cstring> | |
5 | #include <cmath> | |
# | Line 25 | Line 25 | mpiSimulation::mpiSimulation(SimInfo* the_entryPlug) | |
25 | ||
26 | MolToProcMap = new int[entryPlug->n_mol]; | |
27 | MolComponentType = new int[entryPlug->n_mol]; | |
28 | – | |
28 | AtomToProcMap = new int[entryPlug->n_atoms]; | |
29 | ||
30 | mpiSim = this; | |
# | Line 35 | Line 34 | mpiSimulation::~mpiSimulation(){ | |
34 | ||
35 | mpiSimulation::~mpiSimulation(){ | |
36 | ||
37 | + | delete[] MolToProcMap; |
38 | + | delete[] MolComponentType; |
39 | + | delete[] AtomToProcMap; |
40 | + | |
41 | delete mpiPlug; | |
42 | // perhaps we should let fortran know the party is over. | |
43 | ||
# | Line 82 | Line 85 | int* mpiSimulation::divideLabor( void ){ | |
85 | ||
86 | myRandom = new randomSPRNG( baseSeed ); | |
87 | ||
88 | < | a = (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal; |
88 | > | a = 3.0 * (double)mpiPlug->nMolGlobal / (double)mpiPlug->nAtomsGlobal; |
89 | ||
90 | // Initialize things that we'll send out later: | |
91 | for (i = 0; i < mpiPlug->numberProcessors; i++ ) { | |
# | Line 130 | Line 133 | int* mpiSimulation::divideLabor( void ){ | |
133 | // How many atoms does this processor have? | |
134 | ||
135 | old_atoms = AtomsPerProc[which_proc]; | |
133 | – | |
134 | – | // If the processor already had too many atoms, just skip this |
135 | – | // processor and try again. |
136 | – | |
137 | – | if (old_atoms >= nTarget) continue; |
138 | – | |
136 | add_atoms = compStamps[MolComponentType[i]]->getNAtoms(); | |
137 | new_atoms = old_atoms + add_atoms; | |
141 | – | |
142 | – | // If we can add this molecule to this processor without sending |
143 | – | // it above nTarget, then go ahead and do it: |
144 | – | |
145 | – | if (new_atoms <= nTarget) { |
146 | – | MolToProcMap[i] = which_proc; |
147 | – | AtomsPerProc[which_proc] += add_atoms; |
148 | – | for (j = 0 ; j < add_atoms; j++ ) { |
149 | – | atomIndex++; |
150 | – | AtomToProcMap[atomIndex] = which_proc; |
151 | – | } |
152 | – | done = 1; |
153 | – | continue; |
154 | – | } |
138 | ||
139 | // If we've been through this loop too many times, we need | |
140 | // to just give up and assign the molecule to this processor | |
# | Line 169 | Line 152 | int* mpiSimulation::divideLabor( void ){ | |
152 | MolToProcMap[i] = which_proc; | |
153 | AtomsPerProc[which_proc] += add_atoms; | |
154 | for (j = 0 ; j < add_atoms; j++ ) { | |
155 | < | atomIndex++; |
156 | < | AtomToProcMap[atomIndex] = which_proc; |
155 | > | AtomToProcMap[atomIndex] = which_proc; |
156 | > | atomIndex++; |
157 | } | |
158 | done = 1; | |
159 | continue; | |
160 | } | |
161 | + | |
162 | + | // If we can add this molecule to this processor without sending |
163 | + | // it above nTarget, then go ahead and do it: |
164 | + | |
165 | + | if (new_atoms <= nTarget) { |
166 | + | MolToProcMap[i] = which_proc; |
167 | + | AtomsPerProc[which_proc] += add_atoms; |
168 | + | for (j = 0 ; j < add_atoms; j++ ) { |
169 | + | AtomToProcMap[atomIndex] = which_proc; |
170 | + | atomIndex++; |
171 | + | } |
172 | + | done = 1; |
173 | + | continue; |
174 | + | } |
175 | ||
179 | – | // The only situation left is where old_atoms < nTarget, but |
180 | – | // new_atoms > nTarget. We want to accept this with some |
181 | – | // probability that dies off the farther we are from nTarget |
176 | ||
177 | + | // The only situation left is when new_atoms > nTarget. We |
178 | + | // want to accept this with some probability that dies off the |
179 | + | // farther we are from nTarget |
180 | + | |
181 | // roughly: x = new_atoms - nTarget | |
182 | // Pacc(x) = exp(- a * x) | |
183 | < | // where a = 1 / (average atoms per molecule) |
183 | > | // where a = penalty / (average atoms per molecule) |
184 | ||
185 | x = (double) (new_atoms - nTarget); | |
186 | y = myRandom->getRandom(); | |
187 | < | |
188 | < | if (exp(- a * x) > y) { |
187 | > | |
188 | > | if (y < exp(- a * x)) { |
189 | MolToProcMap[i] = which_proc; | |
190 | AtomsPerProc[which_proc] += add_atoms; | |
191 | for (j = 0 ; j < add_atoms; j++ ) { | |
192 | < | atomIndex++; |
193 | < | AtomToProcMap[atomIndex] = which_proc; |
194 | < | } |
192 | > | AtomToProcMap[atomIndex] = which_proc; |
193 | > | atomIndex++; |
194 | > | } |
195 | done = 1; | |
196 | continue; | |
197 | } else { | |
# | Line 205 | Line 203 | int* mpiSimulation::divideLabor( void ){ | |
203 | ||
204 | // Spray out this nonsense to all other processors: | |
205 | ||
206 | < | MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, |
206 | > | MPI::COMM_WORLD.Bcast(MolToProcMap, mpiPlug->nMolGlobal, |
207 | MPI_INT, 0); | |
208 | ||
209 | < | MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, |
209 | > | MPI::COMM_WORLD.Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, |
210 | MPI_INT, 0); | |
211 | ||
212 | < | MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, |
212 | > | MPI::COMM_WORLD.Bcast(MolComponentType, mpiPlug->nMolGlobal, |
213 | MPI_INT, 0); | |
214 | ||
215 | < | MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors, |
215 | > | MPI::COMM_WORLD.Bcast(AtomsPerProc, mpiPlug->numberProcessors, |
216 | MPI_INT, 0); | |
217 | } else { | |
218 | ||
219 | // Listen to your marching orders from processor 0: | |
220 | ||
221 | < | MPI::COMM_WORLD.Bcast(&MolToProcMap, mpiPlug->nMolGlobal, |
221 | > | MPI::COMM_WORLD.Bcast(MolToProcMap, mpiPlug->nMolGlobal, |
222 | MPI_INT, 0); | |
223 | ||
224 | < | MPI::COMM_WORLD.Bcast(&AtomToProcMap, mpiPlug->nAtomsGlobal, |
224 | > | MPI::COMM_WORLD.Bcast(AtomToProcMap, mpiPlug->nAtomsGlobal, |
225 | MPI_INT, 0); | |
226 | ||
227 | < | MPI::COMM_WORLD.Bcast(&MolComponentType, mpiPlug->nMolGlobal, |
227 | > | MPI::COMM_WORLD.Bcast(MolComponentType, mpiPlug->nMolGlobal, |
228 | MPI_INT, 0); | |
229 | ||
230 | < | MPI::COMM_WORLD.Bcast(&AtomsPerProc, mpiPlug->numberProcessors, |
230 | > | MPI::COMM_WORLD.Bcast(AtomsPerProc, mpiPlug->numberProcessors, |
231 | MPI_INT, 0); | |
232 | + | |
233 | + | |
234 | } | |
235 | ||
236 | ||
# | Line 282 | Line 282 | int* mpiSimulation::divideLabor( void ){ | |
282 | local_index = 0; | |
283 | for (i = 0; i < mpiPlug->nAtomsGlobal; i++) { | |
284 | if (AtomToProcMap[i] == mpiPlug->myNode) { | |
285 | – | local_index++; |
285 | globalIndex[local_index] = i; | |
286 | + | local_index++; |
287 | } | |
288 | } | |
289 | < | |
289 | > | |
290 | return globalIndex; | |
291 | } | |
292 | ||
# | Line 296 | Line 296 | void mpiSimulation::mpiRefresh( void ){ | |
296 | int isError, i; | |
297 | int *globalIndex = new int[mpiPlug->myNlocal]; | |
298 | ||
299 | < | for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex(); |
299 | > | // Fortran indexing needs to be increased by 1 in order to get the 2 languages to |
300 | > | // not barf |
301 | ||
302 | + | for(i=0; i<mpiPlug->myNlocal; i++) globalIndex[i] = entryPlug->atoms[i]->getGlobalIndex()+1; |
303 | + | |
304 | ||
305 | isError = 0; | |
306 | setFsimParallel( mpiPlug, &(entryPlug->n_atoms), globalIndex, &isError ); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |