--- trunk/mdtools/mpi_implementation/mpiSimulation.cpp	2002/09/30 19:36:20	122
+++ trunk/mdtools/mpi_implementation/mpiSimulation.cpp	2003/01/30 15:20:21	253
@@ -1,25 +1,230 @@
-#include <mpi++.h>
-#include <mpiSimulation.hpp>
+#include <cstdlib>
+#include <cstring>
+#include <mpi.h>
 
-mpiSimulation::mpiSimulation()
+#include "mpiSimulation.hpp"
+#include "simError.h"
+
+extern "C"{
+  void wrapsimparallelmod_( void (*wrapFunction)(void (*fSub)( mpiSimData*,
+							       int*, int*, 
+							       int*)));
+}
+
+void wrapSimParallel(void (*fSub)(mpiSimData*, int*, int*, int*));
+
+
+mpiSimulation* mpiSim;
+
+mpiSimulation::mpiSimulation(SimInfo* the_entryPlug)
 {
-  int mpi_error;
+  entryPlug = the_entryPlug;
+  mpiPlug = new mpiSimData;
+  
+  mpiPlug->numberProcessors = MPI::COMM_WORLD.Get_size();
+  mpiPlug->myNode = worldRank;
+  
+  mpiSim = this;
+  wrapMe();
 
-  MPI::Init();
+}
 
-  numberProcessors = MPI::Comm::Get_size();
-  myNode = MPI::Comm::Get_rank();
-  MPI::Get_processor_name(processorName,&processorNameLen);
+
+mpiSimulation::~mpiSimulation(){
+  
+  delete mpiPlug;
+  // perhaps we should let fortran know the party is over.
+  
 }
 
-mpiSimulation::mpiInitSimulation(SimInfo* entry_plug)
-{
+void mpiSimulation::wrapMe(){
 
-  // need to get nmol here......
+  wrapsimparallelmod_( wrapSimParallel );
+}
 
 
 
-  myMolStart = nint(float(node)/numberProcessors*entry_plug->n_mol);
-  myMolEnd   = nint(float(node + 1)/numberProcessors*entry_plug->n_mol;);
-  nMolLocal  = myMolEnd - myMolStart + 1
+int* mpiSimulation::divideLabor( void ){
+
+  int* globalIndex;
+
+  int nComponents;
+  MoleculeStamp** compStamps;
+  int* componentsNmol;
+
+  double numerator;
+  double denominator;
+  double precast;
+
+  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;
+
+  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;
+
+  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( i == mpiPlug->myNode ){
+      mpiPlug->myMolStart = molIndex;
+      mpiPlug->myAtomStart = atomIndex;
+    }
+    
+    while( !done ){
+      
+      if( (molIndex-compStart) >= componentsNmol[compIndex] ){
+	compStart = molIndex;
+	compIndex++;
+	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( i == mpiPlug->myNode ){
+      mpiPlug->myMolEnd = (molIndex - 1);
+      mpiPlug->myAtomEnd = (atomIndex - 1);
+      mpiPlug->myNlocal = nLocal;
+      mpiPlug->myMol = molLocal;
+    }
+    
+    numerator = (double)( entryPlug->n_atoms - atomIndex );
+    denominator = (double)( mpiPlug->numberProcessors - (i+1) );
+    precast = numerator / denominator;
+    nTarget = (int)( precast + 0.5 );
+  }
+  
+  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;
+  }
+
+
+  MPI_Allreduce( &nLocal, &testSum, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
+  
+  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();
+    }
+  }
+
+  sprintf( checkPointMsg,
+	   "Successfully divided the molecules among the processors.\n" );
+  MPIcheckPoint();
+
+  // lets create the identity array
+
+  globalIndex = new int[mpiPlug->myNlocal];
+  index = mpiPlug->myAtomStart;
+  for( i=0; i<mpiPlug->myNlocal; i++){
+    globalIndex[i] = index;
+    index++;
+  }
+
+  return globalIndex;
 }
+
+
+void wrapSimParallel(void (*fSub)(mpiSimData*, int*, int*, int*)){
+  
+  mpiSim->setInternal( fSub );
+}
+
+
+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();
+}
+