Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_vector.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
18 #ifdef DEAL_II_WITH_TRILINOS
19 
20 # include <deal.II/base/mpi.h>
22 
27 
28 # include <boost/io/ios_state.hpp>
29 
30 # include <Epetra_Export.h>
31 # include <Epetra_Import.h>
32 # include <Epetra_Vector.h>
33 
34 # include <cmath>
35 # include <memory>
36 
37 
39 
40 namespace TrilinosWrappers
41 {
42 # ifndef DOXYGEN
43  namespace internal
44  {
45  VectorReference::operator TrilinosScalar() const
46  {
47  AssertIndexRange(index, vector.size());
48 
49  // Trilinos allows for vectors to be referenced by the [] or ()
50  // operators but only () checks index bounds. We check these bounds by
51  // ourselves, so we can use []. Note that we can only get local values.
52 
53  const TrilinosWrappers::types::int_type local_index =
54  vector.vector->Map().LID(
55  static_cast<TrilinosWrappers::types::int_type>(index));
56  Assert(local_index >= 0,
58  index,
59  vector.vector->Map().NumMyElements(),
60  vector.vector->Map().MinMyGID(),
61  vector.vector->Map().MaxMyGID()));
62 
63 
64  return (*(vector.vector))[0][local_index];
65  }
66  } // namespace internal
67 # endif
68 
69  namespace MPI
70  {
72  : Subscriptor()
73  , last_action(Zero)
74  , compressed(true)
75  , has_ghosts(false)
76  , vector(new Epetra_FEVector(
77  Epetra_Map(0, 0, 0, Utilities::Trilinos::comm_self())))
78  {}
79 
80 
81 
82  Vector::Vector(const IndexSet &parallel_partitioning,
83  const MPI_Comm communicator)
84  : Vector()
85  {
86  reinit(parallel_partitioning, communicator);
87  }
88 
89 
90 
92  : Vector()
93  {
95  vector = std::make_unique<Epetra_FEVector>(*v.vector);
97  }
98 
99 
100 
101  Vector::Vector(Vector &&v) // NOLINT
102  : Vector()
103  {
104  // initialize a minimal, valid object and swap
105  static_cast<Subscriptor &>(*this) = static_cast<Subscriptor &&>(v);
106  swap(v);
107  }
108 
109 
110 
111  Vector::Vector(const IndexSet &parallel_partitioner,
112  const Vector &v,
113  const MPI_Comm communicator)
114  : Vector()
115  {
116  AssertThrow(parallel_partitioner.size() ==
117  static_cast<size_type>(
119  ExcDimensionMismatch(parallel_partitioner.size(),
121  v.vector->Map())));
122 
123  vector = std::make_unique<Epetra_FEVector>(
124  parallel_partitioner.make_trilinos_map(communicator, true));
125  reinit(v, false, true);
126  }
127 
128 
129 
130  Vector::Vector(const IndexSet &local,
131  const IndexSet &ghost,
132  const MPI_Comm communicator)
133  : Vector()
134  {
135  reinit(local, ghost, communicator, false);
136  }
137 
138 
139 
140  void
142  {
143  // When we clear the vector, reset the pointer and generate an empty
144  // vector.
145  Epetra_Map map(0, 0, Epetra_MpiComm(MPI_COMM_SELF));
146 
147  has_ghosts = false;
148  vector = std::make_unique<Epetra_FEVector>(map);
149  last_action = Zero;
150  }
151 
152 
153 
154  void
155  Vector::reinit(const IndexSet &parallel_partitioner,
156  const MPI_Comm communicator,
157  const bool /*omit_zeroing_entries*/)
158  {
159  nonlocal_vector.reset();
160 
161  const bool overlapping =
162  !parallel_partitioner.is_ascending_and_one_to_one(communicator);
163 
164  Epetra_Map map =
165  parallel_partitioner.make_trilinos_map(communicator, overlapping);
166 
167  vector = std::make_unique<Epetra_FEVector>(map);
168 
169  has_ghosts = vector->Map().UniqueGIDs() == false;
170 
171  // If the IndexSets are overlapping, we don't really know
172  // which process owns what. So we decide that no process
173  // owns anything in that case. In particular asking for
174  // the locally owned elements is not allowed.
175  if (has_ghosts)
176  {
179  }
180  else
181  owned_elements = parallel_partitioner;
182 
183 # ifdef DEBUG
184  const size_type n_elements_global =
186 
187  Assert(has_ghosts || n_elements_global == size(), ExcInternalError());
188 # endif
189 
190  last_action = Zero;
191  }
192 
193 
194 
195  void
197  const bool omit_zeroing_entries,
198  const bool allow_different_maps)
199  {
200  nonlocal_vector.reset();
201 
202  // In case we do not allow to have different maps, this call means that
203  // we have to reset the vector. So clear the vector, initialize our map
204  // with the map in v, and generate the vector.
205  if (allow_different_maps == false)
206  {
207  // check equality for MPI communicators: We can only choose the fast
208  // version in case the underlying Epetra_MpiComm object is the same,
209  // otherwise we might access an MPI_Comm object that has been
210  // deleted
211  const Epetra_MpiComm *my_comm =
212  dynamic_cast<const Epetra_MpiComm *>(&vector->Comm());
213  const Epetra_MpiComm *v_comm =
214  dynamic_cast<const Epetra_MpiComm *>(&v.vector->Comm());
215  const bool same_communicators =
216  my_comm != nullptr && v_comm != nullptr &&
217  my_comm->DataPtr() == v_comm->DataPtr();
218  if (!same_communicators ||
219  vector->Map().SameAs(v.vector->Map()) == false)
220  {
221  vector = std::make_unique<Epetra_FEVector>(v.vector->Map());
223  last_action = Zero;
225  }
226  else if (omit_zeroing_entries == false)
227  {
228  // old and new vectors have exactly the same map, i.e. size and
229  // parallel distribution
230  int ierr;
231  ierr = vector->GlobalAssemble(last_action);
232  (void)ierr;
233  Assert(ierr == 0, ExcTrilinosError(ierr));
234 
235  ierr = vector->PutScalar(0.0);
236  Assert(ierr == 0, ExcTrilinosError(ierr));
237 
238  last_action = Zero;
239  }
240  }
241 
242  // Otherwise, we have to check that the two vectors are already of the
243  // same size, create an object for the data exchange and then insert all
244  // the data. The first assertion is only a check whether the user knows
245  // what they are doing.
246  else
247  {
248  Assert(omit_zeroing_entries == false,
249  ExcMessage(
250  "It is not possible to exchange data with the "
251  "option 'omit_zeroing_entries' set, which would not write "
252  "elements."));
253 
254  AssertThrow(size() == v.size(),
255  ExcDimensionMismatch(size(), v.size()));
256 
257  Epetra_Import data_exchange(vector->Map(), v.vector->Map());
258 
259  const int ierr = vector->Import(*v.vector, data_exchange, Insert);
260  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
261 
262  last_action = Insert;
263  }
264 # ifdef DEBUG
265  const Epetra_MpiComm *comm_ptr =
266  dynamic_cast<const Epetra_MpiComm *>(&(v.vector->Comm()));
267  Assert(comm_ptr != nullptr, ExcInternalError());
268  const size_type n_elements_global =
269  Utilities::MPI::sum(owned_elements.n_elements(), comm_ptr->Comm());
270  Assert(has_ghosts || n_elements_global == size(), ExcInternalError());
271 # endif
272  }
273 
274 
275 
276  void
277  Vector::reinit(const MPI::BlockVector &v, const bool import_data)
278  {
279  nonlocal_vector.reset();
282 
283  // In case we do not allow to have different maps, this call means that
284  // we have to reset the vector. So clear the vector, initialize our map
285  // with the map in v, and generate the vector.
286  if (v.n_blocks() == 0)
287  return;
288 
289  // create a vector that holds all the elements contained in the block
290  // vector. need to manually create an Epetra_Map.
291  size_type n_elements = 0, added_elements = 0, block_offset = 0;
292  for (size_type block = 0; block < v.n_blocks(); ++block)
293  n_elements += v.block(block).vector->Map().NumMyElements();
294  std::vector<TrilinosWrappers::types::int_type> global_ids(n_elements, -1);
295  for (size_type block = 0; block < v.n_blocks(); ++block)
296  {
297  TrilinosWrappers::types::int_type *glob_elements =
299  v.block(block).trilinos_partitioner());
300  size_type vector_size = v.block(block).vector->Map().NumMyElements();
301  for (size_type i = 0; i < vector_size; ++i)
302  global_ids[added_elements++] = glob_elements[i] + block_offset;
303  owned_elements.add_indices(v.block(block).owned_elements,
304  block_offset);
305  block_offset += v.block(block).size();
306  }
307 
308  Assert(n_elements == added_elements, ExcInternalError());
309  Epetra_Map new_map(v.size(),
310  n_elements,
311  global_ids.data(),
312  0,
313  v.block(0).trilinos_partitioner().Comm());
314 
315  auto actual_vec = std::make_unique<Epetra_FEVector>(new_map);
316 
317  TrilinosScalar *entries = (*actual_vec)[0];
318  for (size_type block = 0; block < v.n_blocks(); ++block)
319  {
320  v.block(block).trilinos_vector().ExtractCopy(entries, 0);
321  entries += v.block(block).vector->Map().NumMyElements();
322  }
323 
324  if (import_data == true)
325  {
327  *actual_vec)) == v.size(),
329  *actual_vec),
330  v.size()));
331 
332  Epetra_Import data_exchange(vector->Map(), actual_vec->Map());
333 
334  const int ierr = vector->Import(*actual_vec, data_exchange, Insert);
335  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
336 
337  last_action = Insert;
338  }
339  else
340  vector = std::move(actual_vec);
341 # ifdef DEBUG
342  const Epetra_MpiComm *comm_ptr =
343  dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
344  Assert(comm_ptr != nullptr, ExcInternalError());
345  const size_type n_elements_global =
346  Utilities::MPI::sum(owned_elements.n_elements(), comm_ptr->Comm());
347 
348  Assert(has_ghosts || n_elements_global == size(), ExcInternalError());
349 # endif
350  }
351 
352 
353 
354  void
355  Vector::reinit(const IndexSet &locally_owned_entries,
356  const IndexSet &ghost_entries,
357  const MPI_Comm communicator,
358  const bool vector_writable)
359  {
360  nonlocal_vector.reset();
361  owned_elements = locally_owned_entries;
362  if (vector_writable == false)
363  {
364  IndexSet parallel_partitioner = locally_owned_entries;
365  parallel_partitioner.add_indices(ghost_entries);
366  Epetra_Map map =
367  parallel_partitioner.make_trilinos_map(communicator, true);
368  vector = std::make_unique<Epetra_FEVector>(map);
369  }
370  else
371  {
372  Epetra_Map map =
373  locally_owned_entries.make_trilinos_map(communicator, true);
374  Assert(map.IsOneToOne(),
375  ExcMessage("A writable vector must not have ghost entries in "
376  "its parallel partitioning"));
377 
378  if (vector->Map().SameAs(map) == false)
379  vector = std::make_unique<Epetra_FEVector>(map);
380  else
381  {
382  const int ierr = vector->PutScalar(0.);
383  (void)ierr;
384  Assert(ierr == 0, ExcTrilinosError(ierr));
385  }
386 
387  IndexSet nonlocal_entries(ghost_entries);
388  nonlocal_entries.subtract_set(locally_owned_entries);
389  if (Utilities::MPI::n_mpi_processes(communicator) > 1)
390  {
391  Epetra_Map nonlocal_map =
392  nonlocal_entries.make_trilinos_map(communicator, true);
394  std::make_unique<Epetra_MultiVector>(nonlocal_map, 1);
395  }
396  }
397 
398  has_ghosts = vector->Map().UniqueGIDs() == false;
399 
400  last_action = Zero;
401 
402 # ifdef DEBUG
403  const size_type n_elements_global =
405 
406  Assert(has_ghosts || n_elements_global == size(), ExcInternalError());
407 # endif
408  }
409 
410 
411 
412  void
414  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
415  const bool make_ghosted,
416  const bool vector_writable)
417  {
418  if (make_ghosted)
419  {
420  Assert(partitioner->ghost_indices_initialized(),
421  ExcMessage("You asked to create a ghosted vector, but the "
422  "partitioner does not provide ghost indices."));
423 
424  this->reinit(partitioner->locally_owned_range(),
425  partitioner->ghost_indices(),
426  partitioner->get_mpi_communicator(),
427  vector_writable);
428  }
429  else
430  {
431  this->reinit(partitioner->locally_owned_range(),
432  partitioner->get_mpi_communicator());
433  }
434  }
435 
436 
437 
438  Vector &
440  {
441  Assert(vector.get() != nullptr,
442  ExcMessage("Vector is not constructed properly."));
443 
444  // check equality for MPI communicators to avoid accessing a possibly
445  // invalid MPI_Comm object
446  const Epetra_MpiComm *my_comm =
447  dynamic_cast<const Epetra_MpiComm *>(&vector->Comm());
448  const Epetra_MpiComm *v_comm =
449  dynamic_cast<const Epetra_MpiComm *>(&v.vector->Comm());
450  const bool same_communicators = my_comm != nullptr && v_comm != nullptr &&
451  my_comm->DataPtr() == v_comm->DataPtr();
452  // Need to ask MPI whether the communicators are the same. We would like
453  // to use the following checks but currently we cannot make sure the
454  // memory of my_comm is not stale from some MPI_Comm_free
455  // somewhere. This can happen when a vector lives in GrowingVectorMemory
456  // data structures. Thus, the following code is commented out.
457  //
458  // if (my_comm != nullptr &&
459  // v_comm != nullptr &&
460  // my_comm->DataPtr() != v_comm->DataPtr())
461  // {
462  // int communicators_same = 0;
463  // const int ierr = MPI_Comm_compare (my_comm->GetMpiComm(),
464  // v_comm->GetMpiComm(),
465  // &communicators_same);
466  // AssertThrowMPI(ierr);
467  // if (!(communicators_same == MPI_IDENT ||
468  // communicators_same == MPI_CONGRUENT))
469  // same_communicators = false;
470  // else
471  // same_communicators = true;
472  // }
473 
474  // distinguish three cases. First case: both vectors have the same
475  // layout (just need to copy the local data, not reset the memory and
476  // the underlying Epetra_Map). The third case means that we have to
477  // rebuild the calling vector.
478  if (same_communicators && v.vector->Map().SameAs(vector->Map()))
479  {
480  *vector = *v.vector;
481  if (v.nonlocal_vector.get() != nullptr)
483  std::make_unique<Epetra_MultiVector>(v.nonlocal_vector->Map(), 1);
484  last_action = Zero;
485  }
486  // Second case: vectors have the same global
487  // size, but different parallel layouts (and
488  // one of them a one-to-one mapping). Then we
489  // can call the import/export functionality.
490  else if (size() == v.size() &&
491  (v.vector->Map().UniqueGIDs() || vector->Map().UniqueGIDs()))
492  {
493  reinit(v, false, true);
494  }
495  // Third case: Vectors do not have the same
496  // size.
497  else
498  {
499  vector = std::make_unique<Epetra_FEVector>(*v.vector);
500  last_action = Zero;
503  }
504 
505  if (v.nonlocal_vector.get() != nullptr)
507  std::make_unique<Epetra_MultiVector>(v.nonlocal_vector->Map(), 1);
508 
509  return *this;
510  }
511 
512 
513 
514  Vector &
515  Vector::operator=(Vector &&v) noexcept
516  {
517  static_cast<Subscriptor &>(*this) = static_cast<Subscriptor &&>(v);
518  swap(v);
519  return *this;
520  }
521 
522 
523 
524  template <typename number>
525  Vector &
526  Vector::operator=(const ::Vector<number> &v)
527  {
528  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
529 
530  // this is probably not very efficient but works. in particular, we could
531  // do better if we know that number==TrilinosScalar because then we could
532  // elide the copying of elements
533  //
534  // let's hope this isn't a particularly frequent operation
535  std::pair<size_type, size_type> local_range = this->local_range();
536  for (size_type i = local_range.first; i < local_range.second; ++i)
537  (*vector)[0][i - local_range.first] = v(i);
538 
539  return *this;
540  }
541 
542 
543 
544  void
546  const Vector &v)
547  {
548  Assert(m.trilinos_matrix().Filled() == true,
549  ExcMessage("Matrix is not compressed. "
550  "Cannot find exchange information!"));
551  Assert(v.vector->Map().UniqueGIDs() == true,
552  ExcMessage("The input vector has overlapping data, "
553  "which is not allowed."));
554 
555  if (vector->Map().SameAs(m.trilinos_matrix().ColMap()) == false)
556  vector =
557  std::make_unique<Epetra_FEVector>(m.trilinos_matrix().ColMap());
558 
559  Epetra_Import data_exchange(vector->Map(), v.vector->Map());
560  const int ierr = vector->Import(*v.vector, data_exchange, Insert);
561 
562  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
563 
564  last_action = Insert;
565  }
566 
567 
568  void
570  const VectorOperation::values operation)
571  {
572  Assert(
573  this->size() == rwv.size(),
574  ExcMessage(
575  "Both vectors need to have the same size for import_elements() to work!"));
576  // TODO: a generic import_elements() function should handle any kind of
577  // data layout in ReadWriteVector, but this function is of limited use as
578  // this class will (hopefully) be retired eventually.
581 
582  if (operation == VectorOperation::insert)
583  {
584  for (const auto idx : this->locally_owned_elements())
585  (*this)[idx] = rwv[idx];
586  }
587  else if (operation == VectorOperation::add)
588  {
589  for (const auto idx : this->locally_owned_elements())
590  (*this)[idx] += rwv[idx];
591  }
592  else
593  AssertThrow(false, ExcNotImplemented());
594 
595  this->compress(operation);
596  }
597 
598 
599  void
601  {
602  // Select which mode to send to Trilinos. Note that we use last_action if
603  // available and ignore what the user tells us to detect wrongly mixed
604  // operations. Typically given_last_action is only used on machines that
605  // do not execute an operation (because they have no own cells for
606  // example).
607  Epetra_CombineMode mode = last_action;
608  if (last_action == Zero)
609  {
610  if (given_last_action == VectorOperation::add)
611  mode = Add;
612  else if (given_last_action == VectorOperation::insert)
613  mode = Insert;
614  else
615  Assert(
616  false,
617  ExcMessage(
618  "compress() can only be called with VectorOperation add, insert, or unknown"));
619  }
620  else
621  {
622  Assert(
623  ((last_action == Add) &&
624  (given_last_action == VectorOperation::add)) ||
625  ((last_action == Insert) &&
626  (given_last_action == VectorOperation::insert)),
627  ExcMessage(
628  "The last operation on the Vector and the given last action in the compress() call do not agree!"));
629  }
630 
631 
632 # ifdef DEBUG
633  // check that every process has decided to use the same mode. This will
634  // otherwise result in undefined behavior in the call to
635  // GlobalAssemble().
636  const double double_mode = mode;
637  const Epetra_MpiComm *comm_ptr =
638  dynamic_cast<const Epetra_MpiComm *>(&(trilinos_partitioner().Comm()));
639  Assert(comm_ptr != nullptr, ExcInternalError());
640 
641  const Utilities::MPI::MinMaxAvg result =
642  Utilities::MPI::min_max_avg(double_mode, comm_ptr->GetMpiComm());
643  Assert(result.max == result.min,
644  ExcMessage(
645  "Not all processors agree whether the last operation on "
646  "this vector was an addition or a set operation. This will "
647  "prevent the compress() operation from succeeding."));
648 
649 # endif
650 
651  // Now pass over the information about what we did last to the vector.
652  int ierr = 0;
653  if (nonlocal_vector.get() == nullptr || mode != Add)
654  ierr = vector->GlobalAssemble(mode);
655  else
656  {
657  Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
658  ierr = vector->Export(*nonlocal_vector, exporter, mode);
659  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
660  ierr = nonlocal_vector->PutScalar(0.);
661  }
662  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
663  last_action = Zero;
664 
665  compressed = true;
666  }
667 
668 
669 
671  Vector::operator()(const size_type index) const
672  {
673  // Extract local indices in the vector.
674  TrilinosWrappers::types::int_type trilinos_i = vector->Map().LID(
675  static_cast<TrilinosWrappers::types::int_type>(index));
676  TrilinosScalar value = 0.;
677 
678  // If the element is not present on the current processor, we can't
679  // continue. This is the main difference to the el() function.
680  if (trilinos_i == -1)
681  {
682  Assert(false,
684  vector->Map().NumMyElements(),
685  vector->Map().MinMyGID(),
686  vector->Map().MaxMyGID()));
687  }
688  else
689  value = (*vector)[0][trilinos_i];
690 
691  return value;
692  }
693 
694 
695 
696  void
697  Vector::add(const Vector &v, const bool allow_different_maps)
698  {
699  if (allow_different_maps == false)
700  *this += v;
701  else
702  {
704  AssertThrow(size() == v.size(),
705  ExcDimensionMismatch(size(), v.size()));
706 
707 # if DEAL_II_TRILINOS_VERSION_GTE(11, 11, 0)
708  Epetra_Import data_exchange(vector->Map(), v.vector->Map());
709  int ierr =
710  vector->Import(*v.vector, data_exchange, Epetra_AddLocalAlso);
711  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
712  last_action = Add;
713 # else
714  // In versions older than 11.11 the Import function is broken for
715  // adding Hence, we provide a workaround in this case
716 
717  Epetra_MultiVector dummy(vector->Map(), 1, false);
718  Epetra_Import data_exchange(dummy.Map(), v.vector->Map());
719 
720  int ierr = dummy.Import(*v.vector, data_exchange, Insert);
721  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
722 
723  ierr = vector->Update(1.0, dummy, 1.0);
724  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
725 # endif
726  }
727  }
728 
729 
730 
731  bool
732  Vector::operator==(const Vector &v) const
733  {
734  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
735  if (vector->Map().NumMyElements() != v.vector->Map().NumMyElements())
736  return false;
737 
738  size_type vector_size = vector->Map().NumMyElements();
739  for (size_type i = 0; i < vector_size; ++i)
740  if ((*(v.vector))[0][i] != (*vector)[0][i])
741  return false;
742 
743  return true;
744  }
745 
746 
747 
748  bool
749  Vector::operator!=(const Vector &v) const
750  {
751  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
752 
753  return (!(*this == v));
754  }
755 
756 
757 
758  bool
760  {
761  // get a representation of the vector and
762  // loop over all the elements
763  TrilinosScalar *start_ptr = (*vector)[0];
764  const TrilinosScalar *ptr = start_ptr,
765  *eptr = start_ptr + vector->Map().NumMyElements();
766  unsigned int flag = 0;
767  while (ptr != eptr)
768  {
769  if (*ptr != 0)
770  {
771  flag = 1;
772  break;
773  }
774  ++ptr;
775  }
776 
777  // in parallel, check that the vector
778  // is zero on _all_ processors.
779  const Epetra_MpiComm *mpi_comm =
780  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
781  Assert(mpi_comm != nullptr, ExcInternalError());
782  unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
783  return num_nonzero == 0;
784  }
785 
786 
787 
788  bool
790  {
791  // get a representation of the vector and
792  // loop over all the elements
793  TrilinosScalar *start_ptr = (*vector)[0];
794  const TrilinosScalar *ptr = start_ptr,
795  *eptr = start_ptr + vector->Map().NumMyElements();
796  unsigned int flag = 0;
797  while (ptr != eptr)
798  {
799  if (*ptr < 0.0)
800  {
801  flag = 1;
802  break;
803  }
804  ++ptr;
805  }
806 
807  // in parallel, check that the vector
808  // is zero on _all_ processors.
809  const auto max_n_negative =
811  return max_n_negative == 0;
812  }
813 
814 
815 
816  void
817  Vector::print(std::ostream &out,
818  const unsigned int precision,
819  const bool scientific,
820  const bool across) const
821  {
822  AssertThrow(out.fail() == false, ExcIO());
823  boost::io::ios_flags_saver restore_flags(out);
824 
825 
826  out.precision(precision);
827  if (scientific)
828  out.setf(std::ios::scientific, std::ios::floatfield);
829  else
830  out.setf(std::ios::fixed, std::ios::floatfield);
831 
832  size_type vector_size = vector->Map().NumMyElements();
833  if (size() != vector_size)
834  {
835  auto global_id = [&](const size_type index) {
836  return gid(vector->Map(), index);
837  };
838  out << "size:" << size()
839  << " locally_owned_size:" << vector->Map().NumMyElements() << " :"
840  << std::endl;
841  for (size_type i = 0; i < vector_size; ++i)
842  out << "[" << global_id(i) << "]: " << (*(vector))[0][i]
843  << std::endl;
844  }
845  else
846  {
847  TrilinosScalar *val;
848  int leading_dimension;
849  int ierr = vector->ExtractView(&val, &leading_dimension);
850 
851  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
852  if (across)
853  for (size_type i = 0; i < size(); ++i)
854  out << static_cast<double>(val[i]) << ' ';
855  else
856  for (size_type i = 0; i < size(); ++i)
857  out << static_cast<double>(val[i]) << std::endl;
858  out << std::endl;
859  }
860 
861  AssertThrow(out.fail() == false, ExcIO());
862  }
863 
864 
865 
866  void
868  {
872  std::swap(vector, v.vector);
875  }
876 
877 
878 
879  std::size_t
881  {
882  // TODO[TH]: No accurate memory
883  // consumption for Trilinos vectors
884  // yet. This is a rough approximation with
885  // one index and the value per local
886  // entry.
887  return sizeof(*this) +
888  this->vector->Map().NumMyElements() *
889  (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
890  }
891 
892  // explicit instantiations
893 # ifndef DOXYGEN
894 # include "trilinos_vector.inst"
895 # endif
896  } // namespace MPI
897 } // namespace TrilinosWrappers
898 
900 
901 #endif // DEAL_II_WITH_TRILINOS
virtual size_type size() const override
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
Definition: index_set.cc:929
size_type size() const
Definition: index_set.h:1696
size_type n_elements() const
Definition: index_set.h:1854
void set_size(const size_type size)
Definition: index_set.h:1684
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:840
void clear()
Definition: index_set.h:1672
void subtract_set(const IndexSet &other)
Definition: index_set.cc:288
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
size_type size() const override
const IndexSet & get_stored_elements() const
void compress(VectorOperation::values operation)
const Epetra_BlockMap & trilinos_partitioner() const
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
void import_elements(const LinearAlgebra::ReadWriteVector< double > &rwv, const VectorOperation::values operation)
MPI_Comm get_mpi_communicator() const
reference operator()(const size_type index)
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
std::unique_ptr< Epetra_FEVector > vector
size_type size() const override
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
IndexSet locally_owned_elements() const
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
bool operator!=(const Vector &v) const
bool operator==(const Vector &v) const
std::pair< size_type, size_type > local_range() const
Vector & operator=(const TrilinosScalar s)
::types::global_dof_index size_type
std::size_t memory_consumption() const
const Epetra_CrsMatrix & trilinos_matrix() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
void swap(BlockVector &u, BlockVector &v)
TrilinosWrappers::types::int_type global_length(const Epetra_MultiVector &vector)
int gid(const Epetra_BlockMap &map, int i)
TrilinosWrappers::types::int_type * my_global_elements(const Epetra_BlockMap &map)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
T max(const T &t, const MPI_Comm mpi_communicator)
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm mpi_communicator)
Definition: mpi.cc:123
const Epetra_Comm & comm_self()
double TrilinosScalar
Definition: types.h:175