/[escript]/branches/trilinos_from_5897/dudley/src/NodeFile_gather.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/NodeFile_gather.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 /****************************************************************************  
  *  
  *   Dudley: Mesh: NodeFile  
  *   gathers the NodeFile out from the NodeFile in using the entries  
  *   in index[0:out->numNodes-1] which are between min_index and max_index (exclusive)  
  *   the node index[i]  
  *  
  ****************************************************************************/  
   
17  #include "NodeFile.h"  #include "NodeFile.h"
18    
19    using escript::DataTypes::real_t;
20    
21  namespace dudley {  namespace dudley {
22    
23  void Dudley_NodeFile_gatherEntries(dim_t n, index_t* index, index_t min_index, index_t max_index,  // helper function
24                                     index_t* Id_out, index_t* Id_in,  static void gatherEntries(dim_t n, const index_t* index,
25                                     index_t* Tag_out, index_t* Tag_in,                            index_t min_index, index_t max_index,
26                                     index_t* globalDegreesOfFreedom_out, index_t* globalDegreesOfFreedom_in,                            index_t* Id_out, const index_t* Id_in,
27                                     dim_t numDim, double* Coordinates_out, double* Coordinates_in)                            int* Tag_out, const int* Tag_in,
28                              index_t* globalDegreesOfFreedom_out,
29                              const index_t* globalDegreesOfFreedom_in,
30                              int numDim, real_t* Coordinates_out,
31                              const real_t* Coordinates_in)
32  {  {
33      dim_t i;      const dim_t range = max_index - min_index;
34      index_t k;      const size_t numDim_size = numDim * sizeof(real_t);
35      const index_t range = max_index - min_index;  #pragma omp parallel for
36      const size_t numDim_size = (size_t) numDim * sizeof(double);      for (index_t i = 0; i < n; i++) {
37  #pragma omp parallel for private(i,k) schedule(static)          const index_t k = index[i] - min_index;
38      for (i = 0; i < n; i++)          if (k >= 0 && k < range) {
     {  
         k = index[i] - min_index;  
         if ((k >= 0) && (k < range))  
         {  
39              Id_out[i] = Id_in[k];              Id_out[i] = Id_in[k];
40              Tag_out[i] = Tag_in[k];              Tag_out[i] = Tag_in[k];
41              globalDegreesOfFreedom_out[i] = globalDegreesOfFreedom_in[k];              globalDegreesOfFreedom_out[i] = globalDegreesOfFreedom_in[k];
42              memcpy(&(Coordinates_out[INDEX2(0, i, numDim)]), &(Coordinates_in[INDEX2(0, k, numDim)]), numDim_size);              memcpy(&Coordinates_out[INDEX2(0, i, numDim)],
43                       &Coordinates_in[INDEX2(0, k, numDim)], numDim_size);
44          }          }
45      }      }
46  }  }
47    
48  void Dudley_NodeFile_gather(index_t * index, Dudley_NodeFile * in, Dudley_NodeFile * out)  // helper function
49    static void scatterEntries(dim_t n, const index_t* index,
50                               index_t min_index, index_t max_index,
51                               index_t* Id_out, const index_t* Id_in,
52                               int* Tag_out, const int* Tag_in,
53                               index_t* globalDegreesOfFreedom_out,
54                               const index_t* globalDegreesOfFreedom_in,
55                               int numDim, real_t* Coordinates_out,
56                               const real_t* Coordinates_in)
57  {  {
58      index_t min_id, max_id;      const dim_t range = max_index - min_index;
59      Dudley_NodeFile_setGlobalIdRange(&min_id, &max_id, in);      const size_t numDim_size = numDim * sizeof(real_t);
60      Dudley_NodeFile_gatherEntries(out->numNodes, index, min_id, max_id,  
61                                    out->Id, in->Id,  #pragma omp parallel for
62                                    out->Tag, in->Tag,      for (index_t i = 0; i < n; i++) {
63                                    out->globalDegreesOfFreedom, in->globalDegreesOfFreedom,          const index_t k = index[i] - min_index;
64                                    out->numDim, out->Coordinates, in->Coordinates);          if (k >= 0 && k < range) {
65                Id_out[k] = Id_in[i];
66                Tag_out[k] = Tag_in[i];
67                globalDegreesOfFreedom_out[k] = globalDegreesOfFreedom_in[i];
68                memcpy(&Coordinates_out[INDEX2(0, k, numDim)],
69                       &Coordinates_in[INDEX2(0, i, numDim)], numDim_size);
70            }
71        }
72  }  }
73    
74  void Dudley_NodeFile_gather_global(index_t * index, Dudley_NodeFile * in, Dudley_NodeFile * out)  void NodeFile::gather(const index_t* index, const NodeFile* in)
75  {  {
76      index_t min_id, max_id, undefined_node;      const std::pair<index_t,index_t> idRange(in->getGlobalIdRange());
77      int buffer_rank, *distribution = NULL;      gatherEntries(numNodes, index, idRange.first, idRange.second, Id, in->Id,
78      index_t *Id_buffer = NULL, *Tag_buffer = NULL, *globalDegreesOfFreedom_buffer = NULL;               Tag, in->Tag, globalDegreesOfFreedom, in->globalDegreesOfFreedom,
79      double *Coordinates_buffer = NULL;               numDim, Coordinates, in->Coordinates);
80      dim_t p, buffer_len, n;  }
     char error_msg[100];  
 #ifdef ESYS_MPI  
     int dest, source;  
     MPI_Status status;  
 #endif  
81    
82      /* get the global range of node ids */  void NodeFile::gather_global(const index_t* index, const NodeFile* in)
83      Dudley_NodeFile_setGlobalIdRange(&min_id, &max_id, in);  {
84      undefined_node = min_id - 1;      // get the global range of node IDs
85        const std::pair<index_t,index_t> idRange(in->getGlobalIdRange());
86      distribution = new index_t[in->MPIInfo->size + 1];      const index_t UNDEFINED = idRange.first - 1;
87        std::vector<index_t> distribution(in->MPIInfo->size + 1);
88      /* distribute the range of node ids */  
89      buffer_len = in->MPIInfo->setDistribution(min_id, max_id, distribution);      // distribute the range of node IDs
90      /* allocate buffers */      dim_t buffer_len = MPIInfo->setDistribution(idRange.first, idRange.second,
91      Id_buffer = new  index_t[buffer_len];                                                  &distribution[0]);
92      Tag_buffer = new  index_t[buffer_len];  
93      globalDegreesOfFreedom_buffer = new  index_t[buffer_len];      // allocate buffers
94      Coordinates_buffer = new  double[buffer_len * out->numDim];      index_t* Id_buffer = new index_t[buffer_len];
95      /* fill Id_buffer by the undefined_node marker to check if nodes are defined */      int* Tag_buffer = new int[buffer_len];
96  #pragma omp parallel for private(n) schedule(static)      index_t* globalDegreesOfFreedom_buffer = new index_t[buffer_len];
97      for (n = 0; n < buffer_len; n++)      real_t* Coordinates_buffer = new real_t[buffer_len * numDim];
98          Id_buffer[n] = undefined_node;  
99        // fill Id_buffer by the UNDEFINED marker to check if nodes are
100        // defined
101    #pragma omp parallel for
102        for (index_t n = 0; n < buffer_len; n++)
103            Id_buffer[n] = UNDEFINED;
104    
105      /* fill the buffer by sending portions around in a circle */      // fill the buffer by sending portions around in a circle
106  #ifdef ESYS_MPI  #ifdef ESYS_MPI
107      dest = in->MPIInfo->mod_rank(in->MPIInfo->rank + 1);      MPI_Status status;
108      source = in->MPIInfo->mod_rank(in->MPIInfo->rank - 1);      int dest = MPIInfo->mod_rank(MPIInfo->rank + 1);
109        int source = MPIInfo->mod_rank(MPIInfo->rank - 1);
110  #endif  #endif
111      buffer_rank = in->MPIInfo->rank;      int buffer_rank = MPIInfo->rank;
112      for (p = 0; p < in->MPIInfo->size; ++p)      for (int p = 0; p < MPIInfo->size; ++p) {
     {  
         if (p > 0)  
         {               /* the initial send can be skipped */  
113  #ifdef ESYS_MPI  #ifdef ESYS_MPI
114              MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_INT,          if (p > 0) { // the initial send can be skipped
115                                   dest, in->MPIInfo->counter(), source, in->MPIInfo->counter(),              MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_DIM_T, dest,
116                                   in->MPIInfo->comm, &status);                          MPIInfo->counter(), source, MPIInfo->counter(),
117              MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT,                          MPIInfo->comm, &status);
118                                   dest, in->MPIInfo->counter() + 1, source,              MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT, dest,
119                                   in->MPIInfo->counter() + 1, in->MPIInfo->comm, &status);                          MPIInfo->counter() + 1, source,
120              MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len, MPI_INT, dest,                          MPIInfo->counter() + 1, MPIInfo->comm, &status);
121                                   in->MPIInfo->counter() + 2, source, in->MPIInfo->counter() + 2,              MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len,
122                                   in->MPIInfo->comm, &status);                          MPI_DIM_T, dest, MPIInfo->counter() + 2, source,
123              MPI_Sendrecv_replace(Coordinates_buffer, buffer_len * out->numDim, MPI_DOUBLE, dest,                          MPIInfo->counter() + 2, MPIInfo->comm, &status);
124                                   in->MPIInfo->counter() + 3, source, in->MPIInfo->counter() + 3,              MPI_Sendrecv_replace(Coordinates_buffer, buffer_len * numDim,
125                                   in->MPIInfo->comm, &status);                          MPI_DOUBLE, dest, MPIInfo->counter() + 3, source,
126              in->MPIInfo->incCounter(4);                          MPIInfo->counter() + 3, MPIInfo->comm, &status);
127  #endif              MPIInfo->incCounter(4);
128          }          }
129          buffer_rank = in->MPIInfo->mod_rank(buffer_rank - 1);  #endif
130          Dudley_NodeFile_scatterEntries(in->numNodes, in->Id,          buffer_rank = MPIInfo->mod_rank(buffer_rank - 1);
131                                         distribution[buffer_rank], distribution[buffer_rank + 1],          scatterEntries(in->numNodes, in->Id, distribution[buffer_rank],
132                                         Id_buffer, in->Id,                         distribution[buffer_rank + 1], Id_buffer, in->Id,
133                                         Tag_buffer, in->Tag,                         Tag_buffer, in->Tag, globalDegreesOfFreedom_buffer,
134                                         globalDegreesOfFreedom_buffer, in->globalDegreesOfFreedom,                         in->globalDegreesOfFreedom, numDim, Coordinates_buffer,
135                                         out->numDim, Coordinates_buffer, in->Coordinates);                         in->Coordinates);
136      }      }
137      /* now entries are collected from the buffer again by sending the entries around in a circle */      // now entries are collected from the buffer again by sending the entries
138        // around in a circle
139  #ifdef ESYS_MPI  #ifdef ESYS_MPI
140      dest = in->MPIInfo->mod_rank(in->MPIInfo->rank + 1);      dest = MPIInfo->mod_rank(MPIInfo->rank + 1);
141      source = in->MPIInfo->mod_rank(in->MPIInfo->rank - 1);      source = MPIInfo->mod_rank(MPIInfo->rank - 1);
142  #endif  #endif
143      buffer_rank = in->MPIInfo->rank;      buffer_rank = MPIInfo->rank;
144      for (p = 0; p < in->MPIInfo->size; ++p)      for (int p = 0; p < MPIInfo->size; ++p) {
145      {          gatherEntries(numNodes, index, distribution[buffer_rank],
146          Dudley_NodeFile_gatherEntries(out->numNodes, index,                        distribution[buffer_rank + 1], Id, Id_buffer,
147                                        distribution[buffer_rank], distribution[buffer_rank + 1],                        Tag, Tag_buffer, globalDegreesOfFreedom,
148                                        out->Id, Id_buffer,                        globalDegreesOfFreedom_buffer, numDim,
149                                        out->Tag, Tag_buffer,                        Coordinates, Coordinates_buffer);
                                       out->globalDegreesOfFreedom, globalDegreesOfFreedom_buffer,  
                                       out->numDim, out->Coordinates, Coordinates_buffer);  
         if (p < in->MPIInfo->size - 1)  
         {               /* the last send can be skipped */  
150  #ifdef ESYS_MPI  #ifdef ESYS_MPI
151              MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_INT,          if (p < MPIInfo->size - 1) { // the last send can be skipped
152                                   dest, in->MPIInfo->counter(), source, in->MPIInfo->counter(),              MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_DIM_T, dest,
153                                   in->MPIInfo->comm, &status);                                MPIInfo->counter(), source,
154              MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT,                                MPIInfo->counter(), MPIInfo->comm, &status);
155                                   dest, in->MPIInfo->counter() + 1, source,              MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT, dest,
156                                   in->MPIInfo->counter() + 1, in->MPIInfo->comm, &status);                                MPIInfo->counter() + 1, source,
157              MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len, MPI_INT, dest,                                MPIInfo->counter() + 1, MPIInfo->comm, &status);
158                                   in->MPIInfo->counter() + 2, source, in->MPIInfo->counter() + 2,              MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len,
159                                   in->MPIInfo->comm, &status);                                MPI_DIM_T, dest, MPIInfo->counter() + 2, source,
160              MPI_Sendrecv_replace(Coordinates_buffer, buffer_len * out->numDim, MPI_DOUBLE, dest,                                MPIInfo->counter() + 2, MPIInfo->comm, &status);
161                                   in->MPIInfo->counter() + 3, source, in->MPIInfo->counter() + 3,              MPI_Sendrecv_replace(Coordinates_buffer, buffer_len * numDim,
162                                   in->MPIInfo->comm, &status);                                MPI_DOUBLE, dest, MPIInfo->counter() + 3, source,
163              in->MPIInfo->incCounter(4);                                MPIInfo->counter() + 3, MPIInfo->comm, &status);
164  #endif              MPIInfo->incCounter(4);
         }  
         buffer_rank = in->MPIInfo->mod_rank(buffer_rank - 1);  
     }  
     /* check if all nodes are set: */  
 #pragma omp parallel for private(n) schedule(static)  
     for (n = 0; n < out->numNodes; ++n)  
     {  
         if (out->Id[n] == undefined_node)  
         {  
             sprintf(error_msg,  
                     "Dudley_NodeFile_gather_global: Node id %d at position %d is referenced but is not defined.",  
                     out->Id[n], n);  
             throw DudleyException(error_msg);  
165          }          }
166    #endif
167            buffer_rank = MPIInfo->mod_rank(buffer_rank - 1);
168      }      }
169      delete[] Id_buffer;      delete[] Id_buffer;
170      delete[] Tag_buffer;      delete[] Tag_buffer;
171      delete[] globalDegreesOfFreedom_buffer;      delete[] globalDegreesOfFreedom_buffer;
172      delete[] Coordinates_buffer;      delete[] Coordinates_buffer;
173      delete[] distribution;  #if DOASSERT
174        // check if all nodes are set
175        index_t err = -1;
176    #pragma omp parallel for
177        for (index_t n = 0; n < numNodes; ++n) {
178            if (Id[n] == UNDEFINED) {
179    #pragma omp critical
180                err = n;
181            }
182        }
183        if (err >= 0) {
184            std::stringstream ss;
185            ss << "NodeFile::gather_global: Node id " << Id[err]
186                << " at position " << err << " is referenced but not defined.";
187            throw escript::AssertException(ss.str());
188        }
189    #endif // DOASSERT
190  }  }
191    
192  } // namespace dudley  } // namespace dudley

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26