Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
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
37namespace 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
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.
161 V.get_stored_elements().size()) ||
163 {
166 dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
167 }
168 }
169 else
170 {
172 std::dynamic_pointer_cast<const CommunicationPattern>(
173 communication_pattern);
175 epetra_comm_pattern != nullptr,
177 std::string("The communication pattern is not of type ") +
178 "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
179 }
180
181 Epetra_Import import_map(epetra_comm_pattern->get_epetra_import());
182
183 // The TargetMap and the SourceMap have their roles inverted.
184 Epetra_FEVector source_vector(import_map.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_map, Insert);
190 else if (operation == VectorOperation::add)
191 vector->Export(source_vector, import_map, Add);
192 else if (operation == VectorOperation::max)
193 vector->Export(source_vector, import_map, Epetra_Max);
194 else if (operation == VectorOperation::min)
195 vector->Export(source_vector, import_map, Epetra_Min);
196 else
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 {
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);
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,
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()),
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
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,
525 {
526 this->add(a, V);
527
528 return *this * W;
529 }
530
531
532
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
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
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:1661
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:789
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1697
void compress() const
Definition index_set.h:1669
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
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
virtual double operator*(const VectorSpaceVector< double > &V) const override
virtual void import_elements(const ReadWriteVector< double > &V, VectorOperation::values operation, std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > communication_pattern={}) override
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
void create_epetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
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
void reinit(const IndexSet &parallel_partitioner, const MPI_Comm communicator, const bool omit_zeroing_entries=false)
virtual void sadd(const double s, const double a, const VectorSpaceVector< double > &V) override
const IndexSet & get_stored_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcVectorTypeNotCompatible()
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
static ::ExceptionBase & ExcDifferentParallelPartitioning()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
T sum(const T &t, const MPI_Comm mpi_communicator)