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