double_multi_vector.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#include "double_multi_vector.h"
27#include "matrices.h"
28
29namespace oomph
30{
31 /// The contents of the vector are redistributed to match the new
32 /// distribution. In a non-MPI rebuild this method works, but does nothing.
33 /// \b NOTE 1: The current distribution and the new distribution must have
34 /// the same number of global rows.
35 /// \b NOTE 2: The current distribution and the new distribution must have
36 /// the same Communicator.
38 const LinearAlgebraDistribution* const& dist_pt)
39 {
40#ifdef OOMPH_HAS_MPI
41#ifdef PARANOID
42 if (!Internal_values)
43 {
44 // if this vector does not own the double* values then it cannot be
45 // distributed.
46 // note: this is not stictly necessary - would just need to be careful
47 // with delete[] below.
48 std::ostringstream error_message;
49 error_message
50 << "This multi vector does not own its data (i.e. data has been "
51 << "passed in via set_external_values() and therefore "
52 << "cannot be redistributed";
53 throw OomphLibError(
54 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
55 }
56 // paranoid check that the nrows for both distributions is the
57 // same
58 if (dist_pt->nrow() != this->nrow())
59 {
60 std::ostringstream error_message;
61 error_message << "The number of global rows in the new distribution ("
62 << dist_pt->nrow() << ") is not equal to the number"
63 << " of global rows in the current distribution ("
64 << this->nrow() << ").\n";
65 throw OomphLibError(
66 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
67 }
68 // paranoid check that the current distribution and the new distribution
69 // have the same Communicator
70 OomphCommunicator temp_comm(*dist_pt->communicator_pt());
71 if (!(temp_comm == *this->distribution_pt()->communicator_pt()))
72 {
73 std::ostringstream error_message;
74 error_message << "The new distribution and the current distribution must "
75 << "have the same communicator.";
76 throw OomphLibError(
77 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
78 }
79#endif
80
81 // check the distributions are not the same
82 if (!((*this->distribution_pt()) == *dist_pt))
83 {
84 // Cache the number of vectors
85 const unsigned n_vector = this->Nvector;
86
87 // get the rank and the number of processors
88 int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
89 int nproc = this->distribution_pt()->communicator_pt()->nproc();
90
91 // if both vectors are distributed
92 if (this->distributed() && dist_pt->distributed())
93 {
94 // new nrow_local and first_row data
95 Vector<unsigned> new_first_row_data(nproc);
96 Vector<unsigned> new_nrow_local_data(nproc);
97 Vector<unsigned> current_first_row_data(nproc);
98 Vector<unsigned> current_nrow_local_data(nproc);
99 for (int i = 0; i < nproc; i++)
100 {
101 new_first_row_data[i] = dist_pt->first_row(i);
102 new_nrow_local_data[i] = dist_pt->nrow_local(i);
103 current_first_row_data[i] = this->first_row(i);
104 current_nrow_local_data[i] = this->nrow_local(i);
105 }
106
107 // compute which local rows are expected to be received from each
108 // processor / sent to each processor
109 Vector<unsigned> new_first_row_for_proc(nproc);
110 Vector<unsigned> new_nrow_local_for_proc(nproc);
111 Vector<unsigned> new_first_row_from_proc(nproc);
112 Vector<unsigned> new_nrow_local_from_proc(nproc);
113
114 // for every processor compute first_row and nrow_local that will
115 // will sent and received by this processor
116 for (int p = 0; p < nproc; p++)
117 {
118 // start with data to be sent
119 if ((new_first_row_data[p] < (current_first_row_data[my_rank] +
120 current_nrow_local_data[my_rank])) &&
121 (current_first_row_data[my_rank] <
122 (new_first_row_data[p] + new_nrow_local_data[p])))
123 {
124 new_first_row_for_proc[p] =
125 std::max(current_first_row_data[my_rank], new_first_row_data[p]);
126 new_nrow_local_for_proc[p] =
127 std::min((current_first_row_data[my_rank] +
128 current_nrow_local_data[my_rank]),
129 (new_first_row_data[p] + new_nrow_local_data[p])) -
130 new_first_row_for_proc[p];
131 }
132
133 // and data to be received
134 if ((new_first_row_data[my_rank] <
135 (current_first_row_data[p] + current_nrow_local_data[p])) &&
136 (current_first_row_data[p] <
137 (new_first_row_data[my_rank] + new_nrow_local_data[my_rank])))
138 {
139 new_first_row_from_proc[p] =
140 std::max(current_first_row_data[p], new_first_row_data[my_rank]);
141 new_nrow_local_from_proc[p] =
142 std::min(
143 (current_first_row_data[p] + current_nrow_local_data[p]),
144 (new_first_row_data[my_rank] + new_nrow_local_data[my_rank])) -
145 new_first_row_from_proc[p];
146 }
147 }
148
149 // Storage for the new data
150 double** temp_data = new double*[n_vector];
151 double* contiguous_temp_data =
152 new double[n_vector * new_nrow_local_data[my_rank]];
153 for (unsigned v = 0; v < n_vector; ++v)
154 {
155 temp_data[v] =
156 &contiguous_temp_data[v * new_nrow_local_data[my_rank]];
157 }
158
159 // "send to self" or copy Data that does not need to be sent else where
160 // to temp_data
161 if (new_nrow_local_for_proc[my_rank] != 0)
162 {
163 unsigned j =
164 new_first_row_for_proc[my_rank] - current_first_row_data[my_rank];
165 unsigned k =
166 new_first_row_for_proc[my_rank] - new_first_row_data[my_rank];
167 for (unsigned i = 0; i < new_nrow_local_for_proc[my_rank]; i++)
168 {
169 for (unsigned v = 0; v < n_vector; ++v)
170 {
171 temp_data[v][k + i] = Values[v][j + i];
172 }
173 }
174 }
175
176 // send and receive circularly
177 for (int p = 1; p < nproc; p++)
178 {
179 // next processor to send to
180 unsigned dest_p = (my_rank + p) % nproc;
181
182 // next processor to receive from
183 unsigned source_p = (nproc + my_rank - p) % nproc;
184
185 // send and receive the value
186 MPI_Status status;
187 for (unsigned v = 0; v < n_vector; v++)
188 {
189 MPI_Sendrecv(Values[v] + new_first_row_for_proc[dest_p] -
190 current_first_row_data[my_rank],
191 new_nrow_local_for_proc[dest_p],
192 MPI_DOUBLE,
193 dest_p,
194 1,
195 temp_data[v] + new_first_row_from_proc[source_p] -
196 new_first_row_data[my_rank],
197 new_nrow_local_from_proc[source_p],
198 MPI_DOUBLE,
199 source_p,
200 1,
201 this->distribution_pt()->communicator_pt()->mpi_comm(),
202 &status);
203 }
204 }
205
206 // copy from temp data to Values_pt
207 delete[] Values[0];
208 delete[] Values;
209 Values = temp_data;
210 }
211 // if this vector is distributed but the new distributed is global
212 else if (this->distributed() && !dist_pt->distributed())
213 {
214 // copy existing Values_pt to temp_data
215 unsigned n_local_data = this->nrow_local();
216 double** temp_data = new double*[n_vector];
217 // New continguous data
218 double* contiguous_temp_data = new double[n_vector * n_local_data];
219 for (unsigned v = 0; v < n_vector; ++v)
220 {
221 temp_data[v] = &contiguous_temp_data[v * n_local_data];
222 for (unsigned i = 0; i < n_local_data; i++)
223 {
224 temp_data[v][i] = Values[v][i];
225 }
226 }
227
228 // clear and resize Values_pt
229 delete[] Values[0];
230 double* values = new double[this->nrow() * n_vector];
231 for (unsigned v = 0; v < n_vector; v++)
232 {
233 Values[v] = &values[v * this->nrow()];
234 }
235
236 // create a int vector of first rows
237 int* dist_first_row = new int[nproc];
238 int* dist_nrow_local = new int[nproc];
239 for (int p = 0; p < nproc; p++)
240 {
241 dist_first_row[p] = this->first_row(p);
242 dist_nrow_local[p] = this->nrow_local(p);
243 }
244
245 // gather the local vectors from all processors on all processors
246 int my_local_data(this->nrow_local());
247
248 // Loop over all vectors
249 for (unsigned v = 0; v < n_vector; v++)
250 {
251 MPI_Allgatherv(
252 temp_data[v],
253 my_local_data,
254 MPI_DOUBLE,
255 Values[v],
256 dist_nrow_local,
257 dist_first_row,
258 MPI_DOUBLE,
259 this->distribution_pt()->communicator_pt()->mpi_comm());
260 }
261
262 // update the distribution
263 this->build_distribution(dist_pt);
264
265 // delete the temp_data
266 delete[] temp_data[0];
267 delete[] temp_data;
268
269 // clean up
270 delete[] dist_first_row;
271 delete[] dist_nrow_local;
272 }
273
274 // if this vector is not distrubted but the target vector is
275 else if (!this->distributed() && dist_pt->distributed())
276 {
277 // cache the new nrow_local
278 unsigned nrow_local = dist_pt->nrow_local();
279
280 // and first_row
281 unsigned first_row = dist_pt->first_row();
282
283 const unsigned n_local_data = nrow_local;
284 double** temp_data = new double*[n_vector];
285 double* contiguous_temp_data = new double[n_vector * n_local_data];
286
287 // copy the data
288 for (unsigned v = 0; v < n_vector; v++)
289 {
290 temp_data[v] = &contiguous_temp_data[v * n_local_data];
291 for (unsigned i = 0; i < n_local_data; i++)
292 {
293 temp_data[v][i] = Values[v][first_row + i];
294 }
295 }
296
297 // copy to Values_pt
298 delete[] Values[0];
299 delete[] Values;
300 Values = temp_data;
301
302 // update the distribution
303 this->build_distribution(dist_pt);
304 }
305
306 // copy the Distribution
307 this->build_distribution(dist_pt);
308 }
309#endif
310
311 // Update the doublevector representation
313 }
314
315
316} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
bool distributed() const
distribution is serial or distributed
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
void setup_doublevector_representation()
compute the A-norm using the matrix at matrix_pt
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
Allows are external data to be used by this vector. WARNING: The size of the external data must corre...
unsigned Nvector
The number of vectors.
double ** values()
access function to the underlying values
double ** Values
the local data, need a pointer to a pointer so that the individual vectors can be extracted
bool Internal_values
Boolean flag to indicate whether the vector's data (values_pt) is owned by this vector.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...