From: Daniel Arndt Date: Tue, 30 May 2017 22:25:23 +0000 (+0200) Subject: Remove TrilinosWrappers::VectorBase X-Git-Tag: v9.0.0-rc1~1536^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4581769fc9a03e7f5eabb00a534c88eaa70bebe9;p=dealii.git Remove TrilinosWrappers::VectorBase --- diff --git a/include/deal.II/lac/trilinos_vector.h b/include/deal.II/lac/trilinos_vector.h index 60ffa95c6d..8a8bb874e8 100644 --- a/include/deal.II/lac/trilinos_vector.h +++ b/include/deal.II/lac/trilinos_vector.h @@ -20,36 +20,172 @@ #include #ifdef DEAL_II_WITH_TRILINOS - # include # include # include +# include # include # include -# include # include +# include +# include +# include + 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 - DEAL_II_NAMESPACE_OPEN -// forward declaration -template 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 v(i)=w(i). 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 s. + */ + const VectorReference & + operator= (const TrilinosScalar &s) const; + + /** + * Add s to the referenced element of the vector-> + */ + const VectorReference & + operator+= (const TrilinosScalar &s) const; + + /** + * Subtract s from the referenced element of the vector-> + */ + const VectorReference & + operator-= (const TrilinosScalar &s) const; + + /** + * Multiply the referenced element of the vector by s. + */ + const VectorReference & + operator*= (const TrilinosScalar &s) const; + + /** + * Divide the referenced element of the vector by s. + */ + 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 @@ -75,11 +211,9 @@ namespace TrilinosWrappers /** * 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 { @@ -87,9 +221,19 @@ namespace TrilinosWrappers /** * 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 @@ -159,7 +303,7 @@ namespace TrilinosWrappers * 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; isize(); ++i) * vector(i) = i; @@ -250,18 +394,27 @@ namespace TrilinosWrappers * * @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 C standard libraries + * vector<...> 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 */ //@{ /** @@ -276,116 +429,6 @@ namespace TrilinosWrappers */ 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 v. If omit_zeroing_entries - * is not true, the elements in the vector are initialized with - * zero. If it is set to true, 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, allow_different_maps, - * 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 v - * is a distributed vector. In this case, the variable - * omit_zeroing_entries needs to be set to false, - * 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 - Vector &operator= (const ::dealii::Vector &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::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 @@ -437,7 +480,7 @@ namespace TrilinosWrappers * @ref GlossGhostedVector "vectors with ghost elements" */ Vector (const IndexSet ¶llel_partitioning, - const VectorBase &v, + const Vector &v, const MPI_Comm &communicator = MPI_COMM_WORLD); /** @@ -457,6 +500,50 @@ namespace TrilinosWrappers const dealii::Vector &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 v. If omit_zeroing_entries + * is not true, the elements in the vector are initialized with + * zero. If it is set to true, 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, allow_different_maps, + * 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 v + * is a distributed vector. In this case, the variable + * omit_zeroing_entries needs to be set to false, + * 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 @@ -467,7 +554,6 @@ namespace TrilinosWrappers * 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 @@ -513,85 +599,1562 @@ namespace TrilinosWrappers 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 + * (Add or Insert) 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 + * v=0. + */ + 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 + Vector &operator= (const ::dealii::Vector &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::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 (0,N), otherwise it will be a pair + * (i,i+n), where n=local_size() and + * i 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 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 pth root of the sum of the + * pth 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 + * el() 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 + * el() 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 &indices, + std::vector &values) const; + + /** + * Just as the above, but with pointers. Useful in minimizing copying of + * data around. + */ + template + void extract_subvector_to (ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const; + + /** + * Make the Vector class a bit like the vector<> 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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &values); + + /** + * Take an address where n_elements are stored contiguously and + * add them into the vector. Handles all cases which are not covered by + * the other two add() 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 operator+=. + * + * Though, if the second argument allow_different_maps 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. *this += a*V. + */ + void add (const TrilinosScalar a, + const Vector &V); + + /** + * Multiple addition of scaled vectors, i.e. *this += a*V + b*W. + */ + void add (const TrilinosScalar a, + const Vector &V, + const TrilinosScalar b, + const Vector &W); + + /** + * Scaling and simple vector addition, i.e. *this = s*(*this) + + * V. + */ + void sadd (const TrilinosScalar s, + const Vector &V); + + /** + * Scaling and simple addition, i.e. *this = s*(*this) + a*V. + */ + 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 *this = a*V. + */ + 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 + * swap(u,v) that simply calls u.swap(v), 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 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 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 (r); + + return *this; + } + + + + inline + VectorReference & + VectorReference::operator= (const VectorReference &r) + { + // as above + *this = static_cast (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(*this) * value; + vector.set (1, &index, &new_value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator/= (const TrilinosScalar &value) const + { + TrilinosScalar new_value = static_cast(*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 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 &indices, + std::vector &values) const + { + for (size_type i = 0; i < indices.size(); ++i) + values[i] = operator()(indices[i]); + } + + + + template + 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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &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; iMap().LID(static_cast(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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &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; iMap().LID(static_cast(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(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::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 - Vector::Vector (const IndexSet ¶llel_partitioner, - const dealii::Vector &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(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 + + inline Vector & - Vector::operator= (const ::dealii::Vector &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; iUpdate(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(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(*vector); + } + + + + inline + Epetra_FEVector & + Vector::trilinos_vector () + { + return *vector; + } + + + + inline + const Epetra_Map & + Vector::vector_partitioner () const + { + return static_cast(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(&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 + Vector::Vector (const IndexSet ¶llel_partitioner, + const dealii::Vector &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; iPutScalar(0.); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } return *this; } - + } /* end of namespace MPI */ #endif /* DOXYGEN */ - } /* end of namespace MPI */ } /* end of namespace TrilinosWrappers */ /*@}*/ @@ -630,51 +2193,10 @@ namespace internal } }; - /** - * A helper class used internally in linear_operator.h. Specialization for - * TrilinosWrappers::VectorBase. - */ - template<> - class ReinitHelper - { - public: - template - 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 - 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. @@ -693,5 +2215,4 @@ DEAL_II_NAMESPACE_CLOSE /*---------------------------- trilinos_vector.h ---------------------------*/ -#endif -/*---------------------------- trilinos_vector.h ---------------------------*/ +#endif // dealii__trilinos_vector_h diff --git a/include/deal.II/lac/trilinos_vector_base.h b/include/deal.II/lac/trilinos_vector_base.h index 523db9b0d9..10b5f44f6b 100644 --- a/include/deal.II/lac/trilinos_vector_base.h +++ b/include/deal.II/lac/trilinos_vector_base.h @@ -16,1810 +16,21 @@ #ifndef dealii__trilinos_vector_base_h #define dealii__trilinos_vector_base_h - -#include - -#ifdef DEAL_II_WITH_TRILINOS - -#include -# include -# include -# include -# include - -# include -# include -# include - -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_NAMESPACE_OPEN -// forward declaration -template 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 v(i)=w(i). 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 s. - */ - const VectorReference & - operator = (const TrilinosScalar &s) const; - - /** - * Add s to the referenced element of the vector-> - */ - const VectorReference & - operator += (const TrilinosScalar &s) const; - - /** - * Subtract s from the referenced element of the vector-> - */ - const VectorReference & - operator -= (const TrilinosScalar &s) const; - - /** - * Multiply the referenced element of the vector by s. - */ - const VectorReference & - operator *= (const TrilinosScalar &s) const; - - /** - * Divide the referenced element of the vector by s. - */ - 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 C standard libraries - * vector<...> 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 reinit() 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 - * (Add or Insert) 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 - * v=0. 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 - VectorBase & - operator = (const ::dealii::Vector &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 (0,N), otherwise it will be a pair - * (i,i+n), where n=local_size() and - * i 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 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 pth root of the sum of the - * pth 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 - * el() 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 - * el() 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 &indices, - std::vector &values) const; - - /** - * Just as the above, but with pointers. Useful in minimizing copying of - * data around. - */ - template - void extract_subvector_to (ForwardIterator indices_begin, - const ForwardIterator indices_end, - OutputIterator values_begin) const; - - /** - * Make the Vector class a bit like the vector<> 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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &values); - - /** - * Take an address where n_elements are stored contiguously and - * add them into the vector. Handles all cases which are not covered by - * the other two add() 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 operator +=. - * - * Though, if the second argument allow_different_maps 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. *this += a*V. - */ - void add (const TrilinosScalar a, - const VectorBase &V); - - /** - * Multiple addition of scaled vectors, i.e. *this += a*V + b*W. - */ - void add (const TrilinosScalar a, - const VectorBase &V, - const TrilinosScalar b, - const VectorBase &W); - - /** - * Scaling and simple vector addition, i.e. *this = s*(*this) + - * V. - */ - void sadd (const TrilinosScalar s, - const VectorBase &V); - - /** - * Scaling and simple addition, i.e. *this = s*(*this) + a*V. - */ - 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 *this = a*V. - */ - 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 - * swap(u,v) that simply calls u.swap(v), 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 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 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 (r); - - return *this; - } - - - - inline - VectorReference & - VectorReference::operator = (const VectorReference &r) - { - // as above - *this = static_cast (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(*this) * value; - vector.set (1, &index, &new_value); - return *this; - } - - - - inline - const VectorReference & - VectorReference::operator /= (const TrilinosScalar &value) const - { - TrilinosScalar new_value = static_cast(*this) / value; - vector.set (1, &index, &new_value); - return *this; - } - } - - - - inline - bool - VectorBase::in_local_range (const size_type index) const - { - std::pair 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 &indices, - std::vector &values) const - { - for (size_type i = 0; i < indices.size(); ++i) - values[i] = operator()(indices[i]); - } - - - - template - 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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &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; iMap().LID(static_cast(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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &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; iMap().LID(static_cast(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(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::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(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; iUpdate(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(*vector); - } - - - - inline - Epetra_FEVector & - VectorBase::trilinos_vector () - { - return *vector; - } - - - - inline - const Epetra_Map & - VectorBase::vector_partitioner () const - { - return static_cast(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(&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 diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 5597406a10..3c6dba8d11 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -121,13 +121,12 @@ IF(DEAL_II_WITH_TRILINOS) 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() diff --git a/source/lac/trilinos_vector.cc b/source/lac/trilinos_vector.cc index 446456765c..acab363c1d 100644 --- a/source/lac/trilinos_vector.cc +++ b/source/lac/trilinos_vector.cc @@ -13,19 +13,21 @@ // // --------------------------------------------------------------------- -#include #include #ifdef DEAL_II_WITH_TRILINOS # include +# include # include # include -# include DEAL_II_DISABLE_EXTRA_DIAGNOSTICS # include +# include # include + +# include DEAL_II_ENABLE_EXTRA_DIAGNOSTICS # include @@ -35,13 +37,44 @@ DEAL_II_NAMESPACE_OPEN 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(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 + {} @@ -55,13 +88,13 @@ namespace TrilinosWrappers 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) + {} @@ -78,29 +111,32 @@ namespace TrilinosWrappers 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(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); } @@ -111,6 +147,24 @@ namespace TrilinosWrappers + 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, @@ -150,7 +204,7 @@ namespace TrilinosWrappers void - Vector::reinit (const VectorBase &v, + Vector::reinit (const Vector &v, const bool omit_zeroing_entries, const bool allow_different_maps) { @@ -367,6 +421,9 @@ namespace TrilinosWrappers 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 @@ -443,19 +500,21 @@ namespace TrilinosWrappers + template Vector & - Vector::operator = (const TrilinosWrappers::VectorBase &v) + Vector::operator = (const ::dealii::Vector &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 local_range = this->local_range (); + for (size_type i=local_range.first; i(&(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(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(&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(val[i]) << ' '; + else + for (size_type i=0; i(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 diff --git a/source/lac/trilinos_vector_base.inst.in b/source/lac/trilinos_vector.inst.in similarity index 89% rename from source/lac/trilinos_vector_base.inst.in rename to source/lac/trilinos_vector.inst.in index 2f1c4d3d3c..e94ff81dc9 100644 --- a/source/lac/trilinos_vector_base.inst.in +++ b/source/lac/trilinos_vector.inst.in @@ -18,6 +18,6 @@ for (SCALAR : REAL_SCALARS) { template - VectorBase & - VectorBase::operator = (const ::dealii::Vector &v); + Vector & + Vector::operator = (const ::dealii::Vector &v); } diff --git a/source/lac/trilinos_vector_base.cc b/source/lac/trilinos_vector_base.cc deleted file mode 100644 index ec456bf657..0000000000 --- a/source/lac/trilinos_vector_base.cc +++ /dev/null @@ -1,462 +0,0 @@ -// --------------------------------------------------------------------- -// -// 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 -#include - -#ifdef DEAL_II_WITH_TRILINOS - -# include - -# include - -DEAL_II_DISABLE_EXTRA_DIAGNOSTICS -# include -# include - -# include -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(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 - VectorBase & - VectorBase::operator = (const ::dealii::Vector &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 - local_range = this->local_range (); - for (size_type i=local_range.first; i(&(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(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(&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(val[i]) << ' '; - else - for (size_type i=0; i(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