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