class.
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__communication_pattern_base_h
+#define dealii__communication_pattern_base_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/mpi.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+class IndexSet;
+
+namespace LinearAlgebra
+{
+ /**
+ * CommunicationPattern is an abstract class that is used to define a
+ * communication plan that can be called repeatedly to efficiently obtain
+ * off-processor elements.
+ *
+ * @author Bruno Turcksin, 2015.
+ */
+ class CommunicationPatternBase
+ {
+ public:
+ /**
+ * Reinitialize the communication pattern. The first argument @p
+ * vector_space_vector_index_set is the index set associated to a
+ * VectorSpaceVector object. The second argument @p
+ * read_write_vector_index_set is the index set associated to a
+ * ReadWriteVector object.
+ */
+ virtual void reinit(const IndexSet &vector_space_vector_index_set,
+ const IndexSet &read_write_vector_index_set,
+ const MPI_Comm &communicator) = 0;
+
+ /**
+ * Return a constant reference to the underlying mpi communicator.
+ */
+ virtual const MPI_Comm &get_mpi_communicator() const = 0;
+ };
+
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
/**
* Multiply the entire vector by a fixed factor.
*/
- virtual VectorSpaceVector<Number> &operator*= (const Number factor);
+ virtual Vector<Number> &operator*= (const Number factor) override;
/**
* Divide the entire vector by a fixed factor.
*/
- virtual VectorSpaceVector<Number> &operator/= (const Number factor);
+ virtual Vector<Number> &operator/= (const Number factor) override;
/**
* Add the vector @p V to the present one.
*/
- virtual VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+ virtual Vector<Number> &operator+= (const VectorSpaceVector<Number> &V) override;
/**
* Substract the vector @p V from the present one.
*/
- virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+ virtual Vector<Number> &operator-= (const VectorSpaceVector<Number> &V) override;
/**
* Return the scalar product of two vectors.
*/
- virtual Number operator* (const VectorSpaceVector<Number> &V);
+ virtual Number operator* (const VectorSpaceVector<Number> &V) const override;
/**
* Add @p a to all components. Note that @p a is a scalar not a vector.
*/
- virtual void add(const Number a);
+ virtual void add(const Number a) override;
/**
* Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
*/
- virtual void add(const Number a, const VectorSpaceVector<Number> &V);
+ virtual void add(const Number a, const VectorSpaceVector<Number> &V) override;
/**
* Multiple addition of a multiple of a vector, i.e. <tt>*this +=
* a*V+b*W</tt>.
*/
virtual void add(const Number a, const VectorSpaceVector<Number> &V,
- const Number b, const VectorSpaceVector<Number> &W);
+ const Number b, const VectorSpaceVector<Number> &W) override;
/**
* Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
* s*(*this)+a*V</tt>.
*/
virtual void sadd(const Number s, const Number a,
- const VectorSpaceVector<Number> &V);
+ const VectorSpaceVector<Number> &V) override;
/**
* Scale each element of this vector by the corresponding element in the
* argument. This function is mostly meant to simulate multiplication (and
* immediate re-assignement) by a diagonal scaling matrix.
*/
- virtual void scale(const VectorSpaceVector<Number> &scaling_factors);
+ virtual void scale(const VectorSpaceVector<Number> &scaling_factors) override;
/**
* Assignement <tt>*this = a*V</tt>.
*/
- virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
+ virtual void equ(const Number a, const VectorSpaceVector<Number> &V) override;
/**
* Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
* absolute values of all entries).
*/
- virtual typename VectorSpaceVector<Number>::real_type l1_norm();
+ virtual typename VectorSpaceVector<Number>::real_type l1_norm() const override;
/**
* Return the l<sub>2</sub> norm of the vector (i.e., the square root of
* the sum of the square of all entries among all processors).
*/
- virtual typename VectorSpaceVector<Number>::real_type l2_norm();
+ virtual typename VectorSpaceVector<Number>::real_type l2_norm() const override;
/**
* Return the maximum norm of the vector (i.e., the maximum absolute value
* among all entries and among all processors).
*/
- virtual typename VectorSpaceVector<Number>::real_type linfty_norm();
+ virtual typename VectorSpaceVector<Number>::real_type linfty_norm() const override;
/**
* Perform a combined operation of a vector addition and a subsequent
*/
virtual Number add_and_dot(const Number a,
const VectorSpaceVector<Number> &V,
- const VectorSpaceVector<Number> &W);
+ const VectorSpaceVector<Number> &W) override;
/**
* Return the global size of the vector, equal to the sum of the number of
* locally owned indices among all processors.
*/
- size_type size() const;
+ virtual size_type size() const override;
/**
* Return an index set that describes which elements of this vector are
* vec.locally_owned_elements() == complete_index_set(vec.size())
* @endcode
*/
- const dealii::IndexSet locally_owned_elements() const;
+ virtual dealii::IndexSet locally_owned_elements() const override;
/**
* Prints the vector to the output stream @p out.
virtual void print(std::ostream &out,
const unsigned int precision=3,
const bool scientific=true,
- const bool across=true) const;
+ const bool across=true) const override;
/**
* Write the vector en bloc to a file. This is done in a binary mode, so
* the output is neither readable by humans nor (probably) by other
/**
* Returns the memory consumption of this class in bytes.
*/
- std::size_t memory_consumption() const;
+ virtual std::size_t memory_consumption() const override;
/**
* Attempt to perform an operation between two incompatible vector types.
* large vector.
*/
typename VectorSpaceVector<Number>::real_type l1_norm_recursive(unsigned int i,
- unsigned int j);
+ unsigned int j) const;
/**
* Compute the squared L2 norm in a recursive way by dividing the vector
*/
typename VectorSpaceVector<Number>::real_type l2_norm_squared_recursive(
unsigned int i,
- unsigned int j);
+ unsigned int j) const;
/**
* Serialize the data of this object using boost. This function is
template <typename Number>
inline
- const dealii::IndexSet Vector<Number>::locally_owned_elements() const
+ dealii::IndexSet Vector<Number>::locally_owned_elements() const
{
- return ReadWriteVector<Number>::get_stored_elements();
+ return IndexSet(ReadWriteVector<Number>::get_stored_elements());
}
namespace LinearAlgebra
{
template <typename Number>
- VectorSpaceVector<Number> &Vector<Number>::operator*= (const Number factor)
+ Vector<Number> &Vector<Number>::operator*= (const Number factor)
{
AssertIsFinite(factor);
for (unsigned int i=0; i<this->size(); ++i)
template <typename Number>
- VectorSpaceVector<Number> &Vector<Number>::operator/= (const Number factor)
+ Vector<Number> &Vector<Number>::operator/= (const Number factor)
{
AssertIsFinite(factor);
Assert(factor!=Number(0.), ExcZero());
template <typename Number>
- VectorSpaceVector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
+ Vector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
{
// Check that casting will work.
Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
template <typename Number>
- VectorSpaceVector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
+ Vector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
{
// Check that casting will work.
Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
template <typename Number>
- Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V)
+ Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V) const
{
// Check that casting will work.
Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
template <typename Number>
- typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm()
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm() const
{
Assert (this->size(), ExcEmptyObject());
template <typename Number>
- typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm()
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm() const
{
Assert (this->size(), ExcEmptyObject());
template <typename Number>
- typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm()
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm() const
{
typename ReadWriteVector<Number>::real_type norm = 0.;
for (unsigned int i=0; i<this->size(); ++i)
template <typename Number>
typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(
unsigned int i,
- unsigned int j)
+ unsigned int j) const
{
Assert(j>=i, ExcInternalError());
typename ReadWriteVector<Number>::real_type norm = 0.;
template <typename Number>
typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_squared_recursive(
unsigned int i,
- unsigned int j)
+ unsigned int j) const
{
Assert(j>=i, ExcInternalError());
DEAL_II_NAMESPACE_OPEN
+namespace LinearAlgebra
+{
+ class CommunicationPatternBase;
+}
+
#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
#ifdef DEAL_II_WITH_PETSC
/**
- * Imports all the elements present in the vector's IndexSet from the
- * input vector @p petsc_vec.
- */
- ReadWriteVector<Number> &
- operator= (const PETScWrappers::MPI::Vector &petsc_vec);
+ * Imports all the elements present in the vector's IndexSet from the input
+ * vector @p petsc_vec. VectorOperation::values @p operation is used to decide
+ * if the elements in @p V should be added to the current vector or replace
+ * the current elements. The last parameter can be used if the same
+ * communication pattern is used multiple times. This can be used to improve
+ * performance.
+ */
+ void import(const PETScWrappers::MPI::Vector &petsc_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL);
#endif
#ifdef DEAL_II_WITH_TRILINOS
- /**
- * Imports all the elements present in the vector's IndexSet from the
- * input vector @p trilinos_vec.
- */
- ReadWriteVector<Number> &
- operator= (const TrilinosWrappers::MPI::Vector &trilinos_vec);
-
/**
* Imports all the elements present in the vector's IndexSet from the input
- * vector @p epetra_vec.
- */
- ReadWriteVector<Number> &
- operator= (const EpetraWrappers::Vector &epetra_vec);
+ * vector @p trilinos_vec. VectorOperation::values @p operation is used to
+ * decide if the elements in @p V should be added to the current vector or
+ * replace the current elements. The last parameter can be used if the same
+ * communication pattern is used multiple times. This can be used to improve
+ * performance.
+ */
+ void import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL);
+
+ /**
+ * imports all the elements present in the vector's IndexSet from the input
+ * vector @p epetra_vec. VectorOperation::values @p operation is used to
+ * decide if the elements in @p V should be added to the current vector or
+ * replace the current elements. The last parameter can be used if the same
+ * communication pattern is used multiple times. This can be used to improve
+ * performance.
+ */
+ void import(const EpetraWrappers::Vector &epetra_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL);
#endif
/**
#ifdef DEAL_II_WITH_TRILINOS
/**
* Imports all the elements present in the vector's IndexSet from the input
- * vector @p multivector.
- */
- ReadWriteVector<Number> &import(const Epetra_MultiVector &multivector,
- const IndexSet &locally_owned_elements);
+ * vector @p multivector. This is an helper function and it should not be
+ * used directly.
+ */
+ void import(const Epetra_MultiVector &multivector,
+ const IndexSet &locally_owned_elements,
+ VectorOperation::values operation,
+ const MPI_Comm &mpi_comm,
+ const CommunicationPatternBase *communication_pattern);
#endif
/**
*/
void resize_val (const size_type new_allocated_size);
+ /**
+ *
+ */
+ void create_epetra_comm_pattern(const IndexSet &source_index_set,
+ const MPI_Comm &mpi_comm);
+
/**
* Indices of the elements stored.
*/
IndexSet stored_elements;
+ /**
+ * Indices of the elements stored on a VectorSpaceVector
+ */
+ std::shared_ptr<IndexSet> source_stored_elements;
+
+ /**
+ *
+ */
+ std::shared_ptr<CommunicationPatternBase> comm_pattern;
+
/**
* Pointer to the array of local elements of this vector.
*/
inline
ReadWriteVector<Number>::ReadWriteVector ()
:
- val (NULL)
+ source_stored_elements(nullptr),
+ comm_pattern(nullptr),
+ val(nullptr)
{}
ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
:
Subscriptor(),
- val (NULL)
+ source_stored_elements(nullptr),
+ comm_pattern(nullptr),
+ val(nullptr)
{
this->operator=(v);
}
inline
ReadWriteVector<Number>::ReadWriteVector (const size_type size)
:
- val (NULL)
+ source_stored_elements(nullptr),
+ comm_pattern(nullptr),
+ val(nullptr)
{
reinit (size, false);
}
inline
ReadWriteVector<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
:
- val (NULL)
+ source_stored_elements(nullptr),
+ comm_pattern(nullptr),
+ val(nullptr)
{
reinit (locally_stored_indices);
}
// ---------------------------------------------------------------------
//
-// Copyright (C) 2015, 2016 by the deal.II authors
+// Copyright (C) 2015-2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
#include <deal.II/lac/trilinos_epetra_vector.h>
#ifdef DEAL_II_WITH_TRILINOS
-# include <Epetra_Import.h>
+# include <deal.II/lac/trilinos_epetra_communication_pattern.h>
+# include "Epetra_Import.h"
#endif
DEAL_II_NAMESPACE_OPEN
{
if (new_alloc_size == 0)
{
- if (val != NULL)
+ if (val != nullptr)
free(val);
- val = NULL;
+ val = nullptr;
}
else
{
- if (val != NULL)
+ if (val != nullptr)
free(val);
Utilities::System::posix_memalign ((void **)&val, 64, sizeof(Number)*new_alloc_size);
// set entries to zero if so requested
if (omit_zeroing_entries == false)
this->operator = (Number());
+
+ // reset the communication patter
+ source_stored_elements = nullptr;
+ comm_pattern = nullptr;
}
if (omit_zeroing_entries == false)
this->operator= (Number());
+
+ // reset the communication patter
+ source_stored_elements = nullptr;
+ comm_pattern = nullptr;
}
// initialize to zero
if (omit_zeroing_entries == false)
this->operator= (Number());
+
+ // reset the communication patter
+ source_stored_elements = nullptr;
+ comm_pattern = nullptr;
}
template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const PETScWrappers::MPI::Vector &petsc_vec)
+ void
+ ReadWriteVector<Number>::import(const PETScWrappers::MPI::Vector &petsc_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern)
{
//TODO: this works only if no communication is needed.
Assert(petsc_vec.locally_owned_elements() == stored_elements,
// restore the representation of the vector
ierr = VecRestoreArray (static_cast<const Vec &>(petsc_vec), &start_ptr);
AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- // return a pointer to this object per normal c++ operator overloading
- // semantics
- return *this;
}
#endif
#ifdef DEAL_II_WITH_TRILINOS
template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::import(const Epetra_MultiVector &multivector,
- const IndexSet &locally_owned_elements)
+ void
+ ReadWriteVector<Number>::import(const Epetra_MultiVector &multivector,
+ const IndexSet &source_elements,
+ VectorOperation::values operation,
+ const MPI_Comm &mpi_comm,
+ const CommunicationPatternBase *communication_pattern)
{
- // Copy the local elements
- std::vector<size_type> trilinos_indices, readwrite_indices;
- locally_owned_elements.fill_index_vector(trilinos_indices);
- stored_elements.fill_index_vector(readwrite_indices);
- std::vector<size_type> intersection(trilinos_indices.size());
- std::vector<size_type>::iterator end;
- end = std::set_intersection(trilinos_indices.begin(), trilinos_indices.end(),
- readwrite_indices.begin(), readwrite_indices.end(),
- intersection.begin());
- const size_t intersection_size = end-intersection.begin();
- intersection.resize(intersection_size);
- std::vector<size_t> trilinos_position(intersection_size);
- unsigned int j=0;
- for (size_t i=0; i<intersection_size; ++i)
+ const EpetraWrappers::CommunicationPattern *epetra_comm_pattern = nullptr;
+
+ if (communication_pattern == nullptr)
{
- while (intersection[i]!=trilinos_indices[j])
- ++j;
- trilinos_position[i] = j;
+ if (source_elements == *source_stored_elements)
+ {
+ epetra_comm_pattern =
+ dynamic_cast<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get());
+ if (epetra_comm_pattern == nullptr)
+ create_epetra_comm_pattern(source_elements, mpi_comm);
+ }
+ else
+ create_epetra_comm_pattern(source_elements, mpi_comm);
+
+ epetra_comm_pattern =
+ dynamic_cast<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get());
+ }
+ else
+ {
+ epetra_comm_pattern =
+ dynamic_cast<const EpetraWrappers::CommunicationPattern *> (communication_pattern);
+ AssertThrow(epetra_comm_pattern != nullptr,
+ ExcMessage(std::string("The communication pattern is not of type ") +
+ "LinearAlgebra::EpetraWrappers::CommunicationPattern."));
}
- double *values = multivector.Values();
- for (size_t i=0; i<trilinos_position.size(); ++i)
- val[global_to_local(intersection[i])] = values[trilinos_position[i]];
+ Epetra_Import import(epetra_comm_pattern->get_epetra_import());
-#ifdef DEAL_II_WITH_MPI
- // Copy the off-processor elements if necessary
- if (intersection_size != n_elements())
+ Epetra_FEVector target_vector(import.TargetMap());
+
+ if (operation==VectorOperation::insert)
{
- Epetra_BlockMap source_map = multivector.Map();
- std::vector<size_type> off_proc_indices(readwrite_indices.size());
- // Subtract the elements of readwrite that are in intersection.
- end = std::set_difference(readwrite_indices.begin(), readwrite_indices.end(),
- intersection.begin(), intersection.end(), off_proc_indices.begin());
- off_proc_indices.resize(end-off_proc_indices.begin());
- // Cast size_type to TrilinosWrappers::type::int_type
- TrilinosWrappers::types::int_type *off_proc_array =
- new TrilinosWrappers::types::int_type [off_proc_indices.size()];
- for (size_t i=0; i<off_proc_indices.size(); ++i)
- off_proc_array[i] = static_cast<TrilinosWrappers::types::int_type> (
- off_proc_indices[i]);
- Epetra_Map target_map(off_proc_indices.size(),off_proc_indices.size(),
- off_proc_array,0,source_map.Comm());
- delete [] off_proc_array;
- Epetra_Import import(target_map, source_map);
- Epetra_FEVector target_vector(target_map);
target_vector.Import(multivector, import, Insert);
- values = target_vector.Values();
- for (unsigned int i=0; i<off_proc_indices.size(); ++i)
- val[global_to_local(off_proc_indices[i])] = values[i];
+ double *values = target_vector.Values();
+ for (int i=0; i<target_vector.MyLength(); ++i)
+ val[i] = values[i];
+ }
+ else
+ {
+ target_vector.Import(multivector, import, Add);
+ double *values = target_vector.Values();
+ for (int i=0; i<target_vector.MyLength(); ++i)
+ val[i] += values[i];
}
-#endif
-
- // return a pointer to this object per normal c++ operator overloading
- // semantics
- return *this;
}
-#endif
template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+ void
+ ReadWriteVector<Number>::import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern)
{
- return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements());
+ import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements(),
+ operation, trilinos_vec.get_mpi_communicator(), communication_pattern);
}
+
template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec)
+ void
+ ReadWriteVector<Number>::import(const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern)
{
- return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements());
+ import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements(),
+ operation, trilinos_vec.get_mpi_communicator(), communication_pattern);
}
+#endif
+
template <typename Number>
out.flags (old_flags);
out.precision(old_precision);
}
+
+
+
+ template <typename Number>
+ void
+ ReadWriteVector<Number>::create_epetra_comm_pattern(const IndexSet &source_index_set,
+ const MPI_Comm &mpi_comm)
+ {
+ source_stored_elements = std::make_shared<IndexSet>(source_index_set);
+ comm_pattern.reset(new EpetraWrappers::CommunicationPattern(
+ *source_stored_elements, stored_elements, mpi_comm));
+ }
} // end of namespace LinearAlgebra
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015-2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__trilinos_epetra_communication_pattern_h
+#define dealii__trilinos_epetra_communication_pattern_h
+
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/lac/communication_pattern_base.h>
+#include "Epetra_Import.h"
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ namespace EpetraWrappers
+ {
+ /**
+ * This class implements a wrapper to Trilinos Import.
+ */
+ class CommunicationPattern : public CommunicationPatternBase
+ {
+ public:
+ /**
+ * Reinitialize the communication pattern. The first argument @p
+ * vector_space_vector_index_set is the index set associated to a
+ * VectorSpaceVector object. The second argument @p
+ * read_write_vector_index_set is the index set associated to a
+ * ReadWriteVector object.
+ */
+ CommunicationPattern(const IndexSet &vector_space_vector_index_set,
+ const IndexSet &read_write_vector_index_set,
+ const MPI_Comm &communicator);
+
+ /**
+ * Reinitialize the object.
+ */
+ void reinit(const IndexSet &vector_space_vector_index_set,
+ const IndexSet &read_write_vector_index_set,
+ const MPI_Comm &communicator) override;
+
+ /**
+ * Return the underlying MPI communicator.
+ */
+ const MPI_Comm &get_mpi_communicator() const override;
+
+ /**
+ * Return the underlying Epetra_Import object.
+ */
+ const Epetra_Import &get_epetra_import() const;
+
+ private:
+ /**
+ * Shared pointer to the MPI communicator used.
+ */
+ std_cxx11::shared_ptr<const MPI_Comm> comm;
+
+ /**
+ * Shared pointer to the Epetra_Import object used.
+ */
+ std_cxx11::shared_ptr<Epetra_Import> import;
+ };
+ } // end of namespace EpetraWrappers
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
+
+#endif
namespace EpetraWrappers
{
/**
- * This class implements a wrapper to the the Trilinos distributed vector
+ * This class implements a wrapper to the Trilinos distributed vector
* class Epetra_FEVector. This class is derived from the
* LinearAlgebra::VectorSpaceVector class. Note however that Epetra only
* works with Number = double.
/**
* This constructor takes an IndexSet that defines how to distribute the
* individual components among the MPI processors. Since it also
- * includes information abtou the size of the vector, this is all we
+ * includes information about the size of the vector, this is all we
* need to genereate a %parallel vector.
*/
- Vector(const IndexSet ¶llel_partitioner,
- const MPI_Comm &communicator = MPI_COMM_WORLD);
+ explicit Vector(const IndexSet ¶llel_partitioner,
+ const MPI_Comm &communicator);
/**
* Reinit functionality. This function destroys the old vector content
* filled with zero (false) or left untouched (true).
*/
void reinit (const IndexSet ¶llel_partitioner,
- const MPI_Comm &communicator = MPI_COMM_WORLD,
+ const MPI_Comm &communicator,
const bool omit_zeroing_entries = false);
/**
/**
* Multiply the entire vector by a fixed factor.
*/
- virtual VectorSpaceVector<double> &operator*= (const double factor);
+ virtual Vector &operator*= (const double factor);
/**
* Divide the entire vector by a fixed factor.
*/
- virtual VectorSpaceVector<double> &operator/= (const double factor);
+ virtual Vector &operator/= (const double factor);
/**
* Add the vector @p V to the present one.
*/
- virtual VectorSpaceVector<double> &operator+= (const VectorSpaceVector<double> &V);
+ virtual Vector &operator+= (const VectorSpaceVector<double> &V);
/**
* Substract the vector @p V from the present one.
*/
- virtual VectorSpaceVector<double> &operator-= (const VectorSpaceVector<double> &V);
+ virtual Vector &operator-= (const VectorSpaceVector<double> &V);
/**
* Return the scalar product of two vectors. The vectors need to have the
* same layout.
*/
- virtual double operator* (const VectorSpaceVector<double> &V);
+ virtual double operator* (const VectorSpaceVector<double> &V) const;
/**
* Add @p a to all components. Note that @p is a scalar not a vector.
* Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
* absolute values of all entries among all processors).
*/
- virtual double l1_norm();
+ virtual double l1_norm() const;
/**
* Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
* the sum of the square of all entries among all processors).
*/
- virtual double l2_norm();
+ virtual double l2_norm() const;
/**
* Returns the maximum norm of the vector (i.e., the maximum absolute value
* among all entries and among all processors).
*/
- virtual double linfty_norm();
+ virtual double linfty_norm() const;
/**
* Performs a combined operation of a vector addition and a subsequent
* vec.locally_owned_elements() == complete_index_set(vec.size())
* @endcode
*/
- virtual const ::dealii::IndexSet locally_owned_elements() const;
+ virtual ::dealii::IndexSet locally_owned_elements() const;
/**
* Return a const reference to the underlying Trilinos
*/
virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V) = 0;
+ /**
+ * Imports all the elements present in the vector's IndexSet from the input
+ * vector @p V. VectorOperation::values @p operation is used to decide if
+ * the elements in @p V should be added to the current vector or replace the
+ * current elements. The last parameter can be used if the same
+ * communication pattern is used multiple times. This can be used to improve
+ * performance.
+ */
+ virtual void Import(const ReadWriteVector<Number> &V,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL);
+
/**
* Return the scalar product of two vectors.
*/
- virtual Number operator* (const VectorSpaceVector<Number> &V) = 0;
+ virtual Number operator* (const VectorSpaceVector<Number> &V) const = 0;
/**
* Add @p a to all components. Note that @p a is a scalar not a vector.
* Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
* absolute values of all entries among all processors).
*/
- virtual real_type l1_norm() = 0;
+ virtual real_type l1_norm() const = 0;
/**
* Return the l<sub>2</sub> norm of the vector (i.e., the square root of
* the sum of the square of all entries among all processors).
*/
- virtual real_type l2_norm() = 0;
+ virtual real_type l2_norm() const = 0;
/**
* Return the maximum norm of the vector (i.e., the maximum absolute value
* among all entries and among all processors).
*/
- virtual real_type linfty_norm() = 0;
+ virtual real_type linfty_norm() const = 0;
/**
* Perform a combined operation of a vector addition and a subsequent
* vec.locally_owned_elements() == complete_index_set(vec.size())
* @endcode
*/
- virtual const dealii::IndexSet locally_owned_elements() const = 0;
+ virtual dealii::IndexSet locally_owned_elements() const = 0;
/**
* Print the vector to the output stream @p out.
${_src}
trilinos_block_sparse_matrix.cc
trilinos_block_vector.cc
+ trilinos_epetra_communication_pattern.cc
trilinos_epetra_vector.cc
trilinos_parallel_block_vector.cc
trilinos_precondition.cc
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/trilinos_epetra_communication_pattern.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/index_set.h>
+#include "Epetra_Map.h"
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ namespace EpetraWrappers
+ {
+ CommunicationPattern::CommunicationPattern(const IndexSet &vector_space_vector_index_set,
+ const IndexSet &read_write_vector_index_set,
+ const MPI_Comm &communicator)
+ {
+ reinit(vector_space_vector_index_set, read_write_vector_index_set, communicator);
+ }
+
+
+
+ void CommunicationPattern::reinit(const IndexSet &vector_space_vector_index_set,
+ const IndexSet &read_write_vector_index_set,
+ const MPI_Comm &communicator)
+ {
+ comm = std::make_shared<const MPI_Comm>(communicator);
+
+ Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm);
+ Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm);
+
+ // Target map is read_write_vector_map
+ // Source map is vector_space_vector_map. This map must have uniquely
+ // owned GID.
+ import.reset(new Epetra_Import(read_write_vector_map, vector_space_vector_map));
+ }
+
+
+
+ const MPI_Comm &CommunicationPattern::get_mpi_communicator() const
+ {
+ return *comm;
+ }
+
+
+
+ const Epetra_Import &CommunicationPattern::get_epetra_import() const
+ {
+ return *import;
+ }
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
#include <deal.II/lac/trilinos_epetra_vector.h>
#ifdef DEAL_II_WITH_TRILINOS
+
#ifdef DEAL_II_WITH_MPI
#include <deal.II/base/index_set.h>
- VectorSpaceVector<double> &Vector::operator*= (const double factor)
+ Vector &Vector::operator*= (const double factor)
{
AssertIsFinite(factor);
vector->Scale(factor);
- VectorSpaceVector<double> &Vector::operator/= (const double factor)
+ Vector &Vector::operator/= (const double factor)
{
AssertIsFinite(factor);
Assert(factor!=0., ExcZero());
- VectorSpaceVector<double> &Vector::operator+= (const VectorSpaceVector<double> &V)
+ Vector &Vector::operator+= (const VectorSpaceVector<double> &V)
{
// Check that casting will work.
Assert(dynamic_cast<const Vector *>(&V)!=NULL, ExcVectorTypeNotCompatible());
- VectorSpaceVector<double> &Vector::operator-= (const VectorSpaceVector<double> &V)
+ Vector &Vector::operator-= (const VectorSpaceVector<double> &V)
{
this->add(-1.,V);
- double Vector::operator* (const VectorSpaceVector<double> &V)
+ double Vector::operator* (const VectorSpaceVector<double> &V) const
{
// Check that casting will work.
Assert(dynamic_cast<const Vector *>(&V)!=NULL,
- double Vector::l1_norm()
+ double Vector::l1_norm() const
{
double norm(0.);
int ierr = vector->Norm1(&norm);
- double Vector::l2_norm()
+ double Vector::l2_norm() const
{
double norm(0.);
int ierr = vector->Norm2(&norm);
- double Vector::linfty_norm()
+ double Vector::linfty_norm() const
{
double norm(0.);
int ierr = vector->NormInf(&norm);
- const ::dealii::IndexSet Vector::locally_owned_elements() const
+ ::dealii::IndexSet Vector::locally_owned_elements() const
{
const Epetra_Map *map = dynamic_cast<const Epetra_Map *>(&(vector->Map()));
return IndexSet(*map);
out.setf(std::ios::fixed, std::ios::floatfield);
if (across)
- for (size_type i=0; i<size(); ++i)
+ for (size_type i=0; i<vector->MyLength(); ++i)
out << val[i] << ' ';
else
- for (size_type i=0; i<size(); ++i)
+ for (size_type i=0; i<vector->MyLength(); ++i)
out << val[i] << std::endl;
out << std::endl;
parallel_partitioner_1.compress();
parallel_partitioner_2.compress();
LinearAlgebra::EpetraWrappers::Vector a;
- LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1);
+ LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1, MPI_COMM_WORLD);
LinearAlgebra::EpetraWrappers::Vector c(b);
AssertThrow(a.size()==0, ExcMessage("Vector has the wrong size."));
AssertThrow(b.size()==10, ExcMessage("Vector has the wrong size."));
AssertThrow(c.size()==10, ExcMessage("Vector has the wrong size."));
- a.reinit(parallel_partitioner_2);
+ a.reinit(parallel_partitioner_2, MPI_COMM_WORLD);
AssertThrow(a.size()==10, ExcMessage("Vector has the wrong size."));
AssertThrow(parallel_partitioner_1==b.locally_owned_elements(),
}
parallel_partitioner_1.compress();
parallel_partitioner_2.compress();
- LinearAlgebra::EpetraWrappers::Vector a(parallel_partitioner_1);
- LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1);
- LinearAlgebra::EpetraWrappers::Vector c(parallel_partitioner_2);
+ LinearAlgebra::EpetraWrappers::Vector a(parallel_partitioner_1, MPI_COMM_WORLD);
+ LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1, MPI_COMM_WORLD);
+ LinearAlgebra::EpetraWrappers::Vector c(parallel_partitioner_2, MPI_COMM_WORLD);
IndexSet read_write_index_set(10);
if (rank==0)
projects / dealii.git / commitdiff