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

Contents of /branches/trilinos_from_5897/dudley/src/NodeFile_createDenseLabelings.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 9075 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "NodeFile.h"
18
19 namespace dudley {
20
21 dim_t NodeFile::createDenseDOFLabeling()
22 {
23 const index_t UNSET_ID = -1, SET_ID = 1;
24
25 // get the global range of DOF IDs
26 const std::pair<index_t,index_t> idRange(getGlobalDOFRange());
27
28 // distribute the range of DOF IDs
29 std::vector<index_t> distribution(MPIInfo->size + 1);
30 dim_t bufferLen = MPIInfo->setDistribution(idRange.first, idRange.second,
31 &distribution[0]);
32 const dim_t myDOFs = distribution[MPIInfo->rank + 1] - distribution[MPIInfo->rank];
33
34 index_t* DOF_buffer = new index_t[bufferLen];
35 // fill buffer by the UNSET_ID marker to check if nodes are defined
36 #pragma omp parallel for
37 for (index_t n = 0; n < bufferLen; n++)
38 DOF_buffer[n] = UNSET_ID;
39
40 // fill the buffer by sending portions around in a circle
41 #ifdef ESYS_MPI
42 MPI_Status status;
43 int dest = MPIInfo->mod_rank(MPIInfo->rank + 1);
44 int source = MPIInfo->mod_rank(MPIInfo->rank - 1);
45 #endif
46 int buffer_rank = MPIInfo->rank;
47 for (int p = 0; p < MPIInfo->size; ++p) {
48 #ifdef ESYS_MPI
49 if (p > 0) { // the initial send can be skipped
50 MPI_Sendrecv_replace(DOF_buffer, bufferLen, MPI_DIM_T, dest,
51 MPIInfo->counter(), source, MPIInfo->counter(),
52 MPIInfo->comm, &status);
53 MPIInfo->incCounter();
54 }
55 #endif
56 buffer_rank = MPIInfo->mod_rank(buffer_rank - 1);
57 const index_t dof0 = distribution[buffer_rank];
58 const index_t dof1 = distribution[buffer_rank + 1];
59 #pragma omp parallel for
60 for (index_t n = 0; n < numNodes; n++) {
61 const index_t k = globalDegreesOfFreedom[n];
62 if (dof0 <= k && k < dof1) {
63 DOF_buffer[k - dof0] = SET_ID;
64 }
65 }
66 }
67 // count the entries in the buffer
68 // TODO: OMP parallel
69 dim_t myNewDOFs = 0;
70 for (index_t n = 0; n < myDOFs; ++n) {
71 if (DOF_buffer[n] == SET_ID) {
72 DOF_buffer[n] = myNewDOFs;
73 myNewDOFs++;
74 }
75 }
76
77 std::vector<index_t> loc_offsets(MPIInfo->size);
78 std::vector<index_t> offsets(MPIInfo->size);
79 dim_t new_numGlobalDOFs;
80 bool* set_new_DOF = new bool[numNodes];
81
82 #ifdef ESYS_MPI
83 new_numGlobalDOFs = 0;
84 loc_offsets[MPIInfo->rank] = myNewDOFs;
85 MPI_Allreduce(&loc_offsets[0], &offsets[0], MPIInfo->size, MPI_DIM_T,
86 MPI_SUM, MPIInfo->comm);
87 for (int n = 0; n < MPIInfo->size; ++n) {
88 loc_offsets[n] = new_numGlobalDOFs;
89 new_numGlobalDOFs += offsets[n];
90 }
91 #else
92 new_numGlobalDOFs = myNewDOFs;
93 #endif
94
95 #pragma omp parallel
96 {
97 #pragma omp for
98 for (index_t n = 0; n < myDOFs; ++n)
99 DOF_buffer[n] += loc_offsets[MPIInfo->rank];
100 #pragma omp for
101 for (index_t n = 0; n < numNodes; ++n)
102 set_new_DOF[n] = true;
103 }
104
105 // now entries are collected from the buffer again by sending them around
106 // in a circle
107 #ifdef ESYS_MPI
108 dest = MPIInfo->mod_rank(MPIInfo->rank + 1);
109 source = MPIInfo->mod_rank(MPIInfo->rank - 1);
110 #endif
111 buffer_rank = MPIInfo->rank;
112 for (int p = 0; p < MPIInfo->size; ++p) {
113 const index_t dof0 = distribution[buffer_rank];
114 const index_t dof1 = distribution[buffer_rank + 1];
115 #pragma omp parallel for
116 for (index_t n = 0; n < numNodes; n++) {
117 const index_t k = globalDegreesOfFreedom[n];
118 if (set_new_DOF[n] && dof0 <= k && k < dof1) {
119 globalDegreesOfFreedom[n] = DOF_buffer[k - dof0];
120 set_new_DOF[n] = false;
121 }
122 }
123 #ifdef ESYS_MPI
124 if (p < MPIInfo->size - 1) { // the last send can be skipped
125 MPI_Sendrecv_replace(DOF_buffer, bufferLen, MPI_DIM_T, dest,
126 MPIInfo->counter(), source, MPIInfo->counter(),
127 MPIInfo->comm, &status);
128 MPIInfo->incCounter();
129 }
130 #endif
131 buffer_rank = MPIInfo->mod_rank(buffer_rank - 1);
132 }
133 delete[] DOF_buffer;
134 delete[] set_new_DOF;
135 return new_numGlobalDOFs;
136 }
137
138 dim_t NodeFile::createDenseNodeLabeling(std::vector<index_t>& nodeDistribution,
139 const std::vector<index_t>& dofDistribution)
140 {
141 const index_t UNSET_ID = -1, SET_ID = 1;
142 const index_t myFirstDOF = dofDistribution[MPIInfo->rank];
143 const index_t myLastDOF = dofDistribution[MPIInfo->rank + 1];
144
145 // find the range of node IDs controlled by me
146 index_t min_id = escript::DataTypes::index_t_max();
147 index_t max_id = escript::DataTypes::index_t_min();
148 #pragma omp parallel
149 {
150 index_t loc_min_id = min_id;
151 index_t loc_max_id = max_id;
152 #pragma omp for
153 for (index_t n = 0; n < numNodes; n++) {
154 const index_t dof = globalDegreesOfFreedom[n];
155 if (myFirstDOF <= dof && dof < myLastDOF) {
156 loc_min_id = std::min(loc_min_id, Id[n]);
157 loc_max_id = std::max(loc_max_id, Id[n]);
158 }
159 }
160 #pragma omp critical
161 {
162 min_id = std::min(loc_min_id, min_id);
163 max_id = std::max(loc_max_id, max_id);
164 }
165 }
166 dim_t myBufferLen = (max_id >= min_id ? max_id - min_id + 1 : 0);
167 dim_t bufferLen;
168
169 #ifdef ESYS_MPI
170 MPI_Allreduce(&myBufferLen, &bufferLen, 1, MPI_DIM_T, MPI_MAX,
171 MPIInfo->comm);
172 #else
173 bufferLen = myBufferLen;
174 #endif
175
176 const dim_t headerLen = 2;
177
178 index_t* Node_buffer = new index_t[bufferLen + headerLen];
179 // mark and count the nodes in use
180 #pragma omp parallel
181 {
182 #pragma omp for
183 for (index_t n = 0; n < bufferLen + headerLen; n++)
184 Node_buffer[n] = UNSET_ID;
185 #pragma omp for
186 for (index_t n = 0; n < numNodes; n++) {
187 globalNodesIndex[n] = -1;
188 const index_t dof = globalDegreesOfFreedom[n];
189 if (myFirstDOF <= dof && dof < myLastDOF)
190 Node_buffer[Id[n] - min_id + headerLen] = SET_ID;
191 }
192 }
193 dim_t myNewNumNodes = 0;
194 for (index_t n = 0; n < myBufferLen; n++) {
195 if (Node_buffer[headerLen + n] == SET_ID) {
196 Node_buffer[headerLen + n] = myNewNumNodes;
197 myNewNumNodes++;
198 }
199 }
200 // make the local number of nodes globally available
201 #ifdef ESYS_MPI
202 MPI_Allgather(&myNewNumNodes, 1, MPI_DIM_T, &nodeDistribution[0], 1,
203 MPI_DIM_T, MPIInfo->comm);
204 #else
205 nodeDistribution[0] = myNewNumNodes;
206 #endif
207
208 dim_t globalNumNodes = 0;
209 for (int p = 0; p < MPIInfo->size; ++p) {
210 const dim_t itmp = nodeDistribution[p];
211 nodeDistribution[p] = globalNumNodes;
212 globalNumNodes += itmp;
213 }
214 nodeDistribution[MPIInfo->size] = globalNumNodes;
215
216 // offset node buffer
217 #pragma omp parallel for
218 for (index_t n = 0; n < myBufferLen; n++)
219 Node_buffer[n + headerLen] += nodeDistribution[MPIInfo->rank];
220
221 // now we send this buffer around to assign global node index
222 #ifdef ESYS_MPI
223 int dest = MPIInfo->mod_rank(MPIInfo->rank + 1);
224 int source = MPIInfo->mod_rank(MPIInfo->rank - 1);
225 #endif
226 Node_buffer[0] = min_id;
227 Node_buffer[1] = max_id;
228 int buffer_rank = MPIInfo->rank;
229 for (int p = 0; p < MPIInfo->size; ++p) {
230 const index_t nodeID0 = Node_buffer[0];
231 const index_t nodeID1 = Node_buffer[1];
232 const index_t dof0 = dofDistribution[buffer_rank];
233 const index_t dof1 = dofDistribution[buffer_rank + 1];
234 if (nodeID0 <= nodeID1) {
235 #pragma omp parallel for
236 for (index_t n = 0; n < numNodes; n++) {
237 const index_t dof = globalDegreesOfFreedom[n];
238 const index_t id = Id[n] - nodeID0;
239 if (dof0 <= dof && dof < dof1 && id >= 0 &&
240 id <= nodeID1 - nodeID0)
241 globalNodesIndex[n] = Node_buffer[id + headerLen];
242 }
243 }
244 #ifdef ESYS_MPI
245 if (p < MPIInfo->size - 1) { // the last send can be skipped
246 MPI_Status status;
247 MPI_Sendrecv_replace(Node_buffer, bufferLen + headerLen, MPI_DIM_T,
248 dest, MPIInfo->counter(), source,
249 MPIInfo->counter(), MPIInfo->comm, &status);
250 MPIInfo->incCounter();
251 }
252 #endif
253 buffer_rank = MPIInfo->mod_rank(buffer_rank - 1);
254 }
255 delete[] Node_buffer;
256 return globalNumNodes;
257 }
258
259 } // namespace dudley
260

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/NodeFile_createDenseLabelings.cpp:5567-5588 /branches/lapack2681/finley/src/NodeFile_createDenseLabelings.cpp:2682-2741 /branches/pasowrap/dudley/src/NodeFile_createDenseLabelings.cpp:3661-3674 /branches/py3_attempt2/dudley/src/NodeFile_createDenseLabelings.cpp:3871-3891 /branches/restext/finley/src/NodeFile_createDenseLabelings.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/NodeFile_createDenseLabelings.cpp:3669-3791 /branches/stage3.0/finley/src/NodeFile_createDenseLabelings.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/NodeFile_createDenseLabelings.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/NodeFile_createDenseLabelings.cpp:3517-3974 /release/3.0/finley/src/NodeFile_createDenseLabelings.cpp:2591-2601 /release/4.0/dudley/src/NodeFile_createDenseLabelings.cpp:5380-5406 /trunk/dudley/src/NodeFile_createDenseLabelings.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/NodeFile_createDenseLabelings.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26