| 2 |  | * @file Communicator.hpp | 
| 3 |  | * @author Charles Vardeman <cvardema.at.nd.edu> | 
| 4 |  | * @date 08/18/2010 | 
| 5 | – | * @time 11:56am | 
| 5 |  | * @version 1.0 | 
| 6 |  | * | 
| 7 |  | * @section LICENSE | 
| 41 |  | * | 
| 42 |  | * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | 
| 43 |  | * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | 
| 44 | < | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | 
| 45 | < | * [4]  Vardeman & Gezelter, in progress (2009). | 
| 44 | > | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). | 
| 45 | > | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 46 | > | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 47 |  | */ | 
| 48 |  |  | 
| 49 | < | #ifndef FORCEDECOMPOSITION_COMMUNICATOR_HPP | 
| 50 | < | #define FORCEDECOMPOSITION_COMMUNICATOR_HPP | 
| 49 | > | #ifndef PARALLEL_COMMUNICATOR_HPP | 
| 50 | > | #define PARALLEL_COMMUNICATOR_HPP | 
| 51 |  |  | 
| 52 |  | #include <config.h> | 
| 53 |  | #include <mpi.h> | 
| 54 |  | #include "math/SquareMatrix3.hpp" | 
| 55 |  |  | 
| 56 | + | using namespace std; | 
| 57 |  | namespace OpenMD{ | 
| 58 |  |  | 
| 59 |  | #ifdef IS_MPI | 
| 60 |  |  | 
| 61 | < | enum direction { | 
| 62 | < | Row = 0, | 
| 63 | < | Column = 1 | 
| 61 | > | enum communicatorType { | 
| 62 | > | Global = 0, | 
| 63 | > | Row = 1, | 
| 64 | > | Column = 2 | 
| 65 |  | }; | 
| 66 |  |  | 
| 67 | < | template<typename T> | 
| 68 | < | struct MPITraits | 
| 69 | < | { | 
| 70 | < | static const MPI::Datatype datatype; | 
| 71 | < | static const int dim; | 
| 67 | > | template<class T> | 
| 68 | > | class MPITraits { | 
| 69 | > | public: | 
| 70 | > | static MPI::Datatype Type(); | 
| 71 | > | static int Length() { return 1; }; | 
| 72 |  | }; | 
| 73 |  |  | 
| 74 | < | template<> const MPI::Datatype MPITraits<int>::datatype = MPI_INT; | 
| 75 | < | template<> const int MPITraits<int>::dim = 1; | 
| 76 | < | template<> const MPI::Datatype MPITraits<RealType>::datatype = MPI_REALTYPE; | 
| 77 | < | template<> const int MPITraits<RealType>::dim = 1; | 
| 78 | < | template<> const MPI::Datatype MPITraits<Vector3d>::datatype = MPI_REALTYPE; | 
| 79 | < | template<> const int MPITraits<Vector3d>::dim = 3; | 
| 80 | < | template<> const MPI::Datatype MPITraits<Mat3x3d>::datatype = MPI_REALTYPE; | 
| 81 | < | template<> const int MPITraits<Mat3x3d>::dim = 9; | 
| 74 | > | template<> inline MPI::Datatype MPITraits<int>::Type() { return MPI::INT; } | 
| 75 | > | template<> inline MPI::Datatype MPITraits<RealType>::Type() { return MPI::REALTYPE; } | 
| 76 | > |  | 
| 77 | > | template<class T, unsigned int Dim> | 
| 78 | > | class MPITraits< Vector<T, Dim> > { | 
| 79 | > | public: | 
| 80 | > | static MPI::Datatype Type() { return MPITraits<T>::Type(); } | 
| 81 | > | static int Length() {return Dim;} | 
| 82 | > | }; | 
| 83 | > |  | 
| 84 | > | template<class T> | 
| 85 | > | class MPITraits< Vector3<T> > { | 
| 86 | > | public: | 
| 87 | > | static MPI::Datatype Type() { return MPITraits<T>::Type(); } | 
| 88 | > | static int Length() {return 3;} | 
| 89 | > | }; | 
| 90 | > |  | 
| 91 | > | template<class T, unsigned int R, unsigned int C> | 
| 92 | > | class MPITraits< RectMatrix<T, R, C> > { | 
| 93 | > | public: | 
| 94 | > | static MPI::Datatype Type() { return MPITraits<T>::Type(); } | 
| 95 | > | static int Length() {return R * C;} | 
| 96 | > | }; | 
| 97 | > |  | 
| 98 | > | template<class T> | 
| 99 | > | class MPITraits< SquareMatrix3<T> > { | 
| 100 | > | public: | 
| 101 | > | static MPI::Datatype Type() { return MPITraits<T>::Type(); } | 
| 102 | > | static int Length() {return 9;} | 
| 103 | > | }; | 
| 104 |  |  | 
| 105 | < | template<direction D, typename T> | 
| 105 | > |  | 
| 106 | > | template<communicatorType D> | 
| 107 |  | class Communicator { | 
| 108 |  | public: | 
| 109 |  |  | 
| 110 | < | Communicator<D, T>(int nObjects) { | 
| 110 | > | Communicator<D>() { | 
| 111 |  |  | 
| 112 |  | int nProc = MPI::COMM_WORLD.Get_size(); | 
| 113 |  | int myRank = MPI::COMM_WORLD.Get_rank(); | 
| 114 | < |  | 
| 114 | > |  | 
| 115 |  | int nColumnsMax = (int) sqrt(RealType(nProc)); | 
| 116 |  |  | 
| 117 |  | int nColumns; | 
| 119 |  | if (nProc % i == 0) nColumns = i; | 
| 120 |  | } | 
| 121 |  |  | 
| 122 | < | int nRows = nProc / nColumns; | 
| 122 | > | // int nRows = nProc / nColumns; | 
| 123 |  | rowIndex_ = myRank / nColumns; | 
| 124 |  | columnIndex_ = myRank % nColumns; | 
| 125 |  |  | 
| 126 | < | if (D == Row) { | 
| 126 | > | switch(D) { | 
| 127 | > | case Row : | 
| 128 |  | myComm = MPI::COMM_WORLD.Split(rowIndex_, 0); | 
| 129 | < | } else { | 
| 129 | > | break; | 
| 130 | > | case Column: | 
| 131 |  | myComm = MPI::COMM_WORLD.Split(columnIndex_, 0); | 
| 132 | + | break; | 
| 133 | + | case Global: | 
| 134 | + | myComm = MPI::COMM_WORLD.Split(myRank, 0); | 
| 135 |  | } | 
| 106 | – |  | 
| 107 | – | int nCommProcs = myComm.Get_size(); | 
| 136 |  |  | 
| 137 | < | counts.reserve(nCommProcs); | 
| 138 | < | displacements.reserve(nCommProcs); | 
| 137 | > | } | 
| 138 | > |  | 
| 139 | > | MPI::Intracomm getComm() { return myComm; } | 
| 140 | > |  | 
| 141 | > | private: | 
| 142 | > | int rowIndex_; | 
| 143 | > | int columnIndex_; | 
| 144 | > | MPI::Intracomm myComm; | 
| 145 | > | }; | 
| 146 | > |  | 
| 147 |  |  | 
| 148 | < | planSize_ = MPITraits<T>::dim * nObjects; | 
| 149 | < |  | 
| 148 | > | template<typename T> | 
| 149 | > | class Plan { | 
| 150 | > | public: | 
| 151 | > |  | 
| 152 | > | Plan<T>(MPI::Intracomm comm, int nObjects) : myComm(comm) { | 
| 153 | > | int nCommProcs = myComm.Get_size(); | 
| 154 | > |  | 
| 155 | > | counts.resize(nCommProcs, 0); | 
| 156 | > | displacements.resize(nCommProcs, 0); | 
| 157 | > |  | 
| 158 | > | planSize_ = MPITraits<T>::Length() * nObjects; | 
| 159 | > |  | 
| 160 |  | myComm.Allgather(&planSize_, 1, MPI::INT, &counts[0], 1, MPI::INT); | 
| 161 | < |  | 
| 161 | > |  | 
| 162 |  | displacements[0] = 0; | 
| 163 |  | for (int i = 1; i < nCommProcs; i++) { | 
| 164 |  | displacements[i] = displacements[i-1] + counts[i-1]; | 
| 165 | < | } | 
| 165 | > | } | 
| 166 | > |  | 
| 167 | > | size_ = 0; | 
| 168 | > | for (int i = 0; i < nCommProcs; i++) { | 
| 169 | > | size_ += counts[i]; | 
| 170 | > | } | 
| 171 |  | } | 
| 172 |  |  | 
| 173 | < |  | 
| 174 | < | void gather(std::vector<T>& v1, std::vector<T>& v2) { | 
| 173 | > |  | 
| 174 | > | void gather(vector<T>& v1, vector<T>& v2) { | 
| 175 |  |  | 
| 176 | + | // an assert would be helpful here to make sure the vectors are the | 
| 177 | + | // correct geometry | 
| 178 | + |  | 
| 179 |  | myComm.Allgatherv(&v1[0], | 
| 180 |  | planSize_, | 
| 181 | < | MPITraits<T>::datatype, | 
| 181 | > | MPITraits<T>::Type(), | 
| 182 |  | &v2[0], | 
| 183 |  | &counts[0], | 
| 184 |  | &displacements[0], | 
| 185 | < | MPITraits<T>::datatype); | 
| 186 | < | } | 
| 133 | < |  | 
| 185 | > | MPITraits<T>::Type()); | 
| 186 | > | } | 
| 187 |  |  | 
| 188 | < |  | 
| 189 | < | void scatter(std::vector<T>& v1, std::vector<T>& v2) { | 
| 190 | < |  | 
| 188 | > | void scatter(vector<T>& v1, vector<T>& v2) { | 
| 189 | > | // an assert would be helpful here to make sure the vectors are the | 
| 190 | > | // correct geometry | 
| 191 | > |  | 
| 192 |  | myComm.Reduce_scatter(&v1[0], &v2[0], &counts[0], | 
| 193 | < | MPITraits<T>::datatype, MPI::SUM); | 
| 193 | > | MPITraits<T>::Type(), MPI::SUM); | 
| 194 |  | } | 
| 195 |  |  | 
| 196 | + | int getSize() { | 
| 197 | + | return size_; | 
| 198 | + | } | 
| 199 | + |  | 
| 200 |  | private: | 
| 201 |  | int planSize_;     ///< how many are on local proc | 
| 202 | < | int rowIndex_; | 
| 203 | < | int columnIndex_; | 
| 204 | < | std::vector<int> counts; | 
| 147 | < | std::vector<int> displacements; | 
| 202 | > | int size_; | 
| 203 | > | vector<int> counts; | 
| 204 | > | vector<int> displacements; | 
| 205 |  | MPI::Intracomm myComm; | 
| 206 |  | }; | 
| 207 |  |  |