deal.II version GIT relicensing-2581-gae2745de1b 2025-02-07 22:30: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\}}\)
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 - 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
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 : vector(new Epetra_FEVector(
85 Epetra_Map(0, 0, 0, Utilities::Trilinos::comm_self())))
86 {}
87
88
89
91 : vector(new Epetra_FEVector(V.trilinos_vector()))
92 {}
93
94
95
96 Vector::Vector(const IndexSet &parallel_partitioner,
97 const MPI_Comm communicator)
98 : vector(new Epetra_FEVector(
99 parallel_partitioner.make_trilinos_map(communicator, false)))
100 {}
101
102
103
104 void
105 Vector::reinit(const IndexSet &parallel_partitioner,
106 const MPI_Comm communicator,
107 const bool omit_zeroing_entries)
108 {
109 Epetra_Map input_map =
110 parallel_partitioner.make_trilinos_map(communicator, false);
111 if (vector->Map().SameAs(input_map) == false)
112 vector = std::make_unique<Epetra_FEVector>(input_map);
113 else if (omit_zeroing_entries == false)
114 {
115 const int ierr = vector->PutScalar(0.);
117 (void)ierr;
118 }
119 }
120
121
122
123 void
124 Vector::reinit(const Vector &V, const bool omit_zeroing_entries)
125 {
129 }
130
131
132
133 void
136 const ArrayView<double> &elements) const
137 {
138 AssertDimension(indices.size(), elements.size());
139 const auto &vector = trilinos_vector();
140 const auto &map = vector.Map();
141
142 for (unsigned int i = 0; i < indices.size(); ++i)
143 {
144 AssertIndexRange(indices[i], size());
145 const auto trilinos_i =
146 map.LID(static_cast<TrilinosWrappers::types::int_type>(indices[i]));
147 elements[i] = vector[0][trilinos_i];
148 }
149 }
150
151
152
153 Vector &
155 {
156 // Distinguish three cases:
157 // - First case: both vectors have the same layout.
158 // - Second case: both vectors have the same size but different layout.
159 // - Third case: the vectors have different size.
160 if (vector->Map().SameAs(V.trilinos_vector().Map()))
161 *vector = V.trilinos_vector();
162 else
163 {
164 if (size() == V.size())
165 {
166 Epetra_Import data_exchange(vector->Map(),
167 V.trilinos_vector().Map());
168
169 const int ierr =
172 (void)ierr;
173 }
174 else
175 vector = std::make_unique<Epetra_FEVector>(V.trilinos_vector());
176 }
177
178 return *this;
179 }
180
181
182
183 Vector &
184 Vector::operator=(const double s)
185 {
186 Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector."));
187
188 const int ierr = vector->PutScalar(s);
190 (void)ierr;
191
192 return *this;
193 }
194
195
196
197 void
200 VectorOperation::values operation,
201 const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
202 &communication_pattern)
203 {
204 // If no communication pattern is given, create one. Otherwise, use the
205 // one given.
206 if (communication_pattern == nullptr)
207 {
208 // The first time import is called, a communication pattern is
209 // created. Check if the communication pattern already exists and if
210 // it can be reused.
212 V.get_stored_elements().size()) ||
214 {
217 dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
218 }
219 }
220 else
221 {
223 std::dynamic_pointer_cast<const CommunicationPattern>(
226 epetra_comm_pattern != nullptr,
227 ExcMessage("The communication pattern is not of type "
228 "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
229 }
230
231 Epetra_Import import_map(epetra_comm_pattern->get_epetra_import());
232
233 // The TargetMap and the SourceMap have their roles inverted.
234 Epetra_FEVector source_vector(import_map.TargetMap());
235 double *values = source_vector.Values();
236 std::copy(V.begin(), V.end(), values);
237
246 else
248 }
249
250
251
252 Vector &
253 Vector::operator*=(const double factor)
254 {
255 AssertIsFinite(factor);
256 vector->Scale(factor);
257
258 return *this;
259 }
260
261
262
263 Vector &
264 Vector::operator/=(const double factor)
265 {
266 AssertIsFinite(factor);
267 Assert(factor != 0., ExcZero());
268 *this *= 1. / factor;
269
270 return *this;
271 }
272
273
274
275 Vector &
277 {
278 // If the maps are the same we can Update right away.
279 if (vector->Map().SameAs(V.trilinos_vector().Map()))
280 {
281 const int ierr = vector->Update(1., V.trilinos_vector(), 1.);
283 (void)ierr;
284 }
285 else
286 {
287 Assert(this->size() == V.size(),
288 ExcDimensionMismatch(this->size(), V.size()));
289
290 Epetra_Import data_exchange(vector->Map(), V.trilinos_vector().Map());
291 const int ierr = vector->Import(V.trilinos_vector(),
295 (void)ierr;
296 }
297
298 return *this;
299 }
300
301
302
303 Vector &
305 {
306 this->add(-1., V);
307
308 return *this;
309 }
310
311
312
313 double
314 Vector::operator*(const Vector &V) const
315 {
316 Assert(this->size() == V.size(),
317 ExcDimensionMismatch(this->size(), V.size()));
318 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
320
321 double result(0.);
322 const int ierr = vector->Dot(V.trilinos_vector(), &result);
324 (void)ierr;
325
326 return result;
327 }
328
329
330
331 void
332 Vector::add(const double a)
333 {
335 const unsigned local_size(vector->MyLength());
336 for (unsigned int i = 0; i < local_size; ++i)
337 (*vector)[0][i] += a;
338 }
339
340
341
342 void
343 Vector::add(const double a, const Vector &V)
344 {
346 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
348
349 const int ierr = vector->Update(a, V.trilinos_vector(), 1.);
351 (void)ierr;
352 }
353
354
355
356 void
357 Vector::add(const double a,
358 const Vector &V,
359 const double b,
360 const Vector &W)
361 {
362 Assert(vector->Map().SameAs(V.trilinos_vector().Map()),
364 Assert(vector->Map().SameAs(W.trilinos_vector().Map()),
368
369 const int ierr =
370 vector->Update(a, V.trilinos_vector(), b, W.trilinos_vector(), 1.);
372 (void)ierr;
373 }
374
375
376
377 void
378 Vector::sadd(const double s, const double a, const Vector &V)
379 {
380 *this *= s;
381 Vector tmp(V);
382 tmp *= a;
383 *this += tmp;
384 }
385
386
387
388 void
389 Vector::scale(const Vector &scaling_factors)
390 {
391 Assert(vector->Map().SameAs(scaling_factors.trilinos_vector().Map()),
393
394 const int ierr =
395 vector->Multiply(1.0, scaling_factors.trilinos_vector(), *vector, 0.0);
397 (void)ierr;
398 }
399
400
401
402 void
403 Vector::equ(const double a, const Vector &V)
404 {
405 // If we don't have the same map, copy.
406 if (vector->Map().SameAs(V.trilinos_vector().Map()) == false)
407 this->sadd(0., a, V);
408 else
409 {
410 // Otherwise, just update
411 int ierr = vector->Update(a, V.trilinos_vector(), 0.);
413 (void)ierr;
414 }
415 }
416
417
418
419 bool
421 {
422 // get a representation of the vector and
423 // loop over all the elements
424 double *start_ptr = (*vector)[0];
425 const double *ptr = start_ptr, *eptr = start_ptr + vector->MyLength();
426 unsigned int flag = 0;
427 while (ptr != eptr)
428 {
429 if (*ptr != 0)
430 {
431 flag = 1;
432 break;
433 }
434 ++ptr;
435 }
436
437 // Check that the vector is zero on _all_ processors.
438 const Epetra_MpiComm *mpi_comm =
439 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
440 Assert(mpi_comm != nullptr, ExcInternalError());
441 unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
442
443 return num_nonzero == 0;
444 }
445
446
447
448 double
450 {
451 double mean_value(0.);
452
453 int ierr = vector->MeanValue(&mean_value);
455 (void)ierr;
456
457 return mean_value;
458 }
459
460
461
462 double
464 {
465 double norm(0.);
466 int ierr = vector->Norm1(&norm);
468 (void)ierr;
469
470 return norm;
471 }
472
473
474
475 double
477 {
478 double norm(0.);
479 int ierr = vector->Norm2(&norm);
481 (void)ierr;
482
483 return norm;
484 }
485
486
487
488 double
490 {
491 double norm(0.);
492 int ierr = vector->NormInf(&norm);
494 (void)ierr;
495
496 return norm;
497 }
498
499
500
501 double
502 Vector::add_and_dot(const double a, const Vector &V, const Vector &W)
503 {
504 this->add(a, V);
505
506 return *this * W;
507 }
508
509
510
513 {
514# ifndef DEAL_II_WITH_64BIT_INDICES
515 return vector->GlobalLength();
516# else
517 return vector->GlobalLength64();
518# endif
519 }
520
521
522
525 {
526 return vector->MyLength();
527 }
528
529
530
533 {
534 const Epetra_MpiComm *epetra_comm =
535 dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
536 Assert(epetra_comm != nullptr, ExcInternalError());
537 return epetra_comm->GetMpiComm();
538 }
539
540
541
544 {
545 IndexSet is(size());
546
547 // easy case: local range is contiguous
548 if (vector->Map().LinearMap())
549 {
550# ifndef DEAL_II_WITH_64BIT_INDICES
551 is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID() + 1);
552# else
553 is.add_range(vector->Map().MinMyGID64(),
554 vector->Map().MaxMyGID64() + 1);
555# endif
556 }
557 else if (vector->Map().NumMyElements() > 0)
558 {
559 const size_type n_indices = vector->Map().NumMyElements();
560# ifndef DEAL_II_WITH_64BIT_INDICES
561 unsigned int *vector_indices =
562 reinterpret_cast<unsigned int *>(vector->Map().MyGlobalElements());
563# else
565 reinterpret_cast<size_type *>(vector->Map().MyGlobalElements64());
566# endif
568 }
569 is.compress();
570
571 return is;
572 }
573
574
575 void
577 {}
578
579
580
581 const Epetra_FEVector &
583 {
584 return *vector;
585 }
586
587
588
589 Epetra_FEVector &
591 {
592 return *vector;
593 }
594
595
596
597 void
598 Vector::print(std::ostream &out,
599 const unsigned int precision,
600 const bool scientific,
601 const bool across) const
602 {
603 AssertThrow(out.fail() == false, ExcIO());
604 boost::io::ios_flags_saver restore_flags(out);
605
606 // Get a representation of the vector and loop over all
607 // the elements
608 double *val;
610 int ierr = vector->ExtractView(&val, &leading_dimension);
611
613 (void)ierr;
614 out.precision(precision);
615 if (scientific)
616 out.setf(std::ios::scientific, std::ios::floatfield);
617 else
618 out.setf(std::ios::fixed, std::ios::floatfield);
619
620 if (across)
621 for (int i = 0; i < vector->MyLength(); ++i)
622 out << val[i] << ' ';
623 else
624 for (int i = 0; i < vector->MyLength(); ++i)
625 out << val[i] << std::endl;
626 out << std::endl;
627
628 // restore the representation
629 // of the vector
630 AssertThrow(out.fail() == false, ExcIO());
631 }
632
633
634
635 std::size_t
637 {
638 return sizeof(*this) +
639 vector->MyLength() *
640 (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
641 }
642
643
644
645 void
647 const MPI_Comm mpi_comm)
648 {
651 std::make_shared<CommunicationPattern>(locally_owned_elements(),
653 mpi_comm);
654 }
655 } // namespace EpetraWrappers
656} // namespace LinearAlgebra
657
659
660#endif
std::size_t size() const
Definition array_view.h:684
size_type size() const
Definition index_set.h:1786
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
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)
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, const ArrayView< double > &elements) const override
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:513
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:514
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)