*/
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.
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
*/
void subtract_set (const IndexSet &other);
-
/**
* Fills the given vector with all indices contained in this IndexSet.
*/
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__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
precondition_block.cc
precondition_block_ez.cc
relaxation_block.cc
+ read_write_vector.cc
slepc_solver.cc
slepc_spectral_transformation.cc
solver.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+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);
+ \}
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// 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();
+}