From: Bruno Turcksin Date: Thu, 2 Aug 2018 14:51:24 +0000 (+0000) Subject: Add support for CUDA memory space in LA::distributed::Vector X-Git-Tag: v9.1.0-rc1~789^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fc92fda2bc8177169dc36e049d25d1132a227838;p=dealii.git Add support for CUDA memory space in LA::distributed::Vector Add support for CUDA memory space to LinearAlgebra::distributed::Vector by adding a new template parameter with a default value of Host. --- diff --git a/include/deal.II/base/memory_space.h b/include/deal.II/base/memory_space.h new file mode 100644 index 0000000000..f67b50cf5c --- /dev/null +++ b/include/deal.II/base/memory_space.h @@ -0,0 +1,222 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_memory_space_h +#define dealii_memory_space_h + +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * + */ +namespace MemorySpace +{ + /** + * Structure describing Host memory space. + */ + struct Host + {}; + + + + /** + * Structure describing CUDA memory space. + */ + struct CUDA + {}; + + + +#ifdef DEAL_II_COMPILER_CUDA_AWARE + /** + * Deleter for unique/shared pointer pointing to an array allocated on the + * device. + */ + template + void + deleter_device_array(Number *ptr) noexcept + { + const cudaError_t cuda_error_code = cudaFree(ptr); + (void)cuda_error_code; + AssertNothrow(cuda_error_code == cudaSuccess, + ExcCudaError(cudaGetErrorString(cuda_error_code))); + } +#endif + + + + /** + * Data structure + */ + template + struct MemorySpaceData + { + MemorySpaceData() + { + static_assert(std::is_same::value || + std::is_same::value, + "MemorySpace should be Host or CUDA"); + } + + /** + * Copy the active data (values for Host and values_dev for CUDA) to @p begin. + * If the data is on the device it is moved to the host. + */ + void + copy_to(Number *begin, std::size_t n_elements) + { + (void)begin; + (void)n_elements; + } + + void + /** + * Copy the data in @p begin to the active data of the structure (values for + * Host and values_dev for CUDA). The pointer @p begin must be on the host. + */ + copy_from(Number *begin, std::size_t n_elements) + { + (void)begin; + (void)n_elements; + } + + /** + * Pointer to data on the host. + */ + std::unique_ptr values; + + /** + * Pointer to data on the device. + */ + std::unique_ptr values_dev; + }; + + + + /** + * Swap function similar to std::swap. + */ + template + inline void + swap(MemorySpaceData &, + MemorySpaceData &) + { + static_assert(std::is_same::value || + std::is_same::value, + "MemorySpace should be Host or CUDA"); + } + +#ifndef DOXYGEN + + template + struct MemorySpaceData + { + MemorySpaceData() + : values(nullptr, &free) + {} + + void + copy_to(Number *begin, std::size_t n_elements) + { + std::copy(values.get(), values.get() + n_elements, begin); + } + + void + copy_from(Number *begin, std::size_t n_elements) + { + std::copy(begin, begin + n_elements, values.get()); + } + + std::unique_ptr values; + + // This is not used but it allows to simplify the code until we start using + // CUDA-aware MPI. + std::unique_ptr values_dev; + }; + + + + template + inline void + swap(MemorySpaceData &u, MemorySpaceData &v) + { + std::swap(u.values, v.values); + } + + + +# ifdef DEAL_II_COMPILER_CUDA_AWARE + + template + struct MemorySpaceData + { + MemorySpaceData() + : values(nullptr, &free) + , values_dev(nullptr, deleter_device_array) + {} + + void + copy_to(Number *begin, std::size_t n_elements) + { + const cudaError_t cuda_error_code = + cudaMemcpy(begin, + values_dev.get(), + n_elements * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error_code); + } + + void + copy_from(Number *begin, std::size_t n_elements) + { + const cudaError_t cuda_error_code = + cudaMemcpy(values_dev.get(), + begin, + n_elements * sizeof(Number), + cudaMemcpyHostToDevice); + AssertCuda(cuda_error_code); + } + + std::unique_ptr values; + std::unique_ptr values_dev; + }; + + + + template + inline void + swap(MemorySpaceData &u, MemorySpaceData &v) + { + std::swap(u.values, v.values); + std::swap(u.values_dev, v.values_dev); + } + +# endif + +#endif + +} // namespace MemorySpace + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/base/types.h b/include/deal.II/base/types.h index a6724db133..5d94c7d802 100644 --- a/include/deal.II/base/types.h +++ b/include/deal.II/base/types.h @@ -276,7 +276,6 @@ namespace numbers static_cast(-2); } // namespace numbers - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/lac/la_parallel_block_vector.templates.h b/include/deal.II/lac/la_parallel_block_vector.templates.h index fecbee6055..c1fb7d0f17 100644 --- a/include/deal.II/lac/la_parallel_block_vector.templates.h +++ b/include/deal.II/lac/la_parallel_block_vector.templates.h @@ -573,7 +573,7 @@ namespace LinearAlgebra for (unsigned int i = 0; i < this->n_blocks(); ++i) local_result = std::max(local_result, - -static_cast(this->block(i).all_zero_local())); + (this->block(i).linfty_norm_local() == 0) ? -1 : 0); if (this->block(0).partitioner->n_mpi_processes() > 1) return -Utilities::MPI::max( diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h index 2467590a27..b7ea752085 100644 --- a/include/deal.II/lac/la_parallel_vector.h +++ b/include/deal.II/lac/la_parallel_vector.h @@ -18,6 +18,7 @@ #include +#include #include #include #include @@ -175,13 +176,54 @@ namespace LinearAlgebra * fail in some circumstances. Therefore, it is strongly recommended to * not rely on this class to automatically detect the unsupported case. * - * @author Katharina Kormann, Martin Kronbichler, 2010, 2011, 2016 + *

CUDA support

+ * + * This vector class supports two different memory spaces: Host and CUDA. By + * default, the memory space is Host and all the data are allocated on the + * CPU. When the memory space is CUDA, all the data is allocated on the GPU. + * The operations on the vector are performed on the chosen memory space. * + * From the host, there are two methods to access the elements of the Vector + * when using the CUDA memory space:
    + *
  • use get_values() + * + * Vector vector(local_range, comm); + * double* vector_dev = vector.get_values(); + * std::vector vector_host(local_range.n_elements(), 1.); + * Utilities::CUDA::copy_to_dev(vector_host, vector_dev); + * + *
  • use import() + * + * Vector vector(local_range, comm); + * ReadWriteVector rw_vector(local_range); + * for (auto & val : rw_vector) + * val = 1.; + * vector.import(rw_vector, VectorOperations::insert); + * + *
+ * The import method is a lot safer and will perform an MPI communication if + * necessary. Since an MPI communication may be performed, import needs to + * be called on all the processors. + * + * @note By default, all the ranks will try to access the device 0. This is + * fine is if you have one rank per node and one gpu per node. If you + * have multiple GPUs on one node, we need each process to access a + * different GPU. If each node has the same number of GPUs, this can be done + * as follows: + * int n_devices = 0; cudaGetDeviceCount(&n_devices); int + * device_id = my_rank % n_devices; + * cudaSetDevice(device_id); + * + * @see CUDAWrappers + * + * @author Katharina Kormann, Martin Kronbichler, Bruno Turcksin 2010, 2011, + * 2016, 2018 */ - template + template class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector, public Subscriptor { public: + using memory_space = MemorySpace; using value_type = Number; using pointer = value_type *; using const_pointer = const value_type *; @@ -192,6 +234,11 @@ namespace LinearAlgebra using size_type = types::global_dof_index; using real_type = typename numbers::NumberTraits::real_type; + static_assert( + std::is_same::value || + std::is_same::value, + "MemorySpace should be Host or CUDA"); + /** * @name 1: Basic Object-handling */ @@ -204,7 +251,7 @@ namespace LinearAlgebra /** * Copy constructor. Uses the parallel partitioning of @p in_vector. */ - Vector(const Vector &in_vector); + Vector(const Vector &in_vector); /** * Construct a parallel vector of the given global size without any @@ -270,8 +317,8 @@ namespace LinearAlgebra */ template void - reinit(const Vector &in_vector, - const bool omit_zeroing_entries = false); + reinit(const Vector &in_vector, + const bool omit_zeroing_entries = false); /** * Initialize the vector. The local range is specified by @p @@ -326,7 +373,7 @@ namespace LinearAlgebra * handle memory separately. */ void - swap(Vector &v); + swap(Vector &v); /** * Assigns the vector to the parallel partitioning of the input vector @@ -339,8 +386,8 @@ namespace LinearAlgebra * at all, the vector will also update its ghost values in analogy to * the respective setting the Trilinos and PETSc vectors. */ - Vector & - operator=(const Vector &in_vector); + Vector & + operator=(const Vector &in_vector); /** * Assigns the vector to the parallel partitioning of the input vector @@ -354,8 +401,8 @@ namespace LinearAlgebra * the respective setting the Trilinos and PETSc vectors. */ template - Vector & - operator=(const Vector &in_vector); + Vector & + operator=(const Vector &in_vector); #ifdef DEAL_II_WITH_PETSC /** @@ -369,7 +416,7 @@ namespace LinearAlgebra * ReadWriteVector instead. */ DEAL_II_DEPRECATED - Vector & + Vector & operator=(const PETScWrappers::MPI::Vector &petsc_vec); #endif @@ -386,7 +433,7 @@ namespace LinearAlgebra * ReadWriteVector instead. */ DEAL_II_DEPRECATED - Vector & + Vector & operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec); #endif //@} @@ -564,7 +611,7 @@ namespace LinearAlgebra */ template void - copy_locally_owned_data_from(const Vector &src); + copy_locally_owned_data_from(const Vector &src); //@} @@ -584,25 +631,25 @@ namespace LinearAlgebra /** * Multiply the entire vector by a fixed factor. */ - virtual Vector & + virtual Vector & operator*=(const Number factor) override; /** * Divide the entire vector by a fixed factor. */ - virtual Vector & + virtual Vector & operator/=(const Number factor) override; /** * Add the vector @p V to the present one. */ - virtual Vector & + virtual Vector & operator+=(const VectorSpaceVector &V) override; /** * Subtract the vector @p V from the present one. */ - virtual Vector & + virtual Vector & operator-=(const VectorSpaceVector &V) override; /** @@ -612,6 +659,9 @@ namespace LinearAlgebra * current elements. The last parameter can be used if the same * communication pattern is used multiple times. This can be used to * improve performance. + * + * @note If the MemorySpace is CUDA, the data in the ReadWriteVector will + * be moved to the device. */ virtual void import( @@ -777,7 +827,7 @@ namespace LinearAlgebra * zero, also ghost elements are set to zero, otherwise they remain * unchanged. */ - virtual Vector & + virtual Vector & operator=(const Number s) override; /** @@ -804,7 +854,7 @@ namespace LinearAlgebra * s*(*this)+V. */ void - sadd(const Number s, const Vector &V); + sadd(const Number s, const Vector &V); /** * Scaling and multiple addition. @@ -813,11 +863,11 @@ namespace LinearAlgebra */ DEAL_II_DEPRECATED void - sadd(const Number s, - const Number a, - const Vector &V, - const Number b, - const Vector &W); + sadd(const Number s, + const Number a, + const Vector &V, + const Number b, + const Vector &W); /** * Assignment *this = a*u + b*v. @@ -826,10 +876,10 @@ namespace LinearAlgebra */ DEAL_II_DEPRECATED void - equ(const Number a, - const Vector &u, - const Number b, - const Vector &v); + equ(const Number a, + const Vector &u, + const Number b, + const Vector &v); //@} @@ -904,6 +954,9 @@ namespace LinearAlgebra * of the locally owned elements of this vector. * * It holds that end() - begin() == local_size(). + * + * @note For the CUDA memory space, the iterator points to memory on the + * device. */ iterator begin(); @@ -911,6 +964,9 @@ namespace LinearAlgebra /** * Return constant iterator to the start of the locally owned elements * of the vector. + * + * @note For the CUDA memory space, the iterator points to memory on the + * device. */ const_iterator begin() const; @@ -918,6 +974,9 @@ namespace LinearAlgebra /** * Return an iterator pointing to the element past the end of the array * of locally owned entries. + * + * @note For the CUDA memory space, the iterator points to memory on the + * device. */ iterator end(); @@ -925,6 +984,9 @@ namespace LinearAlgebra /** * Return a constant iterator pointing to the element past the end of * the array of the locally owned entries. + * + * @note For the CUDA memory space, the iterator points to memory on the + * device. */ const_iterator end() const; @@ -992,6 +1054,15 @@ namespace LinearAlgebra Number & local_element(const size_type local_index); + /** + * Return the pointer to the underlying raw array. + * + * @note For the CUDA memory space, the pointer points to memory on the + * device. + */ + Number * + get_values() const; + /** * Instead of getting individual elements of a vector via operator(), * this function allows getting a whole set of elements at once. The @@ -1006,6 +1077,8 @@ namespace LinearAlgebra * @endcode * * @pre The sizes of the @p indices and @p values arrays must be identical. + * + * @note This function is not implemented for CUDA memory space. */ template void @@ -1125,6 +1198,13 @@ namespace LinearAlgebra */ DeclException0(ExcVectorTypeNotCompatible); + /** + * Attempt to perform an operation not implemented on the device. + * + * @ingroup Exceptions + */ + DeclException0(ExcNotAllowedForCuda); + /** * Exception */ @@ -1171,18 +1251,12 @@ namespace LinearAlgebra const Number a, const VectorSpaceVector &V); - /** - * Local part of all_zero(). - */ - bool - all_zero_local() const; - /** * Local part of the inner product of two vectors. */ template Number - inner_product_local(const Vector &V) const; + inner_product_local(const Vector &V) const; /** * Local part of norm_sqr(). @@ -1220,9 +1294,9 @@ namespace LinearAlgebra * the add_and_dot() function. */ Number - add_and_dot_local(const Number a, - const Vector &V, - const Vector &W); + add_and_dot_local(const Number a, + const Vector &V, + const Vector &W); /** * Shared pointer to store the parallel partitioning information. This @@ -1237,13 +1311,9 @@ namespace LinearAlgebra size_type allocated_size; /** - * Pointer to the array of local elements of this vector. - * - * Because we allocate these arrays via Utilities::System::posix_memalign, - * we need to use a custom deleter for this object that does not call - * delete[], but instead calls @p free(). + * Underlying data structure storing the local elements of this vector. */ - std::unique_ptr values; + mutable ::dealii::MemorySpace::MemorySpaceData data; /** * For parallel loops with TBB, this member variable stores the affinity @@ -1309,7 +1379,7 @@ namespace LinearAlgebra /* * Make all other vector types friends. */ - template + template friend class Vector; /** @@ -1325,56 +1395,142 @@ namespace LinearAlgebra #ifndef DOXYGEN + namespace internal + { + template + struct Policy + { + static inline typename Vector::iterator + begin(::dealii::MemorySpace::MemorySpaceData &) + { + return nullptr; + } + + static inline typename Vector::const_iterator + begin( + const ::dealii::MemorySpace::MemorySpaceData &) + { + return nullptr; + } - template + static inline Number * + get_values( + ::dealii::MemorySpace::MemorySpaceData &) + { + return nullptr; + } + }; + + + + template + struct Policy + { + static inline + typename Vector::iterator + begin(::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values.get(); + } + + static inline + typename Vector::const_iterator + begin(const ::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values.get(); + } + + static inline Number * + get_values(::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values.get(); + } + }; + + + + template + struct Policy + { + static inline + typename Vector::iterator + begin(::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values_dev.get(); + } + + static inline + typename Vector::const_iterator + begin(const ::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values_dev.get(); + } + + static inline Number * + get_values(::dealii::MemorySpace:: + MemorySpaceData &data) + { + return data.values_dev.get(); + } + }; + } // namespace internal + + + template inline bool - Vector::has_ghost_elements() const + Vector::has_ghost_elements() const { return vector_is_ghosted; } - template - inline typename Vector::size_type - Vector::size() const + template + inline typename Vector::size_type + Vector::size() const { return partitioner->size(); } - template - inline typename Vector::size_type - Vector::local_size() const + template + inline typename Vector::size_type + Vector::local_size() const { return partitioner->local_size(); } - template - inline std::pair::size_type, - typename Vector::size_type> - Vector::local_range() const + template + inline std::pair::size_type, + typename Vector::size_type> + Vector::local_range() const { return partitioner->local_range(); } - template + template inline bool - Vector::in_local_range(const size_type global_index) const + Vector::in_local_range( + const size_type global_index) const { return partitioner->in_local_range(global_index); } - template + template inline IndexSet - Vector::locally_owned_elements() const + Vector::locally_owned_elements() const { IndexSet is(size()); @@ -1386,73 +1542,79 @@ namespace LinearAlgebra - template - inline typename Vector::size_type - Vector::n_ghost_entries() const + template + inline typename Vector::size_type + Vector::n_ghost_entries() const { return partitioner->n_ghost_indices(); } - template + template inline const IndexSet & - Vector::ghost_elements() const + Vector::ghost_elements() const { return partitioner->ghost_indices(); } - template + template inline bool - Vector::is_ghost_entry(const size_type global_index) const + Vector::is_ghost_entry( + const size_type global_index) const { return partitioner->is_ghost_entry(global_index); } - template - inline typename Vector::iterator - Vector::begin() + template + inline typename Vector::iterator + Vector::begin() { - return values.get(); + return internal::Policy::begin(data); } - template - inline typename Vector::const_iterator - Vector::begin() const + template + inline typename Vector::const_iterator + Vector::begin() const { - return values.get(); + return internal::Policy::begin(data); } - template - inline typename Vector::iterator - Vector::end() + template + inline typename Vector::iterator + Vector::end() { - return values.get() + partitioner->local_size(); + return internal::Policy::begin(data) + + partitioner->local_size(); } - template - inline typename Vector::const_iterator - Vector::end() const + template + inline typename Vector::const_iterator + Vector::end() const { - return values.get() + partitioner->local_size(); + return internal::Policy::begin(data) + + partitioner->local_size(); } - template + template inline Number - Vector::operator()(const size_type global_index) const + Vector::operator()(const size_type global_index) const { + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); Assert( partitioner->in_local_range(global_index) || partitioner->ghost_indices().is_element(global_index), @@ -1465,15 +1627,18 @@ namespace LinearAlgebra vector_is_ghosted == true, ExcMessage("You tried to read a ghost element of this vector, " "but it has not imported its ghost values.")); - return values[partitioner->global_to_local(global_index)]; + return data.values[partitioner->global_to_local(global_index)]; } - template + template inline Number & - Vector::operator()(const size_type global_index) + Vector::operator()(const size_type global_index) { + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); Assert( partitioner->in_local_range(global_index) || partitioner->ghost_indices().is_element(global_index), @@ -1487,31 +1652,37 @@ namespace LinearAlgebra // (then, the compiler picks this method according to the C++ rule book // even if a human would pick the const method when this subsequent use // is just a read) - return values[partitioner->global_to_local(global_index)]; + return data.values[partitioner->global_to_local(global_index)]; } - template - inline Number Vector::operator[](const size_type global_index) const + template + inline Number Vector:: + operator[](const size_type global_index) const { return operator()(global_index); } - template - inline Number &Vector::operator[](const size_type global_index) + template + inline Number &Vector:: + operator[](const size_type global_index) { return operator()(global_index); } - template + template inline Number - Vector::local_element(const size_type local_index) const + Vector::local_element( + const size_type local_index) const { + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); AssertIndexRange(local_index, partitioner->local_size() + partitioner->n_ghost_indices()); @@ -1519,28 +1690,44 @@ namespace LinearAlgebra Assert(local_index < local_size() || vector_is_ghosted == true, ExcMessage("You tried to read a ghost element of this vector, " "but it has not imported its ghost values.")); - return values[local_index]; + + return data.values[local_index]; } - template + template inline Number & - Vector::local_element(const size_type local_index) + Vector::local_element(const size_type local_index) { + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); + AssertIndexRange(local_index, partitioner->local_size() + partitioner->n_ghost_indices()); - return values[local_index]; + + return data.values[local_index]; } - template + template + inline Number * + Vector::get_values() const + { + return internal::Policy::get_values(data); + } + + + + template template inline void - Vector::extract_subvector_to(const std::vector &indices, - std::vector &values) const + Vector::extract_subvector_to( + const std::vector &indices, + std::vector & values) const { for (size_type i = 0; i < indices.size(); ++i) values[i] = operator()(indices[i]); @@ -1548,12 +1735,13 @@ namespace LinearAlgebra - template + template template inline void - Vector::extract_subvector_to(ForwardIterator indices_begin, - const ForwardIterator indices_end, - OutputIterator values_begin) const + Vector::extract_subvector_to( + ForwardIterator indices_begin, + const ForwardIterator indices_end, + OutputIterator values_begin) const { while (indices_begin != indices_end) { @@ -1565,11 +1753,12 @@ namespace LinearAlgebra - template + template template inline void - Vector::add(const std::vector & indices, - const ::dealii::Vector &values) + Vector::add( + const std::vector & indices, + const ::dealii::Vector &values) { AssertDimension(indices.size(), values.size()); for (size_type i = 0; i < indices.size(); ++i) @@ -1584,12 +1773,12 @@ namespace LinearAlgebra - template + template template inline void - Vector::add(const size_type n_elements, - const size_type * indices, - const OtherNumber *values) + Vector::add(const size_type n_elements, + const size_type * indices, + const OtherNumber *values) { for (size_type i = 0; i < n_elements; ++i, ++indices, ++values) { @@ -1603,18 +1792,18 @@ namespace LinearAlgebra - template + template inline const MPI_Comm & - Vector::get_mpi_communicator() const + Vector::get_mpi_communicator() const { return partitioner->get_mpi_communicator(); } - template + template inline const std::shared_ptr & - Vector::get_partitioner() const + Vector::get_partitioner() const { return partitioner; } @@ -1633,10 +1822,10 @@ namespace LinearAlgebra * @relatesalso Vector * @author Katharina Kormann, Martin Kronbichler, 2011 */ -template +template inline void -swap(LinearAlgebra::distributed::Vector &u, - LinearAlgebra::distributed::Vector &v) +swap(LinearAlgebra::distributed::Vector &u, + LinearAlgebra::distributed::Vector &v) { u.swap(v); } @@ -1647,12 +1836,13 @@ swap(LinearAlgebra::distributed::Vector &u, * * @author Uwe Koecher, 2017 */ -template -struct is_serial_vector> +template +struct is_serial_vector> : std::false_type {}; + namespace internal { namespace LinearOperatorImplementation diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h index 6619314f9c..a476d61b9e 100644 --- a/include/deal.II/lac/la_parallel_vector.templates.h +++ b/include/deal.II/lac/la_parallel_vector.templates.h @@ -19,6 +19,8 @@ #include +#include +#include #include #include @@ -40,7 +42,7 @@ namespace LinearAlgebra { // In the import_from_ghosted_array_finish we might need to calculate the // maximal and minimal value for the given number type, which is not - // straight forward for complex numbers. Therefore, comparison of complex + // straightforward for complex numbers. Therefore, comparison of complex // numbers is prohibited and throws an exception. template Number @@ -75,13 +77,304 @@ namespace LinearAlgebra "implemented for complex numbers")); return a; } - } // namespace internal - template + // Resize the underlying array on the host or on the device + template + struct la_parallel_vector_templates_functions + { + static_assert( + std::is_same::value || + std::is_same::value, + "MemorySpace should be Host or CUDA"); + + static void + resize_val(const types::global_dof_index /*new_alloc_size*/, + types::global_dof_index & /*allocated_size*/, + ::dealii::MemorySpace::MemorySpaceData + & /*data*/) + {} + + static void + import(const ::dealii::LinearAlgebra::ReadWriteVector & /*V*/, + ::dealii::VectorOperation::values /*operation*/, + std::shared_ptr + /*communication_pattern*/, + const IndexSet & /*locally_owned_elem*/, + ::dealii::MemorySpace::MemorySpaceData + & /*data*/) + {} + + template + static void + linfty_norm_local( + const ::dealii::MemorySpace::MemorySpaceData + & /*data*/, + const unsigned int /*size*/, + RealType & /*max*/) + {} + }; + + template + struct la_parallel_vector_templates_functions + { + using size_type = types::global_dof_index; + + static void + resize_val(const types::global_dof_index new_alloc_size, + types::global_dof_index & allocated_size, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + if (new_alloc_size > allocated_size) + { + Assert(((allocated_size > 0 && data.values != nullptr) || + data.values == nullptr), + ExcInternalError()); + + Number *new_val; + Utilities::System::posix_memalign( + (void **)&new_val, 64, sizeof(Number) * new_alloc_size); + data.values.reset(new_val); + + allocated_size = new_alloc_size; + } + else if (new_alloc_size == 0) + { + data.values.reset(); + allocated_size = 0; + } + } + + static void + import(const ::dealii::LinearAlgebra::ReadWriteVector &V, + ::dealii::VectorOperation::values operation, + std::shared_ptr + communication_pattern, + const IndexSet &locally_owned_elem, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Assert( + (operation == ::dealii::VectorOperation::add) || + (operation == ::dealii::VectorOperation::insert), + ExcMessage( + "Only VectorOperation::add and VectorOperation::insert are allowed")); + + ::dealii::LinearAlgebra::distributed:: + Vector + tmp_vector(communication_pattern); + + // fill entries from ReadWriteVector into the distributed vector, + // including ghost entries. this is not really efficient right now + // because indices are translated twice, once by nth_index_in_set(i) + // and once for operator() of tmp_vector + const IndexSet &v_stored = V.get_stored_elements(); + for (size_type i = 0; i < v_stored.n_elements(); ++i) + tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i); + + tmp_vector.compress(operation); + + // Copy the local elements of tmp_vector to the right place in val + IndexSet tmp_index_set = tmp_vector.locally_owned_elements(); + if (operation == VectorOperation::add) + { + for (size_type i = 0; i < tmp_index_set.n_elements(); ++i) + { + data.values[locally_owned_elem.index_within_set( + tmp_index_set.nth_index_in_set(i))] += + tmp_vector.local_element(i); + } + } + else + { + for (size_type i = 0; i < tmp_index_set.n_elements(); ++i) + { + data.values[locally_owned_elem.index_within_set( + tmp_index_set.nth_index_in_set(i))] = + tmp_vector.local_element(i); + } + } + } + + template + static void + linfty_norm_local(const ::dealii::MemorySpace::MemorySpaceData< + Number, + ::dealii::MemorySpace::Host> &data, + const unsigned int size, + RealType & max) + { + for (size_type i = 0; i < size; ++i) + max = + std::max(numbers::NumberTraits::abs(data.values[i]), max); + } + }; + +#ifdef DEAL_II_COMPILER_CUDA_AWARE + template + struct la_parallel_vector_templates_functions + { + using size_type = types::global_dof_index; + + static void + resize_val(const types::global_dof_index new_alloc_size, + types::global_dof_index & allocated_size, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + static_assert( + std::is_same::value || + std::is_same::value, + "Number should be float or double for CUDA memory space"); + + if (new_alloc_size > allocated_size) + { + Assert(((allocated_size > 0 && data.values_dev != nullptr) || + data.values_dev == nullptr), + ExcInternalError()); + + Number *new_val_dev; + Utilities::CUDA::malloc(new_val_dev, new_alloc_size); + data.values_dev.reset(new_val_dev); + + allocated_size = new_alloc_size; + } + else if (new_alloc_size == 0) + { + data.values_dev.reset(); + allocated_size = 0; + } + } + + static void + import(const ReadWriteVector &V, + VectorOperation::values operation, + std::shared_ptr + communication_pattern, + const IndexSet &locally_owned_elem, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Assert( + (operation == ::dealii::VectorOperation::add) || + (operation == ::dealii::VectorOperation::insert), + ExcMessage( + "Only VectorOperation::add and VectorOperation::insert are allowed")); + + ::dealii::LinearAlgebra::distributed:: + Vector + tmp_vector(communication_pattern); + + // fill entries from ReadWriteVector into the distributed vector, + // including ghost entries. this is not really efficient right now + // because indices are translated twice, once by nth_index_in_set(i) + // and once for operator() of tmp_vector + const IndexSet & v_stored = V.get_stored_elements(); + const size_type n_elements = v_stored.n_elements(); + std::vector indices(n_elements); + for (size_type i = 0; i < n_elements; ++i) + indices[i] = communication_pattern->global_to_local( + v_stored.nth_index_in_set(i)); + // Move the indices to the device + size_type *indices_dev; + ::dealii::Utilities::CUDA::malloc(indices_dev, n_elements); + ::dealii::Utilities::CUDA::copy_to_dev(indices, indices_dev); + // Move the data to the device + Number *V_dev; + ::dealii::Utilities::CUDA::malloc(V_dev, n_elements); + cudaError_t cuda_error_code = cudaMemcpy(V_dev, + V.begin(), + n_elements * sizeof(Number), + cudaMemcpyHostToDevice); + AssertCuda(cuda_error_code); + + // Set the values in tmp_vector + const int n_blocks = + 1 + (n_elements - 1) / (::dealii::CUDAWrappers::chunk_size * + ::dealii::CUDAWrappers::block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::set_permutated + <<>>( + tmp_vector.begin(), V_dev, indices_dev, n_elements); + + tmp_vector.compress(operation); + + // Copy the local elements of tmp_vector to the right place in val + IndexSet tmp_index_set = tmp_vector.locally_owned_elements(); + const size_type tmp_n_elements = tmp_index_set.n_elements(); + indices.resize(tmp_n_elements); + for (size_type i = 0; i < tmp_n_elements; ++i) + indices[i] = locally_owned_elem.index_within_set( + tmp_index_set.nth_index_in_set(i)); + ::dealii::Utilities::CUDA::free(indices_dev); + ::dealii::Utilities::CUDA::malloc(indices_dev, tmp_n_elements); + ::dealii::Utilities::CUDA::copy_to_dev(indices, indices_dev); + + if (operation == VectorOperation::add) + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_permutated< + Number><<>>( + data.values_dev.get(), + tmp_vector.begin(), + indices_dev, + tmp_n_elements); + else + ::dealii::LinearAlgebra::CUDAWrappers::kernel::set_permutated< + Number><<>>( + data.values_dev.get(), + tmp_vector.begin(), + indices_dev, + tmp_n_elements); + + ::dealii::Utilities::CUDA::free(indices_dev); + ::dealii::Utilities::CUDA::free(V_dev); + } + + template + static void + linfty_norm_local(const ::dealii::MemorySpace::MemorySpaceData< + Number, + ::dealii::MemorySpace::CUDA> &data, + const unsigned int size, + RealType & result) + { + static_assert(std::is_same::value, + "RealType should be the same type as Number"); + + Number * result_device; + cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemset(result_device, Number(), sizeof(Number)); + + const int n_blocks = + 1 + (size - 1) / (::dealii::CUDAWrappers::chunk_size * + ::dealii::CUDAWrappers::block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction< + Number, + ::dealii::LinearAlgebra::CUDAWrappers::kernel::LInfty> + <<>>( + result_device, data.values_dev.get(), size); + + // Copy the result back to the host + error_code = cudaMemcpy(&result, + result_device, + sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + // Free the memory on the device + error_code = cudaFree(result_device); + AssertCuda(error_code); + } + }; +#endif + } // namespace internal + + + template void - Vector::clear_mpi_requests() + Vector::clear_mpi_requests() { #ifdef DEAL_II_WITH_MPI for (size_type j = 0; j < compress_requests.size(); j++) @@ -101,39 +394,23 @@ namespace LinearAlgebra - template + template void - Vector::resize_val(const size_type new_alloc_size) + Vector::resize_val(const size_type new_alloc_size) { - if (new_alloc_size > allocated_size) - { - Assert(((allocated_size > 0 && values != nullptr) || - values == nullptr), - ExcInternalError()); - - Number *new_val; - Utilities::System::posix_memalign((void **)&new_val, - 64, - sizeof(Number) * new_alloc_size); - values.reset(new_val); + internal::la_parallel_vector_templates_functions:: + resize_val(new_alloc_size, allocated_size, data); - allocated_size = new_alloc_size; - } - else if (new_alloc_size == 0) - { - values.reset(); - allocated_size = 0; - } thread_loop_partitioner = std::make_shared<::dealii::parallel::internal::TBBPartitioner>(); } - template + template void - Vector::reinit(const size_type size, - const bool omit_zeroing_entries) + Vector::reinit(const size_type size, + const bool omit_zeroing_entries) { clear_mpi_requests(); @@ -155,11 +432,11 @@ namespace LinearAlgebra - template + template template void - Vector::reinit(const Vector &v, - const bool omit_zeroing_entries) + Vector::reinit(const Vector &v, + const bool omit_zeroing_entries) { clear_mpi_requests(); Assert(v.partitioner.get() != nullptr, ExcNotInitialized()); @@ -192,11 +469,11 @@ namespace LinearAlgebra - template + template void - Vector::reinit(const IndexSet &locally_owned_indices, - const IndexSet &ghost_indices, - const MPI_Comm communicator) + Vector::reinit(const IndexSet &locally_owned_indices, + const IndexSet &ghost_indices, + const MPI_Comm communicator) { // set up parallel partitioner with index sets and communicator std::shared_ptr new_partitioner( @@ -208,10 +485,10 @@ namespace LinearAlgebra - template + template void - Vector::reinit(const IndexSet &locally_owned_indices, - const MPI_Comm communicator) + Vector::reinit(const IndexSet &locally_owned_indices, + const MPI_Comm communicator) { // set up parallel partitioner with index sets and communicator std::shared_ptr new_partitioner( @@ -221,9 +498,9 @@ namespace LinearAlgebra - template + template void - Vector::reinit( + Vector::reinit( const std::shared_ptr &partitioner_in) { clear_mpi_requests(); @@ -249,22 +526,20 @@ namespace LinearAlgebra - template - Vector::Vector() + template + Vector::Vector() : partitioner(new Utilities::MPI::Partitioner()) , allocated_size(0) - , values(nullptr, &free) { reinit(0); } - template - Vector::Vector(const Vector &v) + template + Vector::Vector(const Vector &v) : Subscriptor() , allocated_size(0) - , values(nullptr, &free) , vector_is_ghosted(false) { reinit(v, true); @@ -274,21 +549,19 @@ namespace LinearAlgebra const size_type this_size = local_size(); if (this_size > 0) { - dealii::internal::VectorOperations::Vector_copy - copier(v.values.get(), values.get()); - dealii::internal::VectorOperations::parallel_for( - copier, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::copy( + thread_loop_partitioner, partitioner->local_size(), v.data, data); } } - template - Vector::Vector(const IndexSet &local_range, - const IndexSet &ghost_indices, - const MPI_Comm communicator) + template + Vector::Vector(const IndexSet &local_range, + const IndexSet &ghost_indices, + const MPI_Comm communicator) : allocated_size(0) - , values(nullptr, &free) , vector_is_ghosted(false) { reinit(local_range, ghost_indices, communicator); @@ -296,11 +569,10 @@ namespace LinearAlgebra - template - Vector::Vector(const IndexSet &local_range, - const MPI_Comm communicator) + template + Vector::Vector(const IndexSet &local_range, + const MPI_Comm communicator) : allocated_size(0) - , values(nullptr, &free) , vector_is_ghosted(false) { reinit(local_range, communicator); @@ -308,10 +580,9 @@ namespace LinearAlgebra - template - Vector::Vector(const size_type size) + template + Vector::Vector(const size_type size) : allocated_size(0) - , values(nullptr, &free) , vector_is_ghosted(false) { reinit(size, false); @@ -319,11 +590,10 @@ namespace LinearAlgebra - template - Vector::Vector( + template + Vector::Vector( const std::shared_ptr &partitioner) : allocated_size(0) - , values(nullptr, &free) , vector_is_ghosted(false) { reinit(partitioner); @@ -331,8 +601,8 @@ namespace LinearAlgebra - template - inline Vector::~Vector() + template + inline Vector::~Vector() { try { @@ -344,9 +614,9 @@ namespace LinearAlgebra - template - inline Vector & - Vector::operator=(const Vector &c) + template + inline Vector & + Vector::operator=(const Vector &c) { #ifdef _MSC_VER return this->operator=(c); @@ -357,10 +627,11 @@ namespace LinearAlgebra - template + template template - inline Vector & - Vector::operator=(const Vector &c) + inline Vector & + Vector:: + operator=(const Vector &c) { Assert(c.partitioner.get() != nullptr, ExcNotInitialized()); @@ -411,10 +682,9 @@ namespace LinearAlgebra const size_type this_size = partitioner->local_size(); if (this_size > 0) { - dealii::internal::VectorOperations::Vector_copy - copier(c.values.get(), values.get()); - dealii::internal::VectorOperations::parallel_for( - copier, 0, this_size, thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::copy( + thread_loop_partitioner, this_size, c.data, data); } if (must_update_ghost_values) @@ -426,18 +696,21 @@ namespace LinearAlgebra - template + template template void - Vector::copy_locally_owned_data_from(const Vector &src) + Vector::copy_locally_owned_data_from( + const Vector &src) { AssertDimension(partitioner->local_size(), src.partitioner->local_size()); if (partitioner->local_size() > 0) { - dealii::internal::VectorOperations::Vector_copy - copier(src.values.get(), values.get()); - dealii::internal::VectorOperations::parallel_for( - copier, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::copy( + thread_loop_partitioner, + partitioner->local_size(), + src.data, + data); } } @@ -475,9 +748,10 @@ namespace LinearAlgebra } } // namespace petsc_helpers - template - Vector & - Vector::operator=(const PETScWrappers::MPI::Vector &petsc_vec) + template + Vector & + Vector:: + operator=(const PETScWrappers::MPI::Vector &petsc_vec) { // TODO: We would like to use the same compact infrastructure as for the // Trilinos vector below, but the interface through ReadWriteVector does @@ -517,9 +791,10 @@ namespace LinearAlgebra #ifdef DEAL_II_WITH_TRILINOS - template - Vector & - Vector::operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec) + template + Vector & + Vector:: + operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec) { # ifdef DEAL_II_WITH_MPI IndexSet combined_set = partitioner->locally_owned_range(); @@ -541,9 +816,10 @@ namespace LinearAlgebra - template + template void - Vector::compress(::dealii::VectorOperation::values operation) + Vector::compress( + ::dealii::VectorOperation::values operation) { compress_start(0, operation); compress_finish(operation); @@ -551,9 +827,9 @@ namespace LinearAlgebra - template + template void - Vector::update_ghost_values() const + Vector::update_ghost_values() const { update_ghost_values_start(); update_ghost_values_finish(); @@ -561,23 +837,35 @@ namespace LinearAlgebra - template + template void - Vector::zero_out_ghosts() const + Vector::zero_out_ghosts() const { - if (values != nullptr) - std::fill_n(values.get() + partitioner->local_size(), + if (data.values != nullptr) + std::fill_n(data.values.get() + partitioner->local_size(), partitioner->n_ghost_indices(), Number()); +#ifdef DEAL_II_COMPILER_CUDA_AWARE + if (data.values_dev != nullptr) + { + const cudaError_t cuda_error_code = + cudaMemset(data.values_dev.get() + partitioner->local_size(), + 0, + partitioner->n_ghost_indices() * sizeof(Number)); + AssertCuda(cuda_error_code); + } +#endif + vector_is_ghosted = false; } - template + template void - Vector::compress_start(const unsigned int counter, - ::dealii::VectorOperation::values operation) + Vector::compress_start( + const unsigned int counter, + ::dealii::VectorOperation::values operation) { (void)counter; (void)operation; @@ -593,10 +881,31 @@ namespace LinearAlgebra import_data = std_cxx14::make_unique(partitioner->n_import_indices()); +# ifdef DEAL_II_COMPILER_CUDA_AWARE + // TODO: for now move the data to the host and then move it back to the + // the device. We use values to store the elements because the function + // uses a view of the array and thus we need the data on the host to + // outlive the scope of the function. + if (std::is_same::value) + { + Number *new_val; + Utilities::System::posix_memalign((void **)&new_val, + 64, + sizeof(Number) * allocated_size); + data.values.reset(new_val); + + cudaError_t cuda_error_code = + cudaMemcpy(data.values.get(), + data.values_dev.get(), + allocated_size * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error_code); + } +# endif partitioner->import_from_ghosted_array_start( operation, counter, - ArrayView(values.get() + partitioner->local_size(), + ArrayView(data.values.get() + partitioner->local_size(), partitioner->n_ghost_indices()), ArrayView(import_data.get(), partitioner->n_import_indices()), compress_requests); @@ -605,9 +914,10 @@ namespace LinearAlgebra - template + template void - Vector::compress_finish(::dealii::VectorOperation::values operation) + Vector::compress_finish( + ::dealii::VectorOperation::values operation) { #ifdef DEAL_II_WITH_MPI vector_is_ghosted = false; @@ -625,10 +935,26 @@ namespace LinearAlgebra operation, ArrayView(import_data.get(), partitioner->n_import_indices()), - ArrayView(values.get(), partitioner->local_size()), - ArrayView(values.get() + partitioner->local_size(), + ArrayView(data.values.get(), partitioner->local_size()), + ArrayView(data.values.get() + partitioner->local_size(), partitioner->n_ghost_indices()), compress_requests); + +# ifdef DEAL_II_COMPILER_CUDA_AWARE + // TODO For now, the communication is done on the host, so we need to + // move the data back to the device. + if (std::is_same::value) + { + cudaError_t cuda_error_code = + cudaMemcpy(data.values_dev.get(), + data.values.get(), + allocated_size * sizeof(Number), + cudaMemcpyHostToDevice); + AssertCuda(cuda_error_code); + + data.values.reset(); + } +# endif #else (void)operation; #endif @@ -636,9 +962,10 @@ namespace LinearAlgebra - template + template void - Vector::update_ghost_values_start(const unsigned int counter) const + Vector::update_ghost_values_start( + const unsigned int counter) const { #ifdef DEAL_II_WITH_MPI // nothing to do when we neither have import nor ghost indices. @@ -654,11 +981,34 @@ namespace LinearAlgebra import_data = std_cxx14::make_unique(partitioner->n_import_indices()); +# ifdef DEAL_II_COMPILER_CUDA_AWARE + // TODO: for now move the data to the host and then move it back to the + // the device. We use values to store the elements because the function + // uses a view of the array and thus we need the data on the host to + // outlive the scope of the function. + if (std::is_same::value) + { + Number *new_val; + Utilities::System::posix_memalign((void **)&new_val, + 64, + sizeof(Number) * allocated_size); + + data.values.reset(new_val); + + cudaError_t cuda_error_code = + cudaMemcpy(data.values.get(), + data.values_dev.get(), + allocated_size * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error_code); + } +# endif + partitioner->export_to_ghosted_array_start( counter, - ArrayView(values.get(), partitioner->local_size()), + ArrayView(data.values.get(), partitioner->local_size()), ArrayView(import_data.get(), partitioner->n_import_indices()), - ArrayView(values.get() + partitioner->local_size(), + ArrayView(data.values.get() + partitioner->local_size(), partitioner->n_ghost_indices()), update_ghost_values_requests); @@ -669,9 +1019,9 @@ namespace LinearAlgebra - template + template void - Vector::update_ghost_values_finish() const + Vector::update_ghost_values_finish() const { #ifdef DEAL_II_WITH_MPI // wait for both sends and receives to complete, even though only @@ -685,19 +1035,34 @@ namespace LinearAlgebra Threads::Mutex::ScopedLock lock(mutex); partitioner->export_to_ghosted_array_finish( - ArrayView(values.get() + partitioner->local_size(), + ArrayView(data.values.get() + partitioner->local_size(), partitioner->n_ghost_indices()), update_ghost_values_requests); } +# ifdef DEAL_II_COMPILER_CUDA_AWARE + // TODO For now, the communication is done on the host, so we need to + // move the data back to the device. + if (std::is_same::value) + { + cudaError_t cuda_error_code = + cudaMemcpy(data.values_dev.get() + partitioner->local_size(), + data.values.get() + partitioner->local_size(), + partitioner->n_ghost_indices() * sizeof(Number), + cudaMemcpyHostToDevice); + AssertCuda(cuda_error_code); + + data.values.reset(); + } +# endif #endif vector_is_ghosted = true; } - template + template void - Vector::import( + Vector::import( const ReadWriteVector & V, VectorOperation::values operation, std::shared_ptr communication_pattern) @@ -725,7 +1090,8 @@ namespace LinearAlgebra "Utilities::MPI::Partitioner.")); } Vector tmp_vector(comm_pattern); - std::copy(begin(), end(), tmp_vector.begin()); + + data.copy_to(tmp_vector.begin(), local_size()); // fill entries from ReadWriteVector into the distributed vector, // including ghost entries. this is not really efficient right now @@ -774,12 +1140,12 @@ namespace LinearAlgebra } tmp_vector.compress(operation); - std::copy(tmp_vector.begin(), tmp_vector.end(), begin()); + data.copy_from(tmp_vector.begin(), local_size()); } - template + template void - Vector::swap(Vector &v) + Vector::swap(Vector &v) { #ifdef DEAL_II_WITH_MPI @@ -822,25 +1188,25 @@ namespace LinearAlgebra std::swap(partitioner, v.partitioner); std::swap(thread_loop_partitioner, v.thread_loop_partitioner); std::swap(allocated_size, v.allocated_size); - std::swap(values, v.values); + std::swap(data, v.data); std::swap(import_data, v.import_data); std::swap(vector_is_ghosted, v.vector_is_ghosted); } - template - Vector & - Vector::operator=(const Number s) + template + Vector & + Vector::operator=(const Number s) { const size_type this_size = local_size(); if (this_size > 0) { - dealii::internal::VectorOperations::Vector_set setter( - s, values.get()); - - dealii::internal::VectorOperations::parallel_for( - setter, 0, this_size, thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::set(thread_loop_partitioner, + this_size, + s, + data); } // if we call Vector::operator=0, we want to zero out all the entries @@ -853,36 +1219,37 @@ namespace LinearAlgebra - template + template void - Vector::reinit(const VectorSpaceVector &V, - const bool omit_zeroing_entries) + Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&V) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&V) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &down_V = dynamic_cast &>(V); + const VectorType &down_V = dynamic_cast(V); reinit(down_V, omit_zeroing_entries); } - template - Vector & - Vector::operator+=(const VectorSpaceVector &vv) + template + Vector & + Vector::operator+=(const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_add_v - vector_add(values.get(), v.values.get()); - dealii::internal::VectorOperations::parallel_for( - vector_add, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::add_vector( + thread_loop_partitioner, partitioner->local_size(), v.data, data); if (vector_is_ghosted) update_ghost_values(); @@ -892,21 +1259,21 @@ namespace LinearAlgebra - template - Vector & - Vector::operator-=(const VectorSpaceVector &vv) + template + Vector & + Vector::operator-=(const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_subtract_v - vector_subtract(values.get(), v.values.get()); - dealii::internal::VectorOperations::parallel_for( - vector_subtract, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::subtract_vector( + thread_loop_partitioner, partitioner->local_size(), v.data, data); if (vector_is_ghosted) update_ghost_values(); @@ -916,16 +1283,15 @@ namespace LinearAlgebra - template + template void - Vector::add(const Number a) + Vector::add(const Number a) { AssertIsFinite(a); - dealii::internal::VectorOperations::Vectorization_add_factor - vector_add(values.get(), a); - dealii::internal::VectorOperations::parallel_for( - vector_add, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::add_factor( + thread_loop_partitioner, partitioner->local_size(), a, data); if (vector_is_ghosted) update_ghost_values(); @@ -933,15 +1299,16 @@ namespace LinearAlgebra - template + template void - Vector::add_local(const Number a, - const VectorSpaceVector &vv) + Vector::add_local(const Number a, + const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertIsFinite(a); AssertDimension(local_size(), v.local_size()); @@ -950,17 +1317,17 @@ namespace LinearAlgebra if (a == Number(0.)) return; - dealii::internal::VectorOperations::Vectorization_add_av - vector_add(values.get(), v.values.get(), a); - dealii::internal::VectorOperations::parallel_for( - vector_add, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::add_av( + thread_loop_partitioner, partitioner->local_size(), a, v.data, data); } - template + template void - Vector::add(const Number a, const VectorSpaceVector &vv) + Vector::add(const Number a, + const VectorSpaceVector &vv) { add_local(a, vv); @@ -970,20 +1337,21 @@ namespace LinearAlgebra - template + template void - Vector::add(const Number a, - const VectorSpaceVector &vv, - const Number b, - const VectorSpaceVector &ww) + Vector::add(const Number a, + const VectorSpaceVector &vv, + const Number b, + const VectorSpaceVector &ww) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); - Assert(dynamic_cast *>(&ww) != nullptr, + const VectorType &v = dynamic_cast(vv); + Assert(dynamic_cast(&ww) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &w = dynamic_cast &>(ww); + const VectorType &w = dynamic_cast(ww); AssertIsFinite(a); AssertIsFinite(b); @@ -991,10 +1359,15 @@ namespace LinearAlgebra AssertDimension(local_size(), v.local_size()); AssertDimension(local_size(), w.local_size()); - dealii::internal::VectorOperations::Vectorization_add_avpbw - vector_add(values.get(), v.values.get(), w.values.get(), a, b); - dealii::internal::VectorOperations::parallel_for( - vector_add, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::add_avpbw( + thread_loop_partitioner, + partitioner->local_size(), + a, + b, + v.data, + w.data, + data); if (vector_is_ghosted) update_ghost_values(); @@ -1002,10 +1375,10 @@ namespace LinearAlgebra - template + template void - Vector::add(const std::vector &indices, - const std::vector & values) + Vector::add(const std::vector &indices, + const std::vector & values) { for (std::size_t i = 0; i < indices.size(); ++i) { @@ -1015,17 +1388,17 @@ namespace LinearAlgebra - template + template void - Vector::sadd(const Number x, const Vector &v) + Vector::sadd(const Number x, + const Vector &v) { AssertIsFinite(x); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_sadd_xv - vector_sadd(values.get(), v.values.get(), x); - dealii::internal::VectorOperations::parallel_for( - vector_sadd, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::sadd_xv( + thread_loop_partitioner, partitioner->local_size(), x, v.data, data); if (vector_is_ghosted) update_ghost_values(); @@ -1033,34 +1406,39 @@ namespace LinearAlgebra - template + template void - Vector::sadd_local(const Number x, - const Number a, - const VectorSpaceVector &vv) + Vector::sadd_local(const Number x, + const Number a, + const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert((dynamic_cast(&vv) != nullptr), ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertIsFinite(x); AssertIsFinite(a); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_sadd_xav - vector_sadd(values.get(), v.values.get(), a, x); - dealii::internal::VectorOperations::parallel_for( - vector_sadd, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::sadd_xav( + thread_loop_partitioner, + partitioner->local_size(), + x, + a, + v.data, + data); } - template + template void - Vector::sadd(const Number x, - const Number a, - const VectorSpaceVector &vv) + Vector::sadd(const Number x, + const Number a, + const VectorSpaceVector &vv) { sadd_local(x, a, vv); @@ -1070,13 +1448,13 @@ namespace LinearAlgebra - template + template void - Vector::sadd(const Number x, - const Number a, - const Vector &v, - const Number b, - const Vector &w) + Vector::sadd(const Number x, + const Number a, + const Vector &v, + const Number b, + const Vector &w) { AssertIsFinite(x); AssertIsFinite(a); @@ -1085,10 +1463,16 @@ namespace LinearAlgebra AssertDimension(local_size(), v.local_size()); AssertDimension(local_size(), w.local_size()); - dealii::internal::VectorOperations::Vectorization_sadd_xavbw - vector_sadd(values.get(), v.values.get(), w.values.get(), x, a, b); - dealii::internal::VectorOperations::parallel_for( - vector_sadd, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::sadd_xavbw( + thread_loop_partitioner, + partitioner->local_size(), + x, + a, + b, + v.data, + w.data, + data); if (vector_is_ghosted) update_ghost_values(); @@ -1096,16 +1480,15 @@ namespace LinearAlgebra - template - Vector & - Vector::operator*=(const Number factor) + template + Vector & + Vector::operator*=(const Number factor) { AssertIsFinite(factor); - dealii::internal::VectorOperations::Vectorization_multiply_factor - vector_multiply(values.get(), factor); - dealii::internal::VectorOperations::parallel_for( - vector_multiply, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::multiply_factor( + thread_loop_partitioner, partitioner->local_size(), factor, data); if (vector_is_ghosted) update_ghost_values(); @@ -1115,9 +1498,9 @@ namespace LinearAlgebra - template - Vector & - Vector::operator/=(const Number factor) + template + Vector & + Vector::operator/=(const Number factor) { operator*=(static_cast(1.) / factor); return *this; @@ -1125,21 +1508,23 @@ namespace LinearAlgebra - template + template void - Vector::scale(const VectorSpaceVector &vv) + Vector::scale(const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_scale - vector_scale(values.get(), v.values.get()); - dealii::internal::VectorOperations::parallel_for( - vector_scale, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::scale(thread_loop_partitioner, + local_size(), + v.data, + data); if (vector_is_ghosted) update_ghost_values(); @@ -1147,22 +1532,24 @@ namespace LinearAlgebra - template + template void - Vector::equ(const Number a, const VectorSpaceVector &vv) + Vector::equ(const Number a, + const VectorSpaceVector &vv) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); AssertIsFinite(a); AssertDimension(local_size(), v.local_size()); - dealii::internal::VectorOperations::Vectorization_equ_au - vector_equ(values.get(), v.values.get(), a); - dealii::internal::VectorOperations::parallel_for( - vector_equ, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::equ_au( + thread_loop_partitioner, partitioner->local_size(), a, v.data, data); + if (vector_is_ghosted) update_ghost_values(); @@ -1170,12 +1557,12 @@ namespace LinearAlgebra - template + template void - Vector::equ(const Number a, - const Vector &v, - const Number b, - const Vector &w) + Vector::equ(const Number a, + const Vector &v, + const Number b, + const Vector &w) { AssertIsFinite(a); AssertIsFinite(b); @@ -1183,10 +1570,15 @@ namespace LinearAlgebra AssertDimension(local_size(), v.local_size()); AssertDimension(local_size(), w.local_size()); - dealii::internal::VectorOperations::Vectorization_equ_aubv - vector_equ(values.get(), v.values.get(), w.values.get(), a, b); - dealii::internal::VectorOperations::parallel_for( - vector_equ, 0, partitioner->local_size(), thread_loop_partitioner); + dealii::internal::VectorOperations:: + functions::equ_aubv( + thread_loop_partitioner, + partitioner->local_size(), + a, + b, + v.data, + w.data, + data); if (vector_is_ghosted) update_ghost_values(); @@ -1194,67 +1586,44 @@ namespace LinearAlgebra - template + template bool - Vector::all_zero_local() const + Vector::all_zero() const { - const size_type local_size = partitioner->local_size(); - for (size_type i = 0; i < local_size; ++i) - if (values[i] != Number(0)) - return false; - return true; - } - - - - template - bool - Vector::all_zero() const - { - // use int instead of bool. in order to make global reduction operations - // work also when MPI_Init was not called, only call MPI_Allreduce - // commands when there is more than one processor (note that reinit() - // functions handle this case correctly through the job_supports_mpi() - // query). this is the same in all the reduce functions below - int local_result = -static_cast(all_zero_local()); - if (partitioner->n_mpi_processes() > 1) - return -Utilities::MPI::max(local_result, - partitioner->get_mpi_communicator()); - else - return -local_result; + return (linfty_norm() == 0) ? true : false; } - template + template template Number - Vector::inner_product_local(const Vector &v) const + Vector::inner_product_local( + const Vector &v) const { if (PointerComparison::equal(this, &v)) return norm_sqr_local(); AssertDimension(partitioner->local_size(), v.partitioner->local_size()); - Number sum; - dealii::internal::VectorOperations::Dot dot( - values.get(), v.values.get()); - dealii::internal::VectorOperations::parallel_reduce( - dot, 0, partitioner->local_size(), sum, thread_loop_partitioner); - AssertIsFinite(sum); - - return sum; + return dealii::internal::VectorOperations:: + functions::dot(thread_loop_partitioner, + partitioner->local_size(), + v.data, + data); } - template - Number Vector::operator*(const VectorSpaceVector &vv) const + template + Number Vector:: + operator*(const VectorSpaceVector &vv) const { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert((dynamic_cast(&vv) != nullptr), ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); + const VectorType &v = dynamic_cast(vv); Number local_result = inner_product_local(v); if (partitioner->n_mpi_processes() > 1) @@ -1266,15 +1635,17 @@ namespace LinearAlgebra - template - typename Vector::real_type - Vector::norm_sqr_local() const + template + typename Vector::real_type + Vector::norm_sqr_local() const { - real_type sum; - dealii::internal::VectorOperations::Norm2 norm2( - values.get()); - dealii::internal::VectorOperations::parallel_reduce( - norm2, 0, partitioner->local_size(), sum, thread_loop_partitioner); + real_type sum; + + + dealii::internal::VectorOperations:: + functions::norm_2( + thread_loop_partitioner, partitioner->local_size(), sum, data); + AssertIsFinite(sum); return sum; @@ -1282,28 +1653,27 @@ namespace LinearAlgebra - template + template Number - Vector::mean_value_local() const + Vector::mean_value_local() const { Assert(size() != 0, ExcEmptyObject()); if (partitioner->local_size() == 0) return Number(); - Number sum; - dealii::internal::VectorOperations::MeanValue mean(values.get()); - dealii::internal::VectorOperations::parallel_reduce( - mean, 0, partitioner->local_size(), sum, thread_loop_partitioner); + Number sum = ::dealii::internal::VectorOperations:: + functions::mean_value( + thread_loop_partitioner, partitioner->local_size(), data); return sum / real_type(partitioner->local_size()); } - template + template Number - Vector::mean_value() const + Vector::mean_value() const { Number local_result = mean_value_local(); if (partitioner->n_mpi_processes() > 1) @@ -1317,24 +1687,24 @@ namespace LinearAlgebra - template - typename Vector::real_type - Vector::l1_norm_local() const + template + typename Vector::real_type + Vector::l1_norm_local() const { - real_type sum; - dealii::internal::VectorOperations::Norm1 norm1( - values.get()); - dealii::internal::VectorOperations::parallel_reduce( - norm1, 0, partitioner->local_size(), sum, thread_loop_partitioner); + real_type sum; + + dealii::internal::VectorOperations:: + functions::norm_1( + thread_loop_partitioner, partitioner->local_size(), sum, data); return sum; } - template - typename Vector::real_type - Vector::l1_norm() const + template + typename Vector::real_type + Vector::l1_norm() const { real_type local_result = l1_norm_local(); if (partitioner->n_mpi_processes() > 1) @@ -1346,9 +1716,9 @@ namespace LinearAlgebra - template - typename Vector::real_type - Vector::norm_sqr() const + template + typename Vector::real_type + Vector::norm_sqr() const { real_type local_result = norm_sqr_local(); if (partitioner->n_mpi_processes() > 1) @@ -1360,32 +1730,33 @@ namespace LinearAlgebra - template - typename Vector::real_type - Vector::l2_norm() const + template + typename Vector::real_type + Vector::l2_norm() const { return std::sqrt(norm_sqr()); } - template - typename Vector::real_type - Vector::lp_norm_local(const real_type p) const + template + typename Vector::real_type + Vector::lp_norm_local(const real_type p) const { - real_type sum; - dealii::internal::VectorOperations::NormP normp( - values.get(), p); - dealii::internal::VectorOperations::parallel_reduce( - normp, 0, partitioner->local_size(), sum, thread_loop_partitioner); + real_type sum; + + dealii::internal::VectorOperations:: + functions::norm_p( + thread_loop_partitioner, partitioner->local_size(), sum, p, data); + return std::pow(sum, 1. / p); } - template - typename Vector::real_type - Vector::lp_norm(const real_type p) const + template + typename Vector::real_type + Vector::lp_norm(const real_type p) const { const real_type local_result = lp_norm_local(p); if (partitioner->n_mpi_processes() > 1) @@ -1399,24 +1770,24 @@ namespace LinearAlgebra - template - typename Vector::real_type - Vector::linfty_norm_local() const + template + typename Vector::real_type + Vector::linfty_norm_local() const { real_type max = 0.; const size_type local_size = partitioner->local_size(); - for (size_type i = 0; i < local_size; ++i) - max = std::max(numbers::NumberTraits::abs(values[i]), max); + internal::la_parallel_vector_templates_functions:: + linfty_norm_local(data, local_size, max); return max; } - template - inline typename Vector::real_type - Vector::linfty_norm() const + template + inline typename Vector::real_type + Vector::linfty_norm() const { const real_type local_result = linfty_norm_local(); if (partitioner->n_mpi_processes() > 1) @@ -1428,40 +1799,43 @@ namespace LinearAlgebra - template + template Number - Vector::add_and_dot_local(const Number a, - const Vector &v, - const Vector &w) + Vector::add_and_dot_local( + const Number a, + const Vector &v, + const Vector &w) { const size_type vec_size = partitioner->local_size(); AssertDimension(vec_size, v.local_size()); AssertDimension(vec_size, w.local_size()); - Number sum; - dealii::internal::VectorOperations::AddAndDot adder( - this->values.get(), v.values.get(), w.values.get(), a); - dealii::internal::VectorOperations::parallel_reduce( - adder, 0, vec_size, sum, thread_loop_partitioner); + Number sum = dealii::internal::VectorOperations:: + functions::add_and_dot( + thread_loop_partitioner, vec_size, a, v.data, w.data, data); + AssertIsFinite(sum); + return sum; } - template + template Number - Vector::add_and_dot(const Number a, - const VectorSpaceVector &vv, - const VectorSpaceVector &ww) + Vector::add_and_dot( + const Number a, + const VectorSpaceVector &vv, + const VectorSpaceVector &ww) { // Downcast. Throws an exception if invalid. - Assert(dynamic_cast *>(&vv) != nullptr, + using VectorType = Vector; + Assert((dynamic_cast(&vv) != nullptr), ExcVectorTypeNotCompatible()); - const Vector &v = dynamic_cast &>(vv); - Assert(dynamic_cast *>(&ww) != nullptr, + const VectorType &v = dynamic_cast(vv); + Assert((dynamic_cast(&ww) != nullptr), ExcVectorTypeNotCompatible()); - const Vector &w = dynamic_cast &>(ww); + const VectorType &w = dynamic_cast(ww); Number local_result = add_and_dot_local(a, v, w); if (partitioner->n_mpi_processes() > 1) @@ -1473,9 +1847,9 @@ namespace LinearAlgebra - template + template inline bool - Vector::partitioners_are_compatible( + Vector::partitioners_are_compatible( const Utilities::MPI::Partitioner &part) const { return partitioner->is_compatible(part); @@ -1483,9 +1857,9 @@ namespace LinearAlgebra - template + template inline bool - Vector::partitioners_are_globally_compatible( + Vector::partitioners_are_globally_compatible( const Utilities::MPI::Partitioner &part) const { return partitioner->is_globally_compatible(part); @@ -1493,9 +1867,9 @@ namespace LinearAlgebra - template + template std::size_t - Vector::memory_consumption() const + Vector::memory_consumption() const { std::size_t memory = sizeof(*this); memory += sizeof(Number) * static_cast(allocated_size); @@ -1514,12 +1888,12 @@ namespace LinearAlgebra - template + template void - Vector::print(std::ostream & out, - const unsigned int precision, - const bool scientific, - const bool across) const + Vector::print(std::ostream & out, + const unsigned int precision, + const bool scientific, + const bool across) const { Assert(partitioner.get() != nullptr, ExcInternalError()); AssertThrow(out, ExcIO()); @@ -1544,6 +1918,9 @@ namespace LinearAlgebra } #endif + std::vector stored_elements(allocated_size); + data.copy_to(stored_elements.data(), allocated_size); + out << "Process #" << partitioner->this_mpi_process() << std::endl << "Local range: [" << partitioner->local_range().first << ", " << partitioner->local_range().second @@ -1551,10 +1928,10 @@ namespace LinearAlgebra << "Vector data:" << std::endl; if (across) for (size_type i = 0; i < partitioner->local_size(); ++i) - out << local_element(i) << ' '; + out << stored_elements[i] << ' '; else for (size_type i = 0; i < partitioner->local_size(); ++i) - out << local_element(i) << std::endl; + out << stored_elements[i] << std::endl; out << std::endl; if (vector_is_ghosted) @@ -1563,13 +1940,13 @@ namespace LinearAlgebra if (across) for (size_type i = 0; i < partitioner->n_ghost_indices(); ++i) out << '(' << partitioner->ghost_indices().nth_index_in_set(i) - << '/' << local_element(partitioner->local_size() + i) + << '/' << stored_elements[partitioner->local_size() + i] << ") "; else for (size_type i = 0; i < partitioner->n_ghost_indices(); ++i) out << '(' << partitioner->ghost_indices().nth_index_in_set(i) - << '/' << local_element(partitioner->local_size() + i) << ")" - << std::endl; + << '/' << stored_elements[partitioner->local_size() + i] + << ")" << std::endl; out << std::endl; } out << std::flush; @@ -1596,8 +1973,7 @@ namespace LinearAlgebra out.precision(old_precision); } - } // end of namespace distributed - + } // namespace distributed } // namespace LinearAlgebra diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 9800cb55a6..96f2e06aef 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -20,6 +20,7 @@ #include +#include #include #include #include @@ -42,9 +43,9 @@ namespace LinearAlgebra { namespace distributed { - template + template class Vector; - } + } // namespace distributed } // namespace LinearAlgebra @@ -1915,18 +1916,19 @@ namespace internal } // selection for diagonal matrix around parallel deal.II vector - template + template inline void vector_updates( - const LinearAlgebra::distributed::Vector & src, - const DiagonalMatrix> &jacobi, - const bool start_zero, - const double factor1, - const double factor2, - LinearAlgebra::distributed::Vector &update1, - LinearAlgebra::distributed::Vector &update2, - LinearAlgebra::distributed::Vector &, - LinearAlgebra::distributed::Vector &dst) + const LinearAlgebra::distributed::Vector &src, + const DiagonalMatrix< + LinearAlgebra::distributed::Vector> &jacobi, + const bool start_zero, + const double factor1, + const double factor2, + LinearAlgebra::distributed::Vector & update1, + LinearAlgebra::distributed::Vector & update2, + LinearAlgebra::distributed::Vector &, + LinearAlgebra::distributed::Vector &dst) { VectorUpdater upd(src.begin(), jacobi.get_vector().begin(), @@ -2012,10 +2014,10 @@ namespace internal vector.add(-mean_value); } - template + template void set_initial_guess( - ::dealii::LinearAlgebra::distributed::Vector &vector) + ::dealii::LinearAlgebra::distributed::Vector &vector) { // Choose a high-frequency mode consisting of numbers between 0 and 1 // that is cheap to compute (cheaper than random numbers) but avoids @@ -2222,9 +2224,16 @@ PreconditionChebyshev:: (std::is_same>::value == false && - std::is_same>::value == false)) + ((std::is_same< + VectorType, + LinearAlgebra::distributed::Vector>::value == + false) || + (std::is_same< + VectorType, + LinearAlgebra::distributed::Vector>::value == + false)))) update3.reinit(src, true); const_cast< diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 8cc7b82630..9e67f98178 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -47,9 +47,9 @@ namespace LinearAlgebra class CommunicationPatternBase; namespace distributed { - template + template class Vector; - } + } // namespace distributed } // namespace LinearAlgebra #ifdef DEAL_II_WITH_PETSC @@ -290,9 +290,10 @@ namespace LinearAlgebra * be used if the same communication pattern is used multiple times. This * can be used to improve performance. */ + template void - import(const distributed::Vector &vec, - VectorOperation::values operation, + import(const distributed::Vector &vec, + VectorOperation::values operation, const std::shared_ptr &communication_pattern = std::shared_ptr()); diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index 664d468040..362ac1f1a5 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -40,10 +40,8 @@ # include #endif -#if defined(DEAL_II_WITH_CUDA) +#ifdef DEAL_II_WITH_CUDA # include - -# include #endif DEAL_II_NAMESPACE_OPEN @@ -51,6 +49,163 @@ DEAL_II_NAMESPACE_OPEN namespace LinearAlgebra { + namespace internal + { + // In the import_from_ghosted_array_finish we need to calculate the + // maximal and minimal value for the given number type, which is not + // straight forward for complex numbers. Therefore, comparison of complex + // numbers is prohibited and throws an assert. + template + Number + get_min(const Number a, const Number b) + { + return std::min(a, b); + } + + template + std::complex + get_min(const std::complex a, const std::complex) + { + AssertThrow(false, + ExcMessage("VectorOperation::min not " + "implemented for complex numbers")); + return a; + } + + template + Number + get_max(const Number a, const Number b) + { + return std::max(a, b); + } + + template + std::complex + get_max(const std::complex a, const std::complex) + { + AssertThrow(false, + ExcMessage("VectorOperation::max not " + "implemented for complex numbers")); + return a; + } + + + + template + struct read_write_vector_functions + { + static void + import(const std::shared_ptr + & communication_pattern, + const Number * values, + const ::dealii::VectorOperation::values operation, + ::dealii::LinearAlgebra::ReadWriteVector &rw_vector) + { + (void)communication_pattern; + (void)values; + (void)operation; + (void)rw_vector; + + static_assert( + std::is_same::value || + std::is_same::value, + "MemorySpace should be Host or CUDA"); + } + }; + + + + template + struct read_write_vector_functions + { + using size_type = types::global_dof_index; + + + static void + import(const std::shared_ptr + & communication_pattern, + const Number * values, + const ::dealii::VectorOperation::values operation, + ::dealii::LinearAlgebra::ReadWriteVector &rw_vector) + { + distributed::Vector tmp_vector( + communication_pattern); + + const unsigned int n_elements = communication_pattern->local_size(); + std::copy(values, values + n_elements, tmp_vector.begin()); + tmp_vector.update_ghost_values(); + + const IndexSet &stored = rw_vector.get_stored_elements(); + if (operation == VectorOperation::add) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) += + tmp_vector(stored.nth_index_in_set(i)); + else if (operation == VectorOperation::min) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = + get_min(tmp_vector(stored.nth_index_in_set(i)), + rw_vector.local_element(i)); + else if (operation == VectorOperation::max) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = + get_max(tmp_vector(stored.nth_index_in_set(i)), + rw_vector.local_element(i)); + else + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = tmp_vector(stored.nth_index_in_set(i)); + } + }; + + + +#ifdef DEAL_II_COMPILER_CUDA_AWARE + template + struct read_write_vector_functions + { + using size_type = types::global_dof_index; + + static void + import(const std::shared_ptr + & communication_pattern, + const Number * values, + const ::dealii::VectorOperation::values operation, + ::dealii::LinearAlgebra::ReadWriteVector &rw_vector) + { + distributed::Vector tmp_vector( + communication_pattern); + + const unsigned int n_elements = communication_pattern->local_size(); + cudaError_t cuda_error_code = cudaMemcpy(tmp_vector.begin(), + values, + n_elements * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error_code); + tmp_vector.update_ghost_values(); + + const IndexSet &stored = rw_vector.get_stored_elements(); + if (operation == VectorOperation::add) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) += + tmp_vector(stored.nth_index_in_set(i)); + else if (operation == VectorOperation::min) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = + get_min(tmp_vector(stored.nth_index_in_set(i)), + rw_vector.local_element(i)); + else if (operation == VectorOperation::max) + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = + get_max(tmp_vector(stored.nth_index_in_set(i)), + rw_vector.local_element(i)); + else + for (size_type i = 0; i < stored.n_elements(); ++i) + rw_vector.local_element(i) = tmp_vector(stored.nth_index_in_set(i)); + } + }; +#endif + } // namespace internal + + template void ReadWriteVector::resize_val(const size_type new_alloc_size) @@ -178,10 +333,10 @@ namespace LinearAlgebra ReadWriteVector::apply(const Functor &func) { FunctorTemplate functor(*this, func); - internal::VectorOperations::parallel_for(functor, - 0, - n_elements(), - thread_loop_partitioner); + dealii::internal::VectorOperations::parallel_for(functor, + 0, + n_elements(), + thread_loop_partitioner); } @@ -250,52 +405,14 @@ namespace LinearAlgebra return *this; } - namespace internal - { - // In the import_from_ghosted_array_finish we need to calculate the maximal - // and minimal value for the given number type, which is not straight - // forward for complex numbers. Therefore, comparison of complex numbers is - // prohibited and throws an assert. - template - Number - get_min(const Number a, const Number b) - { - return std::min(a, b); - } - - template - std::complex - get_min(const std::complex a, const std::complex) - { - AssertThrow(false, - ExcMessage("VectorOperation::min not " - "implemented for complex numbers")); - return a; - } - - template - Number - get_max(const Number a, const Number b) - { - return std::max(a, b); - } - template - std::complex - get_max(const std::complex a, const std::complex) - { - AssertThrow(false, - ExcMessage("VectorOperation::max not " - "implemented for complex numbers")); - return a; - } - } // namespace internal template + template void ReadWriteVector::import( - const distributed::Vector &vec, - VectorOperation::values operation, + const distributed::Vector &vec, + VectorOperation::values operation, const std::shared_ptr &communication_pattern) { @@ -318,28 +435,10 @@ namespace LinearAlgebra ExcMessage("The communication pattern is not of type " "Utilities::MPI::Partitioner.")); } - distributed::Vector tmp_vector(comm_pattern); - std::copy(vec.begin(), vec.end(), tmp_vector.begin()); - tmp_vector.update_ghost_values(); - const IndexSet &stored = get_stored_elements(); - if (operation == VectorOperation::add) - for (size_type i = 0; i < stored.n_elements(); ++i) - local_element(i) += tmp_vector(stored.nth_index_in_set(i)); - else if (operation == VectorOperation::min) - for (size_type i = 0; i < stored.n_elements(); ++i) - local_element(i) = - internal::get_min(tmp_vector(stored.nth_index_in_set(i)), - local_element(i)); - else if (operation == VectorOperation::max) - for (size_type i = 0; i < stored.n_elements(); ++i) - local_element(i) = - internal::get_max(tmp_vector(stored.nth_index_in_set(i)), - local_element(i)); - else - for (size_type i = 0; i < stored.n_elements(); ++i) - local_element(i) = tmp_vector(stored.nth_index_in_set(i)); + internal::read_write_vector_functions::import( + comm_pattern, vec.begin(), operation, *this); } @@ -570,7 +669,7 @@ namespace LinearAlgebra -#if defined(DEAL_II_WITH_CUDA) +#ifdef DEAL_II_COMPILER_CUDA_AWARE template void ReadWriteVector::import( diff --git a/include/deal.II/lac/vector_operations_internal.h b/include/deal.II/lac/vector_operations_internal.h index 9fd318574f..9d629d3279 100644 --- a/include/deal.II/lac/vector_operations_internal.h +++ b/include/deal.II/lac/vector_operations_internal.h @@ -17,11 +17,16 @@ #ifndef dealii_vector_operations_internal_h #define dealii_vector_operations_internal_h +#include #include #include #include +#include #include +#include +#include + #include #include @@ -307,7 +312,7 @@ namespace internal template struct Vectorization_add_av { - Vectorization_add_av(Number *val, Number *v_val, Number factor) + Vectorization_add_av(Number *val, const Number *v_val, Number factor) : val(val) , v_val(v_val) , factor(factor) @@ -329,15 +334,18 @@ namespace internal } } - Number *val; - Number *v_val; - Number factor; + Number * val; + const Number *v_val; + Number factor; }; template struct Vectorization_sadd_xav { - Vectorization_sadd_xav(Number *val, Number *v_val, Number a, Number x) + Vectorization_sadd_xav(Number * val, + const Number *v_val, + Number a, + Number x) : val(val) , v_val(v_val) , a(a) @@ -360,16 +368,16 @@ namespace internal } } - Number *val; - Number *v_val; - Number a; - Number x; + Number * val; + const Number *v_val; + Number a; + Number x; }; template struct Vectorization_subtract_v { - Vectorization_subtract_v(Number *val, Number *v_val) + Vectorization_subtract_v(Number *val, const Number *v_val) : val(val) , v_val(v_val) {} @@ -390,8 +398,8 @@ namespace internal } } - Number *val; - Number *v_val; + Number * val; + const Number *v_val; }; template @@ -425,7 +433,7 @@ namespace internal template struct Vectorization_add_v { - Vectorization_add_v(Number *val, Number *v_val) + Vectorization_add_v(Number *val, const Number *v_val) : val(val) , v_val(v_val) {} @@ -446,18 +454,18 @@ namespace internal } } - Number *val; - Number *v_val; + Number * val; + const Number *v_val; }; template struct Vectorization_add_avpbw { - Vectorization_add_avpbw(Number *val, - Number *v_val, - Number *w_val, - Number a, - Number b) + Vectorization_add_avpbw(Number * val, + const Number *v_val, + const Number *w_val, + Number a, + Number b) : val(val) , v_val(v_val) , w_val(w_val) @@ -481,17 +489,17 @@ namespace internal } } - Number *val; - Number *v_val; - Number *w_val; - Number a; - Number b; + Number * val; + const Number *v_val; + const Number *w_val; + Number a; + Number b; }; template struct Vectorization_sadd_xv { - Vectorization_sadd_xv(Number *val, Number *v_val, Number x) + Vectorization_sadd_xv(Number *val, const Number *v_val, Number x) : val(val) , v_val(v_val) , x(x) @@ -513,20 +521,20 @@ namespace internal } } - Number *val; - Number *v_val; - Number x; + Number * val; + const Number *v_val; + Number x; }; template struct Vectorization_sadd_xavbw { - Vectorization_sadd_xavbw(Number *val, - Number *v_val, - Number *w_val, - Number x, - Number a, - Number b) + Vectorization_sadd_xavbw(Number * val, + const Number *v_val, + const Number *w_val, + Number x, + Number a, + Number b) : val(val) , v_val(v_val) , w_val(w_val) @@ -551,18 +559,18 @@ namespace internal } } - Number *val; - Number *v_val; - Number *w_val; - Number x; - Number a; - Number b; + Number * val; + const Number *v_val; + const Number *w_val; + Number x; + Number a; + Number b; }; template struct Vectorization_scale { - Vectorization_scale(Number *val, Number *v_val) + Vectorization_scale(Number *val, const Number *v_val) : val(val) , v_val(v_val) {} @@ -583,14 +591,14 @@ namespace internal } } - Number *val; - Number *v_val; + Number * val; + const Number *v_val; }; template struct Vectorization_equ_au { - Vectorization_equ_au(Number *val, Number *u_val, Number a) + Vectorization_equ_au(Number *val, const Number *u_val, Number a) : val(val) , u_val(u_val) , a(a) @@ -612,19 +620,19 @@ namespace internal } } - Number *val; - Number *u_val; - Number a; + Number * val; + const Number *u_val; + Number a; }; template struct Vectorization_equ_aubv { - Vectorization_equ_aubv(Number *val, - Number *u_val, - Number *v_val, - Number a, - Number b) + Vectorization_equ_aubv(Number * val, + const Number *u_val, + const Number *v_val, + Number a, + Number b) : val(val) , u_val(u_val) , v_val(v_val) @@ -648,11 +656,11 @@ namespace internal } } - Number *val; - Number *u_val; - Number *v_val; - Number a; - Number b; + Number * val; + const Number *u_val; + const Number *v_val; + Number a; + Number b; }; template @@ -1385,6 +1393,1060 @@ namespace internal (void)partitioner; #endif } + + + template + struct functions + { + static void + copy( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + { + static_assert( + std::is_same::value && + std::is_same::value, + "For the CUDA MemorySpace Number and Number2 should be the same type"); + } + + static void + set( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*s*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + add_vector( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + subtract_vector( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + add_factor( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + Number /*a*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + add_av( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*a*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + add_avpbw( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*a*/, + const Number /*b*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*w_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + sadd_xv( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*x*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + sadd_xav( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*x*/, + const Number /*a*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + sadd_xavbw( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*x*/, + const Number /*a*/, + const Number /*b*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*w_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + multiply_factor( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*factor*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + scale( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + equ_au( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*a*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static void + equ_aubv( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*a*/, + const Number /*b*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*w_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static Number + dot( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + { + return Number(); + } + + template + static void + norm_2( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + real_type & /*sum*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static Number + mean_value( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*data*/) + { + return Number(); + } + + template + static void + norm_1(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + real_type & /*sum*/, + Number * /*values*/, + Number * /*values_dev*/) + {} + + template + static void + norm_p( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + real_type & /*sum*/, + real_type /*p*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + {} + + static Number + add_and_dot( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + /*thread_loop_partitioner*/, + const size_type /*size*/, + const Number /*a*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*v_data*/, + const ::dealii::MemorySpace::MemorySpaceData + & /*w_data*/, + ::dealii::MemorySpace::MemorySpaceData & /*data*/) + { + return Number(); + } + }; + + + + template + struct functions + { + static void + copy(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vector_copy copier(v_data.values.get(), + data.values.get()); + parallel_for(copier, 0, size, thread_loop_partitioner); + } + + static void + set(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number s, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vector_set setter(s, data.values.get()); + parallel_for(setter, 0, size, thread_loop_partitioner); + } + + static void + add_vector(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_add_v vector_add(data.values.get(), + v_data.values.get()); + parallel_for(vector_add, 0, size, thread_loop_partitioner); + } + + static void + subtract_vector( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vectorization_subtract_v vector_subtract(data.values.get(), + v_data.values.get()); + parallel_for(vector_subtract, 0, size, thread_loop_partitioner); + } + + static void + add_factor(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + Number a, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_add_factor vector_add(data.values.get(), a); + parallel_for(vector_add, 0, size, thread_loop_partitioner); + } + + static void + add_av(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vectorization_add_av vector_add(data.values.get(), + v_data.values.get(), + a); + parallel_for(vector_add, 0, size, thread_loop_partitioner); + } + + static void + add_avpbw(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_add_avpbw vector_add( + data.values.get(), v_data.values.get(), w_data.values.get(), a, b); + parallel_for(vector_add, 0, size, thread_loop_partitioner); + } + + static void + sadd_xv(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number x, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_sadd_xv vector_sadd(data.values.get(), + v_data.values.get(), + x); + parallel_for(vector_sadd, 0, size, thread_loop_partitioner); + } + + static void + sadd_xav(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number x, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_sadd_xav vector_sadd(data.values.get(), + v_data.values.get(), + a, + x); + parallel_for(vector_sadd, 0, size, thread_loop_partitioner); + } + + static void + sadd_xavbw(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number x, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_sadd_xavbw vector_sadd( + data.values.get(), v_data.values.get(), w_data.values.get(), x, a, b); + parallel_for(vector_sadd, 0, size, thread_loop_partitioner); + } + + static void + multiply_factor( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number factor, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vectorization_multiply_factor vector_multiply(data.values.get(), + factor); + parallel_for(vector_multiply, 0, size, thread_loop_partitioner); + } + + static void + scale(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vectorization_scale vector_scale(data.values.get(), + v_data.values.get()); + parallel_for(vector_scale, 0, size, thread_loop_partitioner); + } + + static void + equ_au(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Vectorization_equ_au vector_equ(data.values.get(), + v_data.values.get(), + a); + parallel_for(vector_equ, 0, size, thread_loop_partitioner); + } + + static void + equ_aubv(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Vectorization_equ_aubv vector_equ( + data.values.get(), v_data.values.get(), w_data.values.get(), a, b); + parallel_for(vector_equ, 0, size, thread_loop_partitioner); + } + + static Number + dot(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Number sum; + dealii::internal::VectorOperations::Dot dot( + data.values.get(), v_data.values.get()); + dealii::internal::VectorOperations::parallel_reduce( + dot, 0, size, sum, thread_loop_partitioner); + AssertIsFinite(sum); + + return sum; + } + + template + static void + norm_2(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + real_type & sum, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Norm2 norm2(data.values.get()); + parallel_reduce(norm2, 0, size, sum, thread_loop_partitioner); + } + + static Number + mean_value(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Number sum; + MeanValue mean(data.values.get()); + parallel_reduce(mean, 0, size, sum, thread_loop_partitioner); + + return sum; + } + + template + static void + norm_1(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + real_type & sum, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Norm1 norm1(data.values.get()); + parallel_reduce(norm1, 0, size, sum, thread_loop_partitioner); + } + + template + static void + norm_p(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + real_type & sum, + real_type p, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + NormP normp(data.values.get(), p); + parallel_reduce(normp, 0, size, sum, thread_loop_partitioner); + } + + static Number + add_and_dot( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Number sum; + AddAndDot adder(data.values.get(), + v_data.values.get(), + w_data.values.get(), + a); + parallel_reduce(adder, 0, size, sum, thread_loop_partitioner); + + return sum; + } + }; + + + +#ifdef DEAL_II_COMPILER_CUDA_AWARE + template + struct functions + { + static const int block_size = + ::dealii::LinearAlgebra::CUDAWrappers::kernel::block_size; + static const int chunk_size = + ::dealii::LinearAlgebra::CUDAWrappers::kernel::chunk_size; + + static void + copy(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + cudaError_t cuda_error_code = cudaMemcpy(data.values_dev.get(), + v_data.values_dev.get(), + size * sizeof(Number), + cudaMemcpyDeviceToDevice); + AssertCuda(cuda_error_code); + } + + static void + set(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number s, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::set + <<>>(data.values_dev.get(), s, size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + add_vector(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV + <<>>(data.values_dev.get(), + 1., + v_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + subtract_vector( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV + <<>>(data.values_dev.get(), + -1., + v_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + add_factor(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + Number a, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_add + <<>>(data.values_dev.get(), a, size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + add_av(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV + <<>>(data.values_dev.get(), + a, + v_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + add_avpbw(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aVbW + <<>>(data.values_dev.get(), + a, + v_data.values_dev.get(), + b, + w_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + sadd_xv(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number x, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd + <<>>( + x, data.values_dev.get(), 1., v_data.values_dev.get(), size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + sadd_xav(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number x, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd + <<>>( + x, data.values_dev.get(), a, v_data.values_dev.get(), size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + sadd_xavbw(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number x, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd + <<>>(x, + data.values_dev.get(), + a, + v_data.values_dev.get(), + b, + w_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + multiply_factor( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number factor, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_scale + <<>>(data.values_dev.get(), factor, size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + scale(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::scale + <<>>(data.values_dev.get(), + v_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + equ_au(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ + <<>>(data.values_dev.get(), + a, + v_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static void + equ_aubv(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number a, + const Number b, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace:: + MemorySpaceData &data) + { + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ + <<>>(data.values_dev.get(), + a, + v_data.values_dev.get(), + b, + w_data.values_dev.get(), + size); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + } + + static Number + dot(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Number * result_device; + cudaError_t error_code = + cudaMalloc(&result_device, size * sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemset(result_device, Number(), sizeof(Number)); + + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::double_vector_reduction< + Number, + ::dealii::LinearAlgebra::CUDAWrappers::kernel::DotProduct> + <<>>(result_device, + data.values_dev.get(), + v_data.values_dev.get(), + static_cast( + size)); + + // Copy the result back to the host + Number result; + error_code = cudaMemcpy(&result, + result_device, + sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + // Free the memory on the device + error_code = cudaFree(result_device); + AssertCuda(error_code); + + AssertIsFinite(result); + + return result; + } + + template + static void + norm_2(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner> + thread_loop_partitioner, + const size_type size, + real_type & sum, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + sum = dot(thread_loop_partitioner, size, data, data); + } + + static Number + mean_value(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const ::dealii::MemorySpace:: + MemorySpaceData &data) + { + Number * result_device; + cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemset(result_device, Number(), sizeof(Number)); + + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction< + Number, + ::dealii::LinearAlgebra::CUDAWrappers::kernel::ElemSum> + <<>>(result_device, + data.values_dev.get(), + size); + + // Copy the result back to the host + Number result; + error_code = cudaMemcpy(&result, + result_device, + sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + // Free the memory on the device + error_code = cudaFree(result_device); + AssertCuda(error_code); + + return result; + } + + template + static void + norm_1(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + real_type & sum, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Number * result_device; + cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemset(result_device, Number(), sizeof(Number)); + + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction< + Number, + ::dealii::LinearAlgebra::CUDAWrappers::kernel::L1Norm> + <<>>(result_device, + data.values_dev.get(), + size); + + // Copy the result back to the host + error_code = cudaMemcpy(&sum, + result_device, + sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + // Free the memory on the device + error_code = cudaFree(result_device); + AssertCuda(error_code); + } + + template + static void + norm_p( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type, + real_type &, + real_type, + ::dealii::MemorySpace::MemorySpaceData &) + { + Assert(false, ExcNotImplemented()); + } + + static Number + add_and_dot( + std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>, + const size_type size, + const Number a, + const ::dealii::MemorySpace:: + MemorySpaceData &v_data, + const ::dealii::MemorySpace:: + MemorySpaceData &w_data, + ::dealii::MemorySpace::MemorySpaceData + &data) + { + Number * res_d; + cudaError_t error_code = cudaMalloc(&res_d, sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemset(res_d, 0., sizeof(Number)); + AssertCuda(error_code); + + const int n_blocks = 1 + (size - 1) / (chunk_size * block_size); + ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_and_dot + <<>>(res_d, + data.values_dev.get(), + v_data.values_dev.get(), + w_data.values_dev.get(), + a, + size); + + Number res; + error_code = + cudaMemcpy(&res, res_d, sizeof(Number), cudaMemcpyDeviceToHost); + AssertCuda(error_code); + error_code = cudaFree(res_d); + + return res; + } + }; +#endif } // namespace VectorOperations } // namespace internal diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 1459d498d8..54d3ff6500 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -45,7 +45,7 @@ namespace LinearAlgebra { namespace distributed { - template + template class Vector; } } // namespace LinearAlgebra diff --git a/source/base/cuda.cu b/source/base/cuda.cu index ac746b62e2..19a7ea64b5 100644 --- a/source/base/cuda.cu +++ b/source/base/cuda.cu @@ -15,9 +15,12 @@ #include #include +#include #include +#include #include +#include DEAL_II_NAMESPACE_OPEN @@ -47,6 +50,13 @@ namespace Utilities dealii::GrowingVectorMemory< LinearAlgebra::CUDAWrappers::Vector>::release_unused_memory(); + dealii::GrowingVectorMemory< + LinearAlgebra::distributed::Vector>:: + release_unused_memory(); + dealii::GrowingVectorMemory< + LinearAlgebra::distributed::Vector>:: + release_unused_memory(); + cusolverStatus_t cusolver_error_code = cusolverDnDestroy(cusolver_dn_handle); AssertCusolver(cusolver_error_code); diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 0abdfcba5d..27729b8e89 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -150,6 +150,8 @@ IF(DEAL_II_WITH_CUDA) cuda_solver_direct.cu cuda_sparse_matrix.cu cuda_vector.cu + la_parallel_vector.cu + read_write_vector.cu ) ENDIF() diff --git a/source/lac/la_parallel_vector.cc b/source/lac/la_parallel_vector.cc index 0ef8028af4..b56ed27e32 100644 --- a/source/lac/la_parallel_vector.cc +++ b/source/lac/la_parallel_vector.cc @@ -29,8 +29,10 @@ namespace LinearAlgebra { namespace distributed { -#define TEMPL_COPY_CONSTRUCTOR(S1, S2) \ - template Vector &Vector::operator=(const Vector &) +#define TEMPL_COPY_CONSTRUCTOR(S1, S2) \ + template Vector \ + &Vector::operator= \ + (const Vector &) TEMPL_COPY_CONSTRUCTOR(double, float); TEMPL_COPY_CONSTRUCTOR(float, double); diff --git a/source/lac/la_parallel_vector.cu b/source/lac/la_parallel_vector.cu new file mode 100644 index 0000000000..352cf2c25c --- /dev/null +++ b/source/lac/la_parallel_vector.cu @@ -0,0 +1,32 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace LinearAlgebra +{ + namespace distributed + { + template class Vector; + template class Vector; + } // namespace distributed +} // namespace LinearAlgebra + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/la_parallel_vector.inst.in b/source/lac/la_parallel_vector.inst.in index 856cdcf4c4..12c85c8a2a 100644 --- a/source/lac/la_parallel_vector.inst.in +++ b/source/lac/la_parallel_vector.inst.in @@ -21,7 +21,7 @@ for (SCALAR : REAL_AND_COMPLEX_SCALARS) \{ namespace distributed \{ - template class Vector; + template class Vector; \} \} } @@ -33,11 +33,15 @@ for (S1 : REAL_AND_COMPLEX_SCALARS; S2 : REAL_SCALARS) namespace distributed \{ template void - Vector::reinit(const Vector &, const bool); + Vector::reinit( + const Vector &, + const bool); template S1 - Vector::inner_product_local(const Vector &) const; + Vector::inner_product_local( + const Vector &) const; template void - Vector::copy_locally_owned_data_from(const Vector &); + Vector::copy_locally_owned_data_from< + S2>(const Vector &); \} \} } @@ -50,11 +54,15 @@ for (S1, S2 : COMPLEX_SCALARS) namespace distributed \{ template void - Vector::reinit(const Vector &, const bool); + Vector::reinit( + const Vector &, + const bool); template S1 - Vector::inner_product_local(const Vector &) const; + Vector::inner_product_local( + const Vector &) const; template void - Vector::copy_locally_owned_data_from(const Vector &); + Vector::copy_locally_owned_data_from< + S2>(const Vector &); \} \} } diff --git a/source/lac/read_write_vector.cu b/source/lac/read_write_vector.cu new file mode 100644 index 0000000000..dd39d52c25 --- /dev/null +++ b/source/lac/read_write_vector.cu @@ -0,0 +1,47 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace LinearAlgebra +{ + template void + ReadWriteVector::import( + const CUDAWrappers::Vector &, + VectorOperation::values, + const std::shared_ptr &); + template void + ReadWriteVector::import( + const distributed::Vector &, + VectorOperation::values, + const std::shared_ptr &); + + template void + ReadWriteVector::import( + const CUDAWrappers::Vector &, + VectorOperation::values, + const std::shared_ptr &); + template void + ReadWriteVector::import( + const distributed::Vector &, + VectorOperation::values, + const std::shared_ptr &); +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE