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