Reference documentation for deal.II version GIT relicensing-468-ge3ce94fd06 2024-04-23 15: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\}}\)
Loading...
Searching...
No Matches
trilinos_epetra_vector.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 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
20# include <deal.II/base/mpi.h>
22
24
25# include <boost/io/ios_state.hpp>
26
27# include <Epetra_Import.h>
28# include <Epetra_Map.h>
29# include <Epetra_MpiComm.h>
30
31# include <memory>
32
33
35
36namespace LinearAlgebra
37{
38 namespace EpetraWrappers
39 {
40# ifndef DOXYGEN
41 namespace internal
42 {
43 VectorReference::operator value_type() const
44 {
45 AssertIndexRange(index, 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
55# ifndef DEAL_II_WITH_64BIT_INDICES
56 Assert(local_index >= 0,
57 ExcAccessToNonLocalElement(index,
58 vector.vector->Map().NumMyElements(),
59 vector.vector->Map().MinMyGID(),
60 vector.vector->Map().MaxMyGID()));
61# else
62 Assert(local_index >= 0,
63 ExcAccessToNonLocalElement(index,
64 vector.vector->Map().NumMyElements(),
65 vector.vector->Map().MinMyGID64(),
66 vector.vector->Map().MaxMyGID64()));
67# endif
68
69 return (*(vector.vector))[0][local_index];
70 }
71 } // namespace internal
72# endif
73
74
75 // Check that the class we declare here satisfies the
76 // vector-space-vector concept. If we catch it here,
77 // any mistake in the vector class declaration would
78 // show up in uses of this class later on as well.
79# ifdef DEAL_II_HAVE_CXX20
81# endif
82
84 : Subscriptor()
85 , vector(new Epetra_FEVector(
86 Epetra_Map(0, 0, 0, Utilities::Trilinos::comm_self())))
87 {}
88
89
90
92 : Subscriptor()
93 , vector(new Epetra_FEVector(V.trilinos_vector()))
94 {}
95
96
97
98 Vector::Vector(const IndexSet &parallel_partitioner,
99 const MPI_Comm communicator)
100 : Subscriptor()
101 , vector(new Epetra_FEVector(
102 parallel_partitioner.make_trilinos_map(communicator, false)))
103 {}
104
105
106
107 void
108 Vector::reinit(const IndexSet &parallel_partitioner,
109 const MPI_Comm communicator,
110 const bool omit_zeroing_entries)
111 {
112 Epetra_Map input_map =
113 parallel_partitioner.make_trilinos_map(communicator, false);
114 if (vector->Map().SameAs(input_map) == false)
115 vector = std::make_unique<Epetra_FEVector>(input_map);
116 else if (omit_zeroing_entries == false)
117 {
118 const int ierr = vector->PutScalar(0.);
119 Assert(ierr == 0, ExcTrilinosError(ierr));
120 (void)ierr;
121 }
122 }
123
124
125
126 void
127 Vector::reinit(const Vector &V, const bool omit_zeroing_entries)
128 {
131 omit_zeroing_entries);
132 }
133
134
135
136 void
139 ArrayView<double> &elements) const
140 {
141 AssertDimension(indices.size(), elements.size());
142 const auto &vector = trilinos_vector();
143 const auto &map = vector.Map();
144
145 for (unsigned int i = 0; i < indices.size(); ++i)
146 {
147 AssertIndexRange(indices[i], size());
148 const auto trilinos_i =
149 map.LID(static_cast<TrilinosWrappers::types::int_type>(indices[i]));
150 elements[i] = vector[0][trilinos_i];
151 }
152 }
153
154
155
156 Vector &
158 {
159 // Distinguish three cases:
160 // - First case: both vectors have the same layout.
161 // - Second case: both vectors have the same size but different layout.
162 // - Third case: the vectors have different size.
163 if (vector->Map().SameAs(V.trilinos_vector().Map()))
164 *vector = V.trilinos_vector();
165 else
166 {
167 if (size() == V.size())
168 {
169 Epetra_Import data_exchange(vector->Map(),
170 V.trilinos_vector().Map());
171
172 const int ierr =
173 vector->Import(V.trilinos_vector(), data_exchange, Insert);
174 Assert(ierr == 0, ExcTrilinosError(ierr));
175 (void)ierr;
176 }
177 else
178 vector = std::make_unique<Epetra_FEVector>(V.trilinos_vector());
179 }
180
181 return *this;
182 }
183
184
185
186 Vector &
187 Vector::operator=(const double s)
188 {
189 Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector."));
190
191 const int ierr = vector->PutScalar(s);
192 Assert(ierr == 0, ExcTrilinosError(ierr));
193 (void)ierr;
194
195 return *this;
196 }
197
198
199
200 void
203 VectorOperation::values operation,
204 const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
205 &communication_pattern)
206 {
207 // If no communication pattern is given, create one. Otherwise, use the
208 // one given.
209 if (communication_pattern == nullptr)
210 {
211 // The first time import is called, a communication pattern is
212 // created. Check if the communication pattern already exists and if
213 // it can be reused.
215 V.get_stored_elements().size()) ||
217 {
220 dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
221 }
222 }
223 else
224 {
226 std::dynamic_pointer_cast<const CommunicationPattern>(
227 communication_pattern);
229 epetra_comm_pattern != nullptr,
231 std::string("The communication pattern is not of type ") +
232 "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
233 }
234
235 Epetra_Import import_map(epetra_comm_pattern->get_epetra_import());
236
237 // The TargetMap and the SourceMap have their roles inverted.
238 Epetra_FEVector source_vector(import_map.TargetMap());
239 double *values = source_vector.Values();
240 std::copy(V.begin(), V.end(), values);
241
242 if (operation == VectorOperation::insert)
243 vector->Export(source_vector, import_map, Insert);
244 else if (operation == VectorOperation::add)
245 vector->Export(source_vector, import_map, Add);
246 else if (operation == VectorOperation::max)
247 vector->Export(source_vector, import_map, Epetra_Max);
248 else if (operation == VectorOperation::min)
249 vector->Export(source_vector, import_map, Epetra_Min);
250 else
252 }
253
254
255
256 Vector &
257 Vector::operator*=(const double factor)
258 {
259 AssertIsFinite(factor);
260 vector->Scale(factor);
261
262 return *this;
263 }
264
265
266
267 Vector &
268 Vector::operator/=(const double factor)
269 {
270 AssertIsFinite(factor);
271 Assert(factor != 0., ExcZero());
272 *this *= 1. / factor;
273
274 return *this;
275 }
276
277
278
279 Vector &
281 {
282 // If the maps are the same we can Update right away.
283 if (vector->Map().SameAs(V.trilinos_vector().Map()))
284 {
285 const int ierr = vector->Update(1., V.trilinos_vector(), 1.);
286 Assert(ierr == 0, ExcTrilinosError(ierr));
287 (void)ierr;
288 }
289 else
290 {
291 Assert(this->size() == V.size(),
292 ExcDimensionMismatch(this->size(), V.size()));
293
294# if DEAL_II_TRILINOS_VERSION_GTE(11, 11, 0)
295 Epetra_Import data_exchange(vector->Map(), V.trilinos_vector().Map());
296 const int ierr = vector->Import(V.trilinos_vector(),
297 data_exchange,
298 Epetra_AddLocalAlso);
299 Assert(ierr == 0, ExcTrilinosError(ierr));
300 (void)ierr;
301# else
302 // In versions older than 11.11 the Import function is broken for
303 // adding Hence, we provide a workaround in this case
304
305 Epetra_MultiVector dummy(vector->Map(), 1, false);
306 Epetra_Import data_exchange(dummy.Map(), V.trilinos_vector().Map());
307
308 int ierr = dummy.Import(V.trilinos_vector(), data_exchange, Insert);
309 Assert(ierr == 0, ExcTrilinosError(ierr));
310
311 ierr = vector->Update(1.0, dummy, 1.0);
312 Assert(ierr == 0, ExcTrilinosError(ierr));
313 (void)ierr;
314# endif
315 }
316
317 return *this;
318 }
319
320
321
322 Vector &
324 {
325 this->add(-1., V);
326
327 return *this;
328 }
329
330
331
332 double
333 Vector::operator*(const Vector &V) const
334 {
335 Assert(this->size() == V.size(),
336 ExcDimensionMismatch(this->size(), V.size()));
337 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
339
340 double result(0.);
341 const int ierr = vector->Dot(V.trilinos_vector(), &result);
342 Assert(ierr == 0, ExcTrilinosError(ierr));
343 (void)ierr;
344
345 return result;
346 }
347
348
349
350 void
351 Vector::add(const double a)
352 {
354 const unsigned local_size(vector->MyLength());
355 for (unsigned int i = 0; i < local_size; ++i)
356 (*vector)[0][i] += a;
357 }
358
359
360
361 void
362 Vector::add(const double a, const Vector &V)
363 {
365 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
367
368 const int ierr = vector->Update(a, V.trilinos_vector(), 1.);
369 Assert(ierr == 0, ExcTrilinosError(ierr));
370 (void)ierr;
371 }
372
373
374
375 void
376 Vector::add(const double a,
377 const Vector &V,
378 const double b,
379 const Vector &W)
380 {
381 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
383 Assert(vector->Map().SameAs(W.trilinos_vector().Map()),
387
388 const int ierr =
389 vector->Update(a, V.trilinos_vector(), b, W.trilinos_vector(), 1.);
390 Assert(ierr == 0, ExcTrilinosError(ierr));
391 (void)ierr;
392 }
393
394
395
396 void
397 Vector::sadd(const double s, const double a, const Vector &V)
398 {
399 *this *= s;
400 Vector tmp(V);
401 tmp *= a;
402 *this += tmp;
403 }
404
405
406
407 void
408 Vector::scale(const Vector &scaling_factors)
409 {
410 Assert(vector->Map().SameAs(scaling_factors.trilinos_vector().Map()),
412
413 const int ierr =
414 vector->Multiply(1.0, scaling_factors.trilinos_vector(), *vector, 0.0);
415 Assert(ierr == 0, ExcTrilinosError(ierr));
416 (void)ierr;
417 }
418
419
420
421 void
422 Vector::equ(const double a, const Vector &V)
423 {
424 // If we don't have the same map, copy.
425 if (vector->Map().SameAs(V.trilinos_vector().Map()) == false)
426 this->sadd(0., a, V);
427 else
428 {
429 // Otherwise, just update
430 int ierr = vector->Update(a, V.trilinos_vector(), 0.);
431 Assert(ierr == 0, ExcTrilinosError(ierr));
432 (void)ierr;
433 }
434 }
435
436
437
438 bool
440 {
441 // get a representation of the vector and
442 // loop over all the elements
443 double *start_ptr = (*vector)[0];
444 const double *ptr = start_ptr, *eptr = start_ptr + vector->MyLength();
445 unsigned int flag = 0;
446 while (ptr != eptr)
447 {
448 if (*ptr != 0)
449 {
450 flag = 1;
451 break;
452 }
453 ++ptr;
454 }
455
456 // Check that the vector is zero on _all_ processors.
457 const Epetra_MpiComm *mpi_comm =
458 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
459 Assert(mpi_comm != nullptr, ExcInternalError());
460 unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
461
462 return num_nonzero == 0;
463 }
464
465
466
467 double
469 {
470 double mean_value(0.);
471
472 int ierr = vector->MeanValue(&mean_value);
473 Assert(ierr == 0, ExcTrilinosError(ierr));
474 (void)ierr;
475
476 return mean_value;
477 }
478
479
480
481 double
483 {
484 double norm(0.);
485 int ierr = vector->Norm1(&norm);
486 Assert(ierr == 0, ExcTrilinosError(ierr));
487 (void)ierr;
488
489 return norm;
490 }
491
492
493
494 double
496 {
497 double norm(0.);
498 int ierr = vector->Norm2(&norm);
499 Assert(ierr == 0, ExcTrilinosError(ierr));
500 (void)ierr;
501
502 return norm;
503 }
504
505
506
507 double
509 {
510 double norm(0.);
511 int ierr = vector->NormInf(&norm);
512 Assert(ierr == 0, ExcTrilinosError(ierr));
513 (void)ierr;
514
515 return norm;
516 }
517
518
519
520 double
521 Vector::add_and_dot(const double a, const Vector &V, const Vector &W)
522 {
523 this->add(a, V);
524
525 return *this * W;
526 }
527
528
529
532 {
533# ifndef DEAL_II_WITH_64BIT_INDICES
534 return vector->GlobalLength();
535# else
536 return vector->GlobalLength64();
537# endif
538 }
539
540
541
544 {
545 return vector->MyLength();
546 }
547
548
549
552 {
553 const Epetra_MpiComm *epetra_comm =
554 dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
555 Assert(epetra_comm != nullptr, ExcInternalError());
556 return epetra_comm->GetMpiComm();
557 }
558
559
560
563 {
564 IndexSet is(size());
565
566 // easy case: local range is contiguous
567 if (vector->Map().LinearMap())
568 {
569# ifndef DEAL_II_WITH_64BIT_INDICES
570 is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID() + 1);
571# else
572 is.add_range(vector->Map().MinMyGID64(),
573 vector->Map().MaxMyGID64() + 1);
574# endif
575 }
576 else if (vector->Map().NumMyElements() > 0)
577 {
578 const size_type n_indices = vector->Map().NumMyElements();
579# ifndef DEAL_II_WITH_64BIT_INDICES
580 unsigned int *vector_indices =
581 reinterpret_cast<unsigned int *>(vector->Map().MyGlobalElements());
582# else
583 size_type *vector_indices =
584 reinterpret_cast<size_type *>(vector->Map().MyGlobalElements64());
585# endif
586 is.add_indices(vector_indices, vector_indices + n_indices);
587 }
588 is.compress();
589
590 return is;
591 }
592
593
594 void
596 {}
597
598
599
600 const Epetra_FEVector &
602 {
603 return *vector;
604 }
605
606
607
608 Epetra_FEVector &
610 {
611 return *vector;
612 }
613
614
615
616 void
617 Vector::print(std::ostream &out,
618 const unsigned int precision,
619 const bool scientific,
620 const bool across) const
621 {
622 AssertThrow(out.fail() == false, ExcIO());
623 boost::io::ios_flags_saver restore_flags(out);
624
625 // Get a representation of the vector and loop over all
626 // the elements
627 double *val;
628 int leading_dimension;
629 int ierr = vector->ExtractView(&val, &leading_dimension);
630
631 Assert(ierr == 0, ExcTrilinosError(ierr));
632 (void)ierr;
633 out.precision(precision);
634 if (scientific)
635 out.setf(std::ios::scientific, std::ios::floatfield);
636 else
637 out.setf(std::ios::fixed, std::ios::floatfield);
638
639 if (across)
640 for (int i = 0; i < vector->MyLength(); ++i)
641 out << val[i] << ' ';
642 else
643 for (int i = 0; i < vector->MyLength(); ++i)
644 out << val[i] << std::endl;
645 out << std::endl;
646
647 // restore the representation
648 // of the vector
649 AssertThrow(out.fail() == false, ExcIO());
650 }
651
652
653
654 std::size_t
656 {
657 return sizeof(*this) +
658 vector->MyLength() *
659 (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
660 }
661
662
663
664 void
666 const MPI_Comm mpi_comm)
667 {
668 source_stored_elements = source_index_set;
670 std::make_shared<CommunicationPattern>(locally_owned_elements(),
671 source_index_set,
672 mpi_comm);
673 }
674 } // namespace EpetraWrappers
675} // namespace LinearAlgebra
676
678
679#endif
std::size_t size() const
Definition array_view.h:684
size_type size() const
Definition index_set.h:1765
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1792
void compress() const
Definition index_set.h:1773
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
void equ(const double a, const Vector &V)
void compress(const VectorOperation::values operation)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
std::unique_ptr< Epetra_FEVector > vector
const Epetra_FEVector & trilinos_vector() const
void import_elements(const ReadWriteVector< double > &V, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
void sadd(const double s, const double a, const Vector &V)
virtual size_type size() const override
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< double > &elements) const override
void scale(const Vector &scaling_factors)
void create_epetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
std::shared_ptr< const CommunicationPattern > epetra_comm_pattern
void reinit(const IndexSet &parallel_partitioner, const MPI_Comm communicator, const bool omit_zeroing_entries=false)
double add_and_dot(const double a, const Vector &V, const Vector &W)
const IndexSet & get_stored_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
static ::ExceptionBase & ExcDifferentParallelPartitioning()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
T sum(const T &t, const MPI_Comm mpi_communicator)