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