]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add ReadWriteVector class. Part of the new vector interface.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 12 Aug 2015 20:53:43 +0000 (15:53 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Thu, 10 Sep 2015 18:56:28 +0000 (13:56 -0500)
12 files changed:
include/deal.II/base/index_set.h
include/deal.II/lac/read_write_vector.h [new file with mode: 0644]
include/deal.II/lac/read_write_vector.templates.h [new file with mode: 0644]
source/lac/CMakeLists.txt
source/lac/read_write_vector.cc [new file with mode: 0644]
source/lac/read_write_vector.inst.in [new file with mode: 0644]
tests/lac/readwritevector_add.cc [new file with mode: 0644]
tests/lac/readwritevector_add.output [new file with mode: 0644]
tests/lac/readwritevector_assignement.cc [new file with mode: 0644]
tests/lac/readwritevector_assignement.output [new file with mode: 0644]
tests/trilinos/readwritevector.cc [new file with mode: 0644]
tests/trilinos/readwritevector.mpirun=2.output [new file with mode: 0644]

index 0eca56231e5abf4477de2e8c4d48f850f25e4358..a82ae2ad9ffbd5549aa7abe6f60008e237bec63d 100644 (file)
@@ -222,6 +222,11 @@ public:
    */
   unsigned int n_intervals () const;
 
+  /**
+   * This function returns the local index of the beginning of the largest range.
+   */
+  unsigned int largest_range_index() const;
+
   /**
    * Compress the internal representation by merging individual elements with
    * contiguous ranges, etc. This function does not have any external effect.
@@ -266,7 +271,6 @@ public:
   IndexSet get_view (const size_type begin,
                      const size_type end) const;
 
-
   /**
    * Removes all elements contained in @p other from this set. In other words,
    * if $x$ is the current object and $o$ the argument, then we compute $x
@@ -274,7 +278,6 @@ public:
    */
   void subtract_set (const IndexSet &other);
 
-
   /**
    * Fills the given vector with all indices contained in this IndexSet.
    */
@@ -1409,6 +1412,18 @@ IndexSet::n_intervals () const
 
 
 
+inline
+unsigned int
+IndexSet::largest_range_index() const
+{
+  compress();
+  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
+
+  return main_range->nth_index_in_set;
+}
+
+
+
 inline
 IndexSet::size_type
 IndexSet::nth_index_in_set (const unsigned int n) const
diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h
new file mode 100644 (file)
index 0000000..f09e011
--- /dev/null
@@ -0,0 +1,715 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__read_write_vector_h
+#define dealii__read_write_vector_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/types.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/lac/vector_view.h>
+
+#include <cstring>
+#include <iomanip>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+#ifdef DEAL_II_WITH_PETSC
+namespace PETScWrappers
+{
+  namespace MPI
+  {
+    class Vector;
+  }
+}
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+namespace TrilinosWrappers
+{
+  namespace MPI
+  {
+    class Vector;
+  }
+}
+#endif
+
+namespace LinearAlgebra
+{
+  /*! @addtogroup Vectors
+   *@{
+   */
+
+  /**
+   * ReadWriteVector stores a subset of elements. It allows to access
+   * individual elements to be read or written. However, it does not allow global
+   * operations such as taking the norm. ReadWriteVector can be used to read and
+   * write elements in vector derived from VectorSpaceVector such as
+   * TrilinosWrappers::MPI::Vector and PETScWrappers::MPI::Vector.
+   *
+   * @author Bruno Turcksin, 2015.
+   */
+  template <typename Number>
+  class ReadWriteVector : public Subscriptor
+  {
+  public:
+    /**
+     * Declare standard types used in all containers. These types parallel
+     * those in the <tt>C++</tt> standard libraries <tt>vector<...></tt>
+     * class.
+     */
+    typedef Number                                            value_type;
+    typedef value_type                                       *pointer;
+    typedef const value_type                                 *const_pointer;
+    typedef value_type                                       *iterator;
+    typedef const value_type                                 *const_iterator;
+    typedef value_type                                       &reference;
+    typedef const value_type                                 &const_reference;
+    typedef types::global_dof_index                           size_type;
+    typedef typename numbers::NumberTraits<Number>::real_type real_type;
+
+    /**
+     * @name 1: Basic Object-handling
+     */
+    //@{
+    /**
+     * Empty constructor.
+     */
+    ReadWriteVector ();
+
+    /**
+     * Copy constructor.
+     */
+    ReadWriteVector (const ReadWriteVector<Number> &in_vector);
+
+    /**
+     * Constructs a vector given the size, the stored elements have their index
+     * in [0,size).
+     */
+    ReadWriteVector (const size_type size);
+
+    /**
+     * Constructs a vector whose stored elements indices are given by the
+     * IndexSet @p locally_stored_indices.
+     */
+    ReadWriteVector (const IndexSet &locally_stored_indices);
+
+    /**
+     * Destructor.
+     */
+    ~ReadWriteVector ();
+
+    /**
+     * Sets the global size of the vector to @p size. The stored elements have
+     * their index in [0,size).
+     *
+     * If the flag @p fast is set to false, the memory will be initialized
+     * with zero, otherwise the memory will be untouched (and the user must
+     * make sure to fill it with reasonable data before using it).
+     */
+    void reinit (const size_type size,
+                 const bool      fast = false);
+
+    /**
+     * Uses the same IndexSet as the one of the input vector @p in_vector and
+     * allocates memory for this vector.
+     *
+     * If the flag @p fast is set to false, the memory will be initialized
+     * with zero, otherwise the memory will be untouched (and the user must
+     * make sure to fill it with reasonable data before using it).
+     */
+    template <typename Number2>
+    void reinit(const ReadWriteVector<Number2> &in_vector,
+                const bool                      fast = false);
+
+    /**
+     * Initializes the vector. The indices are specified by @p
+     *
+     * If the flag @p fast is set to false, the memory will be initialized
+     * with zero, otherwise the memory will be untouched (and the user must
+     * make sure to fill it with reasonable data before using it).
+     * locally_stored_indices.
+     */
+    void reinit (const IndexSet &locally_stored_indices,
+                 const bool      fast = false);
+
+    /**
+     * 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.
+     *
+     * This function is analog to the the @p swap function of all C++
+     * standard containers. Also, there is a global function
+     * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+     * analogy to standard functions.
+     */
+    void swap (ReadWriteVector<Number> &v);
+
+    /**
+     * Copies the data and the IndexSet of the input vector @p in_vector.
+     */
+    ReadWriteVector<Number> &
+    operator= (const ReadWriteVector<Number> &in_vector);
+
+    /**
+     * Copies the data and the IndexSet of the input vector @p in_vector.
+     */
+    template <typename Number2>
+    ReadWriteVector<Number> &
+    operator= (const ReadWriteVector<Number2> &in_vector);
+
+#ifdef DEAL_II_WITH_PETSC
+    /**
+     * Imports all the elements present in the vector's IndexSet from the input
+     * vector @p petsc_vec.
+     */
+    ReadWriteVector<Number> &
+    operator= (const PETScWrappers::MPI::Vector &petsc_vec);
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+    /**
+     * Imports all the elements present in the vector's IndexSet from the input
+     * vector @p trilinos_vec.
+     */
+    ReadWriteVector<Number> &
+    operator= (const TrilinosWrappers::MPI::Vector &trilinos_vec);
+#endif
+
+    /**
+     * Sets all elements of the vector to the scalar @p s. This is operation is
+     * allowed if @p is equal to zero.
+     */
+    ReadWriteVector<Number> &operator = (const Number s);
+
+    /**
+     * Returns the size of the vector, equal to number of elements in
+     * the vector.
+     */
+    size_type size () const;
+
+    /**
+     * Return the IndexSet that stores the indices of the elements stored.
+     */
+    const IndexSet &get_stored_elements () const;
+
+    /**
+     * Make the @p ReadWriteVector class a bit like the <tt>vector<></tt> class of
+     * the C++ standard library by returning iterators to the start and end
+     * of the <i>locally stored</i> elements of this vector.
+     *
+     * It holds that end() - begin() == size().
+     */
+    iterator begin ();
+
+    /**
+     * Returns constant iterator to the start of the locally stored elements
+     * of the vector.
+     */
+    const_iterator begin () const;
+
+    /**
+     * Returns an iterator pointing to the element past the end of the array
+     * of locally stored entries.
+     */
+    iterator end ();
+
+    /**
+     * Returns a constant iterator pointing to the element past the end of
+     * the array of the locally stored entries.
+     */
+    const_iterator end () const;
+    //@}
+
+
+    /**
+     * @name 2: Data-Access
+     */
+    //@{
+
+    /**
+     * Read access to the data in the position corresponding to @p
+     * global_index.
+     */
+    Number operator () (const size_type global_index) const;
+
+    /**
+     * Read and write access to the data in the position corresponding to @p
+     * global_index.
+     */
+    Number &operator () (const size_type global_index);
+
+    /**
+     * Read access to the data in the position corresponding to @p
+     * global_index.
+     *
+     * This function does the same thing as operator().
+     */
+    Number operator [] (const size_type global_index) const;
+
+    /**
+     * Read and write access to the data in the position corresponding to @p
+     * global_index.
+     *
+     * This function does the same thing as operator().
+     */
+    Number &operator [] (const size_type global_index);
+
+    /**
+     * 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.
+     */
+    template <typename Number2>
+    void extract_subvector_to (const std::vector<size_type> &indices,
+                               std::vector<Number2> &values) const;
+
+    /**
+     * Just as the above, but with pointers. Useful in minimizing copying of
+     * data around.
+     */
+    template <typename ForwardIterator, typename OutputIterator>
+    void extract_subvector_to (ForwardIterator          indices_begin,
+                               const ForwardIterator    indices_end,
+                               OutputIterator           values_begin) const;
+
+    /**
+     * Read access to the data field specified by @p local_index. The largest
+     * range of the IndexSet is the first one.
+     *
+     * Performance: Direct array access (fast).
+     */
+    Number local_element (const size_type local_index) const;
+
+    /**
+     * Read and write access to the data field specified by @p local_index. The
+     * largest range of the IndexSet is the first one.
+     *
+     * Performance: Direct array access (fast).
+     */
+    Number &local_element (const size_type local_index);
+    //@}
+
+
+    /**
+     * @name 3: Modification of vectors
+     */
+    //@{
+
+    /**
+     * A collective add operation: This function adds a whole set of values
+     * stored in @p values to the vector components specified by @p indices.
+     */
+    template <typename Number2>
+    void add (const std::vector<size_type>  &indices,
+              const std::vector<Number2> &values);
+
+    /**
+     * This is a second collective add operation. As a difference, this
+     * function takes a ReadWriteVector of values.
+     */
+    template <typename Number2>
+    void add (const std::vector<size_type> &indices,
+              const ReadWriteVector<Number2> &values);
+
+    /**
+     * Take an address where <tt>n_elements</tt> are stored contiguously and
+     * add them into the vector. Handles all cases which are not covered by
+     * the other two <tt>add()</tt> functions above.
+     */
+    template <typename Number2>
+    void add (const size_type  n_elements,
+              const size_type *indices,
+              const Number2   *values);
+
+    /**
+     * Prints the vector to the output stream @p out.
+     */
+    void print (std::ostream       &out,
+                const unsigned int  precision  = 3,
+                const bool          scientific = true) const;
+
+    /**
+     * Returns the memory consumption of this class in bytes.
+     */
+    std::size_t memory_consumption () const;
+    //@}
+
+  protected:
+
+    /**
+     * Return the local position of @p global_index.
+     */
+    unsigned int
+    global_to_local (const types::global_dof_index global_index) const
+    {
+      // the following will throw an exception if the global_index is not
+      // in the remaining_elements
+      return static_cast<unsigned int>(stored_elements.index_within_set(global_index));
+    }
+
+    /**
+     * A helper function that is used to resize the val array.
+     */
+    void resize_val (const size_type new_allocated_size);
+
+    /**
+     * Indices of the elements stored.
+     */
+    IndexSet stored_elements;
+
+    /**
+     * Pointer to the array of local elements of this vector.
+     */
+    // TODO: Could we use VectorizedArray here for storage?
+    Number *val;
+  };
+
+  /*@}*/
+
+
+  /*----------------------- Inline functions ----------------------------------*/
+
+#ifndef DOXYGEN
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number>::ReadWriteVector ()
+    :
+    val (nullptr)
+  {}
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
+    :
+    Subscriptor(),
+    val (nullptr)
+  {
+    this->operator=(v);
+  }
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number>::ReadWriteVector (const size_type size)
+    :
+    val (nullptr)
+  {
+    reinit (size, false);
+  }
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
+    :
+    val (nullptr)
+  {
+    reinit (locally_stored_indices);
+  }
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number>::~ReadWriteVector ()
+  {
+    resize_val(0);
+  }
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number> &
+  ReadWriteVector<Number>::operator= (const ReadWriteVector<Number> &in_vector)
+  {
+    resize_val(in_vector.size());
+    stored_elements = in_vector.get_stored_elements();
+    std::copy(in_vector.begin(),in_vector.end(),begin());
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  ReadWriteVector<Number> &
+  ReadWriteVector<Number>::operator= (const ReadWriteVector<Number2> &in_vector)
+  {
+    resize_val(in_vector.size());
+    stored_elements = in_vector.get_stored_elements();
+    std::copy(in_vector.begin(),in_vector.end(),begin());
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::size_type
+  ReadWriteVector<Number>::size () const
+  {
+    return stored_elements.n_elements();
+  }
+
+
+
+  template <typename Number>
+  inline
+  const IndexSet &
+  ReadWriteVector<Number>::get_stored_elements () const
+  {
+    return stored_elements;
+  }
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::iterator
+  ReadWriteVector<Number>::begin ()
+  {
+    return &val[0];
+  }
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::const_iterator
+  ReadWriteVector<Number>::begin () const
+  {
+    return &val[0];
+  }
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::iterator
+  ReadWriteVector<Number>::end ()
+  {
+    return &val[this->size()];
+  }
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::const_iterator
+  ReadWriteVector<Number>::end () const
+  {
+    return &val[this->size()];
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number
+  ReadWriteVector<Number>::operator() (const size_type global_index) const
+  {
+    return val[global_to_local(global_index)];
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number &
+  ReadWriteVector<Number>::operator() (const size_type global_index)
+  {
+    return val[global_to_local (global_index)];
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number
+  ReadWriteVector<Number>::operator[] (const size_type global_index) const
+  {
+    return operator()(global_index);
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number &
+  ReadWriteVector<Number>::operator[] (const size_type global_index)
+  {
+    return operator()(global_index);
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  void ReadWriteVector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+                                                      std::vector<Number2> &values) const
+  {
+    for (size_type i = 0; i < indices.size(); ++i)
+      values[i] = operator()(indices[i]);
+  }
+
+
+
+  template <typename Number>
+  template <typename ForwardIterator, typename OutputIterator>
+  inline
+  void ReadWriteVector<Number>::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++;
+      }
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number
+  ReadWriteVector<Number>::local_element (const size_type local_index) const
+  {
+    AssertIndexRange (local_index, this->size());
+
+    return val[local_index];
+  }
+
+
+
+  template <typename Number>
+  inline
+  Number &
+  ReadWriteVector<Number>::local_element (const size_type local_index)
+  {
+    AssertIndexRange (local_index, this->size());
+
+    return val[local_index];
+  }
+
+
+
+  template <typename Number>
+  inline
+  ReadWriteVector<Number> &
+  ReadWriteVector<Number>::operator= (const Number s)
+  {
+    Assert(s==static_cast<Number>(0), ExcMessage("Only 0 can be assigned to a vector."));
+    std::fill(begin(),end(),s);
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  void
+  ReadWriteVector<Number>::add (const std::vector<size_type> &indices,
+                                const std::vector<Number2>   &values)
+  {
+    AssertDimension (indices.size(), values.size());
+    add (indices.size(), &indices[0], &values[0]);
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  void
+  ReadWriteVector<Number>::add (const std::vector<size_type>   &indices,
+                                const ReadWriteVector<Number2> &values)
+  {
+    for (size_type i=0; i<indices.size(); ++i)
+      {
+        Assert (numbers::is_finite(values[i]),
+                ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+        this->operator()(indices[i]) += values[indices[i]];
+      }
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  void
+  ReadWriteVector<Number>::add (const size_type    n_indices,
+                                const size_type   *indices,
+                                const Number2     *values)
+  {
+    for (size_type i=0; i<n_indices; ++i)
+      {
+        Assert (numbers::is_finite(values[i]),
+                ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+        this->operator()(indices[i]) += values[i];
+      }
+  }
+
+#endif  // ifndef DOXYGEN
+
+} // end of namespace LinearAlgebra
+
+
+
+
+/**
+ * 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 Vector
+ * @author Katharina Kormann, Martin Kronbichler, 2011
+ */
+template <typename Number>
+inline
+void swap (LinearAlgebra::ReadWriteVector<Number> &u,
+           LinearAlgebra::ReadWriteVector<Number> &v)
+{
+  u.swap (v);
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h
new file mode 100644 (file)
index 0000000..0cc0d90
--- /dev/null
@@ -0,0 +1,286 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__parallel_vector_templates_h
+#define dealii__parallel_vector_templates_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+#  include <Epetra_Import.h>
+#endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace LinearAlgebra
+{
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::resize_val (const size_type new_alloc_size)
+  {
+    if (new_alloc_size == 0)
+      {
+        if (val != nullptr)
+          free(val);
+        val = nullptr;
+      }
+    else
+      {
+        if (val != nullptr)
+          free(val);
+
+        Utilities::System::posix_memalign ((void **)&val, 64, sizeof(Number)*new_alloc_size);
+      }
+  }
+
+
+
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::reinit (const size_type size,
+                                   const bool      fast)
+  {
+    // check whether we need to reallocate
+    resize_val (size);
+
+    stored_elements = complete_index_set(size);
+
+    // set entries to zero if so requested
+    if (fast == false)
+      this->operator = (Number());
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  void
+  ReadWriteVector<Number>::reinit (const ReadWriteVector<Number2> &v,
+                                   const bool             fast)
+  {
+    resize_val (v.size());
+
+    stored_elements = v.get_stored_elements();
+
+    if (fast == false)
+      this->operator= (Number());
+  }
+
+
+
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::reinit (const IndexSet &locally_stored_indices,
+                                   const bool      fast)
+  {
+    stored_elements = locally_stored_indices;
+
+    // set vector size and allocate memory
+    resize_val(stored_elements.n_elements());
+
+    // initialize to zero
+    if (fast == false)
+      this->operator= (Number());
+  }
+
+
+
+#ifdef DEAL_II_WITH_PETSC
+  namespace internal
+  {
+    template <typename PETSC_Number, typename Number>
+    void copy_petsc_vector (const PETSC_Number *petsc_start_ptr,
+                            const PETSC_Number *petsc_end_ptr,
+                            Number *ptr)
+    {
+      std::copy(petsc_start_ptr, petsc_end_ptr, ptr);
+    }
+
+    template <typename PETSC_Number, typename Number>
+    void copy_petsc_vector (const std::complex<PETSC_Number> *petsc_start_ptr,
+                            const std::complex<PETSC_Number> *petsc_end_ptr,
+                            std::complex<Number> *ptr)
+    {
+      std::copy(petsc_start_ptr, petsc_end_ptr, ptr);
+    }
+
+    template <typename PETSC_Number, typename Number>
+    void copy_petsc_vector (const std::complex<PETSC_Number> *petsc_start_ptr,
+                            const std::complex<PETSC_Number> *petsc_end_ptr,
+                            Number *ptr)
+    {
+      AssertThrow(false, ExcMessage("Tried to copy complex -> real"));
+    }
+  }
+
+
+
+  template <typename Number>
+  ReadWriteVector<Number> &
+  ReadWriteVector<Number>::operator = (const PETScWrappers::MPI::Vector &petsc_vec)
+  {
+    //TODO: this works only if no communication is needed.
+    Assert(petsc_vec.locally_owned_elements() == stored_elements,
+           StandardExceptions::ExcInvalidState());
+
+    // get a representation of the vector and copy it
+    PetscScalar *start_ptr;
+    int ierr = VecGetArray (static_cast<const Vec &>(petsc_vec), &start_ptr);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    const size_type vec_size = local_size();
+    internal::copy_petsc_vector (start_ptr, start_ptr + vec_size, begin());
+
+    // restore the representation of the vector
+    ierr = VecRestoreArray (static_cast<const Vec &>(petsc_vec), &start_ptr);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    // return a pointer to this object per normal c++ operator overloading
+    // semantics
+    return *this;
+  }
+#endif
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template <typename Number>
+  ReadWriteVector<Number> &
+  ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+  {
+    // Copy the local elements
+    std::vector<size_type> trilinos_indices, readwrite_indices;
+    trilinos_vec.locally_owned_elements().fill_index_vector(trilinos_indices);
+    stored_elements.fill_index_vector(readwrite_indices);
+    std::vector<size_type> intersection(trilinos_indices.size());
+    std::vector<size_type>::iterator end;
+    end = std::set_intersection(trilinos_indices.begin(), trilinos_indices.end(),
+                                readwrite_indices.begin(), readwrite_indices.end(),
+                                intersection.begin());
+    const size_t intersection_size = end-intersection.begin();
+    intersection.resize(intersection_size);
+    std::vector<size_t> trilinos_position(intersection_size);
+    unsigned int j=0;
+    for (size_t i=0; i<intersection_size; ++i)
+      {
+        while (intersection[i]!=trilinos_indices[j])
+          ++j;
+        trilinos_position[i] = j;
+      }
+    double *values = trilinos_vec.trilinos_vector().Values();
+    for (size_t i=0; i<trilinos_position.size(); ++i)
+      val[global_to_local(intersection[i])] = values[trilinos_position[i]];
+
+
+#ifdef DEAL_II_WITH_MPI
+    // Copy the off-processor elements if necessary
+    if (intersection_size != size())
+      {
+        Epetra_BlockMap source_map = trilinos_vec.trilinos_vector().Map();
+        std::vector<size_type> off_proc_indices(readwrite_indices.size());
+        // Subtract the elements of readwrite that are in intersection.
+        end = std::set_difference(readwrite_indices.begin(), readwrite_indices.end(),
+                                  intersection.begin(), intersection.end(), off_proc_indices.begin());
+        off_proc_indices.resize(end-off_proc_indices.begin());
+        // Cast size_type to TrilinosWrappers::type::int_type
+        TrilinosWrappers::types::int_type *off_proc_array =
+          new TrilinosWrappers::types::int_type [off_proc_indices.size()];
+        for (size_t i=0; i<off_proc_indices.size(); ++i)
+          off_proc_array[i] = static_cast<TrilinosWrappers::types::int_type> (
+                                off_proc_indices[i]);
+        Epetra_Map target_map(off_proc_indices.size(),off_proc_indices.size(),
+                              off_proc_array,0,source_map.Comm());
+        delete [] off_proc_array;
+        Epetra_Import import(target_map, source_map);
+        Epetra_FEVector target_vector(target_map);
+        target_vector.Import(trilinos_vec.trilinos_vector(), import, Insert);
+        values = target_vector.Values();
+        for (unsigned int i=0; i<off_proc_indices.size(); ++i)
+          val[global_to_local(off_proc_indices[i])] = values[i];
+      }
+#endif
+
+    // return a pointer to this object per normal c++ operator overloading
+    // semantics
+    return *this;
+  }
+#endif
+
+
+
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::swap (ReadWriteVector<Number> &v)
+  {
+    std::swap(stored_elements, v.stored_elements);
+    std::swap(val, v.val);
+  }
+
+
+
+  template <typename Number>
+  std::size_t
+  ReadWriteVector<Number>::memory_consumption () const
+  {
+    std::size_t memory = sizeof(*this);
+    memory += sizeof (Number) * static_cast<std::size_t>(this->size());
+
+    memory += stored_elements.memory_consumption();
+
+    return memory;
+  }
+
+
+
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::print (std::ostream      &out,
+                                  const unsigned int precision,
+                                  const bool         scientific) const
+  {
+    AssertThrow (out, ExcIO());
+    std::ios::fmtflags old_flags = out.flags();
+    unsigned int old_precision = out.precision (precision);
+
+    out.precision (precision);
+    if (scientific)
+      out.setf (std::ios::scientific, std::ios::floatfield);
+    else
+      out.setf (std::ios::fixed, std::ios::floatfield);
+
+    out << "IndexSet: ";
+    stored_elements.print(out);
+    out << std::endl;
+    for (unsigned int i=0; i<this->size(); ++i)
+      out << val[i] << std::endl;
+    out << std::flush;
+
+
+    AssertThrow (out, ExcIO());
+    // reset output format
+    out.flags (old_flags);
+    out.precision(old_precision);
+  }
+} // end of namespace LinearAlgebra
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index b0726f5a622231234eabdfda9f97c609f7a39ad9..6b06ee08cd2d5ae5721a061588ad42e123a2e960 100644 (file)
@@ -46,6 +46,7 @@ SET(_src
   precondition_block.cc
   precondition_block_ez.cc
   relaxation_block.cc
+  read_write_vector.cc
   slepc_solver.cc
   slepc_spectral_transformation.cc
   solver.cc
@@ -88,6 +89,7 @@ SET(_inst
   parallel_vector.inst.in
   precondition_block.inst.in
   relaxation_block.inst.in
+  read_write_vector.inst.in
   solver.inst.in
   sparse_matrix_ez.inst.in
   sparse_matrix.inst.in
diff --git a/source/lac/read_write_vector.cc b/source/lac/read_write_vector.cc
new file mode 100644 (file)
index 0000000..d92a39f
--- /dev/null
@@ -0,0 +1,23 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/read_write_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#include "read_write_vector.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/lac/read_write_vector.inst.in b/source/lac/read_write_vector.inst.in
new file mode 100644 (file)
index 0000000..eceace6
--- /dev/null
@@ -0,0 +1,34 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (SCALAR : REAL_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    template class ReadWriteVector<SCALAR>;
+  \}
+}
+
+for (S1, S2 : REAL_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    template void ReadWriteVector<S1>::reinit<S2> (const ReadWriteVector<S2>&,
+                                          const bool);
+  \}
+}
+
diff --git a/tests/lac/readwritevector_add.cc b/tests/lac/readwritevector_add.cc
new file mode 100644 (file)
index 0000000..6aa43be
--- /dev/null
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check reinit and add
+
+#include "../tests.h"
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/read_write_vector.h>
+
+#include <fstream>
+#include <vector>
+
+void test()
+{
+  unsigned int size(10);
+  std::vector<types::global_dof_index> indices(size);
+  std::vector<float> float_values(size);
+  std::vector<double> double_values(size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      indices[i] = i;
+      float_values[i] = 2.*i;
+      double_values[i] = 3.*i;
+    }
+
+  LinearAlgebra::ReadWriteVector<float> float_vector;
+  LinearAlgebra::ReadWriteVector<double> double_vector_1;
+  LinearAlgebra::ReadWriteVector<double> double_vector_2;
+  IndexSet is(50);
+  is.add_range(0,10);
+  is.add_index(46);
+  is.add_range(11,25);
+  is.compress();
+
+  float_vector.reinit(size);
+  double_vector_1.reinit(float_vector);
+  double_vector_2.reinit(is);
+
+  IndexSet is_copy(double_vector_2.get_stored_elements());
+  is.print(std::cout);
+  is_copy.print(std::cout);
+
+  float_vector.add(indices,float_values);
+  double_vector_1.add(indices,float_vector);
+  double_vector_2.add(size,&indices[0],&double_values[0]);
+  double_vector_1.swap(double_vector_2);
+  LinearAlgebra::ReadWriteVector<double>::iterator val_1 = double_vector_1.begin();
+  LinearAlgebra::ReadWriteVector<double>::iterator end_1 = double_vector_1.end();
+  LinearAlgebra::ReadWriteVector<double>::iterator val_2 = double_vector_2.begin();
+  LinearAlgebra::ReadWriteVector<double>::iterator end_2 = double_vector_2.end();
+  deallog << "double_vector_1" << std::endl;
+  for (; val_1<end_1; ++val_1)
+    deallog << *val_1 <<std::endl;
+  deallog << "double_vector_2" << std::endl;
+  for (; val_2<end_2; ++val_2)
+    deallog << *val_2 <<std::endl;
+  double_vector_1.extract_subvector_to(indices,float_values);
+  deallog << "subvector" << std::endl;
+  for (unsigned int i=0; i<indices.size(); ++i)
+    deallog << float_values[i] << std::endl;
+}
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  test();
+
+  return 0;
+}
diff --git a/tests/lac/readwritevector_add.output b/tests/lac/readwritevector_add.output
new file mode 100644 (file)
index 0000000..e6155be
--- /dev/null
@@ -0,0 +1,49 @@
+
+DEAL::double_vector_1
+DEAL::0.00000
+DEAL::3.00000
+DEAL::6.00000
+DEAL::9.00000
+DEAL::12.0000
+DEAL::15.0000
+DEAL::18.0000
+DEAL::21.0000
+DEAL::24.0000
+DEAL::27.0000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::0.00000
+DEAL::double_vector_2
+DEAL::0.00000
+DEAL::2.00000
+DEAL::4.00000
+DEAL::6.00000
+DEAL::8.00000
+DEAL::10.0000
+DEAL::12.0000
+DEAL::14.0000
+DEAL::16.0000
+DEAL::18.0000
+DEAL::subvector
+DEAL::0.00000
+DEAL::3.00000
+DEAL::6.00000
+DEAL::9.00000
+DEAL::12.0000
+DEAL::15.0000
+DEAL::18.0000
+DEAL::21.0000
+DEAL::24.0000
+DEAL::27.0000
diff --git a/tests/lac/readwritevector_assignement.cc b/tests/lac/readwritevector_assignement.cc
new file mode 100644 (file)
index 0000000..893e021
--- /dev/null
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check constructor and operator=
+
+#include "../tests.h"
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/read_write_vector.h>
+
+#include <fstream>
+
+
+void test()
+{
+  unsigned int double_size = 2;
+  unsigned int float_size = 10;
+  IndexSet is(50);
+  is.add_range(0,2);
+  is.add_index(46);
+  is.add_range(10,15);
+  LinearAlgebra::ReadWriteVector<double> double_vector(is);
+  LinearAlgebra::ReadWriteVector<float> float_vector(float_size);
+  deallog << "double_size " << double_vector.size() <<std::endl;
+  deallog << "float_size " << float_vector.size() <<std::endl;
+
+  double_vector = 0.;
+  for (unsigned int i=0; i<double_vector.size(); ++i)
+    double_vector.local_element(i) += i;
+  for (unsigned int i=0; i<float_vector.size(); ++i)
+    float_vector[i] = i;
+
+  double_vector.print(deallog.get_file_stream());
+  float_vector.print(deallog.get_file_stream());
+
+  float_vector = double_vector;
+  float_vector.print(deallog.get_file_stream());
+
+  LinearAlgebra::ReadWriteVector<double> double_vector2(double_vector);
+  double_vector2.print(deallog.get_file_stream());
+
+  for (unsigned int i=0; i<double_vector.size(); ++i)
+    double_vector2.local_element(i) += i;
+  double_vector = double_vector2;
+  double_vector.print(deallog.get_file_stream());
+}
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  test();
+
+  return 0;
+}
diff --git a/tests/lac/readwritevector_assignement.output b/tests/lac/readwritevector_assignement.output
new file mode 100644 (file)
index 0000000..2e8a7ec
--- /dev/null
@@ -0,0 +1,55 @@
+
+DEAL::double_size 8
+DEAL::float_size 10
+IndexSet: {[0,1], [10,14], 46}
+
+0.000e+00
+1.000e+00
+2.000e+00
+3.000e+00
+4.000e+00
+5.000e+00
+6.000e+00
+7.000e+00
+IndexSet: {[0,9]}
+
+0.000e+00
+1.000e+00
+2.000e+00
+3.000e+00
+4.000e+00
+5.000e+00
+6.000e+00
+7.000e+00
+8.000e+00
+9.000e+00
+IndexSet: {[0,1], [10,14], 46}
+
+0.000e+00
+1.000e+00
+2.000e+00
+3.000e+00
+4.000e+00
+5.000e+00
+6.000e+00
+7.000e+00
+IndexSet: {[0,1], [10,14], 46}
+
+0.000e+00
+1.000e+00
+2.000e+00
+3.000e+00
+4.000e+00
+5.000e+00
+6.000e+00
+7.000e+00
+IndexSet: {[0,1], [10,14], 46}
+
+0.000e+00
+2.000e+00
+4.000e+00
+6.000e+00
+8.000e+00
+1.000e+01
+1.200e+01
+1.400e+01
diff --git a/tests/trilinos/readwritevector.cc b/tests/trilinos/readwritevector.cc
new file mode 100644 (file)
index 0000000..2a06fb3
--- /dev/null
@@ -0,0 +1,92 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// test operator =
+
+#include "../tests.h"
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <fstream>
+#include <vector>
+
+void test()
+{
+  IndexSet is(8);
+  unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  if (rank==0)
+    is.add_range(0,4);
+  if (rank==1)
+    is.add_range(4,8);
+  is.compress();
+  TrilinosWrappers::MPI::Vector tril_vector(is);
+  Vector<double> tmp(8);
+  for (unsigned int i=0; i<8; ++i)
+    tmp[i] = i;
+  tril_vector = tmp;
+
+  tril_vector.compress(VectorOperation::insert);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
+  IndexSet readwrite_is(8);
+  if (rank==0)
+    {
+      readwrite_is.add_range(0,2);
+      readwrite_is.add_range(6,8);
+    }
+  if (rank==1)
+    readwrite_is.add_range(2,6);
+  readwrite_is.compress();
+  LinearAlgebra::ReadWriteVector<double> readwrite(readwrite_is);
+  readwrite = tril_vector;
+  if (rank==0)
+    {
+      std::vector<double> comp(4);
+      comp[0] = 0.;
+      comp[1] = 1.;
+      comp[2] = 6.;
+      comp[3] = 7.;
+      for (unsigned int i=0; i<4; ++i)
+        AssertThrow(readwrite.local_element(i)==comp[i],
+                    ExcMessage("Element not copied correctly"));
+    }
+  MPI_Barrier(MPI_COMM_WORLD);
+  if (rank==1)
+    {
+      std::vector<double> comp(4);
+      comp[0] = 2.;
+      comp[1] = 3.;
+      comp[2] = 4.;
+      comp[3] = 5.;
+      for (unsigned int i=0; i<4; ++i)
+        AssertThrow(readwrite.local_element(i)==comp[i],
+                    ExcMessage("Element not copied correctly"));
+    }
+}
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
+  test();
+}
diff --git a/tests/trilinos/readwritevector.mpirun=2.output b/tests/trilinos/readwritevector.mpirun=2.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.