Reference documentation for deal.II version 8.1.0
parallel_vector.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: parallel_vector.templates.h 31301 2013-10-18 12:19:18Z kronbichler @f$
3 //
4 // Copyright (C) 2011 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__parallel_vector_templates_h
18 #define __deal2__parallel_vector_templates_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/lac/parallel_vector.h>
23 
25 
26 
27 namespace parallel
28 {
29  namespace distributed
30  {
31 
32  template <typename Number>
33  void
35  {
36 #ifdef DEAL_II_WITH_MPI
37  for (size_type j=0; j<compress_requests.size(); j++)
38  MPI_Request_free(&compress_requests[j]);
39  compress_requests.clear();
40  for (size_type j=0; j<update_ghost_values_requests.size(); j++)
41  MPI_Request_free(&update_ghost_values_requests[j]);
42  update_ghost_values_requests.clear();
43 #endif
44  }
45 
46 
47 
48  template <typename Number>
49  void
50  Vector<Number>::resize_val (const size_type new_alloc_size)
51  {
52  if (new_alloc_size > allocated_size)
53  {
54  Assert (((allocated_size > 0 && val != 0) ||
55  val == 0), ExcInternalError());
56  if (val != 0)
57  delete [] val;
58  val = new Number[new_alloc_size];
59  allocated_size = new_alloc_size;
60  }
61  else if (new_alloc_size == 0)
62  {
63  if (val != 0)
64  delete [] val;
65  val = 0;
66  allocated_size = 0;
67  }
68  }
69 
70 
71 
72  template <typename Number>
73  void
74  Vector<Number>::reinit (const size_type size,
75  const bool fast)
76  {
77  clear_mpi_requests();
78  // check whether we need to reallocate
79  resize_val (size);
80 
81  // reset vector view
82  vector_view.reinit (size, val);
83 
84  // delete previous content in import data
85  if (import_data != 0)
86  delete[] import_data;
87  import_data = 0;
88 
89  // set partitioner to serial version
90  partitioner.reset (new Utilities::MPI::Partitioner (size));
91 
92  // set entries to zero if so requested
93  if (fast == false)
94  this->operator = (Number());
95 
96  vector_is_ghosted = false;
97  }
98 
99 
100 
101  template <typename Number>
102  template <typename Number2>
103  void
105  const bool fast)
106  {
107  clear_mpi_requests();
108  Assert (v.partitioner.get() != 0, ExcNotInitialized());
109 
110  // check whether the partitioners are
111  // different (check only if the are allocated
112  // differently, not if the actual data is
113  // different)
114  if (partitioner.get() != v.partitioner.get())
115  {
116  partitioner = v.partitioner;
117  const size_type new_allocated_size = partitioner->local_size() +
118  partitioner->n_ghost_indices();
119  resize_val (new_allocated_size);
120  vector_view.reinit (partitioner->local_size(), val);
121  }
122  else
123  Assert (vector_view.size() == partitioner->local_size(),
124  ExcInternalError());
125 
126  if (fast == false)
127  this->operator= (Number());
128 
129  if (import_data != 0)
130  {
131  delete [] import_data;
132 
133  // do not reallocate import_data directly, but only upon request. It
134  // is only used as temporary storage for compress() and
135  // update_ghost_values, and we might have vectors where we never
136  // call these methods and hence do not need to have the storage.
137  import_data = 0;
138  }
139 
140  vector_is_ghosted = false;
141  }
142 
143 
144 
145  template <typename Number>
146  void
147  Vector<Number>::reinit (const IndexSet &locally_owned_indices,
148  const IndexSet &ghost_indices,
149  const MPI_Comm communicator)
150  {
151  // set up parallel partitioner with index sets and communicator
152  std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> new_partitioner
153  (new Utilities::MPI::Partitioner (locally_owned_indices,
154  ghost_indices, communicator));
155  reinit (new_partitioner);
156  }
157 
158 
159 
160  template <typename Number>
161  void
162  Vector<Number>::reinit (const IndexSet &locally_owned_indices,
163  const MPI_Comm communicator)
164  {
165  // set up parallel partitioner with index sets and communicator
166  IndexSet ghost_indices(locally_owned_indices.size());
167  std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> new_partitioner
168  (new Utilities::MPI::Partitioner (locally_owned_indices,
169  ghost_indices, communicator));
170  reinit (new_partitioner);
171  }
172 
173 
174 
175  template <typename Number>
176  void
177  Vector<Number>::reinit (const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_in)
178  {
179  clear_mpi_requests();
180  partitioner = partitioner_in;
181 
182  // set vector size and allocate memory
183  const size_type new_allocated_size = partitioner->local_size() +
184  partitioner->n_ghost_indices();
185  resize_val (new_allocated_size);
186  vector_view.reinit (partitioner->local_size(), val);
187 
188  // initialize to zero
189  this->operator= (Number());
190 
191  if (import_data != 0)
192  {
193  delete [] import_data;
194 
195  // do not reallocate import_data directly, but only upon request. It
196  // is only used as temporary storage for compress() and
197  // update_ghost_values, and we might have vectors where we never
198  // call these methods and hence do not need to have the storage.
199  import_data = 0;
200  }
201 
202  vector_is_ghosted = false;
203  }
204 
205 
206 
207  template <typename Number>
208  void
210  const bool call_update_ghost_values)
211  {
212  AssertDimension (local_range().first, c.local_range().first);
213  AssertDimension (local_range().second, c.local_range().second);
214  AssertDimension (vector_view.size(), c.vector_view.size());
215  vector_view = c.vector_view;
216  if (call_update_ghost_values == true)
217  update_ghost_values();
218  else
219  vector_is_ghosted = false;
220  }
221 
222 
223 
224  template <typename Number>
225  void
226  Vector<Number>::compress_start (const unsigned int counter,
227  ::VectorOperation::values operation)
228  {
229  Assert (vector_is_ghosted == false,
230  ExcMessage ("Cannot call compress() on a ghosted vector"));
231 
232 #ifdef DEAL_II_WITH_MPI
233  // nothing to do for insert (only need to zero ghost entries in
234  // compress_finish()). in debug mode we want to check consistency
235  // of the inserted data, therefore the communication is still
236  // initialized. Having different code in debug and optimized mode is
237  // somewhat dangerous, but it really saves communication so it seems
238  // still worthwhile.
239 #ifndef DEBUG
240  if (operation == VectorOperation::insert)
241  return;
242 #endif
243 
244  const Utilities::MPI::Partitioner &part = *partitioner;
245 
246  // nothing to do when we neither have import
247  // nor ghost indices.
248  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
249  return;
250 
251  // make this function thread safe
252  Threads::Mutex::ScopedLock lock (mutex);
253 
254  const size_type n_import_targets = part.import_targets().size();
255  const size_type n_ghost_targets = part.ghost_targets().size();
256 
257  // Need to send and receive the data. Use non-blocking communication,
258  // where it is generally less overhead to first initiate the receive and
259  // then actually send the data
260  if (compress_requests.size() == 0)
261  {
262  // set channels in different range from update_ghost_values channels
263  const unsigned int channel = counter + 400;
264  unsigned int current_index_start = 0;
265  compress_requests.resize (n_import_targets + n_ghost_targets);
266 
267  // allocate import_data in case it is not set up yet
268  if (import_data == 0)
269  import_data = new Number[part.n_import_indices()];
270  for (size_type i=0; i<n_import_targets; i++)
271  {
272  MPI_Recv_init (&import_data[current_index_start],
273  part.import_targets()[i].second*sizeof(Number),
274  MPI_BYTE,
275  part.import_targets()[i].first,
276  part.import_targets()[i].first +
277  part.n_mpi_processes()*channel,
278  part.get_communicator(),
279  &compress_requests[i]);
280  current_index_start += part.import_targets()[i].second;
281  }
282  AssertDimension(current_index_start, part.n_import_indices());
283 
284  Assert (part.local_size() == vector_view.size(), ExcInternalError());
285  current_index_start = part.local_size();
286  for (size_type i=0; i<n_ghost_targets; i++)
287  {
288  MPI_Send_init (&this->val[current_index_start],
289  part.ghost_targets()[i].second*sizeof(Number),
290  MPI_BYTE,
291  part.ghost_targets()[i].first,
292  part.this_mpi_process() +
293  part.n_mpi_processes()*channel,
294  part.get_communicator(),
295  &compress_requests[n_import_targets+i]);
296  current_index_start += part.ghost_targets()[i].second;
297  }
298  AssertDimension (current_index_start,
299  part.local_size()+part.n_ghost_indices());
300  }
301 
302  AssertDimension(n_import_targets + n_ghost_targets,
303  compress_requests.size());
304  if (compress_requests.size() > 0)
305  {
306  int ierr = MPI_Startall(compress_requests.size(),&compress_requests[0]);
307  Assert (ierr == MPI_SUCCESS, ExcInternalError());
308  }
309 
310 #else // ifdef DEAL_II_WITH_MPI
311  (void)counter;
312  (void)operation;
313 #endif
314  }
315 
316 
317 
318  template <typename Number>
319  void
320  Vector<Number>::compress_finish (::VectorOperation::values operation)
321  {
322 #ifdef DEAL_II_WITH_MPI
323 
324  // in optimized mode, no communication was started, so leave the
325  // function directly (and only clear ghosts)
326 #ifndef DEBUG
327  if (operation == VectorOperation::insert)
328  {
329  zero_out_ghosts();
330  return;
331  }
332 #endif
333 
334  const Utilities::MPI::Partitioner &part = *partitioner;
335 
336  // nothing to do when we neither have import nor ghost indices.
337  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
338  return;
339 
340  // make this function thread safe
341  Threads::Mutex::ScopedLock lock (mutex);
342 
343  const size_type n_import_targets = part.import_targets().size();
344  const size_type n_ghost_targets = part.ghost_targets().size();
345 
346  if (operation != ::VectorOperation::insert)
347  AssertDimension (n_ghost_targets+n_import_targets,
348  compress_requests.size());
349 
350  // first wait for the receive to complete
351  if (compress_requests.size() > 0 && n_import_targets > 0)
352  {
353  int ierr = MPI_Waitall (n_import_targets, &compress_requests[0],
354  MPI_STATUSES_IGNORE);
355  Assert (ierr == MPI_SUCCESS, ExcInternalError());
356 
357  Number *read_position = import_data;
358  std::vector<std::pair<size_type, size_type> >::const_iterator
359  my_imports = part.import_indices().begin();
360 
361  // If the operation is no insertion, add the imported data to the
362  // local values. For insert, nothing is done here (but in debug mode
363  // we assert that the specified value is either zero or matches with
364  // the ones already present
365  if (operation != ::VectorOperation::insert)
366  for ( ; my_imports!=part.import_indices().end(); ++my_imports)
367  for (size_type j=my_imports->first; j<my_imports->second; j++)
368  local_element(j) += *read_position++;
369  else
370  for ( ; my_imports!=part.import_indices().end(); ++my_imports)
371  for (size_type j=my_imports->first; j<my_imports->second;
372  j++, read_position++)
373  Assert(*read_position == 0. ||
374  std::abs(local_element(j) - *read_position) <
375  std::abs(local_element(j)) * 1000. *
376  std::numeric_limits<Number>::epsilon(),
377  ExcMessage("Inserted elements do not match."));
378  AssertDimension(read_position-import_data,part.n_import_indices());
379  }
380 
381  if (compress_requests.size() > 0 && n_ghost_targets > 0)
382  {
383  int ierr = MPI_Waitall (n_ghost_targets,
384  &compress_requests[n_import_targets],
385  MPI_STATUSES_IGNORE);
386  Assert (ierr == MPI_SUCCESS, ExcInternalError());
387  }
388  else
389  AssertDimension (part.n_ghost_indices(), 0);
390 
391  zero_out_ghosts ();
392 #else
393  (void)operation;
394 #endif
395  }
396 
397 
398 
399  template <typename Number>
400  void
401  Vector<Number>::update_ghost_values_start (const unsigned int counter) const
402  {
403 #ifdef DEAL_II_WITH_MPI
404  const Utilities::MPI::Partitioner &part = *partitioner;
405 
406  // nothing to do when we neither have import nor ghost indices.
407  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
408  return;
409 
410  // make this function thread safe
411  Threads::Mutex::ScopedLock lock (mutex);
412 
413  const size_type n_import_targets = part.import_targets().size();
414  const size_type n_ghost_targets = part.ghost_targets().size();
415 
416  // Need to send and receive the data. Use non-blocking communication,
417  // where it is generally less overhead to first initiate the receive and
418  // then actually send the data
419  if (update_ghost_values_requests.size() == 0)
420  {
421  Assert (part.local_size() == vector_view.size(),
422  ExcInternalError());
423  size_type current_index_start = part.local_size();
424  update_ghost_values_requests.resize (n_import_targets+n_ghost_targets);
425  for (size_type i=0; i<n_ghost_targets; i++)
426  {
427  // allow writing into ghost indices even though we are in a
428  // const function
429  MPI_Recv_init (const_cast<Number *>(&val[current_index_start]),
430  part.ghost_targets()[i].second*sizeof(Number),
431  MPI_BYTE,
432  part.ghost_targets()[i].first,
433  part.ghost_targets()[i].first +
434  counter*part.n_mpi_processes(),
435  part.get_communicator(),
436  &update_ghost_values_requests[i]);
437  current_index_start += part.ghost_targets()[i].second;
438  }
439  AssertDimension (current_index_start,
440  part.local_size()+part.n_ghost_indices());
441 
442  // allocate import_data in case it is not set up yet
443  if (import_data == 0 && part.n_import_indices() > 0)
444  import_data = new Number[part.n_import_indices()];
445  current_index_start = 0;
446  for (size_type i=0; i<n_import_targets; i++)
447  {
448  MPI_Send_init (&import_data[current_index_start],
449  part.import_targets()[i].second*sizeof(Number),
450  MPI_BYTE, part.import_targets()[i].first,
451  part.this_mpi_process() +
452  part.n_mpi_processes()*counter,
453  part.get_communicator(),
454  &update_ghost_values_requests[n_ghost_targets+i]);
455  current_index_start += part.import_targets()[i].second;
456  }
457  AssertDimension (current_index_start, part.n_import_indices());
458  }
459 
460  // copy the data that is actually to be send to the import_data field
461  if (part.n_import_indices() > 0)
462  {
463  Assert (import_data != 0, ExcInternalError());
464  Number *write_position = import_data;
465  std::vector<std::pair<size_type, size_type> >::const_iterator
466  my_imports = part.import_indices().begin();
467  for ( ; my_imports!=part.import_indices().end(); ++my_imports)
468  for (size_type j=my_imports->first; j<my_imports->second; j++)
469  *write_position++ = local_element(j);
470  }
471 
472  AssertDimension (n_import_targets+n_ghost_targets,
473  update_ghost_values_requests.size());
474  if (update_ghost_values_requests.size() > 0)
475  {
476  int ierr = MPI_Startall(update_ghost_values_requests.size(),
477  &update_ghost_values_requests[0]);
478  Assert (ierr == MPI_SUCCESS, ExcInternalError());
479  }
480 #else
481  (void)counter;
482 #endif
483  }
484 
485 
486 
487  template <typename Number>
488  void
490  {
491 #ifdef DEAL_II_WITH_MPI
492  // wait for both sends and receives to complete, even though only
493  // receives are really necessary. this gives (much) better performance
494  AssertDimension (partitioner->ghost_targets().size() +
495  partitioner->import_targets().size(),
496  update_ghost_values_requests.size());
497  if (update_ghost_values_requests.size() > 0)
498  {
499  // make this function thread safe
500  Threads::Mutex::ScopedLock lock (mutex);
501 
502  int ierr = MPI_Waitall (update_ghost_values_requests.size(),
503  &update_ghost_values_requests[0],
504  MPI_STATUSES_IGNORE);
505  Assert (ierr == MPI_SUCCESS, ExcInternalError());
506  }
507 #endif
508  vector_is_ghosted = true;
509  }
510 
511 
512 
513  template <typename Number>
514  void
516  {
517 #ifdef DEAL_II_WITH_MPI
518 
519 #ifdef DEBUG
520  // make sure that there are not outstanding requests from updating ghost
521  // values or compress
522  int flag = 1;
523  int ierr = MPI_Testall (update_ghost_values_requests.size(),
524  &update_ghost_values_requests[0],
525  &flag, MPI_STATUSES_IGNORE);
526  Assert (ierr == MPI_SUCCESS, ExcInternalError());
527  Assert (flag == 1,
528  ExcMessage("MPI found unfinished update_ghost_values() requests"
529  "when calling swap, which is not allowed"));
530  ierr = MPI_Testall (compress_requests.size(), &compress_requests[0],
531  &flag, MPI_STATUSES_IGNORE);
532  Assert (ierr == MPI_SUCCESS, ExcInternalError());
533  Assert (flag == 1,
534  ExcMessage("MPI found unfinished compress() requests "
535  "when calling swap, which is not allowed"));
536 #endif
537 
538  std::swap (compress_requests, v.compress_requests);
539  std::swap (update_ghost_values_requests, v.update_ghost_values_requests);
540 #endif
541 
542  std::swap (partitioner, v.partitioner);
543  std::swap (allocated_size, v.allocated_size);
544  std::swap (val, v.val);
545  std::swap (import_data, v.import_data);
546  std::swap (vector_is_ghosted, v.vector_is_ghosted);
547 
548  // vector view cannot be swapped so reset it manually (without touching
549  // the vector elements)
550  vector_view.reinit (partitioner->local_size(), val);
551  v.vector_view.reinit (v.partitioner->local_size(), v.val);
552  }
553 
554 
555 
556  template <typename Number>
557  std::size_t
559  {
560  std::size_t memory = sizeof(*this);
561  memory += sizeof (Number) * static_cast<std::size_t>(allocated_size);
562 
563  // if the partitioner is shared between more processors, just count a
564  // fraction of that memory, since we're not actually using more memory
565  // for it.
566  if (partitioner.use_count() > 0)
567  memory += partitioner->memory_consumption()/partitioner.use_count()+1;
568  if (import_data != 0)
569  memory += (static_cast<std::size_t>(partitioner->n_import_indices())*
570  sizeof(Number));
571  return memory;
572  }
573 
574 
575 
576  template <typename Number>
577  void
578  Vector<Number>::print (std::ostream &out,
579  const unsigned int precision,
580  const bool scientific,
581  const bool across) const
582  {
583  Assert (partitioner.get() !=0, ExcInternalError());
584  AssertThrow (out, ExcIO());
585  std::ios::fmtflags old_flags = out.flags();
586  unsigned int old_precision = out.precision (precision);
587 
588  out.precision (precision);
589  if (scientific)
590  out.setf (std::ios::scientific, std::ios::floatfield);
591  else
592  out.setf (std::ios::fixed, std::ios::floatfield);
593 
594  // to make the vector write out all the information in order, use as
595  // many barriers as there are processors and start writing when it's our
596  // turn
597 #ifdef DEAL_II_WITH_MPI
598  if (partitioner->n_mpi_processes() > 1)
599  for (unsigned int i=0; i<partitioner->this_mpi_process(); i++)
600  MPI_Barrier (partitioner->get_communicator());
601 #endif
602 
603  out << "Process #" << partitioner->this_mpi_process() << std::endl
604  << "Local range: [" << partitioner->local_range().first << ", "
605  << partitioner->local_range().second << "), global size: "
606  << partitioner->size() << std::endl
607  << "Vector data:" << std::endl;
608  if (across)
609  for (size_type i=0; i<partitioner->local_size(); ++i)
610  out << local_element(i) << ' ';
611  else
612  for (size_type i=0; i<partitioner->local_size(); ++i)
613  out << local_element(i) << std::endl;
614  out << std::endl;
615 
616  if (vector_is_ghosted)
617  {
618  out << "Ghost entries (global index / value):" << std::endl;
619  if (across)
620  for (size_type i=0; i<partitioner->n_ghost_indices(); ++i)
621  out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
622  << '/' << local_element(partitioner->local_size()+i) << ") ";
623  else
624  for (size_type i=0; i<partitioner->n_ghost_indices(); ++i)
625  out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
626  << '/' << local_element(partitioner->local_size()+i) << ")"
627  << std::endl;
628  out << std::endl;
629  }
630  out << std::flush;
631 
632 #ifdef DEAL_II_WITH_MPI
633  if (partitioner->n_mpi_processes() > 1)
634  {
635  MPI_Barrier (partitioner->get_communicator());
636 
637  for (unsigned int i=partitioner->this_mpi_process()+1;
638  i<partitioner->n_mpi_processes(); i++)
639  MPI_Barrier (partitioner->get_communicator());
640  }
641 #endif
642 
643  AssertThrow (out, ExcIO());
644  // reset output format
645  out.flags (old_flags);
646  out.precision(old_precision);
647  }
648 
649  } // end of namespace distributed
650 
651 } // end of namespace parallel
652 
653 
654 DEAL_II_NAMESPACE_CLOSE
655 
656 #endif
unsigned int n_ghost_indices() const
const std::vector< std::pair< types::global_dof_index, types::global_dof_index > > & import_indices() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void compress_start(const unsigned int communication_channel=0,::VectorOperation::values operation=VectorOperation::add)
unsigned int this_mpi_process() const
types::global_dof_index size() const
Definition: index_set.h:685
void compress_finish(::VectorOperation::values operation)
std_cxx1x::shared_ptr< const Utilities::MPI::Partitioner > partitioner
::ExceptionBase & ExcMessage(std::string arg1)
void reinit(const size_type size, const bool fast=false)
void resize_val(const size_type new_allocated_size)
const MPI_Comm & get_communicator() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
const std::vector< std::pair< unsigned int, types::global_dof_index > > & ghost_targets() const
::ExceptionBase & ExcIO()
const std::vector< std::pair< unsigned int, types::global_dof_index > > & import_targets() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
#define Assert(cond, exc)
Definition: exceptions.h:299
void update_ghost_values_start(const unsigned int communication_channel=0) const
unsigned int n_import_indices() const
unsigned int n_mpi_processes() const
::ExceptionBase & ExcNotInitialized()
void copy_from(const Vector< Number > &in_vector, const bool call_update_ghost_values=false)
::ExceptionBase & ExcInternalError()
unsigned int local_size() const