| 46 |  | * [4]  Vardeman & Gezelter, in progress (2009). | 
| 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 Row, unsigned int Col> | 
| 92 | > | class MPITraits< RectMatrix<T, Row, Col> > { | 
| 93 | > | public: | 
| 94 | > | static MPI::Datatype Type() { return MPITraits<T>::Type(); } | 
| 95 | > | static int Length() {return Row * Col;} | 
| 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; | 
| 123 |  | rowIndex_ = myRank / nColumns; | 
| 124 |  | columnIndex_ = myRank % nColumns; | 
| 125 |  |  | 
| 126 | < | if (D == Row) { | 
| 126 | > | cerr << "rowIndex = " << rowIndex_ << "\t colIndex = " << columnIndex_ << "\n"; | 
| 127 | > | switch(D) { | 
| 128 | > | case Row : | 
| 129 |  | myComm = MPI::COMM_WORLD.Split(rowIndex_, 0); | 
| 130 | < | } else { | 
| 130 | > | break; | 
| 131 | > | case Column: | 
| 132 |  | myComm = MPI::COMM_WORLD.Split(columnIndex_, 0); | 
| 133 | + | break; | 
| 134 | + | case Global: | 
| 135 | + | myComm = MPI::COMM_WORLD.Split(myRank, 0); | 
| 136 |  | } | 
| 106 | – |  | 
| 107 | – | int nCommProcs = myComm.Get_size(); | 
| 137 |  |  | 
| 138 | < | counts.reserve(nCommProcs); | 
| 139 | < | displacements.reserve(nCommProcs); | 
| 138 | > | } | 
| 139 | > |  | 
| 140 | > | MPI::Intracomm getComm() { return myComm; } | 
| 141 | > |  | 
| 142 | > | private: | 
| 143 | > | int rowIndex_; | 
| 144 | > | int columnIndex_; | 
| 145 | > | MPI::Intracomm myComm; | 
| 146 | > | }; | 
| 147 | > |  | 
| 148 |  |  | 
| 149 | < | planSize_ = MPITraits<T>::dim * nObjects; | 
| 150 | < |  | 
| 149 | > | template<typename T> | 
| 150 | > | class Plan { | 
| 151 | > | public: | 
| 152 | > |  | 
| 153 | > | Plan<T>(MPI::Intracomm comm, int nObjects) { | 
| 154 | > | myComm = comm; | 
| 155 | > | int nCommProcs = myComm.Get_size(); | 
| 156 | > |  | 
| 157 | > | counts.resize(nCommProcs, 0); | 
| 158 | > | displacements.resize(nCommProcs, 0); | 
| 159 | > |  | 
| 160 | > | planSize_ = MPITraits<T>::Length() * nObjects; | 
| 161 | > |  | 
| 162 |  | myComm.Allgather(&planSize_, 1, MPI::INT, &counts[0], 1, MPI::INT); | 
| 163 | < |  | 
| 116 | < |  | 
| 163 | > |  | 
| 164 |  | displacements[0] = 0; | 
| 165 |  | for (int i = 1; i < nCommProcs; i++) { | 
| 166 |  | displacements[i] = displacements[i-1] + counts[i-1]; | 
| 120 | – | size_ += count[i-1]; | 
| 167 |  | } | 
| 168 |  |  | 
| 169 |  | size_ = 0; | 
| 172 |  | } | 
| 173 |  | } | 
| 174 |  |  | 
| 175 | < |  | 
| 176 | < | void gather(std::vector<T>& v1, std::vector<T>& v2) { | 
| 175 | > |  | 
| 176 | > | void gather(vector<T>& v1, vector<T>& v2) { | 
| 177 |  |  | 
| 178 | + | // an assert would be helpful here to make sure the vectors are the | 
| 179 | + | // correct geometry | 
| 180 | + |  | 
| 181 |  | myComm.Allgatherv(&v1[0], | 
| 182 |  | planSize_, | 
| 183 | < | MPITraits<T>::datatype, | 
| 183 | > | MPITraits<T>::Type(), | 
| 184 |  | &v2[0], | 
| 185 |  | &counts[0], | 
| 186 |  | &displacements[0], | 
| 187 | < | MPITraits<T>::datatype); | 
| 188 | < | } | 
| 140 | < |  | 
| 187 | > | MPITraits<T>::Type()); | 
| 188 | > | } | 
| 189 |  |  | 
| 190 | < |  | 
| 191 | < | void scatter(std::vector<T>& v1, std::vector<T>& v2) { | 
| 192 | < |  | 
| 190 | > | void scatter(vector<T>& v1, vector<T>& v2) { | 
| 191 | > | // an assert would be helpful here to make sure the vectors are the | 
| 192 | > | // correct geometry | 
| 193 | > |  | 
| 194 |  | myComm.Reduce_scatter(&v1[0], &v2[0], &counts[0], | 
| 195 | < | MPITraits<T>::datatype, MPI::SUM); | 
| 195 | > | MPITraits<T>::Type(), MPI::SUM); | 
| 196 |  | } | 
| 197 | < |  | 
| 197 | > |  | 
| 198 |  | int getSize() { | 
| 199 |  | return size_; | 
| 200 |  | } | 
| 201 |  |  | 
| 202 |  | private: | 
| 203 |  | int planSize_;     ///< how many are on local proc | 
| 155 | – | int rowIndex_; | 
| 156 | – | int columnIndex_; | 
| 204 |  | int size_; | 
| 205 | < | std::vector<int> counts; | 
| 206 | < | std::vector<int> displacements; | 
| 205 | > | vector<int> counts; | 
| 206 | > | vector<int> displacements; | 
| 207 |  | MPI::Intracomm myComm; | 
| 208 |  | }; | 
| 209 |  |  |