#include <deal.II/base/config.h>
#ifdef DEAL_II_WITH_TRILINOS
-
# include <deal.II/base/subscriptor.h>
# include <deal.II/base/index_set.h>
# include <deal.II/base/utilities.h>
+# include <deal.II/base/mpi.h>
# include <deal.II/lac/exceptions.h>
# include <deal.II/lac/vector.h>
-# include <deal.II/lac/trilinos_vector_base.h>
# include <deal.II/lac/vector_type_traits.h>
+# include <vector>
+# include <utility>
+# include <memory>
+
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+# include "Epetra_ConfigDefs.h"
+# ifdef DEAL_II_WITH_MPI // only if MPI is installed
+# include "mpi.h"
+# include "Epetra_MpiComm.h"
+# else
+# include "Epetra_SerialComm.h"
+# endif
+# include "Epetra_FEVector.h"
# include "Epetra_Map.h"
# include "Epetra_LocalMap.h"
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
-#include <memory>
-
DEAL_II_NAMESPACE_OPEN
-// forward declaration
-template <typename> class Vector;
-
/**
* @addtogroup TrilinosWrappers
* @{
*/
+
+/**
+ * A namespace in which wrapper classes for Trilinos objects reside.
+ *
+ * @ingroup TrilinosWrappers
+ */
namespace TrilinosWrappers
{
class SparseMatrix;
+ /**
+ * @cond internal
+ */
+
+ /**
+ * A namespace for internal implementation details of the TrilinosWrapper
+ * members.
+ *
+ * @ingroup TrilinosWrappers
+ */
+ namespace internal
+ {
+ /**
+ * Declare type for container size.
+ */
+ typedef dealii::types::global_dof_index size_type;
+
+ /**
+ * This class implements a wrapper for accessing the Trilinos vector in
+ * the same way as we access deal.II objects: it is initialized with a
+ * vector and an element within it, and has a conversion operatorto
+ * extract the scalar value of this element. It also has a variety of
+ * assignment operatorfor writing to this one element.
+ *
+ * @ingroup TrilinosWrappers
+ */
+ class VectorReference
+ {
+ private:
+ /**
+ * Constructor. It is made private so as to only allow the actual vector
+ * class to create it.
+ */
+ VectorReference (MPI::Vector &vector,
+ const size_type index);
+
+ public:
+
+ /**
+ * This looks like a copy operator, but does something different than
+ * usual. In particular, it does not copy the member variables of this
+ * reference. Rather, it handles the situation where we have two vectors
+ * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
+ * both left and right hand side of the assignment have data type
+ * VectorReference, but what we really mean is to assign the vector
+ * elements represented by the two references. This operatorimplements
+ * this operation. Note also that this allows us to make the assignment
+ * operatorconst.
+ */
+ const VectorReference &
+ operator= (const VectorReference &r) const;
+
+ /**
+ * Same as above but for non-const reference objects.
+ */
+ VectorReference &
+ operator= (const VectorReference &r);
+
+ /**
+ * Set the referenced element of the vector to <tt>s</tt>.
+ */
+ const VectorReference &
+ operator= (const TrilinosScalar &s) const;
+
+ /**
+ * Add <tt>s</tt> to the referenced element of the vector->
+ */
+ const VectorReference &
+ operator+= (const TrilinosScalar &s) const;
+
+ /**
+ * Subtract <tt>s</tt> from the referenced element of the vector->
+ */
+ const VectorReference &
+ operator-= (const TrilinosScalar &s) const;
+
+ /**
+ * Multiply the referenced element of the vector by <tt>s</tt>.
+ */
+ const VectorReference &
+ operator*= (const TrilinosScalar &s) const;
+
+ /**
+ * Divide the referenced element of the vector by <tt>s</tt>.
+ */
+ const VectorReference &
+ operator/= (const TrilinosScalar &s) const;
+
+ /**
+ * Convert the reference to an actual value, i.e. return the value of
+ * the referenced element of the vector.
+ */
+ operator TrilinosScalar () const;
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occurred while calling a Trilinos function");
+
+ private:
+ /**
+ * Point to the vector we are referencing.
+ */
+ MPI::Vector &vector;
+
+ /**
+ * Index of the referenced element of the vector.
+ */
+ const size_type index;
+
+ /**
+ * Make the vector class a friend, so that it can create objects of the
+ * present type.
+ */
+ friend class ::dealii::TrilinosWrappers::MPI::Vector;
+ };
+ }
+ /**
+ * @endcond
+ */
+
namespace
{
#ifndef DEAL_II_WITH_64BIT_INDICES
/**
* Namespace for Trilinos vector classes that work in parallel over MPI.
- * This namespace is restricted to vectors only, whereas matrices are always
- * MPI based when run on more than one processor.
*
* @ingroup TrilinosWrappers
- * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ * @author Martin Kronbichler, Wolfgang Bangerth, Daniel Arndt, 2008, 2017
*/
namespace MPI
{
/**
* This class implements a wrapper to use the Trilinos distributed vector
- * class Epetra_FEVector. This class is derived from the
- * TrilinosWrappers::VectorBase class and provides all functionality
- * included there.
+ * class Epetra_FEVector, the (parallel) partitioning of which
+ * is governed by an Epetra_Map.
+ * The Epetra_FEVector is precisely the kind of vector
+ * we deal with all the time - we probably get it from some assembly
+ * process, where also entries not locally owned might need to written and
+ * hence need to be forwarded to the owner. The only requirement for this
+ * class to work is that Trilinos is installed with the same compiler as is
+ * used for compilation of deal.II.
+ *
+ * The interface of this class is modeled after the existing Vector class in
+ * deal.II. It has almost the same member functions, and is often
+ * exchangeable. However, since Trilinos only supports a single scalar type
+ * (double), it is not templated, and only works with that type.
*
* Note that Trilinos only guarantees that operations do what you expect
* if the function @p GlobalAssemble has been called after vector assembly
* code that may for example look like this:
*
* @code
- * TrilinosWrappers::VectorBase vector;
+ * TrilinosWrappers::MPI::Vector vector;
* // do some write operations on the vector
* for (size_type i=0; i<vector->size(); ++i)
* vector(i) = i;
*
* @ingroup TrilinosWrappers
* @ingroup Vectors
- * @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009
+ * @author Martin Kronbichler, Wolfgang Bangerth, Daniel Arndt,
+ * 2008, 2009, 2017
*/
- class Vector : public VectorBase
+ class Vector : public Subscriptor
{
public:
/**
- * Declare type for container size.
+ * Declare some of the standard types used in all containers. These types
+ * parallel those in the <tt>C</tt> standard libraries
+ * <tt>vector<...></tt> class.
*/
+ typedef TrilinosScalar value_type;
+ typedef TrilinosScalar real_type;
typedef dealii::types::global_dof_index size_type;
+ typedef value_type *iterator;
+ typedef const value_type *const_iterator;
+ typedef internal::VectorReference reference;
+ typedef const internal::VectorReference const_reference;
/**
- * @name Basic constructors and initialization.
+ * @name 1: Basic Object-handling
*/
//@{
/**
*/
Vector (const Vector &v);
- /**
- * Move constructor. Creates a new vector by stealing the internal data
- * of the vector @p v.
- */
- Vector (Vector &&v);
-
- /**
- * Destructor.
- */
- ~Vector ();
-
- /**
- * Reinit functionality. This function sets the calling vector to the
- * dimension and the parallel distribution of the input vector, but does
- * not copy the elements in <tt>v</tt>. If <tt>omit_zeroing_entries</tt>
- * is not <tt>true</tt>, the elements in the vector are initialized with
- * zero. If it is set to <tt>true</tt>, the vector entries are in an
- * unspecified state and the user has to set all elements. In the
- * current implementation, this method does not touch the vector entries
- * in case the vector layout is unchanged from before, otherwise entries
- * are set to zero. Note that this behavior might change between
- * releases without notification.
- *
- * This function has a third argument, <tt>allow_different_maps</tt>,
- * that allows for an exchange of data between two equal-sized vectors
- * (but being distributed differently among the processors). A trivial
- * application of this function is to generate a replication of a whole
- * vector on each machine, when the calling vector is built according to
- * the localized vector class TrilinosWrappers::VectorBase, and <tt>v</tt>
- * is a distributed vector. In this case, the variable
- * <tt>omit_zeroing_entries</tt> needs to be set to <tt>false</tt>,
- * since it does not make sense to exchange data between differently
- * parallelized vectors without touching the elements.
- */
- void reinit (const VectorBase &v,
- const bool omit_zeroing_entries = false,
- const bool allow_different_maps = false);
-
- /**
- * Create vector by merging components from a block vector.
- */
- void reinit (const BlockVector &v,
- const bool import_data = false);
-
- /**
- * Set all components of the vector to the given number @p s. Simply
- * pass this down to the base class, but we still need to declare this
- * function to make the example given in the discussion about making the
- * constructor explicit work.
- */
- Vector &operator= (const TrilinosScalar s);
-
- /**
- * Copy the given vector. Resize the present vector if necessary. In
- * this case, also the Epetra_Map that designs the parallel partitioning
- * is taken from the input vector.
- */
- Vector &operator= (const Vector &v);
-
- /**
- * Move the given vector. This operator replaces the present vector with
- * @p v by efficiently swapping the internal data structures.
- */
- Vector &operator= (Vector &&v);
-
- /**
- * Copy operator from a given localized vector (present on all
- * processes) in TrilinosWrappers format to the current distributed
- * vector. This function assumes that the calling vector (left hand
- * object) already is of the same size as the right hand side vector.
- * Otherwise, an exception will be thrown.
- */
- Vector &operator= (const ::dealii::TrilinosWrappers::VectorBase &v);
-
- /**
- * Another copy function. This one takes a deal.II vector and copies it
- * into a TrilinosWrapper vector. Note that since we do not provide any
- * Epetra_map that tells about the partitioning of the vector among the
- * MPI processes, the size of the TrilinosWrapper vector has to be the
- * same as the size of the input vector. In order to change the map, use
- * the reinit(const Epetra_Map &input_map) function.
- */
- template <typename Number>
- Vector &operator= (const ::dealii::Vector<Number> &v);
-
- /**
- * This reinit function is meant to be used for parallel calculations
- * where some non-local data has to be used. The typical situation where
- * one needs this function is the call of the
- * FEValues<dim>::get_function_values function (or of some derivatives)
- * in parallel. Since it is usually faster to retrieve the data in
- * advance, this function can be called before the assembly forks out to
- * the different processors. What this function does is the following:
- * It takes the information in the columns of the given matrix and looks
- * which data couples between the different processors. That data is
- * then queried from the input vector. Note that you should not write to
- * the resulting vector any more, since the some data can be stored
- * several times on different processors, leading to unpredictable
- * results. In particular, such a vector cannot be used for matrix-
- * vector products as for example done during the solution of linear
- * systems.
- */
- void import_nonlocal_data_for_fe
- (const dealii::TrilinosWrappers::SparseMatrix &matrix,
- const Vector &vector);
-//@}
- /**
- * @name Initialization with an IndexSet
- */
-//@{
/**
* This constructor takes an IndexSet that defines how to distribute the
* individual components among the MPI processors. Since it also
* @ref GlossGhostedVector "vectors with ghost elements"
*/
Vector (const IndexSet ¶llel_partitioning,
- const VectorBase &v,
+ const Vector &v,
const MPI_Comm &communicator = MPI_COMM_WORLD);
/**
const dealii::Vector<Number> &v,
const MPI_Comm &communicator = MPI_COMM_WORLD);
+ /**
+ * Move constructor. Creates a new vector by stealing the internal data
+ * of the vector @p v.
+ */
+ Vector (Vector &&v);
+
+ /**
+ * Destructor.
+ */
+ ~Vector ();
+
+ /**
+ * Release all memory and return to a state just like after having called
+ * the default constructor.
+ */
+ void clear ();
+
+ /**
+ * Reinit functionality. This function sets the calling vector to the
+ * dimension and the parallel distribution of the input vector, but does
+ * not copy the elements in <tt>v</tt>. If <tt>omit_zeroing_entries</tt>
+ * is not <tt>true</tt>, the elements in the vector are initialized with
+ * zero. If it is set to <tt>true</tt>, the vector entries are in an
+ * unspecified state and the user has to set all elements. In the
+ * current implementation, this method does not touch the vector entries
+ * in case the vector layout is unchanged from before, otherwise entries
+ * are set to zero. Note that this behavior might change between
+ * releases without notification.
+ *
+ * This function has a third argument, <tt>allow_different_maps</tt>,
+ * that allows for an exchange of data between two equal-sized vectors
+ * (but being distributed differently among the processors). A trivial
+ * application of this function is to generate a replication of a whole
+ * vector on each machine, when the calling vector is built with a map
+ * consisting of all indices on each process, and <tt>v</tt>
+ * is a distributed vector. In this case, the variable
+ * <tt>omit_zeroing_entries</tt> needs to be set to <tt>false</tt>,
+ * since it does not make sense to exchange data between differently
+ * parallelized vectors without touching the elements.
+ */
+ void reinit (const Vector &v,
+ const bool omit_zeroing_entries = false,
+ const bool allow_different_maps = false);
+
/**
* Reinit functionality. This function destroys the old vector content
* and generates a new one based on the input partitioning. The flag
* the entries to zero, but this might change between releases without
* notification.
*
- *
* Depending on whether the @p parallel_partitioning argument uniquely
* subdivides elements among processors or not, the resulting vector may
* or may not have ghost elements. See the general documentation of this
const IndexSet &ghost_entries,
const MPI_Comm &communicator = MPI_COMM_WORLD,
const bool vector_writable = false);
-//@}
- };
-
+ /**
+ * Create vector by merging components from a block vector.
+ */
+ void reinit (const BlockVector &v,
+ const bool import_data = false);
+ /**
+ * Compress the underlying representation of the Trilinos object, i.e.
+ * flush the buffers of the vector object if it has any. This function is
+ * necessary after writing into a vector element-by-element and before
+ * anything else can be done on it.
+ *
+ * The (defaulted) argument can be used to specify the compress mode
+ * (<code>Add</code> or <code>Insert</code>) in case the vector has not
+ * been written to since the last time this function was called. The
+ * argument is ignored if the vector has been added or written to since
+ * the last time compress() was called.
+ *
+ * See
+ * @ref GlossCompress "Compressing distributed objects"
+ * for more information.
+ */
+ void compress (::dealii::VectorOperation::values operation);
-// ------------------- inline and template functions --------------
+ /**
+ * Set all components of the vector to the given number @p s. Simply
+ * pass this down to the base class, but we still need to declare this
+ * function to make the example given in the discussion about making the
+ * constructor explicit work.
+ * the constructor explicit work.
+ *
+ * Since the semantics of assigning a scalar to a vector are not
+ * immediately clear, this operatorcan only be used if you want
+ * to set the entire vector to zero. This allows the intuitive notation
+ * <tt>v=0</tt>.
+ */
+ Vector &operator= (const TrilinosScalar s);
+ /**
+ * Copy the given vector. Resize the present vector if necessary. In
+ * this case, also the Epetra_Map that designs the parallel partitioning
+ * is taken from the input vector.
+ */
+ Vector &operator= (const Vector &v);
- /**
- * Global function @p swap which overloads the default implementation of
- * the C++ standard library which uses a temporary object. The function
- * simply exchanges the data of the two vectors.
- *
- * @relates TrilinosWrappers::MPI::Vector
- * @author Martin Kronbichler, Wolfgang Bangerth, 2008
- */
- inline
- void swap (Vector &u, Vector &v)
- {
- u.swap (v);
- }
+ /**
+ * Move the given vector. This operatorreplaces the present vector with
+ * @p v by efficiently swapping the internal data structures.
+ */
+ Vector &operator= (Vector &&v);
+ /**
+ * Another copy function. This one takes a deal.II vector and copies it
+ * into a TrilinosWrapper vector. Note that since we do not provide any
+ * Epetra_map that tells about the partitioning of the vector among the
+ * MPI processes, the size of the TrilinosWrapper vector has to be the
+ * same as the size of the input vector.
+ */
+ template <typename Number>
+ Vector &operator= (const ::dealii::Vector<Number> &v);
-#ifndef DOXYGEN
+ /**
+ * This reinit function is meant to be used for parallel calculations
+ * where some non-local data has to be used. The typical situation where
+ * one needs this function is the call of the
+ * FEValues<dim>::get_function_values function (or of some derivatives)
+ * in parallel. Since it is usually faster to retrieve the data in
+ * advance, this function can be called before the assembly forks out to
+ * the different processors. What this function does is the following:
+ * It takes the information in the columns of the given matrix and looks
+ * which data couples between the different processors. That data is
+ * then queried from the input vector. Note that you should not write to
+ * the resulting vector any more, since the some data can be stored
+ * several times on different processors, leading to unpredictable
+ * results. In particular, such a vector cannot be used for matrix-
+ * vector products as for example done during the solution of linear
+ * systems.
+ */
+ void import_nonlocal_data_for_fe
+ (const dealii::TrilinosWrappers::SparseMatrix &matrix,
+ const Vector &vector);
+
+ /**
+ * Test for equality. This function assumes that the present vector and
+ * the one to compare with have the same size already, since comparing
+ * vectors of different sizes makes not much sense anyway.
+ */
+ bool operator== (const Vector &v) const;
+
+ /**
+ * Test for inequality. This function assumes that the present vector and
+ * the one to compare with have the same size already, since comparing
+ * vectors of different sizes makes not much sense anyway.
+ */
+ bool operator!= (const Vector &v) const;
+
+ /**
+ * Return the global dimension of the vector.
+ */
+ size_type size () const;
+
+ /**
+ * Return the local dimension of the vector, i.e. the number of elements
+ * stored on the present MPI process. For sequential vectors, this number
+ * is the same as size(), but for parallel vectors it may be smaller.
+ *
+ * To figure out which elements exactly are stored locally, use
+ * local_range().
+ *
+ * If the vector contains ghost elements, they are included in this
+ * number.
+ */
+ size_type local_size () const;
+
+ /**
+ * Return a pair of indices indicating which elements of this vector are
+ * stored locally. The first number is the index of the first element
+ * stored, the second the index of the one past the last one that is
+ * stored locally. If this is a sequential vector, then the result will be
+ * the pair <code>(0,N)</code>, otherwise it will be a pair
+ * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
+ * <code>i</code> is the first element of the vector stored on this
+ * processor, corresponding to the half open interval $[i,i+n)$
+ *
+ * @note The description above is true most of the time, but not always.
+ * In particular, Trilinos vectors need not store contiguous ranges of
+ * elements such as $[i,i+n)$. Rather, it can store vectors where the
+ * elements are distributed in an arbitrary way across all processors and
+ * each processor simply stores a particular subset, not necessarily
+ * contiguous. In this case, this function clearly makes no sense since it
+ * could, at best, return a range that includes all elements that are
+ * stored locally. Thus, the function only succeeds if the locally stored
+ * range is indeed contiguous. It will trigger an assertion if the local
+ * portion of the vector is not contiguous.
+ */
+ std::pair<size_type, size_type> local_range () const;
+
+ /**
+ * Return whether @p index is in the local range or not, see also
+ * local_range().
+ *
+ * @note The same limitation for the applicability of this function
+ * applies as listed in the documentation of local_range().
+ */
+ bool in_local_range (const size_type index) const;
+
+ /**
+ * Return an index set that describes which elements of this vector are
+ * owned by the current processor. Note that this index set does not
+ * include elements this vector may store locally as ghost elements but
+ * that are in fact owned by another processor. As a consequence, the
+ * index sets returned on different processors if this is a distributed
+ * vector will form disjoint sets that add up to the complete index set.
+ * Obviously, if a vector is created on only one processor, then the
+ * result would satisfy
+ * @code
+ * vec.locally_owned_elements() == complete_index_set (vec.size())
+ * @endcode
+ */
+ IndexSet locally_owned_elements () const;
+
+ /**
+ * Return if the vector contains ghost elements. This answer is true if
+ * there are ghost elements on at least one process.
+ *
+ * @see
+ * @ref GlossGhostedVector "vectors with ghost elements"
+ */
+ bool has_ghost_elements() const;
+
+ /**
+ * This function only exists for compatibility with the @p
+ * LinearAlgebra::distributed::Vector class and does nothing: this class
+ * implements ghost value updates in a different way that is a better fit
+ * with the underlying Trilinos vector object.
+ */
+ void update_ghost_values () const;
+
+ /**
+ * Return the scalar (inner) product of two vectors. The vectors must have
+ * the same size.
+ */
+ TrilinosScalar operator* (const Vector &vec) const;
+
+ /**
+ * Return square of the $l_2$-norm.
+ */
+ real_type norm_sqr () const;
+
+ /**
+ * Mean value of the elements of this vector.
+ */
+ TrilinosScalar mean_value () const;
+
+ /**
+ * Compute the minimal value of the elements of this vector.
+ */
+ TrilinosScalar min () const;
+
+ /**
+ * Compute the maximal value of the elements of this vector.
+ */
+ TrilinosScalar max () const;
+
+ /**
+ * $l_1$-norm of the vector. The sum of the absolute values.
+ */
+ real_type l1_norm () const;
+
+ /**
+ * $l_2$-norm of the vector. The square root of the sum of the squares of
+ * the elements.
+ */
+ real_type l2_norm () const;
+
+ /**
+ * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
+ * <i>p</i>th powers of the absolute values of the elements.
+ */
+ real_type lp_norm (const TrilinosScalar p) const;
+
+ /**
+ * Maximum absolute value of the elements.
+ */
+ real_type linfty_norm () const;
+
+ /**
+ * Performs a combined operation of a vector addition and a subsequent
+ * inner product, returning the value of the inner product. In other
+ * words, the result of this function is the same as if the user called
+ * @code
+ * this->add(a, V);
+ * return_value = *this * W;
+ * @endcode
+ *
+ * The reason this function exists is for compatibility with deal.II's own
+ * vector classes which can implement this functionality with less memory
+ * transfer. However, for Trilinos vectors such a combined operation is
+ * not natively supported and thus the cost is completely equivalent as
+ * calling the two methods separately.
+ */
+ TrilinosScalar add_and_dot (const TrilinosScalar a,
+ const Vector &V,
+ const Vector &W);
+
+ /**
+ * Return whether the vector contains only elements with value zero. This
+ * is a collective operation. This function is expensive, because
+ * potentially all elements have to be checked.
+ */
+ bool all_zero () const;
+
+ /**
+ * Return @p true if the vector has no negative entries, i.e. all entries
+ * are zero or positive. This function is used, for example, to check
+ * whether refinement indicators are really all positive (or zero).
+ */
+ bool is_non_negative () const;
+ //@}
+
+
+ /**
+ * @name 2: Data-Access
+ */
+ //@{
+
+ /**
+ * Provide access to a given element, both read and write.
+ *
+ * When using a vector distributed with MPI, this operation only makes
+ * sense for elements that are actually present on the calling processor.
+ * Otherwise, an exception is thrown. This is different from the
+ * <code>el()</code> function below that always succeeds (but returns zero
+ * on non-local elements).
+ */
+ reference
+ operator() (const size_type index);
+
+ /**
+ * Provide read-only access to an element.
+ *
+ * When using a vector distributed with MPI, this operation only makes
+ * sense for elements that are actually present on the calling processor.
+ * Otherwise, an exception is thrown. This is different from the
+ * <code>el()</code> function below that always succeeds (but returns zero
+ * on non-local elements).
+ */
+ TrilinosScalar
+ operator() (const size_type index) const;
+
+ /**
+ * Provide access to a given element, both read and write.
+ *
+ * Exactly the same as operator().
+ */
+ reference
+ operator[] (const size_type index);
+
+ /**
+ * Provide read-only access to an element.
+ *
+ * Exactly the same as operator().
+ */
+ TrilinosScalar
+ operator[] (const size_type index) const;
+
+ /**
+ * A collective get operation: instead of getting individual elements of a
+ * vector, this function allows to get a whole set of elements at once.
+ * The indices of the elements to be read are stated in the first
+ * argument, the corresponding values are returned in the second.
+ */
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<TrilinosScalar> &values) const;
+
+ /**
+ * Just as the above, but with pointers. Useful in minimizing copying of
+ * data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
+
+ /**
+ * Make the Vector class a bit like the <tt>vector<></tt> class of the C++
+ * standard library by returning iterators to the start and end of the
+ * locally owned elements of this vector. The ordering of local elements
+ * corresponds to the one given by the global indices in case the vector
+ * is constructed from an IndexSet or other methods in deal.II (note that
+ * an Epetra_Map can contain elements in arbitrary orders, though).
+ *
+ * It holds that end() - begin() == local_size().
+ */
+ iterator begin ();
+
+ /**
+ * Return constant iterator to the start of the locally owned elements of
+ * the vector.
+ */
+ const_iterator begin () const;
+
+ /**
+ * Return an iterator pointing to the element past the end of the array of
+ * locally owned entries.
+ */
+ iterator end ();
+
+ /**
+ * Return a constant iterator pointing to the element past the end of the
+ * array of the locally owned entries.
+ */
+ const_iterator end () const;
+
+ //@}
+
+
+ /**
+ * @name 3: Modification of vectors
+ */
+ //@{
+
+ /**
+ * A collective set operation: instead of setting individual elements of a
+ * vector, this function allows to set a whole set of elements at once.
+ * The indices of the elements to be set are stated in the first argument,
+ * the corresponding values in the second.
+ */
+ void set (const std::vector<size_type> &indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * This is a second collective set operation. As a difference, this
+ * function takes a deal.II vector of values.
+ */
+ void set (const std::vector<size_type> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values);
+
+ /**
+ * This collective set operation is of lower level and can handle anything
+ * else — the only thing you have to provide is an address where all
+ * the indices are stored and the number of elements to be set.
+ */
+ void set (const size_type n_elements,
+ const size_type *indices,
+ const TrilinosScalar *values);
+
+ /**
+ * A collective add operation: This function adds a whole set of values
+ * stored in @p values to the vector components specified by @p indices.
+ */
+ void add (const std::vector<size_type> &indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * This is a second collective add operation. As a difference, this
+ * function takes a deal.II vector of values.
+ */
+ void add (const std::vector<size_type> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values);
+
+ /**
+ * Take an address where <tt>n_elements</tt> are stored contiguously and
+ * add them into the vector. Handles all cases which are not covered by
+ * the other two <tt>add()</tt> functions above.
+ */
+ void add (const size_type n_elements,
+ const size_type *indices,
+ const TrilinosScalar *values);
+
+ /**
+ * Multiply the entire vector by a fixed factor.
+ */
+ Vector &operator*= (const TrilinosScalar factor);
+
+ /**
+ * Divide the entire vector by a fixed factor.
+ */
+ Vector &operator/= (const TrilinosScalar factor);
+
+ /**
+ * Add the given vector to the present one.
+ */
+ Vector &operator+= (const Vector &V);
+
+ /**
+ * Subtract the given vector from the present one.
+ */
+ Vector &operator-= (const Vector &V);
+
+ /**
+ * Addition of @p s to all components. Note that @p s is a scalar and not
+ * a vector.
+ */
+ void add (const TrilinosScalar s);
+
+ /**
+ * Simple vector addition, equal to the <tt>operator+=</tt>.
+ *
+ * Though, if the second argument <tt>allow_different_maps</tt> is set,
+ * then it is possible to add data from a vector that uses a different
+ * map, i.e., a vector whose elements are split across processors
+ * differently. This may include vectors with ghost elements, for example.
+ * In general, however, adding vectors with a different element-to-
+ * processor map requires communicating data among processors and,
+ * consequently, is a slower operation than when using vectors using the
+ * same map.
+ */
+ void add (const Vector &V,
+ const bool allow_different_maps = false);
+
+ /**
+ * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+ */
+ void add (const TrilinosScalar a,
+ const Vector &V);
+
+ /**
+ * Multiple addition of scaled vectors, i.e. <tt>*this += a*V + b*W</tt>.
+ */
+ void add (const TrilinosScalar a,
+ const Vector &V,
+ const TrilinosScalar b,
+ const Vector &W);
+
+ /**
+ * Scaling and simple vector addition, i.e. <tt>*this = s*(*this) +
+ * V</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const Vector &V);
+
+ /**
+ * Scaling and simple addition, i.e. <tt>*this = s*(*this) + a*V</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const Vector &V);
+
+ /**
+ * Scale each element of this vector by the corresponding element in the
+ * argument. This function is mostly meant to simulate multiplication (and
+ * immediate re-assignment) by a diagonal scaling matrix.
+ */
+ void scale (const Vector &scaling_factors);
+
+ /**
+ * Assignment <tt>*this = a*V</tt>.
+ */
+ void equ (const TrilinosScalar a,
+ const Vector &V);
+ //@}
+
+ /**
+ * @name 4: Mixed stuff
+ */
+ //@{
+
+ /**
+ * Return a const reference to the underlying Trilinos Epetra_MultiVector
+ * class.
+ */
+ const Epetra_MultiVector &trilinos_vector () const;
+
+ /**
+ * Return a (modifyable) reference to the underlying Trilinos
+ * Epetra_FEVector class.
+ */
+ Epetra_FEVector &trilinos_vector ();
+
+ /**
+ * Return a const reference to the underlying Trilinos Epetra_Map that
+ * sets the parallel partitioning of the vector.
+ */
+ const Epetra_Map &vector_partitioner () const;
+
+ /**
+ * Print to a stream. @p precision denotes the desired precision with
+ * which values shall be printed, @p scientific whether scientific
+ * notation shall be used. If @p across is @p true then the vector is
+ * printed in a line, while if @p false then the elements are printed on a
+ * separate line each.
+ */
+ void print (std::ostream &out,
+ const unsigned int precision = 3,
+ const bool scientific = true,
+ const bool across = true) const;
+
+ /**
+ * Swap the contents of this vector and the other vector @p v. One could
+ * do this operation with a temporary variable and copying over the data
+ * elements, but this function is significantly more efficient since it
+ * only swaps the pointers to the data of the two vectors and therefore
+ * does not need to allocate temporary storage and move data around. Note
+ * that the vectors need to be of the same size and base on the same map.
+ *
+ * This function is analogous to the the @p swap function of all C++
+ * standard containers. Also, there is a global function
+ * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+ * analogy to standard functions.
+ */
+ void swap (Vector &v);
+
+ /**
+ * Estimate for the memory consumption in bytes.
+ */
+ std::size_t memory_consumption () const;
+
+ /**
+ * Return a reference to the MPI communicator object in use with this
+ * object.
+ */
+ const MPI_Comm &get_mpi_communicator () const;
+ //@}
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcDifferentParallelPartitioning);
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occurred while calling a Trilinos function");
+
+ /**
+ * Exception
+ */
+ DeclException4 (ExcAccessToNonLocalElement,
+ size_type, size_type, size_type, size_type,
+ << "You tried to access element " << arg1
+ << " of a distributed vector, but this element is not stored "
+ << "on the current processor. Note: There are "
+ << arg2 << " elements stored "
+ << "on the current processor from within the range "
+ << arg3 << " through " << arg4
+ << " but Trilinos vectors need not store contiguous "
+ << "ranges on each processor, and not every element in "
+ << "this range may in fact be stored locally.");
+
+ private:
+ /**
+ * Trilinos doesn't allow to mix additions to matrix entries and
+ * overwriting them (to make synchronisation of parallel computations
+ * simpler). The way we do it is to, for each access operation, store
+ * whether it is an insertion or an addition. If the previous one was of
+ * different type, then we first have to flush the Trilinos buffers;
+ * otherwise, we can simply go on. Luckily, Trilinos has an object for
+ * this which does already all the parallel communications in such a case,
+ * so we simply use their model, which stores whether the last operation
+ * was an addition or an insertion.
+ */
+ Epetra_CombineMode last_action;
+
+ /**
+ * A boolean variable to hold information on whether the vector is
+ * compressed or not.
+ */
+ bool compressed;
+
+ /**
+ * Whether this vector has ghost elements. This is true on all processors
+ * even if only one of them has any ghost elements.
+ */
+ bool has_ghosts;
+
+ /**
+ * Pointer to the actual Epetra vector object. This may represent a vector
+ * that is in fact distributed among multiple processors. The object
+ * requires an existing Epetra_Map for storing data when setting it up.
+ */
+ std::shared_ptr<Epetra_FEVector> vector;
+
+ /**
+ * A vector object in Trilinos to be used for collecting the non-local
+ * elements if the vector was constructed with an additional IndexSet
+ * describing ghost elements.
+ */
+ std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
+
+ /**
+ * An IndexSet storing the indices this vector owns exclusively.
+ */
+ IndexSet owned_elements;
+
+ /**
+ * Make the reference class a friend.
+ */
+ friend class internal::VectorReference;
+ friend class Vector;
+ friend class MPI::Vector;
+ };
+
+
+
+
+// ------------------- inline and template functions --------------
+
+
+ /**
+ * Global function @p swap which overloads the default implementation of
+ * the C++ standard library which uses a temporary object. The function
+ * simply exchanges the data of the two vectors.
+ *
+ * @relates TrilinosWrappers::MPI::Vector
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+ inline
+ void swap (Vector &u, Vector &v)
+ {
+ u.swap (v);
+ }
+ }
+
+#ifndef DOXYGEN
+
+ namespace internal
+ {
+ inline
+ VectorReference::VectorReference (MPI::Vector &vector,
+ const size_type index)
+ :
+ vector (vector),
+ index (index)
+ {}
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator= (const VectorReference &r) const
+ {
+ // as explained in the class
+ // documentation, this is not the copy
+ // operator. so simply pass on to the
+ // "correct" assignment operator
+ *this = static_cast<TrilinosScalar> (r);
+
+ return *this;
+ }
+
+
+
+ inline
+ VectorReference &
+ VectorReference::operator= (const VectorReference &r)
+ {
+ // as above
+ *this = static_cast<TrilinosScalar> (r);
+
+ return *this;
+ }
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator= (const TrilinosScalar &value) const
+ {
+ vector.set (1, &index, &value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator+= (const TrilinosScalar &value) const
+ {
+ vector.add (1, &index, &value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator-= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = -value;
+ vector.add (1, &index, &new_value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator*= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator/= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+ }
+
+ namespace MPI
+ {
+ inline
+ bool
+ Vector::in_local_range (const size_type index) const
+ {
+ std::pair<size_type, size_type> range = local_range();
+
+ return ((index >= range.first) && (index < range.second));
+ }
+
+
+
+ inline
+ IndexSet
+ Vector::locally_owned_elements() const
+ {
+ Assert(owned_elements.size()==size(),
+ ExcMessage("The locally owned elements have not been properly initialized!"
+ " This happens for example if this object has been initialized"
+ " with exactly one overlapping IndexSet."));
+ return owned_elements;
+ }
+
+
+
+ inline
+ bool
+ Vector::has_ghost_elements() const
+ {
+ return has_ghosts;
+ }
+
+
+
+ inline
+ void
+ Vector::update_ghost_values () const
+ {}
+
+
+
+ inline
+ internal::VectorReference
+ Vector::operator() (const size_type index)
+ {
+ return internal::VectorReference (*this, index);
+ }
+
+
+
+ inline
+ internal::VectorReference
+ Vector::operator[] (const size_type index)
+ {
+ return operator() (index);
+ }
+
+
+
+ inline
+ TrilinosScalar
+ Vector::operator[] (const size_type index) const
+ {
+ return operator() (index);
+ }
+
+
+
+ inline
+ void Vector::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<TrilinosScalar> &values) const
+ {
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+ }
+
+
+
+ template <typename ForwardIterator, typename OutputIterator>
+ inline
+ void Vector::extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+ {
+ while (indices_begin != indices_end)
+ {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++;
+ values_begin++;
+ }
+ }
+
+
+
+ inline
+ Vector::iterator
+ Vector::begin()
+ {
+ return (*vector)[0];
+ }
+
+
+
+ inline
+ Vector::iterator
+ Vector::end()
+ {
+ return (*vector)[0]+local_size();
+ }
+
+
+
+ inline
+ Vector::const_iterator
+ Vector::begin() const
+ {
+ return (*vector)[0];
+ }
+
+
+
+ inline
+ Vector::const_iterator
+ Vector::end() const
+ {
+ return (*vector)[0]+local_size();
+ }
+
+
+
+ inline
+ void
+ Vector::set (const std::vector<size_type> &indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ set (indices.size(), &indices[0], &values[0]);
+ }
+
+
+
+ inline
+ void
+ Vector::set (const std::vector<size_type> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ set (indices.size(), &indices[0], values.begin());
+ }
+
+
+
+ inline
+ void
+ Vector::set (const size_type n_elements,
+ const size_type *indices,
+ const TrilinosScalar *values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ if (last_action == Add)
+ {
+ const int ierr = vector->GlobalAssemble(Add);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+ if (last_action != Insert)
+ last_action = Insert;
+
+ for (size_type i=0; i<n_elements; ++i)
+ {
+ const size_type row = indices[i];
+ const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+ if (local_row != -1)
+ (*vector)[0][local_row] = values[i];
+ else
+ {
+ const int ierr = vector->ReplaceGlobalValues (1,
+ (const TrilinosWrappers::types::int_type *)(&row),
+ &values[i]);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ compressed = false;
+ }
+ // in set operation, do not use the pre-allocated vector for nonlocal
+ // entries even if it exists. This is to ensure that we really only
+ // set the elements touched by the set() method and not all contained
+ // in the nonlocal entries vector (there is no way to distinguish them
+ // on the receiving processor)
+ }
+ }
+
+
+
+ inline
+ void
+ Vector::add (const std::vector<size_type> &indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ add (indices.size(), &indices[0], &values[0]);
+ }
+
+
+
+ inline
+ void
+ Vector::add (const std::vector<size_type> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ add (indices.size(), &indices[0], values.begin());
+ }
+
+
+
+ inline
+ void
+ Vector::add (const size_type n_elements,
+ const size_type *indices,
+ const TrilinosScalar *values)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ if (last_action != Add)
+ {
+ if (last_action == Insert)
+ {
+ const int ierr = vector->GlobalAssemble(Insert);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+ last_action = Add;
+ }
+
+ for (size_type i=0; i<n_elements; ++i)
+ {
+ const size_type row = indices[i];
+ const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+ if (local_row != -1)
+ (*vector)[0][local_row] += values[i];
+ else if (nonlocal_vector.get() == nullptr)
+ {
+ const int ierr = vector->SumIntoGlobalValues (1,
+ (const TrilinosWrappers::types::int_type *)(&row),
+ &values[i]);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ compressed = false;
+ }
+ else
+ {
+ // use pre-allocated vector for non-local entries if it exists for
+ // addition operation
+ const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
+ Assert(my_row != -1,
+ ExcMessage("Attempted to write into off-processor vector entry "
+ "that has not be specified as being writable upon "
+ "initialization"));
+ (*nonlocal_vector)[0][my_row] += values[i];
+ compressed = false;
+ }
+ }
+ }
+
+
+
+ inline
+ Vector::size_type
+ Vector::size () const
+ {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+ return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
+#else
+ return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
+#endif
+ }
+
+
+
+ inline
+ Vector::size_type
+ Vector::local_size () const
+ {
+ return (size_type) vector->Map().NumMyElements();
+ }
+
+
+
+ inline
+ std::pair<Vector::size_type, Vector::size_type>
+ Vector::local_range () const
+ {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+ const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
+ const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
+#else
+ const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
+ const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
+#endif
+
+ Assert (end-begin == vector->Map().NumMyElements(),
+ ExcMessage ("This function only makes sense if the elements that this "
+ "vector stores on the current processor form a contiguous range. "
+ "This does not appear to be the case for the current vector."));
+
+ return std::make_pair (begin, end);
+ }
- template <typename number>
- Vector::Vector (const IndexSet ¶llel_partitioner,
- const dealii::Vector<number> &v,
- const MPI_Comm &communicator)
+
+
+ inline
+ TrilinosScalar
+ Vector::operator* (const Vector &vec) const
{
- *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
- v);
- owned_elements = parallel_partitioner;
+ Assert (vector->Map().SameAs(vec.vector->Map()),
+ ExcDifferentParallelPartitioning());
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ TrilinosScalar result;
+
+ const int ierr = vector->Dot(*(vec.vector), &result);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return result;
+ }
+
+
+
+ inline
+ Vector::real_type
+ Vector::norm_sqr () const
+ {
+ const TrilinosScalar d = l2_norm();
+ return d*d;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ Vector::mean_value () const
+ {
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ TrilinosScalar mean;
+ const int ierr = vector->MeanValue (&mean);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return mean;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ Vector::min () const
+ {
+ TrilinosScalar min_value;
+ const int ierr = vector->MinValue (&min_value);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return min_value;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ Vector::max () const
+ {
+ TrilinosScalar max_value;
+ const int ierr = vector->MaxValue (&max_value);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return max_value;
+ }
+
+
+
+ inline
+ Vector::real_type
+ Vector::l1_norm () const
+ {
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ TrilinosScalar d;
+ const int ierr = vector->Norm1 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ inline
+ Vector::real_type
+ Vector::l2_norm () const
+ {
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ TrilinosScalar d;
+ const int ierr = vector->Norm2 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ inline
+ Vector::real_type
+ Vector::lp_norm (const TrilinosScalar p) const
+ {
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
+ TrilinosScalar norm = 0;
+ TrilinosScalar sum=0;
+ const size_type n_local = local_size();
+
+ // loop over all the elements because
+ // Trilinos does not support lp norms
+ for (size_type i=0; i<n_local; i++)
+ sum += std::pow(std::fabs((*vector)[0][i]), p);
+
+ norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+
+ return norm;
+ }
+
+
+
+ inline
+ Vector::real_type
+ Vector::linfty_norm () const
+ {
+ // while we disallow the other
+ // norm operations on ghosted
+ // vectors, this particular norm
+ // is safe to run even in the
+ // presence of ghost elements
+ TrilinosScalar d;
+ const int ierr = vector->NormInf (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ Vector::add_and_dot (const TrilinosScalar a,
+ const Vector &V,
+ const Vector &W)
+ {
+ this->add(a, V);
+ return *this * W;
}
+ // inline also scalar products, vector
+ // additions etc. since they are all
+ // representable by a single Trilinos
+ // call. This reduces the overhead of the
+ // wrapper class.
inline
Vector &
- Vector::operator= (const TrilinosScalar s)
+ Vector::operator*= (const TrilinosScalar a)
+ {
+ AssertIsFinite(a);
+
+ const int ierr = vector->Scale(a);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ Vector &
+ Vector::operator/= (const TrilinosScalar a)
+ {
+ AssertIsFinite(a);
+
+ const TrilinosScalar factor = 1./a;
+
+ AssertIsFinite(factor);
+
+ const int ierr = vector->Scale(factor);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ Vector &
+ Vector::operator+= (const Vector &v)
{
- VectorBase::operator= (s);
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (vector->Map().SameAs(v.vector->Map()),
+ ExcDifferentParallelPartitioning());
+
+ const int ierr = vector->Update (1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
return *this;
}
- template <typename Number>
+
+ inline
Vector &
- Vector::operator= (const ::dealii::Vector<Number> &v)
+ Vector::operator-= (const Vector &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (vector->Map().SameAs(v.vector->Map()),
+ ExcDifferentParallelPartitioning());
+
+ const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ void
+ Vector::add (const TrilinosScalar s)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ AssertIsFinite(s);
+
+ size_type n_local = local_size();
+ for (size_type i=0; i<n_local; i++)
+ (*vector)[0][i] += s;
+ }
+
+
+
+ inline
+ void
+ Vector::add (const TrilinosScalar a,
+ const Vector &v)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (local_size() == v.local_size(),
+ ExcDimensionMismatch(local_size(), v.local_size()));
+
+ AssertIsFinite(a);
+
+ const int ierr = vector->Update(a, *(v.vector), 1.);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ Vector::add (const TrilinosScalar a,
+ const Vector &v,
+ const TrilinosScalar b,
+ const Vector &w)
{
- if (size() != v.size())
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (local_size() == v.local_size(),
+ ExcDimensionMismatch(local_size(), v.local_size()));
+ Assert (local_size() == w.local_size(),
+ ExcDimensionMismatch(local_size(), w.local_size()));
+
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ Vector::sadd (const TrilinosScalar s,
+ const Vector &v)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
+
+ AssertIsFinite(s);
+
+ // We assume that the vectors have the same Map
+ // if the local size is the same and if the vectors are not ghosted
+ if (local_size() == v.local_size() && !v.has_ghost_elements())
+ {
+ Assert (this->vector->Map().SameAs(v.vector->Map())==true,
+ ExcDifferentParallelPartitioning());
+ const int ierr = vector->Update(1., *(v.vector), s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+ else
+ {
+ (*this) *= s;
+ this->add(v, true);
+ }
+ }
+
+
+
+ inline
+ void
+ Vector::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const Vector &v)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
+ AssertIsFinite(s);
+ AssertIsFinite(a);
+
+ // We assume that the vectors have the same Map
+ // if the local size is the same and if the vectors are not ghosted
+ if (local_size() == v.local_size() && !v.has_ghost_elements())
+ {
+ Assert (this->vector->Map().SameAs(v.vector->Map())==true,
+ ExcDifferentParallelPartitioning());
+ const int ierr = vector->Update(a, *(v.vector), s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+ else
+ {
+ (*this) *= s;
+ Vector tmp = v;
+ tmp *= a;
+ this->add(tmp, true);
+ }
+ }
+
+
+
+ inline
+ void
+ Vector::scale (const Vector &factors)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ Assert (local_size() == factors.local_size(),
+ ExcDimensionMismatch(local_size(), factors.local_size()));
+
+ const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ Vector::equ (const TrilinosScalar a,
+ const Vector &v)
+ {
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ AssertIsFinite(a);
+
+ // If we don't have the same map, copy.
+ if (vector->Map().SameAs(v.vector->Map())==false)
+ {
+ this->sadd(0., a, v);
+ }
+ else
{
- vector.reset (new Epetra_FEVector(Epetra_Map
- (static_cast<TrilinosWrappers::types::int_type>(v.size()), 0,
+ // Otherwise, just update
+ int ierr = vector->Update(a, *v.vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ last_action = Zero;
+ }
+
+ }
+
+
+
+ inline
+ const Epetra_MultiVector &
+ Vector::trilinos_vector () const
+ {
+ return static_cast<const Epetra_MultiVector &>(*vector);
+ }
+
+
+
+ inline
+ Epetra_FEVector &
+ Vector::trilinos_vector ()
+ {
+ return *vector;
+ }
+
+
+
+ inline
+ const Epetra_Map &
+ Vector::vector_partitioner () const
+ {
+ return static_cast<const Epetra_Map &>(vector->Map());
+ }
+
+
+
+ inline
+ const MPI_Comm &
+ Vector::get_mpi_communicator () const
+ {
+ static MPI_Comm comm;
+
#ifdef DEAL_II_WITH_MPI
- Epetra_MpiComm(MPI_COMM_SELF)
+
+ const Epetra_MpiComm *mpi_comm
+ = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
+ comm = mpi_comm->Comm();
+
#else
- Epetra_SerialComm()
+
+ comm = MPI_COMM_SELF;
+
#endif
- )));
- }
- const Epetra_Map &map = vector_partitioner();
- const int size = map.NumMyElements();
+ return comm;
+ }
+
+ template <typename number>
+ Vector::Vector (const IndexSet ¶llel_partitioner,
+ const dealii::Vector<number> &v,
+ const MPI_Comm &communicator)
+ {
+ *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
+ v);
+ owned_elements = parallel_partitioner;
+ }
+
+
+
+ inline
+ Vector &
+ Vector::operator= (const TrilinosScalar s)
+ {
+ AssertIsFinite(s);
+
+ int ierr = vector->PutScalar(s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- // Need to copy out values, since the deal.II might not use doubles, so
- // that a direct access is not possible.
- for (int i=0; i<size; ++i)
- (*vector)[0][i] = v(gid(map,i));
+ if (nonlocal_vector.get() != nullptr)
+ {
+ ierr = nonlocal_vector->PutScalar(0.);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
return *this;
}
-
+ } /* end of namespace MPI */
#endif /* DOXYGEN */
- } /* end of namespace MPI */
} /* end of namespace TrilinosWrappers */
/*@}*/
}
};
- /**
- * A helper class used internally in linear_operator.h. Specialization for
- * TrilinosWrappers::VectorBase.
- */
- template<>
- class ReinitHelper<TrilinosWrappers::VectorBase>
- {
- public:
- template <typename Matrix>
- static
- void reinit_range_vector (const Matrix &matrix,
- TrilinosWrappers::VectorBase &v,
- bool omit_zeroing_entries)
- {
- v.reinit(matrix.locally_owned_range_indices(),
- matrix.get_mpi_communicator(),
- omit_zeroing_entries);
- }
-
- template <typename Matrix>
- static
- void reinit_domain_vector(const Matrix &matrix,
- TrilinosWrappers::VectorBase &v,
- bool omit_zeroing_entries)
- {
- v.reinit(matrix.locally_owned_domain_indices(),
- matrix.get_mpi_communicator(),
- omit_zeroing_entries);
- }
- };
-
} /* namespace LinearOperator */
} /* namespace internal */
-/**
- * Declare dealii::TrilinosWrappers::VectorBase as serial vector.
- *
- * @author Uwe Koecher, 2017
- */
-template <>
-struct is_serial_vector< TrilinosWrappers::VectorBase > : std::true_type
-{
-};
-
/**
* Declare dealii::TrilinosWrappers::MPI::Vector as distributed vector.
/*---------------------------- trilinos_vector.h ---------------------------*/
-#endif
-/*---------------------------- trilinos_vector.h ---------------------------*/
+#endif // dealii__trilinos_vector_h
#ifndef dealii__trilinos_vector_base_h
#define dealii__trilinos_vector_base_h
-
-#include <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_TRILINOS
-
-#include <deal.II/base/utilities.h>
-# include <deal.II/base/subscriptor.h>
-# include <deal.II/lac/exceptions.h>
-# include <deal.II/lac/vector.h>
-# include <deal.II/base/mpi.h>
-
-# include <vector>
-# include <utility>
-# include <memory>
-
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-# include "Epetra_ConfigDefs.h"
-# ifdef DEAL_II_WITH_MPI // only if MPI is installed
-# include "mpi.h"
-# include "Epetra_MpiComm.h"
-# else
-# include "Epetra_SerialComm.h"
-# endif
-# include "Epetra_FEVector.h"
-DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+#include <deal.II/lac/trilinos_vector.h>
DEAL_II_NAMESPACE_OPEN
-// forward declaration
-template <typename number> class Vector;
-
-
-/**
- * @addtogroup TrilinosWrappers
- * @{
- */
-
-/**
- * A namespace in which wrapper classes for Trilinos objects reside.
- *
- * @ingroup TrilinosWrappers
- */
+#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
- // forward declaration
- class VectorBase;
- namespace MPI
- {
- class Vector;
- }
-
- /**
- * @cond internal
- */
-
- /**
- * A namespace for internal implementation details of the TrilinosWrapper
- * members.
- *
- * @ingroup TrilinosWrappers
- */
- namespace internal
- {
- /**
- * Declare type for container size.
- */
- typedef dealii::types::global_dof_index size_type;
-
- /**
- * This class implements a wrapper for accessing the Trilinos vector in
- * the same way as we access deal.II objects: it is initialized with a
- * vector and an element within it, and has a conversion operator to
- * extract the scalar value of this element. It also has a variety of
- * assignment operator for writing to this one element.
- *
- * @ingroup TrilinosWrappers
- */
- class VectorReference
- {
- private:
- /**
- * Constructor. It is made private so as to only allow the actual vector
- * class to create it.
- */
- VectorReference (VectorBase &vector,
- const size_type index);
-
- public:
-
- /**
- * This looks like a copy operator, but does something different than
- * usual. In particular, it does not copy the member variables of this
- * reference. Rather, it handles the situation where we have two vectors
- * @p v and @p w, and assign elements like in <tt>v(i)=w(i)</tt>. Here,
- * both left and right hand side of the assignment have data type
- * VectorReference, but what we really mean is to assign the vector
- * elements represented by the two references. This operator implements
- * this operation. Note also that this allows us to make the assignment
- * operator const.
- */
- const VectorReference &
- operator = (const VectorReference &r) const;
-
- /**
- * Same as above but for non-const reference objects.
- */
- VectorReference &
- operator = (const VectorReference &r);
-
- /**
- * Set the referenced element of the vector to <tt>s</tt>.
- */
- const VectorReference &
- operator = (const TrilinosScalar &s) const;
-
- /**
- * Add <tt>s</tt> to the referenced element of the vector->
- */
- const VectorReference &
- operator += (const TrilinosScalar &s) const;
-
- /**
- * Subtract <tt>s</tt> from the referenced element of the vector->
- */
- const VectorReference &
- operator -= (const TrilinosScalar &s) const;
-
- /**
- * Multiply the referenced element of the vector by <tt>s</tt>.
- */
- const VectorReference &
- operator *= (const TrilinosScalar &s) const;
-
- /**
- * Divide the referenced element of the vector by <tt>s</tt>.
- */
- const VectorReference &
- operator /= (const TrilinosScalar &s) const;
-
- /**
- * Convert the reference to an actual value, i.e. return the value of
- * the referenced element of the vector.
- */
- operator TrilinosScalar () const;
-
- /**
- * Exception
- */
- DeclException1 (ExcTrilinosError,
- int,
- << "An error with error number " << arg1
- << " occurred while calling a Trilinos function");
-
- private:
- /**
- * Point to the vector we are referencing.
- */
- VectorBase &vector;
-
- /**
- * Index of the referenced element of the vector.
- */
- const size_type index;
-
- /**
- * Make the vector class a friend, so that it can create objects of the
- * present type.
- */
- friend class ::dealii::TrilinosWrappers::VectorBase;
- };
- }
- /**
- * @endcond
- */
-
-
- /**
- * Base class for the the distributed memory vector MPI::Vector that is
- * designed for use in either serial implementations or as a localized copy
- * on each processor. The implementation of this class is based on the
- * Trilinos vector class Epetra_FEVector, the (parallel) partitioning of which
- * is governed by an Epetra_Map. This means that the vector type is generic
- * and can be done in this base class, while the definition of the partition
- * map (and hence, the constructor and reinit function) will have to be done
- * in the derived classes. The Epetra_FEVector is precisely the kind of vector
- * we deal with all the time - we probably get it from some assembly
- * process, where also entries not locally owned might need to written and
- * hence need to be forwarded to the owner. The only requirement for this
- * class to work is that Trilinos is installed with the same compiler as is
- * used for compilation of deal.II.
- *
- * The interface of this class is modeled after the existing Vector class in
- * deal.II. It has almost the same member functions, and is often
- * exchangeable. However, since Trilinos only supports a single scalar type
- * (double), it is not templated, and only works with that type.
- *
- * Note that Trilinos only guarantees that operations do what you expect if
- * the function @p GlobalAssemble has been called after vector assembly in
- * order to distribute the data. Therefore, you need to call
- * Vector::compress() before you actually use the vectors.
- *
- * @ingroup TrilinosWrappers
- * @ingroup Vectors
- * @author Martin Kronbichler, 2008
- */
- class VectorBase : public Subscriptor
- {
- public:
- /**
- * Declare some of the standard types used in all containers. These types
- * parallel those in the <tt>C</tt> standard libraries
- * <tt>vector<...></tt> class.
- */
- typedef TrilinosScalar value_type;
- typedef TrilinosScalar real_type;
- typedef dealii::types::global_dof_index size_type;
- typedef value_type *iterator;
- typedef const value_type *const_iterator;
- typedef internal::VectorReference reference;
- typedef const internal::VectorReference const_reference;
-
- /**
- * @name 1: Basic Object-handling
- */
- //@{
-
- /**
- * Default constructor that generates an empty (zero size) vector. The
- * function <tt>reinit()</tt> will have to give the vector the correct
- * size and distribution among processes in case of an MPI run.
- */
- VectorBase ();
-
- /**
- * Copy constructor. Sets the dimension to that of the given vector, and
- * copies all the elements.
- */
- VectorBase (const VectorBase &v);
-
- /**
- * Destructor
- */
- virtual ~VectorBase ();
-
- /**
- * Release all memory and return to a state just like after having called
- * the default constructor.
- */
- void clear ();
-
- /**
- * Reinit functionality, sets the dimension and possibly the parallel
- * partitioning (Epetra_Map) of the calling vector to the settings of the
- * input vector.
- */
- void reinit (const VectorBase &v,
- const bool omit_zeroing_entries = false);
-
- /**
- * Compress the underlying representation of the Trilinos object, i.e.
- * flush the buffers of the vector object if it has any. This function is
- * necessary after writing into a vector element-by-element and before
- * anything else can be done on it.
- *
- * The (defaulted) argument can be used to specify the compress mode
- * (<code>Add</code> or <code>Insert</code>) in case the vector has not
- * been written to since the last time this function was called. The
- * argument is ignored if the vector has been added or written to since
- * the last time compress() was called.
- *
- * See
- * @ref GlossCompress "Compressing distributed objects"
- * for more information.
- */
- void compress (::dealii::VectorOperation::values operation);
-
- /**
- * Set all components of the vector to the given number @p s. Simply pass
- * this down to the Trilinos Epetra object, but we still need to declare
- * this function to make the example given in the discussion about making
- * the constructor explicit work.
- *
- * Since the semantics of assigning a scalar to a vector are not
- * immediately clear, this operator should really only be used if you want
- * to set the entire vector to zero. This allows the intuitive notation
- * <tt>v=0</tt>. Assigning other values is deprecated and may be
- * disallowed in the future.
- */
- VectorBase &
- operator = (const TrilinosScalar s);
-
- /**
- * Copy function. This function takes a VectorBase vector and copies all
- * the elements. The target vector will have the same parallel
- * distribution as the calling vector.
- */
- VectorBase &
- operator = (const VectorBase &v);
-
- /**
- * Another copy function. This one takes a deal.II vector and copies it
- * into a TrilinosWrapper vector. Note that since we do not provide any
- * Epetra_map that tells about the partitioning of the vector among the
- * MPI processes, the size of the TrilinosWrapper vector has to be the
- * same as the size of the input vector. In order to change the map, use
- * the reinit(const Epetra_Map &input_map) function.
- */
- template <typename Number>
- VectorBase &
- operator = (const ::dealii::Vector<Number> &v);
-
- /**
- * Test for equality. This function assumes that the present vector and
- * the one to compare with have the same size already, since comparing
- * vectors of different sizes makes not much sense anyway.
- */
- bool operator == (const VectorBase &v) const;
-
- /**
- * Test for inequality. This function assumes that the present vector and
- * the one to compare with have the same size already, since comparing
- * vectors of different sizes makes not much sense anyway.
- */
- bool operator != (const VectorBase &v) const;
-
- /**
- * Return the global dimension of the vector.
- */
- size_type size () const;
-
- /**
- * Return the local dimension of the vector, i.e. the number of elements
- * stored on the present MPI process. For sequential vectors, this number
- * is the same as size(), but for parallel vectors it may be smaller.
- *
- * To figure out which elements exactly are stored locally, use
- * local_range().
- *
- * If the vector contains ghost elements, they are included in this
- * number.
- */
- size_type local_size () const;
-
- /**
- * Return a pair of indices indicating which elements of this vector are
- * stored locally. The first number is the index of the first element
- * stored, the second the index of the one past the last one that is
- * stored locally. If this is a sequential vector, then the result will be
- * the pair <code>(0,N)</code>, otherwise it will be a pair
- * <code>(i,i+n)</code>, where <code>n=local_size()</code> and
- * <code>i</code> is the first element of the vector stored on this
- * processor, corresponding to the half open interval $[i,i+n)$
- *
- * @note The description above is true most of the time, but not always.
- * In particular, Trilinos vectors need not store contiguous ranges of
- * elements such as $[i,i+n)$. Rather, it can store vectors where the
- * elements are distributed in an arbitrary way across all processors and
- * each processor simply stores a particular subset, not necessarily
- * contiguous. In this case, this function clearly makes no sense since it
- * could, at best, return a range that includes all elements that are
- * stored locally. Thus, the function only succeeds if the locally stored
- * range is indeed contiguous. It will trigger an assertion if the local
- * portion of the vector is not contiguous.
- */
- std::pair<size_type, size_type> local_range () const;
-
- /**
- * Return whether @p index is in the local range or not, see also
- * local_range().
- *
- * @note The same limitation for the applicability of this function
- * applies as listed in the documentation of local_range().
- */
- bool in_local_range (const size_type index) const;
-
- /**
- * Return an index set that describes which elements of this vector are
- * owned by the current processor. Note that this index set does not
- * include elements this vector may store locally as ghost elements but
- * that are in fact owned by another processor. As a consequence, the
- * index sets returned on different processors if this is a distributed
- * vector will form disjoint sets that add up to the complete index set.
- * Obviously, if a vector is created on only one processor, then the
- * result would satisfy
- * @code
- * vec.locally_owned_elements() == complete_index_set (vec.size())
- * @endcode
- */
- IndexSet locally_owned_elements () const;
-
- /**
- * Return if the vector contains ghost elements. This answer is true if
- * there are ghost elements on at least one process.
- *
- * @see
- * @ref GlossGhostedVector "vectors with ghost elements"
- */
- bool has_ghost_elements() const;
-
- /**
- * This function only exists for compatibility with the @p
- * LinearAlgebra::distributed::Vector class and does nothing: this class
- * implements ghost value updates in a different way that is a better fit
- * with the underlying Trilinos vector object.
- */
- void update_ghost_values () const;
-
- /**
- * Return the scalar (inner) product of two vectors. The vectors must have
- * the same size.
- */
- TrilinosScalar operator * (const VectorBase &vec) const;
-
- /**
- * Return square of the $l_2$-norm.
- */
- real_type norm_sqr () const;
-
- /**
- * Mean value of the elements of this vector.
- */
- TrilinosScalar mean_value () const;
-
- /**
- * Compute the minimal value of the elements of this vector.
- */
- TrilinosScalar min () const;
-
- /**
- * Compute the maximal value of the elements of this vector.
- */
- TrilinosScalar max () const;
-
- /**
- * $l_1$-norm of the vector. The sum of the absolute values.
- */
- real_type l1_norm () const;
-
- /**
- * $l_2$-norm of the vector. The square root of the sum of the squares of
- * the elements.
- */
- real_type l2_norm () const;
-
- /**
- * $l_p$-norm of the vector. The <i>p</i>th root of the sum of the
- * <i>p</i>th powers of the absolute values of the elements.
- */
- real_type lp_norm (const TrilinosScalar p) const;
-
- /**
- * Maximum absolute value of the elements.
- */
- real_type linfty_norm () const;
-
- /**
- * Performs a combined operation of a vector addition and a subsequent
- * inner product, returning the value of the inner product. In other
- * words, the result of this function is the same as if the user called
- * @code
- * this->add(a, V);
- * return_value = *this * W;
- * @endcode
- *
- * The reason this function exists is for compatibility with deal.II's own
- * vector classes which can implement this functionality with less memory
- * transfer. However, for Trilinos vectors such a combined operation is
- * not natively supported and thus the cost is completely equivalent as
- * calling the two methods separately.
- */
- TrilinosScalar add_and_dot (const TrilinosScalar a,
- const VectorBase &V,
- const VectorBase &W);
-
- /**
- * Return whether the vector contains only elements with value zero. This
- * is a collective operation. This function is expensive, because
- * potentially all elements have to be checked.
- */
- bool all_zero () const;
-
- /**
- * Return @p true if the vector has no negative entries, i.e. all entries
- * are zero or positive. This function is used, for example, to check
- * whether refinement indicators are really all positive (or zero).
- */
- bool is_non_negative () const;
- //@}
-
-
- /**
- * @name 2: Data-Access
- */
- //@{
-
- /**
- * Provide access to a given element, both read and write.
- *
- * When using a vector distributed with MPI, this operation only makes
- * sense for elements that are actually present on the calling processor.
- * Otherwise, an exception is thrown. This is different from the
- * <code>el()</code> function below that always succeeds (but returns zero
- * on non-local elements).
- */
- reference
- operator () (const size_type index);
-
- /**
- * Provide read-only access to an element.
- *
- * When using a vector distributed with MPI, this operation only makes
- * sense for elements that are actually present on the calling processor.
- * Otherwise, an exception is thrown. This is different from the
- * <code>el()</code> function below that always succeeds (but returns zero
- * on non-local elements).
- */
- TrilinosScalar
- operator () (const size_type index) const;
-
- /**
- * Provide access to a given element, both read and write.
- *
- * Exactly the same as operator().
- */
- reference
- operator [] (const size_type index);
-
- /**
- * Provide read-only access to an element.
- *
- * Exactly the same as operator().
- */
- TrilinosScalar
- operator [] (const size_type index) const;
-
- /**
- * A collective get operation: instead of getting individual elements of a
- * vector, this function allows to get a whole set of elements at once.
- * The indices of the elements to be read are stated in the first
- * argument, the corresponding values are returned in the second.
- */
- void extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<TrilinosScalar> &values) const;
-
- /**
- * Just as the above, but with pointers. Useful in minimizing copying of
- * data around.
- */
- template <typename ForwardIterator, typename OutputIterator>
- void extract_subvector_to (ForwardIterator indices_begin,
- const ForwardIterator indices_end,
- OutputIterator values_begin) const;
-
- /**
- * Make the Vector class a bit like the <tt>vector<></tt> class of the C++
- * standard library by returning iterators to the start and end of the
- * locally owned elements of this vector. The ordering of local elements
- * corresponds to the one given by the global indices in case the vector
- * is constructed from an IndexSet or other methods in deal.II (note that
- * an Epetra_Map can contain elements in arbitrary orders, though).
- *
- * It holds that end() - begin() == local_size().
- */
- iterator begin ();
-
- /**
- * Return constant iterator to the start of the locally owned elements of
- * the vector.
- */
- const_iterator begin () const;
-
- /**
- * Return an iterator pointing to the element past the end of the array of
- * locally owned entries.
- */
- iterator end ();
-
- /**
- * Return a constant iterator pointing to the element past the end of the
- * array of the locally owned entries.
- */
- const_iterator end () const;
-
- //@}
-
-
- /**
- * @name 3: Modification of vectors
- */
- //@{
-
- /**
- * A collective set operation: instead of setting individual elements of a
- * vector, this function allows to set a whole set of elements at once.
- * The indices of the elements to be set are stated in the first argument,
- * the corresponding values in the second.
- */
- void set (const std::vector<size_type> &indices,
- const std::vector<TrilinosScalar> &values);
-
- /**
- * This is a second collective set operation. As a difference, this
- * function takes a deal.II vector of values.
- */
- void set (const std::vector<size_type> &indices,
- const ::dealii::Vector<TrilinosScalar> &values);
-
- /**
- * This collective set operation is of lower level and can handle anything
- * else — the only thing you have to provide is an address where all
- * the indices are stored and the number of elements to be set.
- */
- void set (const size_type n_elements,
- const size_type *indices,
- const TrilinosScalar *values);
-
- /**
- * A collective add operation: This function adds a whole set of values
- * stored in @p values to the vector components specified by @p indices.
- */
- void add (const std::vector<size_type> &indices,
- const std::vector<TrilinosScalar> &values);
-
- /**
- * This is a second collective add operation. As a difference, this
- * function takes a deal.II vector of values.
- */
- void add (const std::vector<size_type> &indices,
- const ::dealii::Vector<TrilinosScalar> &values);
-
- /**
- * Take an address where <tt>n_elements</tt> are stored contiguously and
- * add them into the vector. Handles all cases which are not covered by
- * the other two <tt>add()</tt> functions above.
- */
- void add (const size_type n_elements,
- const size_type *indices,
- const TrilinosScalar *values);
-
- /**
- * Multiply the entire vector by a fixed factor.
- */
- VectorBase &operator *= (const TrilinosScalar factor);
-
- /**
- * Divide the entire vector by a fixed factor.
- */
- VectorBase &operator /= (const TrilinosScalar factor);
-
- /**
- * Add the given vector to the present one.
- */
- VectorBase &operator += (const VectorBase &V);
-
- /**
- * Subtract the given vector from the present one.
- */
- VectorBase &operator -= (const VectorBase &V);
-
- /**
- * Addition of @p s to all components. Note that @p s is a scalar and not
- * a vector.
- */
- void add (const TrilinosScalar s);
-
- /**
- * Simple vector addition, equal to the <tt>operator +=</tt>.
- *
- * Though, if the second argument <tt>allow_different_maps</tt> is set,
- * then it is possible to add data from a vector that uses a different
- * map, i.e., a vector whose elements are split across processors
- * differently. This may include vectors with ghost elements, for example.
- * In general, however, adding vectors with a different element-to-
- * processor map requires communicating data among processors and,
- * consequently, is a slower operation than when using vectors using the
- * same map.
- */
- void add (const VectorBase &V,
- const bool allow_different_maps = false);
-
- /**
- * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
- */
- void add (const TrilinosScalar a,
- const VectorBase &V);
-
- /**
- * Multiple addition of scaled vectors, i.e. <tt>*this += a*V + b*W</tt>.
- */
- void add (const TrilinosScalar a,
- const VectorBase &V,
- const TrilinosScalar b,
- const VectorBase &W);
-
- /**
- * Scaling and simple vector addition, i.e. <tt>*this = s*(*this) +
- * V</tt>.
- */
- void sadd (const TrilinosScalar s,
- const VectorBase &V);
-
- /**
- * Scaling and simple addition, i.e. <tt>*this = s*(*this) + a*V</tt>.
- */
- void sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const VectorBase &V);
-
- /**
- * Scale each element of this vector by the corresponding element in the
- * argument. This function is mostly meant to simulate multiplication (and
- * immediate re-assignment) by a diagonal scaling matrix.
- */
- void scale (const VectorBase &scaling_factors);
-
- /**
- * Assignment <tt>*this = a*V</tt>.
- */
- void equ (const TrilinosScalar a,
- const VectorBase &V);
- //@}
-
-
- /**
- * @name 4: Mixed stuff
- */
- //@{
-
- /**
- * Return a const reference to the underlying Trilinos Epetra_MultiVector
- * class.
- */
- const Epetra_MultiVector &trilinos_vector () const;
-
- /**
- * Return a (modifyable) reference to the underlying Trilinos
- * Epetra_FEVector class.
- */
- Epetra_FEVector &trilinos_vector ();
-
- /**
- * Return a const reference to the underlying Trilinos Epetra_Map that
- * sets the parallel partitioning of the vector.
- */
- const Epetra_Map &vector_partitioner () const;
-
- /**
- * Print to a stream. @p precision denotes the desired precision with
- * which values shall be printed, @p scientific whether scientific
- * notation shall be used. If @p across is @p true then the vector is
- * printed in a line, while if @p false then the elements are printed on a
- * separate line each.
- */
- void print (std::ostream &out,
- const unsigned int precision = 3,
- const bool scientific = true,
- const bool across = true) const;
-
- /**
- * Swap the contents of this vector and the other vector @p v. One could
- * do this operation with a temporary variable and copying over the data
- * elements, but this function is significantly more efficient since it
- * only swaps the pointers to the data of the two vectors and therefore
- * does not need to allocate temporary storage and move data around. Note
- * that the vectors need to be of the same size and base on the same map.
- *
- * This function is analogous to the the @p swap function of all C++
- * standard containers. Also, there is a global function
- * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
- * analogy to standard functions.
- */
- void swap (VectorBase &v);
-
- /**
- * Estimate for the memory consumption in bytes.
- */
- std::size_t memory_consumption () const;
-
- /**
- * Return a reference to the MPI communicator object in use with this
- * object.
- */
- const MPI_Comm &get_mpi_communicator () const;
- //@}
-
- /**
- * Exception
- */
- DeclException0 (ExcDifferentParallelPartitioning);
-
- /**
- * Exception
- */
- DeclException1 (ExcTrilinosError,
- int,
- << "An error with error number " << arg1
- << " occurred while calling a Trilinos function");
-
- /**
- * Exception
- */
- DeclException4 (ExcAccessToNonLocalElement,
- size_type, size_type, size_type, size_type,
- << "You tried to access element " << arg1
- << " of a distributed vector, but this element is not stored "
- << "on the current processor. Note: There are "
- << arg2 << " elements stored "
- << "on the current processor from within the range "
- << arg3 << " through " << arg4
- << " but Trilinos vectors need not store contiguous "
- << "ranges on each processor, and not every element in "
- << "this range may in fact be stored locally.");
-
-
- private:
- /**
- * Trilinos doesn't allow to mix additions to matrix entries and
- * overwriting them (to make synchronisation of parallel computations
- * simpler). The way we do it is to, for each access operation, store
- * whether it is an insertion or an addition. If the previous one was of
- * different type, then we first have to flush the Trilinos buffers;
- * otherwise, we can simply go on. Luckily, Trilinos has an object for
- * this which does already all the parallel communications in such a case,
- * so we simply use their model, which stores whether the last operation
- * was an addition or an insertion.
- */
- Epetra_CombineMode last_action;
-
- /**
- * A boolean variable to hold information on whether the vector is
- * compressed or not.
- */
- bool compressed;
-
- /**
- * Whether this vector has ghost elements. This is true on all processors
- * even if only one of them has any ghost elements.
- */
- bool has_ghosts;
-
- /**
- * Pointer to the actual Epetra vector object. This may represent a vector
- * that is in fact distributed among multiple processors. The object
- * requires an existing Epetra_Map for storing data when setting it up.
- */
- std::shared_ptr<Epetra_FEVector> vector;
-
- /**
- * A vector object in Trilinos to be used for collecting the non-local
- * elements if the vector was constructed with an additional IndexSet
- * describing ghost elements.
- */
- std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
-
- /**
- * An IndexSet storing the indices this vector owns exclusively.
- */
- IndexSet owned_elements;
-
- /**
- * Make the reference class a friend.
- */
- friend class internal::VectorReference;
- friend class Vector;
- friend class MPI::Vector;
- };
-
-
-
-
-// ------------------- inline and template functions --------------
-
- /**
- * Global function swap which overloads the default implementation of the C
- * standard library which uses a temporary object. The function simply
- * exchanges the data of the two vectors.
- *
- * @relates TrilinosWrappers::VectorBase
- * @author Martin Kronbichler, Wolfgang Bangerth, 2008
- */
- inline
- void swap (VectorBase &u, VectorBase &v)
- {
- u.swap (v);
- }
-
-
-#ifndef DOXYGEN
-
- namespace internal
- {
- inline
- VectorReference::VectorReference (VectorBase &vector,
- const size_type index)
- :
- vector (vector),
- index (index)
- {}
-
-
- inline
- const VectorReference &
- VectorReference::operator = (const VectorReference &r) const
- {
- // as explained in the class
- // documentation, this is not the copy
- // operator. so simply pass on to the
- // "correct" assignment operator
- *this = static_cast<TrilinosScalar> (r);
-
- return *this;
- }
-
-
-
- inline
- VectorReference &
- VectorReference::operator = (const VectorReference &r)
- {
- // as above
- *this = static_cast<TrilinosScalar> (r);
-
- return *this;
- }
-
-
- inline
- const VectorReference &
- VectorReference::operator = (const TrilinosScalar &value) const
- {
- vector.set (1, &index, &value);
- return *this;
- }
-
-
-
- inline
- const VectorReference &
- VectorReference::operator += (const TrilinosScalar &value) const
- {
- vector.add (1, &index, &value);
- return *this;
- }
-
-
-
- inline
- const VectorReference &
- VectorReference::operator -= (const TrilinosScalar &value) const
- {
- TrilinosScalar new_value = -value;
- vector.add (1, &index, &new_value);
- return *this;
- }
-
-
-
- inline
- const VectorReference &
- VectorReference::operator *= (const TrilinosScalar &value) const
- {
- TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
- vector.set (1, &index, &new_value);
- return *this;
- }
-
-
-
- inline
- const VectorReference &
- VectorReference::operator /= (const TrilinosScalar &value) const
- {
- TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
- vector.set (1, &index, &new_value);
- return *this;
- }
- }
-
-
-
- inline
- bool
- VectorBase::in_local_range (const size_type index) const
- {
- std::pair<size_type, size_type> range = local_range();
-
- return ((index >= range.first) && (index < range.second));
- }
-
-
-
- inline
- IndexSet
- VectorBase::locally_owned_elements() const
- {
- Assert(owned_elements.size()==size(),
- ExcMessage("The locally owned elements have not been properly initialized!"
- " This happens for example if this object has been initialized"
- " with exactly one overlapping IndexSet."));
- return owned_elements;
- }
-
-
-
- inline
- bool
- VectorBase::has_ghost_elements() const
- {
- return has_ghosts;
- }
-
-
-
- inline
- void
- VectorBase::update_ghost_values () const
- {}
-
-
-
- inline
- internal::VectorReference
- VectorBase::operator () (const size_type index)
- {
- return internal::VectorReference (*this, index);
- }
-
-
-
- inline
- internal::VectorReference
- VectorBase::operator [] (const size_type index)
- {
- return operator() (index);
- }
-
-
- inline
- TrilinosScalar
- VectorBase::operator [] (const size_type index) const
- {
- return operator() (index);
- }
-
-
-
- inline
- void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<TrilinosScalar> &values) const
- {
- for (size_type i = 0; i < indices.size(); ++i)
- values[i] = operator()(indices[i]);
- }
-
-
-
- template <typename ForwardIterator, typename OutputIterator>
- inline
- void VectorBase::extract_subvector_to (ForwardIterator indices_begin,
- const ForwardIterator indices_end,
- OutputIterator values_begin) const
- {
- while (indices_begin != indices_end)
- {
- *values_begin = operator()(*indices_begin);
- indices_begin++;
- values_begin++;
- }
- }
-
-
-
- inline
- VectorBase::iterator
- VectorBase::begin()
- {
- return (*vector)[0];
- }
-
-
-
- inline
- VectorBase::iterator
- VectorBase::end()
- {
- return (*vector)[0]+local_size();
- }
-
-
-
- inline
- VectorBase::const_iterator
- VectorBase::begin() const
- {
- return (*vector)[0];
- }
-
-
-
- inline
- VectorBase::const_iterator
- VectorBase::end() const
- {
- return (*vector)[0]+local_size();
- }
-
-
-
- inline
- void
- VectorBase::reinit (const VectorBase &v,
- const bool omit_zeroing_entries)
- {
- Assert (vector.get() != nullptr,
- ExcMessage("Vector has not been constructed properly."));
-
- if (omit_zeroing_entries == false ||
- vector_partitioner().SameAs(v.vector_partitioner())==false)
- vector.reset (new Epetra_FEVector(*v.vector));
-
- if (v.nonlocal_vector.get() != nullptr)
- nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1));
-
- owned_elements = v.owned_elements;
- }
-
-
-
- inline
- VectorBase &
- VectorBase::operator = (const TrilinosScalar s)
- {
- AssertIsFinite(s);
-
- int ierr = vector->PutScalar(s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- if (nonlocal_vector.get() != nullptr)
- {
- ierr = nonlocal_vector->PutScalar(0.);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
- return *this;
- }
-
-
-
- inline
- void
- VectorBase::set (const std::vector<size_type> &indices,
- const std::vector<TrilinosScalar> &values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- set (indices.size(), &indices[0], &values[0]);
- }
-
-
-
- inline
- void
- VectorBase::set (const std::vector<size_type> &indices,
- const ::dealii::Vector<TrilinosScalar> &values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- set (indices.size(), &indices[0], values.begin());
- }
-
-
-
- inline
- void
- VectorBase::set (const size_type n_elements,
- const size_type *indices,
- const TrilinosScalar *values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- if (last_action == Add)
- {
- const int ierr = vector->GlobalAssemble(Add);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
- if (last_action != Insert)
- last_action = Insert;
-
- for (size_type i=0; i<n_elements; ++i)
- {
- const size_type row = indices[i];
- const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
- if (local_row != -1)
- (*vector)[0][local_row] = values[i];
- else
- {
- const int ierr = vector->ReplaceGlobalValues (1,
- (const TrilinosWrappers::types::int_type *)(&row),
- &values[i]);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- compressed = false;
- }
- // in set operation, do not use the pre-allocated vector for nonlocal
- // entries even if it exists. This is to ensure that we really only
- // set the elements touched by the set() method and not all contained
- // in the nonlocal entries vector (there is no way to distinguish them
- // on the receiving processor)
- }
- }
-
-
-
- inline
- void
- VectorBase::add (const std::vector<size_type> &indices,
- const std::vector<TrilinosScalar> &values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- add (indices.size(), &indices[0], &values[0]);
- }
-
-
-
- inline
- void
- VectorBase::add (const std::vector<size_type> &indices,
- const ::dealii::Vector<TrilinosScalar> &values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- add (indices.size(), &indices[0], values.begin());
- }
-
-
-
- inline
- void
- VectorBase::add (const size_type n_elements,
- const size_type *indices,
- const TrilinosScalar *values)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- if (last_action != Add)
- {
- if (last_action == Insert)
- {
- const int ierr = vector->GlobalAssemble(Insert);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
- last_action = Add;
- }
-
- for (size_type i=0; i<n_elements; ++i)
- {
- const size_type row = indices[i];
- const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
- if (local_row != -1)
- (*vector)[0][local_row] += values[i];
- else if (nonlocal_vector.get() == nullptr)
- {
- const int ierr = vector->SumIntoGlobalValues (1,
- (const TrilinosWrappers::types::int_type *)(&row),
- &values[i]);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- compressed = false;
- }
- else
- {
- // use pre-allocated vector for non-local entries if it exists for
- // addition operation
- const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
- Assert(my_row != -1,
- ExcMessage("Attempted to write into off-processor vector entry "
- "that has not be specified as being writable upon "
- "initialization"));
- (*nonlocal_vector)[0][my_row] += values[i];
- compressed = false;
- }
- }
- }
-
-
-
- inline
- VectorBase::size_type
- VectorBase::size () const
- {
-#ifndef DEAL_II_WITH_64BIT_INDICES
- return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
-#else
- return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
-#endif
- }
-
-
-
- inline
- VectorBase::size_type
- VectorBase::local_size () const
- {
- return (size_type) vector->Map().NumMyElements();
- }
-
-
-
- inline
- std::pair<VectorBase::size_type, VectorBase::size_type>
- VectorBase::local_range () const
- {
-#ifndef DEAL_II_WITH_64BIT_INDICES
- const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
- const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
-#else
- const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
- const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
-#endif
-
- Assert (end-begin == vector->Map().NumMyElements(),
- ExcMessage ("This function only makes sense if the elements that this "
- "vector stores on the current processor form a contiguous range. "
- "This does not appear to be the case for the current vector."));
-
- return std::make_pair (begin, end);
- }
-
-
-
- inline
- TrilinosScalar
- VectorBase::operator * (const VectorBase &vec) const
- {
- Assert (vector->Map().SameAs(vec.vector->Map()),
- ExcDifferentParallelPartitioning());
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- TrilinosScalar result;
-
- const int ierr = vector->Dot(*(vec.vector), &result);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return result;
- }
-
-
-
- inline
- VectorBase::real_type
- VectorBase::norm_sqr () const
- {
- const TrilinosScalar d = l2_norm();
- return d*d;
- }
-
-
-
- inline
- TrilinosScalar
- VectorBase::mean_value () const
- {
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- TrilinosScalar mean;
- const int ierr = vector->MeanValue (&mean);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return mean;
- }
-
-
-
- inline
- TrilinosScalar
- VectorBase::min () const
- {
- TrilinosScalar min_value;
- const int ierr = vector->MinValue (&min_value);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return min_value;
- }
-
-
-
- inline
- TrilinosScalar
- VectorBase::max () const
- {
- TrilinosScalar max_value;
- const int ierr = vector->MaxValue (&max_value);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return max_value;
- }
-
-
-
- inline
- VectorBase::real_type
- VectorBase::l1_norm () const
- {
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- TrilinosScalar d;
- const int ierr = vector->Norm1 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
-
- inline
- VectorBase::real_type
- VectorBase::l2_norm () const
- {
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- TrilinosScalar d;
- const int ierr = vector->Norm2 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
-
- inline
- VectorBase::real_type
- VectorBase::lp_norm (const TrilinosScalar p) const
- {
- Assert (!has_ghost_elements(), ExcGhostsPresent());
-
- TrilinosScalar norm = 0;
- TrilinosScalar sum=0;
- const size_type n_local = local_size();
-
- // loop over all the elements because
- // Trilinos does not support lp norms
- for (size_type i=0; i<n_local; i++)
- sum += std::pow(std::fabs((*vector)[0][i]), p);
-
- norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
-
- return norm;
- }
-
-
-
- inline
- VectorBase::real_type
- VectorBase::linfty_norm () const
- {
- // while we disallow the other
- // norm operations on ghosted
- // vectors, this particular norm
- // is safe to run even in the
- // presence of ghost elements
- TrilinosScalar d;
- const int ierr = vector->NormInf (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
-
- inline
- TrilinosScalar
- VectorBase::add_and_dot (const TrilinosScalar a,
- const VectorBase &V,
- const VectorBase &W)
- {
- this->add(a, V);
- return *this * W;
- }
-
-
-
- // inline also scalar products, vector
- // additions etc. since they are all
- // representable by a single Trilinos
- // call. This reduces the overhead of the
- // wrapper class.
- inline
- VectorBase &
- VectorBase::operator *= (const TrilinosScalar a)
- {
- AssertIsFinite(a);
-
- const int ierr = vector->Scale(a);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- inline
- VectorBase &
- VectorBase::operator /= (const TrilinosScalar a)
- {
- AssertIsFinite(a);
-
- const TrilinosScalar factor = 1./a;
-
- AssertIsFinite(factor);
-
- const int ierr = vector->Scale(factor);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- inline
- VectorBase &
- VectorBase::operator += (const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- Assert (vector->Map().SameAs(v.vector->Map()),
- ExcDifferentParallelPartitioning());
-
- const int ierr = vector->Update (1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- inline
- VectorBase &
- VectorBase::operator -= (const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- Assert (vector->Map().SameAs(v.vector->Map()),
- ExcDifferentParallelPartitioning());
-
- const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- inline
- void
- VectorBase::add (const TrilinosScalar s)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- AssertIsFinite(s);
-
- size_type n_local = local_size();
- for (size_type i=0; i<n_local; i++)
- (*vector)[0][i] += s;
- }
-
-
-
- inline
- void
- VectorBase::add (const TrilinosScalar a,
- const VectorBase &v)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (local_size() == v.local_size(),
- ExcDimensionMismatch(local_size(), v.local_size()));
-
- AssertIsFinite(a);
-
- const int ierr = vector->Update(a, *(v.vector), 1.);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- inline
- void
- VectorBase::add (const TrilinosScalar a,
- const VectorBase &v,
- const TrilinosScalar b,
- const VectorBase &w)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (local_size() == v.local_size(),
- ExcDimensionMismatch(local_size(), v.local_size()));
- Assert (local_size() == w.local_size(),
- ExcDimensionMismatch(local_size(), w.local_size()));
-
- AssertIsFinite(a);
- AssertIsFinite(b);
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- inline
- void
- VectorBase::sadd (const TrilinosScalar s,
- const VectorBase &v)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
-
- AssertIsFinite(s);
-
- // We assume that the vectors have the same Map
- // if the local size is the same and if the vectors are not ghosted
- if (local_size() == v.local_size() && !v.has_ghost_elements())
- {
- Assert (this->vector->Map().SameAs(v.vector->Map())==true,
- ExcDifferentParallelPartitioning());
- const int ierr = vector->Update(1., *(v.vector), s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
- else
- {
- (*this) *= s;
- this->add(v, true);
- }
- }
-
-
-
- inline
- void
- VectorBase::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const VectorBase &v)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
- AssertIsFinite(s);
- AssertIsFinite(a);
-
- // We assume that the vectors have the same Map
- // if the local size is the same and if the vectors are not ghosted
- if (local_size() == v.local_size() && !v.has_ghost_elements())
- {
- Assert (this->vector->Map().SameAs(v.vector->Map())==true,
- ExcDifferentParallelPartitioning());
- const int ierr = vector->Update(a, *(v.vector), s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
- else
- {
- (*this) *= s;
- VectorBase tmp = v;
- tmp *= a;
- this->add(tmp, true);
- }
- }
-
-
-
- inline
- void
- VectorBase::scale (const VectorBase &factors)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- Assert (local_size() == factors.local_size(),
- ExcDimensionMismatch(local_size(), factors.local_size()));
-
- const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- inline
- void
- VectorBase::equ (const TrilinosScalar a,
- const VectorBase &v)
- {
- // if we have ghost values, do not allow
- // writing to this vector at all.
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- AssertIsFinite(a);
-
- // If we don't have the same map, copy.
- if (vector->Map().SameAs(v.vector->Map())==false)
- {
- this->sadd(0., a, v);
- }
- else
- {
- // Otherwise, just update
- int ierr = vector->Update(a, *v.vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Zero;
- }
-
- }
-
-
-
- inline
- const Epetra_MultiVector &
- VectorBase::trilinos_vector () const
- {
- return static_cast<const Epetra_MultiVector &>(*vector);
- }
-
-
-
- inline
- Epetra_FEVector &
- VectorBase::trilinos_vector ()
- {
- return *vector;
- }
-
-
-
- inline
- const Epetra_Map &
- VectorBase::vector_partitioner () const
- {
- return static_cast<const Epetra_Map &>(vector->Map());
- }
-
-
-
- inline
- const MPI_Comm &
- VectorBase::get_mpi_communicator () const
- {
- static MPI_Comm comm;
-
-#ifdef DEAL_II_WITH_MPI
-
- const Epetra_MpiComm *mpi_comm
- = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
- comm = mpi_comm->Comm();
-
-#else
-
- comm = MPI_COMM_SELF;
-
-#endif
-
- return comm;
- }
-
-
-
-#endif // DOXYGEN
-
+ // All functionality from the former base class VectorBase was moved to
+ // MPI::Vector as this ended up to be the only derived class.
+ // Hence, using the class VectorBase is deprecated and we only provide
+ // an alias for backward compatibility.
+ using VectorBase = MPI::Vector;
}
-
-/*@}*/
+#endif
DEAL_II_NAMESPACE_CLOSE
-#endif // DEAL_II_WITH_TRILINOS
-
-/*---------------------------- trilinos_vector_base.h ---------------------------*/
-
-#endif
-/*---------------------------- trilinos_vector_base.h ---------------------------*/
+#endif // dealii__trilinos_vector_base_h
trilinos_solver.cc
trilinos_sparse_matrix.cc
trilinos_sparsity_pattern.cc
- trilinos_vector_base.cc
trilinos_vector.cc
)
SET(_inst
${_inst}
trilinos_sparse_matrix.inst.in
- trilinos_vector_base.inst.in
+ trilinos_vector.inst.in
)
ENDIF()
//
// ---------------------------------------------------------------------
-#include <deal.II/lac/trilinos_index_access.h>
#include <deal.II/lac/trilinos_vector.h>
#ifdef DEAL_II_WITH_TRILINOS
# include <deal.II/base/mpi.h>
+# include <deal.II/lac/trilinos_sparse_matrix.h>
# include <deal.II/lac/trilinos_parallel_block_vector.h>
# include <deal.II/lac/trilinos_index_access.h>
-# include <deal.II/lac/trilinos_sparse_matrix.h>
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
# include <Epetra_Import.h>
+# include <Epetra_Export.h>
# include <Epetra_Vector.h>
+
+# include <boost/io/ios_state.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
# include <cmath>
namespace TrilinosWrappers
{
- namespace MPI
+ namespace internal
{
- Vector::Vector ()
+ VectorReference::operator TrilinosScalar () const
{
- last_action = Zero;
- vector.reset(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self())));
+ Assert (index < vector.size(),
+ ExcIndexRange (index, 0, vector.size()));
+
+ // Trilinos allows for vectors to be referenced by the [] or ()
+ // operators but only () checks index bounds. We check these bounds by
+ // ourselves, so we can use []. Note that we can only get local values.
+
+ const TrilinosWrappers::types::int_type local_index =
+ vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
+ Assert (local_index >= 0,
+ MPI::Vector::ExcAccessToNonLocalElement (index, vector.local_size(),
+ vector.vector->Map().MinMyGID(),
+ vector.vector->Map().MaxMyGID()));
+
+
+ return (*(vector.vector))[0][local_index];
}
+ }
+
+ namespace MPI
+ {
+ Vector::Vector ()
+ :
+ last_action (Zero),
+ compressed (true),
+ has_ghosts (false),
+#ifdef DEAL_II_WITH_MPI
+ vector(new Epetra_FEVector(
+ Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_SELF))))
+#else
+ vector(new Epetra_FEVector(
+ Epetra_Map(0,0,Epetra_SerialComm())))
+#endif
+ {}
Vector::Vector (const Vector &v)
:
- VectorBase()
- {
- last_action = Zero;
- vector.reset (new Epetra_FEVector(*v.vector));
- has_ghosts = v.has_ghosts;
- owned_elements = v.owned_elements;
- }
+ Subscriptor(),
+ last_action (Zero),
+ compressed (true),
+ has_ghosts (v.has_ghosts),
+ vector(new Epetra_FEVector(*v.vector)),
+ owned_elements(v.owned_elements)
+ {}
Vector::Vector (const IndexSet ¶llel_partitioner,
- const VectorBase &v,
+ const Vector &v,
const MPI_Comm &communicator)
:
- VectorBase()
+ last_action (Zero),
+ compressed (true),
+ has_ghosts (false)
{
AssertThrow (parallel_partitioner.size() ==
static_cast<size_type>(TrilinosWrappers::n_global_elements(v.vector->Map())),
ExcDimensionMismatch (parallel_partitioner.size(),
TrilinosWrappers::n_global_elements(v.vector->Map())));
- last_action = Zero;
-
vector.reset (new Epetra_FEVector
(parallel_partitioner.make_trilinos_map(communicator,
true)));
reinit (v, false, true);
}
+
+
Vector::Vector (const IndexSet &local,
const IndexSet &ghost,
const MPI_Comm &communicator)
:
- VectorBase()
+ last_action (Zero),
+ compressed (true)
{
reinit(local, ghost, communicator, false);
}
+ void
+ Vector::clear ()
+ {
+ // When we clear the vector, reset the pointer and generate an empty
+ // vector.
+#ifdef DEAL_II_WITH_MPI
+ Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
+#else
+ Epetra_Map map (0, 0, Epetra_SerialComm());
+#endif
+
+ has_ghosts = false;
+ vector.reset (new Epetra_FEVector(map));
+ last_action = Zero;
+ }
+
+
+
void
Vector::reinit (const IndexSet ¶llel_partitioner,
const MPI_Comm &communicator,
void
- Vector::reinit (const VectorBase &v,
+ Vector::reinit (const Vector &v,
const bool omit_zeroing_entries,
const bool allow_different_maps)
{
Vector &
Vector::operator = (const Vector &v)
{
+ Assert (vector.get() != nullptr,
+ ExcMessage("Vector is not constructed properly."));
+
// check equality for MPI communicators to avoid accessing a possibly
// invalid MPI_Comm object
#ifdef DEAL_II_WITH_MPI
+ template <typename number>
Vector &
- Vector::operator = (const TrilinosWrappers::VectorBase &v)
+ Vector::operator = (const ::dealii::Vector<number> &v)
{
- nonlocal_vector.reset();
-
- Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
-
- Epetra_Import data_exchange (vector->Map(), v.vector->Map());
- const int ierr = vector->Import(*v.vector, data_exchange, Insert);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
- last_action = Insert;
+ // this is probably not very efficient but works. in particular, we could
+ // do better if we know that number==TrilinosScalar because then we could
+ // elide the copying of elements
+ //
+ // let's hope this isn't a particularly frequent operation
+ std::pair<size_type, size_type> local_range = this->local_range ();
+ for (size_type i=local_range.first; i<local_range.second; ++i)
+ (*vector)[0][i-local_range.first] = v(i);
return *this;
}
last_action = Insert;
}
+
+
+ void
+ Vector::compress (::dealii::VectorOperation::values given_last_action)
+ {
+ //Select which mode to send to Trilinos. Note that we use last_action if
+ //available and ignore what the user tells us to detect wrongly mixed
+ //operations. Typically given_last_action is only used on machines that do
+ //not execute an operation (because they have no own cells for example).
+ Epetra_CombineMode mode = last_action;
+ if (last_action == Zero)
+ {
+ if (given_last_action==::dealii::VectorOperation::add)
+ mode = Add;
+ else if (given_last_action==::dealii::VectorOperation::insert)
+ mode = Insert;
+ }
+
+#ifdef DEBUG
+# ifdef DEAL_II_WITH_MPI
+ // check that every process has decided to use the same mode. This will
+ // otherwise result in undefined behaviour in the call to
+ // GlobalAssemble().
+ double double_mode = mode;
+ const Epetra_MpiComm *comm_ptr
+ = dynamic_cast<const Epetra_MpiComm *>(&(vector_partitioner().Comm()));
+ Assert (comm_ptr != nullptr, ExcInternalError());
+ Utilities::MPI::MinMaxAvg result
+ = Utilities::MPI::min_max_avg (double_mode, comm_ptr->GetMpiComm());
+ Assert(result.max-result.min<1e-5,
+ ExcMessage ("Not all processors agree whether the last operation on "
+ "this vector was an addition or a set operation. This will "
+ "prevent the compress() operation from succeeding."));
+
+# endif
+#endif
+
+ // Now pass over the information about what we did last to the vector.
+ int ierr = 0;
+ if (nonlocal_vector.get() == nullptr || mode != Add)
+ ierr = vector->GlobalAssemble(mode);
+ else
+ {
+ Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
+ ierr = vector->Export(*nonlocal_vector, exporter, mode);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ ierr = nonlocal_vector->PutScalar(0.);
+ }
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ last_action = Zero;
+
+ compressed = true;
+ }
+
+
+
+ TrilinosScalar
+ Vector::operator () (const size_type index) const
+ {
+ // Extract local indices in the vector.
+ TrilinosWrappers::types::int_type trilinos_i =
+ vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
+ TrilinosScalar value = 0.;
+
+ // If the element is not present on the current processor, we can't
+ // continue. This is the main difference to the el() function.
+ if (trilinos_i == -1)
+ {
+ Assert (false, ExcAccessToNonLocalElement(index, local_size(),
+ vector->Map().MinMyGID(),
+ vector->Map().MaxMyGID()));
+ }
+ else
+ value = (*vector)[0][trilinos_i];
+
+ return value;
+ }
+
+
+
+ void
+ Vector::add (const Vector &v,
+ const bool allow_different_maps)
+ {
+ if (allow_different_maps == false)
+ *this += v;
+ else
+ {
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+ AssertThrow (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
+
+#if DEAL_II_TRILINOS_VERSION_GTE(11,11,0)
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+ int ierr = vector->Import(*v.vector, data_exchange, Epetra_AddLocalAlso);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ last_action = Add;
+#else
+ // In versions older than 11.11 the Import function is broken for adding
+ // Hence, we provide a workaround in this case
+
+ Epetra_MultiVector dummy(vector->Map(), 1, false);
+ Epetra_Import data_exchange (dummy.Map(), v.vector->Map());
+
+ int ierr = dummy.Import(*v.vector, data_exchange, Insert);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = vector->Update (1.0, dummy, 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+#endif
+ }
+ }
+
+
+
+ bool
+ Vector::operator == (const Vector &v) const
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ if (local_size() != v.local_size())
+ return false;
+
+ size_type i;
+ for (i=0; i<local_size(); i++)
+ if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
+
+ return true;
+ }
+
+
+
+ bool
+ Vector::operator != (const Vector &v) const
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ return (!(*this==v));
+ }
+
+
+
+ bool
+ Vector::all_zero () const
+ {
+ // get a representation of the vector and
+ // loop over all the elements
+ TrilinosScalar *start_ptr = (*vector)[0];
+ const TrilinosScalar *ptr = start_ptr,
+ *eptr = start_ptr + local_size();
+ unsigned int flag = 0;
+ while (ptr != eptr)
+ {
+ if (*ptr != 0)
+ {
+ flag = 1;
+ break;
+ }
+ ++ptr;
+ }
+
+#ifdef DEAL_II_WITH_MPI
+ // in parallel, check that the vector
+ // is zero on _all_ processors.
+ const Epetra_MpiComm *mpi_comm
+ = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
+ Assert(mpi_comm != nullptr, ExcInternalError());
+ unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
+ return num_nonzero == 0;
+#else
+ return flag == 0;
+#endif
+
+ }
+
+
+
+ bool
+ Vector::is_non_negative () const
+ {
+#ifdef DEAL_II_WITH_MPI
+ // if this vector is a parallel one, then
+ // we need to communicate to determine
+ // the answer to the current
+ // function. this still has to be
+ // implemented
+ AssertThrow(local_size() == size(), ExcNotImplemented());
+#endif
+ // get a representation of the vector and
+ // loop over all the elements
+ TrilinosScalar *start_ptr;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // TODO: This
+ // won't work in parallel like
+ // this. Find out a better way to
+ // this in that case.
+ const TrilinosScalar *ptr = start_ptr,
+ *eptr = start_ptr + size();
+ bool flag = true;
+ while (ptr != eptr)
+ {
+ if (*ptr < 0.0)
+ {
+ flag = false;
+ break;
+ }
+ ++ptr;
+ }
+
+ return flag;
+ }
+
+
+
+ void
+ Vector::print (std::ostream &out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool across) const
+ {
+ AssertThrow (out, ExcIO());
+ boost::io::ios_flags_saver restore_flags(out);
+
+ // get a representation of the
+ // vector and loop over all
+ // the elements TODO: up to
+ // now only local data printed
+ // out! Find a way to neatly
+ // output distributed data...
+ TrilinosScalar *val;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&val, &leading_dimension);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ out.precision (precision);
+ if (scientific)
+ out.setf (std::ios::scientific, std::ios::floatfield);
+ else
+ out.setf (std::ios::fixed, std::ios::floatfield);
+
+ if (across)
+ for (size_type i=0; i<size(); ++i)
+ out << static_cast<double>(val[i]) << ' ';
+ else
+ for (size_type i=0; i<size(); ++i)
+ out << static_cast<double>(val[i]) << std::endl;
+ out << std::endl;
+
+ // restore the representation
+ // of the vector
+ AssertThrow (out, ExcIO());
+ }
+
+
+
+ void
+ Vector::swap (Vector &v)
+ {
+ std::swap(last_action, v.last_action);
+ std::swap(compressed, v.compressed);
+ std::swap(vector, v.vector);
+ }
+
+
+
+ std::size_t
+ Vector::memory_consumption () const
+ {
+ //TODO[TH]: No accurate memory
+ //consumption for Trilinos vectors
+ //yet. This is a rough approximation with
+ //one index and the value per local
+ //entry.
+ return sizeof(*this)
+ + this->local_size()*( sizeof(double)+
+ sizeof(TrilinosWrappers::types::int_type) );
+ }
} /* end of namespace MPI */
} /* end of namespace TrilinosWrappers */
+namespace TrilinosWrappers
+{
+ namespace MPI
+ {
+# include "trilinos_vector.inst"
+ }
+}
+
DEAL_II_NAMESPACE_CLOSE
#endif // DEAL_II_WITH_TRILINOS
for (SCALAR : REAL_SCALARS)
{
template
- VectorBase &
- VectorBase::operator = (const ::dealii::Vector<SCALAR> &v);
+ Vector &
+ Vector::operator = (const ::dealii::Vector<SCALAR> &v);
}
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2008 - 2017 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/base/memory_consumption.h>
-#include <deal.II/lac/trilinos_vector_base.h>
-
-#ifdef DEAL_II_WITH_TRILINOS
-
-# include <deal.II/base/mpi.h>
-
-# include <cmath>
-
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-# include <Epetra_Import.h>
-# include <Epetra_Export.h>
-
-# include <boost/io/ios_state.hpp>
-DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
-
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace TrilinosWrappers
-{
- namespace internal
- {
- VectorReference::operator TrilinosScalar () const
- {
- Assert (index < vector.size(),
- ExcIndexRange (index, 0, vector.size()));
-
- // Trilinos allows for vectors to be referenced by the [] or ()
- // operators but only () checks index bounds. We check these bounds by
- // ourselves, so we can use []. Note that we can only get local values.
-
- const TrilinosWrappers::types::int_type local_index =
- vector.vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
- Assert (local_index >= 0,
- VectorBase::ExcAccessToNonLocalElement (index, vector.local_size(),
- vector.vector->Map().MinMyGID(),
- vector.vector->Map().MaxMyGID()));
-
-
- return (*(vector.vector))[0][local_index];
- }
- }
-
-
-
- VectorBase::VectorBase ()
- :
- last_action (Zero),
- compressed (true),
- has_ghosts (false),
-#ifdef DEAL_II_WITH_MPI
- vector(new Epetra_FEVector(
- Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_SELF))))
-#else
- vector(new Epetra_FEVector(
- Epetra_Map(0,0,Epetra_SerialComm())))
-#endif
- {}
-
-
-
- VectorBase::VectorBase (const VectorBase &v)
- :
- Subscriptor(),
- last_action (Zero),
- compressed (true),
- has_ghosts (v.has_ghosts),
- vector(new Epetra_FEVector(*v.vector)),
- owned_elements(complete_index_set(v.size()))
- {}
-
-
-
- VectorBase::~VectorBase ()
- {}
-
-
-
- void
- VectorBase::clear ()
- {
- // When we clear the vector, reset the pointer and generate an empty
- // vector.
-#ifdef DEAL_II_WITH_MPI
- Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_SELF));
-#else
- Epetra_Map map (0, 0, Epetra_SerialComm());
-#endif
-
- has_ghosts = false;
- vector.reset (new Epetra_FEVector(map));
- last_action = Zero;
- }
-
-
-
- VectorBase &
- VectorBase::operator = (const VectorBase &v)
- {
- Assert (vector.get() != nullptr,
- ExcMessage("Vector is not constructed properly."));
-
- if (local_range() != v.local_range())
- {
- last_action = Zero;
- vector.reset (new Epetra_FEVector(*v.vector));
- has_ghosts = v.has_ghosts;
- owned_elements = v.owned_elements;
- }
- else
- {
- Assert (vector->Map().SameAs(v.vector->Map()) == true,
- ExcMessage ("The Epetra maps in the assignment operator ="
- " do not match, even though the local_range "
- " seems to be the same. Check vector setup!"));
- int ierr;
- ierr = vector->GlobalAssemble(last_action);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- ierr = vector->Update(1.0, *v.vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Zero;
- }
-
- return *this;
- }
-
-
-
- template <typename number>
- VectorBase &
- VectorBase::operator = (const ::dealii::Vector<number> &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- // this is probably not very efficient
- // but works. in particular, we could do
- // better if we know that
- // number==TrilinosScalar because then we
- // could elide the copying of elements
- //
- // let's hope this isn't a
- // particularly frequent operation
- std::pair<size_type, size_type>
- local_range = this->local_range ();
- for (size_type i=local_range.first; i<local_range.second; ++i)
- (*vector)[0][i-local_range.first] = v(i);
-
- return *this;
- }
-
-
-
- void
- VectorBase::compress (::dealii::VectorOperation::values given_last_action)
- {
- //Select which mode to send to Trilinos. Note that we use last_action if
- //available and ignore what the user tells us to detect wrongly mixed
- //operations. Typically given_last_action is only used on machines that do
- //not execute an operation (because they have no own cells for example).
- Epetra_CombineMode mode = last_action;
- if (last_action == Zero)
- {
- if (given_last_action==::dealii::VectorOperation::add)
- mode = Add;
- else if (given_last_action==::dealii::VectorOperation::insert)
- mode = Insert;
- }
-
-#ifdef DEBUG
-# ifdef DEAL_II_WITH_MPI
- // check that every process has decided to use the same mode. This will
- // otherwise result in undefined behaviour in the call to
- // GlobalAssemble().
- double double_mode = mode;
- const Epetra_MpiComm *comm_ptr
- = dynamic_cast<const Epetra_MpiComm *>(&(vector_partitioner().Comm()));
- Assert (comm_ptr != nullptr, ExcInternalError());
- Utilities::MPI::MinMaxAvg result
- = Utilities::MPI::min_max_avg (double_mode, comm_ptr->GetMpiComm());
- Assert(result.max-result.min<1e-5,
- ExcMessage ("Not all processors agree whether the last operation on "
- "this vector was an addition or a set operation. This will "
- "prevent the compress() operation from succeeding."));
-
-# endif
-#endif
-
- // Now pass over the information about what we did last to the vector.
- int ierr = 0;
- if (nonlocal_vector.get() == nullptr || mode != Add)
- ierr = vector->GlobalAssemble(mode);
- else
- {
- Epetra_Export exporter(nonlocal_vector->Map(), vector->Map());
- ierr = vector->Export(*nonlocal_vector, exporter, mode);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = nonlocal_vector->PutScalar(0.);
- }
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- last_action = Zero;
-
- compressed = true;
- }
-
-
-
- TrilinosScalar
- VectorBase::operator () (const size_type index) const
- {
- // Extract local indices in the vector.
- TrilinosWrappers::types::int_type trilinos_i =
- vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(index));
- TrilinosScalar value = 0.;
-
- // If the element is not present on the current processor, we can't
- // continue. This is the main difference to the el() function.
- if (trilinos_i == -1)
- {
- Assert (false, ExcAccessToNonLocalElement(index, local_size(),
- vector->Map().MinMyGID(),
- vector->Map().MaxMyGID()));
- }
- else
- value = (*vector)[0][trilinos_i];
-
- return value;
- }
-
-
-
- void
- VectorBase::add (const VectorBase &v,
- const bool allow_different_maps)
- {
- if (allow_different_maps == false)
- *this += v;
- else
- {
- Assert (!has_ghost_elements(), ExcGhostsPresent());
- AssertThrow (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
-
-#if DEAL_II_TRILINOS_VERSION_GTE(11,11,0)
- Epetra_Import data_exchange (vector->Map(), v.vector->Map());
- int ierr = vector->Import(*v.vector, data_exchange, Epetra_AddLocalAlso);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- last_action = Add;
-#else
- // In versions older than 11.11 the Import function is broken for adding
- // Hence, we provide a workaround in this case
-
- Epetra_MultiVector dummy(vector->Map(), 1, false);
- Epetra_Import data_exchange (dummy.Map(), v.vector->Map());
-
- int ierr = dummy.Import(*v.vector, data_exchange, Insert);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- ierr = vector->Update (1.0, dummy, 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-#endif
- }
- }
-
-
-
- bool
- VectorBase::operator == (const VectorBase &v) const
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- if (local_size() != v.local_size())
- return false;
-
- size_type i;
- for (i=0; i<local_size(); i++)
- if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
-
- return true;
- }
-
-
-
- bool
- VectorBase::operator != (const VectorBase &v) const
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- return (!(*this==v));
- }
-
-
-
- bool
- VectorBase::all_zero () const
- {
- // get a representation of the vector and
- // loop over all the elements
- TrilinosScalar *start_ptr = (*vector)[0];
- const TrilinosScalar *ptr = start_ptr,
- *eptr = start_ptr + local_size();
- unsigned int flag = 0;
- while (ptr != eptr)
- {
- if (*ptr != 0)
- {
- flag = 1;
- break;
- }
- ++ptr;
- }
-
-#ifdef DEAL_II_WITH_MPI
- // in parallel, check that the vector
- // is zero on _all_ processors.
- const Epetra_MpiComm *mpi_comm
- = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
- Assert(mpi_comm != nullptr, ExcInternalError());
- unsigned int num_nonzero = Utilities::MPI::sum(flag, mpi_comm->Comm());
- return num_nonzero == 0;
-#else
- return flag == 0;
-#endif
-
- }
-
-
-
- bool
- VectorBase::is_non_negative () const
- {
-#ifdef DEAL_II_WITH_MPI
- // if this vector is a parallel one, then
- // we need to communicate to determine
- // the answer to the current
- // function. this still has to be
- // implemented
- AssertThrow(local_size() == size(), ExcNotImplemented());
-#endif
- // get a representation of the vector and
- // loop over all the elements
- TrilinosScalar *start_ptr;
- int leading_dimension;
- int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- // TODO: This
- // won't work in parallel like
- // this. Find out a better way to
- // this in that case.
- const TrilinosScalar *ptr = start_ptr,
- *eptr = start_ptr + size();
- bool flag = true;
- while (ptr != eptr)
- {
- if (*ptr < 0.0)
- {
- flag = false;
- break;
- }
- ++ptr;
- }
-
- return flag;
- }
-
-
-
- void
- VectorBase::print (std::ostream &out,
- const unsigned int precision,
- const bool scientific,
- const bool across) const
- {
- AssertThrow (out, ExcIO());
- boost::io::ios_flags_saver restore_flags(out);
-
- // get a representation of the
- // vector and loop over all
- // the elements TODO: up to
- // now only local data printed
- // out! Find a way to neatly
- // output distributed data...
- TrilinosScalar *val;
- int leading_dimension;
- int ierr = vector->ExtractView (&val, &leading_dimension);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- out.precision (precision);
- if (scientific)
- out.setf (std::ios::scientific, std::ios::floatfield);
- else
- out.setf (std::ios::fixed, std::ios::floatfield);
-
- if (across)
- for (size_type i=0; i<size(); ++i)
- out << static_cast<double>(val[i]) << ' ';
- else
- for (size_type i=0; i<size(); ++i)
- out << static_cast<double>(val[i]) << std::endl;
- out << std::endl;
-
- // restore the representation
- // of the vector
- AssertThrow (out, ExcIO());
- }
-
-
-
- void
- VectorBase::swap (VectorBase &v)
- {
- std::swap(last_action, v.last_action);
- std::swap(compressed, v.compressed);
- std::swap(vector, v.vector);
- }
-
-
-
- std::size_t
- VectorBase::memory_consumption () const
- {
- //TODO[TH]: No accurate memory
- //consumption for Trilinos vectors
- //yet. This is a rough approximation with
- //one index and the value per local
- //entry.
- return sizeof(*this)
- + this->local_size()*( sizeof(double)+
- sizeof(TrilinosWrappers::types::int_type) );
- }
-
-} /* end of namespace TrilinosWrappers */
-
-
-namespace TrilinosWrappers
-{
-#include "trilinos_vector_base.inst"
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_WITH_TRILINOS
projects / dealii.git / commitdiff