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