Reference documentation for deal.II version GIT relicensing-43-g3fc2fae7c3 2024-03-03 21:10: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 // Check that the class we declare here satisfies the
41 // vector-space-vector concept. If we catch it here,
42 // any mistake in the vector class declaration would
43 // show up in uses of this class later on as well.
44# ifdef DEAL_II_HAVE_CXX20
46# endif
47
49 : Subscriptor()
50 , vector(new Epetra_FEVector(
51 Epetra_Map(0, 0, 0, Utilities::Trilinos::comm_self())))
52 {}
53
54
55
57 : Subscriptor()
58 , vector(new Epetra_FEVector(V.trilinos_vector()))
59 {}
60
61
62
63 Vector::Vector(const IndexSet &parallel_partitioner,
64 const MPI_Comm communicator)
65 : Subscriptor()
66 , vector(new Epetra_FEVector(
67 parallel_partitioner.make_trilinos_map(communicator, false)))
68 {}
69
70
71
72 void
73 Vector::reinit(const IndexSet &parallel_partitioner,
74 const MPI_Comm communicator,
75 const bool omit_zeroing_entries)
76 {
77 Epetra_Map input_map =
78 parallel_partitioner.make_trilinos_map(communicator, false);
79 if (vector->Map().SameAs(input_map) == false)
80 vector = std::make_unique<Epetra_FEVector>(input_map);
81 else if (omit_zeroing_entries == false)
82 {
83 const int ierr = vector->PutScalar(0.);
84 Assert(ierr == 0, ExcTrilinosError(ierr));
85 (void)ierr;
86 }
87 }
88
89
90
91 void
92 Vector::reinit(const Vector &V, const bool omit_zeroing_entries)
93 {
96 omit_zeroing_entries);
97 }
98
99
100
101 void
104 ArrayView<double> &elements) const
105 {
106 AssertDimension(indices.size(), elements.size());
107 const auto &vector = trilinos_vector();
108 const auto &map = vector.Map();
109
110 for (unsigned int i = 0; i < indices.size(); ++i)
111 {
112 AssertIndexRange(indices[i], size());
113 const auto trilinos_i =
114 map.LID(static_cast<TrilinosWrappers::types::int_type>(indices[i]));
115 elements[i] = vector[0][trilinos_i];
116 }
117 }
118
119
120
121 Vector &
123 {
124 // Distinguish three cases:
125 // - First case: both vectors have the same layout.
126 // - Second case: both vectors have the same size but different layout.
127 // - Third case: the vectors have different size.
128 if (vector->Map().SameAs(V.trilinos_vector().Map()))
129 *vector = V.trilinos_vector();
130 else
131 {
132 if (size() == V.size())
133 {
134 Epetra_Import data_exchange(vector->Map(),
135 V.trilinos_vector().Map());
136
137 const int ierr =
138 vector->Import(V.trilinos_vector(), data_exchange, Insert);
139 Assert(ierr == 0, ExcTrilinosError(ierr));
140 (void)ierr;
141 }
142 else
143 vector = std::make_unique<Epetra_FEVector>(V.trilinos_vector());
144 }
145
146 return *this;
147 }
148
149
150
151 Vector &
152 Vector::operator=(const double s)
153 {
154 Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector."));
155
156 const int ierr = vector->PutScalar(s);
157 Assert(ierr == 0, ExcTrilinosError(ierr));
158 (void)ierr;
159
160 return *this;
161 }
162
163
164
165 void
168 VectorOperation::values operation,
169 const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
170 &communication_pattern)
171 {
172 // If no communication pattern is given, create one. Otherwise, use the
173 // one given.
174 if (communication_pattern == nullptr)
175 {
176 // The first time import is called, a communication pattern is
177 // created. Check if the communication pattern already exists and if
178 // it can be reused.
180 V.get_stored_elements().size()) ||
182 {
185 dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
186 }
187 }
188 else
189 {
191 std::dynamic_pointer_cast<const CommunicationPattern>(
192 communication_pattern);
194 epetra_comm_pattern != nullptr,
196 std::string("The communication pattern is not of type ") +
197 "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
198 }
199
200 Epetra_Import import_map(epetra_comm_pattern->get_epetra_import());
201
202 // The TargetMap and the SourceMap have their roles inverted.
203 Epetra_FEVector source_vector(import_map.TargetMap());
204 double *values = source_vector.Values();
205 std::copy(V.begin(), V.end(), values);
206
207 if (operation == VectorOperation::insert)
208 vector->Export(source_vector, import_map, Insert);
209 else if (operation == VectorOperation::add)
210 vector->Export(source_vector, import_map, Add);
211 else if (operation == VectorOperation::max)
212 vector->Export(source_vector, import_map, Epetra_Max);
213 else if (operation == VectorOperation::min)
214 vector->Export(source_vector, import_map, Epetra_Min);
215 else
217 }
218
219
220
221 Vector &
222 Vector::operator*=(const double factor)
223 {
224 AssertIsFinite(factor);
225 vector->Scale(factor);
226
227 return *this;
228 }
229
230
231
232 Vector &
233 Vector::operator/=(const double factor)
234 {
235 AssertIsFinite(factor);
236 Assert(factor != 0., ExcZero());
237 *this *= 1. / factor;
238
239 return *this;
240 }
241
242
243
244 Vector &
246 {
247 // If the maps are the same we can Update right away.
248 if (vector->Map().SameAs(V.trilinos_vector().Map()))
249 {
250 const int ierr = vector->Update(1., V.trilinos_vector(), 1.);
251 Assert(ierr == 0, ExcTrilinosError(ierr));
252 (void)ierr;
253 }
254 else
255 {
256 Assert(this->size() == V.size(),
257 ExcDimensionMismatch(this->size(), V.size()));
258
259# if DEAL_II_TRILINOS_VERSION_GTE(11, 11, 0)
260 Epetra_Import data_exchange(vector->Map(), V.trilinos_vector().Map());
261 const int ierr = vector->Import(V.trilinos_vector(),
262 data_exchange,
263 Epetra_AddLocalAlso);
264 Assert(ierr == 0, ExcTrilinosError(ierr));
265 (void)ierr;
266# else
267 // In versions older than 11.11 the Import function is broken for
268 // adding Hence, we provide a workaround in this case
269
270 Epetra_MultiVector dummy(vector->Map(), 1, false);
271 Epetra_Import data_exchange(dummy.Map(), V.trilinos_vector().Map());
272
273 int ierr = dummy.Import(V.trilinos_vector(), data_exchange, Insert);
274 Assert(ierr == 0, ExcTrilinosError(ierr));
275
276 ierr = vector->Update(1.0, dummy, 1.0);
277 Assert(ierr == 0, ExcTrilinosError(ierr));
278 (void)ierr;
279# endif
280 }
281
282 return *this;
283 }
284
285
286
287 Vector &
289 {
290 this->add(-1., V);
291
292 return *this;
293 }
294
295
296
297 double
298 Vector::operator*(const Vector &V) const
299 {
300 Assert(this->size() == V.size(),
301 ExcDimensionMismatch(this->size(), V.size()));
302 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
304
305 double result(0.);
306 const int ierr = vector->Dot(V.trilinos_vector(), &result);
307 Assert(ierr == 0, ExcTrilinosError(ierr));
308 (void)ierr;
309
310 return result;
311 }
312
313
314
315 void
316 Vector::add(const double a)
317 {
319 const unsigned local_size(vector->MyLength());
320 for (unsigned int i = 0; i < local_size; ++i)
321 (*vector)[0][i] += a;
322 }
323
324
325
326 void
327 Vector::add(const double a, const Vector &V)
328 {
330 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
332
333 const int ierr = vector->Update(a, V.trilinos_vector(), 1.);
334 Assert(ierr == 0, ExcTrilinosError(ierr));
335 (void)ierr;
336 }
337
338
339
340 void
341 Vector::add(const double a,
342 const Vector &V,
343 const double b,
344 const Vector &W)
345 {
346 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
348 Assert(vector->Map().SameAs(W.trilinos_vector().Map()),
352
353 const int ierr =
354 vector->Update(a, V.trilinos_vector(), b, W.trilinos_vector(), 1.);
355 Assert(ierr == 0, ExcTrilinosError(ierr));
356 (void)ierr;
357 }
358
359
360
361 void
362 Vector::sadd(const double s, const double a, const Vector &V)
363 {
364 *this *= s;
365 Vector tmp(V);
366 tmp *= a;
367 *this += tmp;
368 }
369
370
371
372 void
373 Vector::scale(const Vector &scaling_factors)
374 {
375 Assert(vector->Map().SameAs(scaling_factors.trilinos_vector().Map()),
377
378 const int ierr =
379 vector->Multiply(1.0, scaling_factors.trilinos_vector(), *vector, 0.0);
380 Assert(ierr == 0, ExcTrilinosError(ierr));
381 (void)ierr;
382 }
383
384
385
386 void
387 Vector::equ(const double a, const Vector &V)
388 {
389 // If we don't have the same map, copy.
390 if (vector->Map().SameAs(V.trilinos_vector().Map()) == false)
391 this->sadd(0., a, V);
392 else
393 {
394 // Otherwise, just update
395 int ierr = vector->Update(a, V.trilinos_vector(), 0.);
396 Assert(ierr == 0, ExcTrilinosError(ierr));
397 (void)ierr;
398 }
399 }
400
401
402
403 bool
405 {
406 // get a representation of the vector and
407 // loop over all the elements
408 double *start_ptr = (*vector)[0];
409 const double *ptr = start_ptr, *eptr = start_ptr + vector->MyLength();
410 unsigned int flag = 0;
411 while (ptr != eptr)
412 {
413 if (*ptr != 0)
414 {
415 flag = 1;
416 break;
417 }
418 ++ptr;
419 }
420
421 // Check that the vector is zero on _all_ processors.
422 const Epetra_MpiComm *mpi_comm =
423 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
424 Assert(mpi_comm != nullptr, ExcInternalError());
425 unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
426
427 return num_nonzero == 0;
428 }
429
430
431
432 double
434 {
435 double mean_value(0.);
436
437 int ierr = vector->MeanValue(&mean_value);
438 Assert(ierr == 0, ExcTrilinosError(ierr));
439 (void)ierr;
440
441 return mean_value;
442 }
443
444
445
446 double
448 {
449 double norm(0.);
450 int ierr = vector->Norm1(&norm);
451 Assert(ierr == 0, ExcTrilinosError(ierr));
452 (void)ierr;
453
454 return norm;
455 }
456
457
458
459 double
461 {
462 double norm(0.);
463 int ierr = vector->Norm2(&norm);
464 Assert(ierr == 0, ExcTrilinosError(ierr));
465 (void)ierr;
466
467 return norm;
468 }
469
470
471
472 double
474 {
475 double norm(0.);
476 int ierr = vector->NormInf(&norm);
477 Assert(ierr == 0, ExcTrilinosError(ierr));
478 (void)ierr;
479
480 return norm;
481 }
482
483
484
485 double
486 Vector::add_and_dot(const double a, const Vector &V, const Vector &W)
487 {
488 this->add(a, V);
489
490 return *this * W;
491 }
492
493
494
497 {
498# ifndef DEAL_II_WITH_64BIT_INDICES
499 return vector->GlobalLength();
500# else
501 return vector->GlobalLength64();
502# endif
503 }
504
505
506
509 {
510 return vector->MyLength();
511 }
512
513
514
517 {
518 const Epetra_MpiComm *epetra_comm =
519 dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
520 Assert(epetra_comm != nullptr, ExcInternalError());
521 return epetra_comm->GetMpiComm();
522 }
523
524
525
528 {
529 IndexSet is(size());
530
531 // easy case: local range is contiguous
532 if (vector->Map().LinearMap())
533 {
534# ifndef DEAL_II_WITH_64BIT_INDICES
535 is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID() + 1);
536# else
537 is.add_range(vector->Map().MinMyGID64(),
538 vector->Map().MaxMyGID64() + 1);
539# endif
540 }
541 else if (vector->Map().NumMyElements() > 0)
542 {
543 const size_type n_indices = vector->Map().NumMyElements();
544# ifndef DEAL_II_WITH_64BIT_INDICES
545 unsigned int *vector_indices =
546 reinterpret_cast<unsigned int *>(vector->Map().MyGlobalElements());
547# else
548 size_type *vector_indices =
549 reinterpret_cast<size_type *>(vector->Map().MyGlobalElements64());
550# endif
551 is.add_indices(vector_indices, vector_indices + n_indices);
552 }
553 is.compress();
554
555 return is;
556 }
557
558
559 void
561 {}
562
563
564
565 const Epetra_FEVector &
567 {
568 return *vector;
569 }
570
571
572
573 Epetra_FEVector &
575 {
576 return *vector;
577 }
578
579
580
581 void
582 Vector::print(std::ostream &out,
583 const unsigned int precision,
584 const bool scientific,
585 const bool across) const
586 {
587 AssertThrow(out.fail() == false, ExcIO());
588 boost::io::ios_flags_saver restore_flags(out);
589
590 // Get a representation of the vector and loop over all
591 // the elements
592 double *val;
593 int leading_dimension;
594 int ierr = vector->ExtractView(&val, &leading_dimension);
595
596 Assert(ierr == 0, ExcTrilinosError(ierr));
597 (void)ierr;
598 out.precision(precision);
599 if (scientific)
600 out.setf(std::ios::scientific, std::ios::floatfield);
601 else
602 out.setf(std::ios::fixed, std::ios::floatfield);
603
604 if (across)
605 for (int i = 0; i < vector->MyLength(); ++i)
606 out << val[i] << ' ';
607 else
608 for (int i = 0; i < vector->MyLength(); ++i)
609 out << val[i] << std::endl;
610 out << std::endl;
611
612 // restore the representation
613 // of the vector
614 AssertThrow(out.fail() == false, ExcIO());
615 }
616
617
618
619 std::size_t
621 {
622 return sizeof(*this) +
623 vector->MyLength() *
624 (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
625 }
626
627
628
629 void
631 const MPI_Comm mpi_comm)
632 {
633 source_stored_elements = source_index_set;
635 std::make_shared<CommunicationPattern>(locally_owned_elements(),
636 source_index_set,
637 mpi_comm);
638 }
639 } // namespace EpetraWrappers
640} // namespace LinearAlgebra
641
643
644#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)