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