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