Reference documentation for deal.II version Git 8596a7cd07 2020-12-04 07:30:43 +0100
\(\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 - 2019 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 # ifdef DEAL_II_WITH_MPI
21 
22 # include <deal.II/base/index_set.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 CommunicationPatternBase> communication_pattern)
151  {
152  // If no communication pattern is given, create one. Otherwise, use the
153  // one given.
154  if (communication_pattern == nullptr)
155  {
156  // The first time import is called, a communication pattern is
157  // created. Check if the communication pattern already exists and if
158  // it can be reused.
159  if ((source_stored_elements.size() !=
160  V.get_stored_elements().size()) ||
162  {
165  dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
166  }
167  }
168  else
169  {
171  std::dynamic_pointer_cast<const CommunicationPattern>(
172  communication_pattern);
173  AssertThrow(
174  epetra_comm_pattern != nullptr,
175  ExcMessage(
176  std::string("The communication pattern is not of type ") +
177  "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
178  }
179 
180  Epetra_Import import(epetra_comm_pattern->get_epetra_import());
181 
182  // The TargetMap and the SourceMap have their roles inverted.
183  Epetra_FEVector source_vector(import.TargetMap());
184  double * values = source_vector.Values();
185  std::copy(V.begin(), V.end(), values);
186 
187  if (operation == VectorOperation::insert)
188  vector->Export(source_vector, import, Insert);
189  else if (operation == VectorOperation::add)
190  vector->Export(source_vector, import, Add);
191  else if (operation == VectorOperation::max)
192  vector->Export(source_vector, import, Epetra_Max);
193  else if (operation == VectorOperation::min)
194  vector->Export(source_vector, import, Epetra_Min);
195  else
196  AssertThrow(false, ExcNotImplemented());
197  }
198 
199 
200 
201  Vector &
202  Vector::operator*=(const double factor)
203  {
204  AssertIsFinite(factor);
205  vector->Scale(factor);
206 
207  return *this;
208  }
209 
210 
211 
212  Vector &
213  Vector::operator/=(const double factor)
214  {
215  AssertIsFinite(factor);
216  Assert(factor != 0., ExcZero());
217  *this *= 1. / factor;
218 
219  return *this;
220  }
221 
222 
223 
224  Vector &
226  {
227  // Check that casting will work.
228  Assert(dynamic_cast<const Vector *>(&V) != nullptr,
230 
231  // Downcast V. If fails, throws an exception.
232  const Vector &down_V = dynamic_cast<const Vector &>(V);
233  // If the maps are the same we can Update right away.
234  if (vector->Map().SameAs(down_V.trilinos_vector().Map()))
235  {
236  const int ierr = vector->Update(1., down_V.trilinos_vector(), 1.);
237  Assert(ierr == 0, ExcTrilinosError(ierr));
238  (void)ierr;
239  }
240  else
241  {
242  Assert(this->size() == down_V.size(),
243  ExcDimensionMismatch(this->size(), down_V.size()));
244 
245 # if DEAL_II_TRILINOS_VERSION_GTE(11, 11, 0)
246  Epetra_Import data_exchange(vector->Map(),
247  down_V.trilinos_vector().Map());
248  const int ierr = vector->Import(down_V.trilinos_vector(),
249  data_exchange,
250  Epetra_AddLocalAlso);
251  Assert(ierr == 0, ExcTrilinosError(ierr));
252  (void)ierr;
253 # else
254  // In versions older than 11.11 the Import function is broken for
255  // adding Hence, we provide a workaround in this case
256 
257  Epetra_MultiVector dummy(vector->Map(), 1, false);
258  Epetra_Import data_exchange(dummy.Map(),
259  down_V.trilinos_vector().Map());
260 
261  int ierr =
262  dummy.Import(down_V.trilinos_vector(), data_exchange, Insert);
263  Assert(ierr == 0, ExcTrilinosError(ierr));
264 
265  ierr = vector->Update(1.0, dummy, 1.0);
266  Assert(ierr == 0, ExcTrilinosError(ierr));
267  (void)ierr;
268 # endif
269  }
270 
271  return *this;
272  }
273 
274 
275 
276  Vector &
278  {
279  this->add(-1., V);
280 
281  return *this;
282  }
283 
284 
285 
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 
543  MPI_Comm
545  {
546  const Epetra_MpiComm *epetra_comm =
547  dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
548  Assert(epetra_comm != nullptr, ExcInternalError());
549  return epetra_comm->GetMpiComm();
550  }
551 
552 
553 
554  ::IndexSet
556  {
557  IndexSet is(size());
558 
559  // easy case: local range is contiguous
560  if (vector->Map().LinearMap())
561  {
562 # ifndef DEAL_II_WITH_64BIT_INDICES
563  is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID() + 1);
564 # else
565  is.add_range(vector->Map().MinMyGID64(),
566  vector->Map().MaxMyGID64() + 1);
567 # endif
568  }
569  else if (vector->Map().NumMyElements() > 0)
570  {
571  const size_type n_indices = vector->Map().NumMyElements();
572 # ifndef DEAL_II_WITH_64BIT_INDICES
573  unsigned int *vector_indices =
574  reinterpret_cast<unsigned int *>(vector->Map().MyGlobalElements());
575 # else
576  size_type *vector_indices =
577  reinterpret_cast<size_type *>(vector->Map().MyGlobalElements64());
578 # endif
579  is.add_indices(vector_indices, vector_indices + n_indices);
580  }
581  is.compress();
582 
583  return is;
584  }
585 
586 
587 
588  const Epetra_FEVector &
590  {
591  return *vector;
592  }
593 
594 
595 
596  Epetra_FEVector &
598  {
599  return *vector;
600  }
601 
602 
603 
604  void
605  Vector::print(std::ostream & out,
606  const unsigned int precision,
607  const bool scientific,
608  const bool across) const
609  {
610  AssertThrow(out, ExcIO());
611  boost::io::ios_flags_saver restore_flags(out);
612 
613  // Get a representation of the vector and loop over all
614  // the elements
615  double *val;
616  int leading_dimension;
617  int ierr = vector->ExtractView(&val, &leading_dimension);
618 
619  Assert(ierr == 0, ExcTrilinosError(ierr));
620  (void)ierr;
621  out.precision(precision);
622  if (scientific)
623  out.setf(std::ios::scientific, std::ios::floatfield);
624  else
625  out.setf(std::ios::fixed, std::ios::floatfield);
626 
627  if (across)
628  for (int i = 0; i < vector->MyLength(); ++i)
629  out << val[i] << ' ';
630  else
631  for (int i = 0; i < vector->MyLength(); ++i)
632  out << val[i] << std::endl;
633  out << std::endl;
634 
635  // restore the representation
636  // of the vector
637  AssertThrow(out, ExcIO());
638  }
639 
640 
641 
642  std::size_t
644  {
645  return sizeof(*this) +
646  vector->MyLength() *
647  (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
648  }
649 
650 
651 
652  void
654  const MPI_Comm &mpi_comm)
655  {
656  source_stored_elements = source_index_set;
658  std::make_shared<CommunicationPattern>(locally_owned_elements(),
659  source_index_set,
660  mpi_comm);
661  }
662  } // namespace EpetraWrappers
663 } // namespace LinearAlgebra
664 
666 
667 # endif
668 
669 #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:1701
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:1576
virtual std::size_t memory_consumption() const override
std::unique_ptr< Epetra_FEVector > vector
const Epetra_Comm & comm_self()
Definition: utilities.cc:1114
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
size_type size() const
Definition: index_set.h:1632
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:1466
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:372
virtual void import(const ReadWriteVector< double > &V, VectorOperation::values operation, std::shared_ptr< const CommunicationPatternBase > communication_pattern=std::shared_ptr< const CommunicationPatternBase >()) override
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:1671
virtual Vector & operator+=(const VectorSpaceVector< double > &V) override
Definition: cuda.h:32
void compress() const
Definition: index_set.h:1640
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:371
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:1722
virtual double linfty_norm() const override
const IndexSet & get_stored_elements() const
static ::ExceptionBase & ExcInternalError()
virtual double add_and_dot(const double a, const VectorSpaceVector< double > &V, const VectorSpaceVector< double > &W) override