--- branches/development/src/parallel/Communicator.hpp	2011/03/18 19:31:52	1544
+++ trunk/src/parallel/Communicator.hpp	2012/09/10 18:38:44	1796
@@ -43,7 +43,8 @@
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 
 #ifndef PARALLEL_COMMUNICATOR_HPP
@@ -53,40 +54,65 @@
 #include <mpi.h>
 #include "math/SquareMatrix3.hpp"
 
+using namespace std;
 namespace OpenMD{
   
 #ifdef IS_MPI
 
-  enum direction {
-    Row = 0,
-    Column = 1
+  enum communicatorType {
+    Global = 0,
+    Row = 1,
+    Column = 2
   };
     
-  template<typename T> 
-  struct MPITraits
-  {
-    static const MPI::Datatype datatype;
-    static const int dim;
+  template<class T>
+  class MPITraits {
+  public:
+    static MPI::Datatype Type();
+    static int Length() { return 1; };
   };
   
-  template<> const MPI::Datatype MPITraits<int>::datatype = MPI_INT;
-  template<> const int MPITraits<int>::dim = 1;
-  template<> const MPI::Datatype MPITraits<RealType>::datatype = MPI_REALTYPE;
-  template<> const int MPITraits<RealType>::dim = 1;
-  template<> const MPI::Datatype MPITraits<Vector3d>::datatype = MPI_REALTYPE;
-  template<> const int MPITraits<Vector3d>::dim = 3;
-  template<> const MPI::Datatype MPITraits<Mat3x3d>::datatype = MPI_REALTYPE;
-  template<> const int MPITraits<Mat3x3d>::dim = 9;
+  template<> inline MPI::Datatype MPITraits<int>::Type() { return MPI::INT; }
+  template<> inline MPI::Datatype MPITraits<RealType>::Type() { return MPI::REALTYPE; }
+
+  template<class T, unsigned int Dim>
+  class MPITraits< Vector<T, Dim> > {
+  public:
+    static MPI::Datatype Type() { return MPITraits<T>::Type(); }
+    static int Length() {return Dim;}
+  };
+
+  template<class T>
+  class MPITraits< Vector3<T> > {
+  public:
+    static MPI::Datatype Type() { return MPITraits<T>::Type(); }
+    static int Length() {return 3;}
+  };
+
+  template<class T, unsigned int Row, unsigned int Col>
+  class MPITraits< RectMatrix<T, Row, Col> > {
+  public:
+    static MPI::Datatype Type() { return MPITraits<T>::Type(); }
+    static int Length() {return Row * Col;}
+  };
+
+  template<class T>
+  class MPITraits< SquareMatrix3<T> > {
+  public:
+    static MPI::Datatype Type() { return MPITraits<T>::Type(); }
+    static int Length() {return 9;}
+  };
   
-  template<direction D, typename T>
+  
+  template<communicatorType D>
   class Communicator { 
   public: 
     
-    Communicator<D, T>(int nObjects) {
+    Communicator<D>() {
       
       int nProc = MPI::COMM_WORLD.Get_size();
       int myRank = MPI::COMM_WORLD.Get_rank();
-
+      
       int nColumnsMax = (int) sqrt(RealType(nProc));
 
       int nColumns;
@@ -94,29 +120,49 @@ namespace OpenMD{
         if (nProc % i == 0) nColumns = i;        
       }
         
-      int nRows = nProc / nColumns;
       rowIndex_ = myRank / nColumns;      
       columnIndex_ = myRank % nColumns;
 
-      if (D == Row) {
+      switch(D) {
+      case Row :
         myComm = MPI::COMM_WORLD.Split(rowIndex_, 0);
-      } else {
+        break;
+      case Column:
         myComm = MPI::COMM_WORLD.Split(columnIndex_, 0);
+        break;
+      case Global:
+        myComm = MPI::COMM_WORLD.Split(myRank, 0);
       }
-         
-      int nCommProcs = myComm.Get_size();
 
-      counts.reserve(nCommProcs);
-      displacements.reserve(nCommProcs);
+    }
+    
+    MPI::Intracomm getComm() { return myComm; }
+    
+  private:
+    int rowIndex_;
+    int columnIndex_;
+    MPI::Intracomm myComm;
+  };
+  
 
-      planSize_ = MPITraits<T>::dim * nObjects; 
-
+  template<typename T>
+  class Plan {
+  public:
+    
+    Plan<T>(MPI::Intracomm comm, int nObjects) {
+      myComm = comm;
+      int nCommProcs = myComm.Get_size();
+      
+      counts.resize(nCommProcs, 0);
+      displacements.resize(nCommProcs, 0);
+      
+      planSize_ = MPITraits<T>::Length() * nObjects; 
+      
       myComm.Allgather(&planSize_, 1, MPI::INT, &counts[0], 1, MPI::INT);
-
+      
       displacements[0] = 0;
       for (int i = 1; i < nCommProcs; i++) {
         displacements[i] = displacements[i-1] + counts[i-1];
-        size_ += counts[i-1];
       }
 
       size_ = 0;
@@ -125,37 +171,38 @@ namespace OpenMD{
       }
     }
 
-
-    void gather(std::vector<T>& v1, std::vector<T>& v2) {
+    
+    void gather(vector<T>& v1, vector<T>& v2) {
       
+      // an assert would be helpful here to make sure the vectors are the
+      // correct geometry
+      
       myComm.Allgatherv(&v1[0], 
                         planSize_, 
-                        MPITraits<T>::datatype, 
+                        MPITraits<T>::Type(), 
                         &v2[0], 
                         &counts[0], 
                         &displacements[0], 
-                        MPITraits<T>::datatype);      
-    }
-
+                        MPITraits<T>::Type());      
+    }       
     
-   
-    void scatter(std::vector<T>& v1, std::vector<T>& v2) {
-
+    void scatter(vector<T>& v1, vector<T>& v2) {
+      // an assert would be helpful here to make sure the vectors are the
+      // correct geometry
+            
       myComm.Reduce_scatter(&v1[0], &v2[0], &counts[0], 
-                            MPITraits<T>::datatype, MPI::SUM);
+                            MPITraits<T>::Type(), MPI::SUM);
     }
-
+    
     int getSize() {
       return size_;
     }
     
   private:
     int planSize_;     ///< how many are on local proc
-    int rowIndex_;
-    int columnIndex_;
     int size_;
-    std::vector<int> counts;
-    std::vector<int> displacements;
+    vector<int> counts;
+    vector<int> displacements;
     MPI::Intracomm myComm;
   };