From: Wolfgang Bangerth Date: Wed, 12 Aug 2015 20:53:43 +0000 (-0500) Subject: Add ReadWriteVector class. Part of the new vector interface. X-Git-Tag: v8.4.0-rc2~379^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8389d5574cf0fcaa3ed2617946c11c4cd879fc28;p=dealii.git Add ReadWriteVector class. Part of the new vector interface. --- diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h index 0eca56231e..a82ae2ad9f 100644 --- a/include/deal.II/base/index_set.h +++ b/include/deal.II/base/index_set.h @@ -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::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 index 0000000000..f09e0112c5 --- /dev/null +++ b/include/deal.II/lac/read_write_vector.h @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +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 + class ReadWriteVector : public Subscriptor + { + public: + /** + * Declare standard types used in all containers. These types parallel + * those in the C++ standard libraries vector<...> + * 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::real_type real_type; + + /** + * @name 1: Basic Object-handling + */ + //@{ + /** + * Empty constructor. + */ + ReadWriteVector (); + + /** + * Copy constructor. + */ + ReadWriteVector (const ReadWriteVector &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 + void reinit(const ReadWriteVector &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 + * swap(u,v) that simply calls u.swap(v), again in + * analogy to standard functions. + */ + void swap (ReadWriteVector &v); + + /** + * Copies the data and the IndexSet of the input vector @p in_vector. + */ + ReadWriteVector & + operator= (const ReadWriteVector &in_vector); + + /** + * Copies the data and the IndexSet of the input vector @p in_vector. + */ + template + ReadWriteVector & + operator= (const ReadWriteVector &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 & + 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 & + 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 &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 vector<> class of + * the C++ standard library by returning iterators to the start and end + * of the locally stored 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 + 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; + + /** + * 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 + void add (const std::vector &indices, + const std::vector &values); + + /** + * This is a second collective add operation. As a difference, this + * function takes a ReadWriteVector of values. + */ + template + void add (const std::vector &indices, + const ReadWriteVector &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. + */ + template + 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(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 + inline + ReadWriteVector::ReadWriteVector () + : + val (nullptr) + {} + + + + template + inline + ReadWriteVector::ReadWriteVector (const ReadWriteVector &v) + : + Subscriptor(), + val (nullptr) + { + this->operator=(v); + } + + + + template + inline + ReadWriteVector::ReadWriteVector (const size_type size) + : + val (nullptr) + { + reinit (size, false); + } + + + + template + inline + ReadWriteVector::ReadWriteVector (const IndexSet &locally_stored_indices) + : + val (nullptr) + { + reinit (locally_stored_indices); + } + + + + template + inline + ReadWriteVector::~ReadWriteVector () + { + resize_val(0); + } + + + + template + inline + ReadWriteVector & + ReadWriteVector::operator= (const ReadWriteVector &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 + template + inline + ReadWriteVector & + ReadWriteVector::operator= (const ReadWriteVector &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 + inline + typename ReadWriteVector::size_type + ReadWriteVector::size () const + { + return stored_elements.n_elements(); + } + + + + template + inline + const IndexSet & + ReadWriteVector::get_stored_elements () const + { + return stored_elements; + } + + + + template + inline + typename ReadWriteVector::iterator + ReadWriteVector::begin () + { + return &val[0]; + } + + + + template + inline + typename ReadWriteVector::const_iterator + ReadWriteVector::begin () const + { + return &val[0]; + } + + + + template + inline + typename ReadWriteVector::iterator + ReadWriteVector::end () + { + return &val[this->size()]; + } + + + + template + inline + typename ReadWriteVector::const_iterator + ReadWriteVector::end () const + { + return &val[this->size()]; + } + + + + template + inline + Number + ReadWriteVector::operator() (const size_type global_index) const + { + return val[global_to_local(global_index)]; + } + + + + template + inline + Number & + ReadWriteVector::operator() (const size_type global_index) + { + return val[global_to_local (global_index)]; + } + + + + template + inline + Number + ReadWriteVector::operator[] (const size_type global_index) const + { + return operator()(global_index); + } + + + + template + inline + Number & + ReadWriteVector::operator[] (const size_type global_index) + { + return operator()(global_index); + } + + + + template + template + inline + void ReadWriteVector::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 + template + inline + void ReadWriteVector::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 + inline + Number + ReadWriteVector::local_element (const size_type local_index) const + { + AssertIndexRange (local_index, this->size()); + + return val[local_index]; + } + + + + template + inline + Number & + ReadWriteVector::local_element (const size_type local_index) + { + AssertIndexRange (local_index, this->size()); + + return val[local_index]; + } + + + + template + inline + ReadWriteVector & + ReadWriteVector::operator= (const Number s) + { + Assert(s==static_cast(0), ExcMessage("Only 0 can be assigned to a vector.")); + std::fill(begin(),end(),s); + + return *this; + } + + + + template + template + inline + void + ReadWriteVector::add (const std::vector &indices, + const std::vector &values) + { + AssertDimension (indices.size(), values.size()); + add (indices.size(), &indices[0], &values[0]); + } + + + + template + template + inline + void + ReadWriteVector::add (const std::vector &indices, + const ReadWriteVector &values) + { + for (size_type i=0; ioperator()(indices[i]) += values[indices[i]]; + } + } + + + + template + template + inline + void + ReadWriteVector::add (const size_type n_indices, + const size_type *indices, + const Number2 *values) + { + for (size_type i=0; ioperator()(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 +inline +void swap (LinearAlgebra::ReadWriteVector &u, + LinearAlgebra::ReadWriteVector &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 index 0000000000..0cc0d90d02 --- /dev/null +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -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 +#include +#include +#include + +#ifdef DEAL_II_WITH_TRILINOS +# include +#endif + +DEAL_II_NAMESPACE_OPEN + + +namespace LinearAlgebra +{ + template + void + ReadWriteVector::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 + void + ReadWriteVector::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 + template + void + ReadWriteVector::reinit (const ReadWriteVector &v, + const bool fast) + { + resize_val (v.size()); + + stored_elements = v.get_stored_elements(); + + if (fast == false) + this->operator= (Number()); + } + + + + template + void + ReadWriteVector::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 + 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 + void copy_petsc_vector (const std::complex *petsc_start_ptr, + const std::complex *petsc_end_ptr, + std::complex *ptr) + { + std::copy(petsc_start_ptr, petsc_end_ptr, ptr); + } + + template + void copy_petsc_vector (const std::complex *petsc_start_ptr, + const std::complex *petsc_end_ptr, + Number *ptr) + { + AssertThrow(false, ExcMessage("Tried to copy complex -> real")); + } + } + + + + template + ReadWriteVector & + ReadWriteVector::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(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(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 + ReadWriteVector & + ReadWriteVector::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec) + { + // Copy the local elements + std::vector trilinos_indices, readwrite_indices; + trilinos_vec.locally_owned_elements().fill_index_vector(trilinos_indices); + stored_elements.fill_index_vector(readwrite_indices); + std::vector intersection(trilinos_indices.size()); + std::vector::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 trilinos_position(intersection_size); + unsigned int j=0; + for (size_t i=0; i 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[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 + void + ReadWriteVector::swap (ReadWriteVector &v) + { + std::swap(stored_elements, v.stored_elements); + std::swap(val, v.val); + } + + + + template + std::size_t + ReadWriteVector::memory_consumption () const + { + std::size_t memory = sizeof(*this); + memory += sizeof (Number) * static_cast(this->size()); + + memory += stored_elements.memory_consumption(); + + return memory; + } + + + + template + void + ReadWriteVector::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; isize(); ++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 diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index b0726f5a62..6b06ee08cd 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -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 index 0000000000..d92a39f040 --- /dev/null +++ b/source/lac/read_write_vector.cc @@ -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 +#include + +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 index 0000000000..eceace6fd4 --- /dev/null +++ b/source/lac/read_write_vector.inst.in @@ -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; + \} +} + +for (S1, S2 : REAL_SCALARS) +{ + namespace LinearAlgebra + \{ + template void ReadWriteVector::reinit (const ReadWriteVector&, + const bool); + \} +} + diff --git a/tests/lac/readwritevector_add.cc b/tests/lac/readwritevector_add.cc new file mode 100644 index 0000000000..6aa43beeff --- /dev/null +++ b/tests/lac/readwritevector_add.cc @@ -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 +#include +#include + +#include +#include + +void test() +{ + unsigned int size(10); + std::vector indices(size); + std::vector float_values(size); + std::vector double_values(size); + for (unsigned int i=0; i float_vector; + LinearAlgebra::ReadWriteVector double_vector_1; + LinearAlgebra::ReadWriteVector 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::iterator val_1 = double_vector_1.begin(); + LinearAlgebra::ReadWriteVector::iterator end_1 = double_vector_1.end(); + LinearAlgebra::ReadWriteVector::iterator val_2 = double_vector_2.begin(); + LinearAlgebra::ReadWriteVector::iterator end_2 = double_vector_2.end(); + deallog << "double_vector_1" << std::endl; + for (; val_1 +#include +#include + +#include + + +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_vector(is); + LinearAlgebra::ReadWriteVector float_vector(float_size); + deallog << "double_size " << double_vector.size() < double_vector2(double_vector); + double_vector2.print(deallog.get_file_stream()); + + for (unsigned int i=0; i +#include +#include +#include +#include + +#include +#include + +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 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 readwrite(readwrite_is); + readwrite = tril_vector; + if (rank==0) + { + std::vector 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 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 index 0000000000..8b13789179 --- /dev/null +++ b/tests/trilinos/readwritevector.mpirun=2.output @@ -0,0 +1 @@ +