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