--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/cuda.h>
+#include <deal.II/base/exceptions.h>
+
+#include <memory>
+
+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 <typename Number>
+ 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 <typename Number, typename MemorySpace>
+ struct MemorySpaceData
+ {
+ MemorySpaceData()
+ {
+ static_assert(std::is_same<MemorySpace, Host>::value ||
+ std::is_same<MemorySpace, CUDA>::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<Number[], decltype(&free)> values;
+
+ /**
+ * Pointer to data on the device.
+ */
+ std::unique_ptr<Number[]> values_dev;
+ };
+
+
+
+ /**
+ * Swap function similar to std::swap.
+ */
+ template <typename Number, typename MemorySpace>
+ inline void
+ swap(MemorySpaceData<Number, MemorySpace> &,
+ MemorySpaceData<Number, MemorySpace> &)
+ {
+ static_assert(std::is_same<MemorySpace, Host>::value ||
+ std::is_same<MemorySpace, CUDA>::value,
+ "MemorySpace should be Host or CUDA");
+ }
+
+#ifndef DOXYGEN
+
+ template <typename Number>
+ struct MemorySpaceData<Number, Host>
+ {
+ 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<Number[], decltype(&free)> values;
+
+ // This is not used but it allows to simplify the code until we start using
+ // CUDA-aware MPI.
+ std::unique_ptr<Number[]> values_dev;
+ };
+
+
+
+ template <typename Number>
+ inline void
+ swap(MemorySpaceData<Number, Host> &u, MemorySpaceData<Number, Host> &v)
+ {
+ std::swap(u.values, v.values);
+ }
+
+
+
+# ifdef DEAL_II_COMPILER_CUDA_AWARE
+
+ template <typename Number>
+ struct MemorySpaceData<Number, CUDA>
+ {
+ MemorySpaceData()
+ : values(nullptr, &free)
+ , values_dev(nullptr, deleter_device_array<Number>)
+ {}
+
+ 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<Number[], decltype(&free)> values;
+ std::unique_ptr<Number[], void (*)(Number *)> values_dev;
+ };
+
+
+
+ template <typename Number>
+ inline void
+ swap(MemorySpaceData<Number, CUDA> &u, MemorySpaceData<Number, CUDA> &v)
+ {
+ std::swap(u.values, v.values);
+ std::swap(u.values_dev, v.values_dev);
+ }
+
+# endif
+
+#endif
+
+} // namespace MemorySpace
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
static_cast<types::subdomain_id>(-2);
} // namespace numbers
-
DEAL_II_NAMESPACE_CLOSE
#endif
for (unsigned int i = 0; i < this->n_blocks(); ++i)
local_result =
std::max(local_result,
- -static_cast<int>(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(
#include <deal.II/base/config.h>
+#include <deal.II/base/memory_space.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/partitioner.h>
* 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
+ * <h4>CUDA support</h4>
+ *
+ * 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: <ul>
+ * <li> use get_values()
+ * <code>
+ * Vector<double, MemorySpace::CUDA> vector(local_range, comm);
+ * double* vector_dev = vector.get_values();
+ * std::vector<double> vector_host(local_range.n_elements(), 1.);
+ * Utilities::CUDA::copy_to_dev(vector_host, vector_dev);
+ * </code>
+ * <li> use import()
+ * <code>
+ * Vector<double, MemorySpace::CUDA> vector(local_range, comm);
+ * ReadWriteVector<double> rw_vector(local_range);
+ * for (auto & val : rw_vector)
+ * val = 1.;
+ * vector.import(rw_vector, VectorOperations::insert);
+ * </code>
+ * </ul>
+ * 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:
+ * <code> int n_devices = 0; cudaGetDeviceCount(&n_devices); int
+ * device_id = my_rank % n_devices;
+ * cudaSetDevice(device_id);
+ * </code>
+ * @see CUDAWrappers
+ *
+ * @author Katharina Kormann, Martin Kronbichler, Bruno Turcksin 2010, 2011,
+ * 2016, 2018
*/
- template <typename Number>
+ template <typename Number, typename MemorySpace = MemorySpace::Host>
class Vector : public ::dealii::LinearAlgebra::VectorSpaceVector<Number>,
public Subscriptor
{
public:
+ using memory_space = MemorySpace;
using value_type = Number;
using pointer = value_type *;
using const_pointer = const value_type *;
using size_type = types::global_dof_index;
using real_type = typename numbers::NumberTraits<Number>::real_type;
+ static_assert(
+ std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value ||
+ std::is_same<MemorySpace, ::dealii::MemorySpace::CUDA>::value,
+ "MemorySpace should be Host or CUDA");
+
/**
* @name 1: Basic Object-handling
*/
/**
* Copy constructor. Uses the parallel partitioning of @p in_vector.
*/
- Vector(const Vector<Number> &in_vector);
+ Vector(const Vector<Number, MemorySpace> &in_vector);
/**
* Construct a parallel vector of the given global size without any
*/
template <typename Number2>
void
- reinit(const Vector<Number2> &in_vector,
- const bool omit_zeroing_entries = false);
+ reinit(const Vector<Number2, MemorySpace> &in_vector,
+ const bool omit_zeroing_entries = false);
/**
* Initialize the vector. The local range is specified by @p
* handle memory separately.
*/
void
- swap(Vector<Number> &v);
+ swap(Vector<Number, MemorySpace> &v);
/**
* Assigns the vector to the parallel partitioning of the input vector
* at all, the vector will also update its ghost values in analogy to
* the respective setting the Trilinos and PETSc vectors.
*/
- Vector<Number> &
- operator=(const Vector<Number> &in_vector);
+ Vector<Number, MemorySpace> &
+ operator=(const Vector<Number, MemorySpace> &in_vector);
/**
* Assigns the vector to the parallel partitioning of the input vector
* the respective setting the Trilinos and PETSc vectors.
*/
template <typename Number2>
- Vector<Number> &
- operator=(const Vector<Number2> &in_vector);
+ Vector<Number, MemorySpace> &
+ operator=(const Vector<Number2, MemorySpace> &in_vector);
#ifdef DEAL_II_WITH_PETSC
/**
* ReadWriteVector instead.
*/
DEAL_II_DEPRECATED
- Vector<Number> &
+ Vector<Number, MemorySpace> &
operator=(const PETScWrappers::MPI::Vector &petsc_vec);
#endif
* ReadWriteVector instead.
*/
DEAL_II_DEPRECATED
- Vector<Number> &
+ Vector<Number, MemorySpace> &
operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec);
#endif
//@}
*/
template <typename Number2>
void
- copy_locally_owned_data_from(const Vector<Number2> &src);
+ copy_locally_owned_data_from(const Vector<Number2, MemorySpace> &src);
//@}
/**
* Multiply the entire vector by a fixed factor.
*/
- virtual Vector<Number> &
+ virtual Vector<Number, MemorySpace> &
operator*=(const Number factor) override;
/**
* Divide the entire vector by a fixed factor.
*/
- virtual Vector<Number> &
+ virtual Vector<Number, MemorySpace> &
operator/=(const Number factor) override;
/**
* Add the vector @p V to the present one.
*/
- virtual Vector<Number> &
+ virtual Vector<Number, MemorySpace> &
operator+=(const VectorSpaceVector<Number> &V) override;
/**
* Subtract the vector @p V from the present one.
*/
- virtual Vector<Number> &
+ virtual Vector<Number, MemorySpace> &
operator-=(const VectorSpaceVector<Number> &V) override;
/**
* 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(
* zero, also ghost elements are set to zero, otherwise they remain
* unchanged.
*/
- virtual Vector<Number> &
+ virtual Vector<Number, MemorySpace> &
operator=(const Number s) override;
/**
* s*(*this)+V</tt>.
*/
void
- sadd(const Number s, const Vector<Number> &V);
+ sadd(const Number s, const Vector<Number, MemorySpace> &V);
/**
* Scaling and multiple addition.
*/
DEAL_II_DEPRECATED
void
- sadd(const Number s,
- const Number a,
- const Vector<Number> &V,
- const Number b,
- const Vector<Number> &W);
+ sadd(const Number s,
+ const Number a,
+ const Vector<Number, MemorySpace> &V,
+ const Number b,
+ const Vector<Number, MemorySpace> &W);
/**
* Assignment <tt>*this = a*u + b*v</tt>.
*/
DEAL_II_DEPRECATED
void
- equ(const Number a,
- const Vector<Number> &u,
- const Number b,
- const Vector<Number> &v);
+ equ(const Number a,
+ const Vector<Number, MemorySpace> &u,
+ const Number b,
+ const Vector<Number, MemorySpace> &v);
//@}
* of the <i>locally owned</i> 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();
/**
* 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;
/**
* 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();
/**
* 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;
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
* @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 <typename OtherNumber>
void
*/
DeclException0(ExcVectorTypeNotCompatible);
+ /**
+ * Attempt to perform an operation not implemented on the device.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException0(ExcNotAllowedForCuda);
+
/**
* Exception
*/
const Number a,
const VectorSpaceVector<Number> &V);
- /**
- * Local part of all_zero().
- */
- bool
- all_zero_local() const;
-
/**
* Local part of the inner product of two vectors.
*/
template <typename Number2>
Number
- inner_product_local(const Vector<Number2> &V) const;
+ inner_product_local(const Vector<Number2, MemorySpace> &V) const;
/**
* Local part of norm_sqr().
* the add_and_dot() function.
*/
Number
- add_and_dot_local(const Number a,
- const Vector<Number> &V,
- const Vector<Number> &W);
+ add_and_dot_local(const Number a,
+ const Vector<Number, MemorySpace> &V,
+ const Vector<Number, MemorySpace> &W);
/**
* Shared pointer to store the parallel partitioning information. This
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
- * <code>delete[]</code>, but instead calls @p free().
+ * Underlying data structure storing the local elements of this vector.
*/
- std::unique_ptr<Number[], decltype(&free)> values;
+ mutable ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> data;
/**
* For parallel loops with TBB, this member variable stores the affinity
/*
* Make all other vector types friends.
*/
- template <typename Number2>
+ template <typename Number2, typename MemorySpace2>
friend class Vector;
/**
#ifndef DOXYGEN
+ namespace internal
+ {
+ template <typename Number, typename MemorySpace>
+ struct Policy
+ {
+ static inline typename Vector<Number, MemorySpace>::iterator
+ begin(::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
+ {
+ return nullptr;
+ }
+
+ static inline typename Vector<Number, MemorySpace>::const_iterator
+ begin(
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
+ {
+ return nullptr;
+ }
- template <typename Number>
+ static inline Number *
+ get_values(
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> &)
+ {
+ return nullptr;
+ }
+ };
+
+
+
+ template <typename Number>
+ struct Policy<Number, ::dealii::MemorySpace::Host>
+ {
+ static inline
+ typename Vector<Number, ::dealii::MemorySpace::Host>::iterator
+ begin(::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ return data.values.get();
+ }
+
+ static inline
+ typename Vector<Number, ::dealii::MemorySpace::Host>::const_iterator
+ begin(const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ return data.values.get();
+ }
+
+ static inline Number *
+ get_values(::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ return data.values.get();
+ }
+ };
+
+
+
+ template <typename Number>
+ struct Policy<Number, ::dealii::MemorySpace::CUDA>
+ {
+ static inline
+ typename Vector<Number, ::dealii::MemorySpace::CUDA>::iterator
+ begin(::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ return data.values_dev.get();
+ }
+
+ static inline
+ typename Vector<Number, ::dealii::MemorySpace::CUDA>::const_iterator
+ begin(const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ return data.values_dev.get();
+ }
+
+ static inline Number *
+ get_values(::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ return data.values_dev.get();
+ }
+ };
+ } // namespace internal
+
+
+ template <typename Number, typename MemorySpace>
inline bool
- Vector<Number>::has_ghost_elements() const
+ Vector<Number, MemorySpace>::has_ghost_elements() const
{
return vector_is_ghosted;
}
- template <typename Number>
- inline typename Vector<Number>::size_type
- Vector<Number>::size() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::size_type
+ Vector<Number, MemorySpace>::size() const
{
return partitioner->size();
}
- template <typename Number>
- inline typename Vector<Number>::size_type
- Vector<Number>::local_size() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::size_type
+ Vector<Number, MemorySpace>::local_size() const
{
return partitioner->local_size();
}
- template <typename Number>
- inline std::pair<typename Vector<Number>::size_type,
- typename Vector<Number>::size_type>
- Vector<Number>::local_range() const
+ template <typename Number, typename MemorySpace>
+ inline std::pair<typename Vector<Number, MemorySpace>::size_type,
+ typename Vector<Number, MemorySpace>::size_type>
+ Vector<Number, MemorySpace>::local_range() const
{
return partitioner->local_range();
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline bool
- Vector<Number>::in_local_range(const size_type global_index) const
+ Vector<Number, MemorySpace>::in_local_range(
+ const size_type global_index) const
{
return partitioner->in_local_range(global_index);
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline IndexSet
- Vector<Number>::locally_owned_elements() const
+ Vector<Number, MemorySpace>::locally_owned_elements() const
{
IndexSet is(size());
- template <typename Number>
- inline typename Vector<Number>::size_type
- Vector<Number>::n_ghost_entries() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::size_type
+ Vector<Number, MemorySpace>::n_ghost_entries() const
{
return partitioner->n_ghost_indices();
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline const IndexSet &
- Vector<Number>::ghost_elements() const
+ Vector<Number, MemorySpace>::ghost_elements() const
{
return partitioner->ghost_indices();
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline bool
- Vector<Number>::is_ghost_entry(const size_type global_index) const
+ Vector<Number, MemorySpace>::is_ghost_entry(
+ const size_type global_index) const
{
return partitioner->is_ghost_entry(global_index);
}
- template <typename Number>
- inline typename Vector<Number>::iterator
- Vector<Number>::begin()
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::iterator
+ Vector<Number, MemorySpace>::begin()
{
- return values.get();
+ return internal::Policy<Number, MemorySpace>::begin(data);
}
- template <typename Number>
- inline typename Vector<Number>::const_iterator
- Vector<Number>::begin() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::const_iterator
+ Vector<Number, MemorySpace>::begin() const
{
- return values.get();
+ return internal::Policy<Number, MemorySpace>::begin(data);
}
- template <typename Number>
- inline typename Vector<Number>::iterator
- Vector<Number>::end()
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::iterator
+ Vector<Number, MemorySpace>::end()
{
- return values.get() + partitioner->local_size();
+ return internal::Policy<Number, MemorySpace>::begin(data) +
+ partitioner->local_size();
}
- template <typename Number>
- inline typename Vector<Number>::const_iterator
- Vector<Number>::end() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::const_iterator
+ Vector<Number, MemorySpace>::end() const
{
- return values.get() + partitioner->local_size();
+ return internal::Policy<Number, MemorySpace>::begin(data) +
+ partitioner->local_size();
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline Number
- Vector<Number>::operator()(const size_type global_index) const
+ Vector<Number, MemorySpace>::operator()(const size_type global_index) const
{
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::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),
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 <typename Number>
+ template <typename Number, typename MemorySpace>
inline Number &
- Vector<Number>::operator()(const size_type global_index)
+ Vector<Number, MemorySpace>::operator()(const size_type global_index)
{
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::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),
// (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 <typename Number>
- inline Number Vector<Number>::operator[](const size_type global_index) const
+ template <typename Number, typename MemorySpace>
+ inline Number Vector<Number, MemorySpace>::
+ operator[](const size_type global_index) const
{
return operator()(global_index);
}
- template <typename Number>
- inline Number &Vector<Number>::operator[](const size_type global_index)
+ template <typename Number, typename MemorySpace>
+ inline Number &Vector<Number, MemorySpace>::
+ operator[](const size_type global_index)
{
return operator()(global_index);
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline Number
- Vector<Number>::local_element(const size_type local_index) const
+ Vector<Number, MemorySpace>::local_element(
+ const size_type local_index) const
{
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value),
+ ExcMessage(
+ "This function is only implemented for the Host memory space"));
AssertIndexRange(local_index,
partitioner->local_size() +
partitioner->n_ghost_indices());
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 <typename Number>
+ template <typename Number, typename MemorySpace>
inline Number &
- Vector<Number>::local_element(const size_type local_index)
+ Vector<Number, MemorySpace>::local_element(const size_type local_index)
{
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::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 <typename Number>
+ template <typename Number, typename MemorySpace>
+ inline Number *
+ Vector<Number, MemorySpace>::get_values() const
+ {
+ return internal::Policy<Number, MemorySpace>::get_values(data);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
template <typename OtherNumber>
inline void
- Vector<Number>::extract_subvector_to(const std::vector<size_type> &indices,
- std::vector<OtherNumber> &values) const
+ Vector<Number, MemorySpace>::extract_subvector_to(
+ const std::vector<size_type> &indices,
+ std::vector<OtherNumber> & values) const
{
for (size_type i = 0; i < indices.size(); ++i)
values[i] = operator()(indices[i]);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename ForwardIterator, typename OutputIterator>
inline void
- Vector<Number>::extract_subvector_to(ForwardIterator indices_begin,
- const ForwardIterator indices_end,
- OutputIterator values_begin) const
+ Vector<Number, MemorySpace>::extract_subvector_to(
+ ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
{
while (indices_begin != indices_end)
{
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename OtherNumber>
inline void
- Vector<Number>::add(const std::vector<size_type> & indices,
- const ::dealii::Vector<OtherNumber> &values)
+ Vector<Number, MemorySpace>::add(
+ const std::vector<size_type> & indices,
+ const ::dealii::Vector<OtherNumber> &values)
{
AssertDimension(indices.size(), values.size());
for (size_type i = 0; i < indices.size(); ++i)
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename OtherNumber>
inline void
- Vector<Number>::add(const size_type n_elements,
- const size_type * indices,
- const OtherNumber *values)
+ Vector<Number, MemorySpace>::add(const size_type n_elements,
+ const size_type * indices,
+ const OtherNumber *values)
{
for (size_type i = 0; i < n_elements; ++i, ++indices, ++values)
{
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline const MPI_Comm &
- Vector<Number>::get_mpi_communicator() const
+ Vector<Number, MemorySpace>::get_mpi_communicator() const
{
return partitioner->get_mpi_communicator();
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
- Vector<Number>::get_partitioner() const
+ Vector<Number, MemorySpace>::get_partitioner() const
{
return partitioner;
}
* @relatesalso Vector
* @author Katharina Kormann, Martin Kronbichler, 2011
*/
-template <typename Number>
+template <typename Number, typename MemorySpace>
inline void
-swap(LinearAlgebra::distributed::Vector<Number> &u,
- LinearAlgebra::distributed::Vector<Number> &v)
+swap(LinearAlgebra::distributed::Vector<Number, MemorySpace> &u,
+ LinearAlgebra::distributed::Vector<Number, MemorySpace> &v)
{
u.swap(v);
}
*
* @author Uwe Koecher, 2017
*/
-template <typename Number>
-struct is_serial_vector<LinearAlgebra::distributed::Vector<Number>>
+template <typename Number, typename MemorySpace>
+struct is_serial_vector<LinearAlgebra::distributed::Vector<Number, MemorySpace>>
: std::false_type
{};
+
namespace internal
{
namespace LinearOperatorImplementation
#include <deal.II/base/config.h>
+#include <deal.II/base/cuda.h>
+#include <deal.II/base/cuda_size.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/lac/exceptions.h>
{
// 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 <typename Number>
Number
"implemented for complex numbers"));
return a;
}
- } // namespace internal
- template <typename Number>
+ // Resize the underlying array on the host or on the device
+ template <typename Number, typename MemorySpace>
+ struct la_parallel_vector_templates_functions
+ {
+ static_assert(
+ std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value ||
+ std::is_same<MemorySpace, ::dealii::MemorySpace::CUDA>::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<Number, MemorySpace>
+ & /*data*/)
+ {}
+
+ static void
+ import(const ::dealii::LinearAlgebra::ReadWriteVector<Number> & /*V*/,
+ ::dealii::VectorOperation::values /*operation*/,
+ std::shared_ptr<const ::dealii::Utilities::MPI::Partitioner>
+ /*communication_pattern*/,
+ const IndexSet & /*locally_owned_elem*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*data*/)
+ {}
+
+ template <typename RealType>
+ static void
+ linfty_norm_local(
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*data*/,
+ const unsigned int /*size*/,
+ RealType & /*max*/)
+ {}
+ };
+
+ template <typename Number>
+ struct la_parallel_vector_templates_functions<Number,
+ ::dealii::MemorySpace::Host>
+ {
+ 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<Number, ::dealii::MemorySpace::Host> &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<Number> &V,
+ ::dealii::VectorOperation::values operation,
+ std::shared_ptr<const ::dealii::Utilities::MPI::Partitioner>
+ communication_pattern,
+ const IndexSet &locally_owned_elem,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Assert(
+ (operation == ::dealii::VectorOperation::add) ||
+ (operation == ::dealii::VectorOperation::insert),
+ ExcMessage(
+ "Only VectorOperation::add and VectorOperation::insert are allowed"));
+
+ ::dealii::LinearAlgebra::distributed::
+ Vector<Number, ::dealii::MemorySpace::Host>
+ 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 <typename RealType>
+ 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<Number>::abs(data.values[i]), max);
+ }
+ };
+
+#ifdef DEAL_II_COMPILER_CUDA_AWARE
+ template <typename Number>
+ struct la_parallel_vector_templates_functions<Number,
+ ::dealii::MemorySpace::CUDA>
+ {
+ 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<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ static_assert(
+ std::is_same<Number, float>::value ||
+ std::is_same<Number, double>::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<Number> &V,
+ VectorOperation::values operation,
+ std::shared_ptr<const Utilities::MPI::Partitioner>
+ communication_pattern,
+ const IndexSet &locally_owned_elem,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ Assert(
+ (operation == ::dealii::VectorOperation::add) ||
+ (operation == ::dealii::VectorOperation::insert),
+ ExcMessage(
+ "Only VectorOperation::add and VectorOperation::insert are allowed"));
+
+ ::dealii::LinearAlgebra::distributed::
+ Vector<Number, ::dealii::MemorySpace::CUDA>
+ 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<size_type> 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<Number>
+ <<<n_blocks, ::dealii::CUDAWrappers::block_size>>>(
+ 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><<<n_blocks, ::dealii::CUDAWrappers::block_size>>>(
+ data.values_dev.get(),
+ tmp_vector.begin(),
+ indices_dev,
+ tmp_n_elements);
+ else
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::set_permutated<
+ Number><<<n_blocks, ::dealii::CUDAWrappers::block_size>>>(
+ 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 <typename RealType>
+ 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<Number, RealType>::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<Number>>
+ <<<dim3(n_blocks, 1), dim3(::dealii::CUDAWrappers::block_size)>>>(
+ 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 <typename Number, typename MemorySpace>
void
- Vector<Number>::clear_mpi_requests()
+ Vector<Number, MemorySpace>::clear_mpi_requests()
{
#ifdef DEAL_II_WITH_MPI
for (size_type j = 0; j < compress_requests.size(); j++)
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::resize_val(const size_type new_alloc_size)
+ Vector<Number, MemorySpace>::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<Number, MemorySpace>::
+ 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 <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::reinit(const size_type size,
- const bool omit_zeroing_entries)
+ Vector<Number, MemorySpace>::reinit(const size_type size,
+ const bool omit_zeroing_entries)
{
clear_mpi_requests();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename Number2>
void
- Vector<Number>::reinit(const Vector<Number2> &v,
- const bool omit_zeroing_entries)
+ Vector<Number, MemorySpace>::reinit(const Vector<Number2, MemorySpace> &v,
+ const bool omit_zeroing_entries)
{
clear_mpi_requests();
Assert(v.partitioner.get() != nullptr, ExcNotInitialized());
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::reinit(const IndexSet &locally_owned_indices,
- const IndexSet &ghost_indices,
- const MPI_Comm communicator)
+ Vector<Number, MemorySpace>::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<const Utilities::MPI::Partitioner> new_partitioner(
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::reinit(const IndexSet &locally_owned_indices,
- const MPI_Comm communicator)
+ Vector<Number, MemorySpace>::reinit(const IndexSet &locally_owned_indices,
+ const MPI_Comm communicator)
{
// set up parallel partitioner with index sets and communicator
std::shared_ptr<const Utilities::MPI::Partitioner> new_partitioner(
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::reinit(
+ Vector<Number, MemorySpace>::reinit(
const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_in)
{
clear_mpi_requests();
- template <typename Number>
- Vector<Number>::Vector()
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::Vector()
: partitioner(new Utilities::MPI::Partitioner())
, allocated_size(0)
- , values(nullptr, &free)
{
reinit(0);
}
- template <typename Number>
- Vector<Number>::Vector(const Vector<Number> &v)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::Vector(const Vector<Number, MemorySpace> &v)
: Subscriptor()
, allocated_size(0)
- , values(nullptr, &free)
, vector_is_ghosted(false)
{
reinit(v, true);
const size_type this_size = local_size();
if (this_size > 0)
{
- dealii::internal::VectorOperations::Vector_copy<Number, Number>
- copier(v.values.get(), values.get());
- dealii::internal::VectorOperations::parallel_for(
- copier, 0, partitioner->local_size(), thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::copy(
+ thread_loop_partitioner, partitioner->local_size(), v.data, data);
}
}
- template <typename Number>
- Vector<Number>::Vector(const IndexSet &local_range,
- const IndexSet &ghost_indices,
- const MPI_Comm communicator)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::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);
- template <typename Number>
- Vector<Number>::Vector(const IndexSet &local_range,
- const MPI_Comm communicator)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::Vector(const IndexSet &local_range,
+ const MPI_Comm communicator)
: allocated_size(0)
- , values(nullptr, &free)
, vector_is_ghosted(false)
{
reinit(local_range, communicator);
- template <typename Number>
- Vector<Number>::Vector(const size_type size)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::Vector(const size_type size)
: allocated_size(0)
- , values(nullptr, &free)
, vector_is_ghosted(false)
{
reinit(size, false);
- template <typename Number>
- Vector<Number>::Vector(
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace>::Vector(
const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner)
: allocated_size(0)
- , values(nullptr, &free)
, vector_is_ghosted(false)
{
reinit(partitioner);
- template <typename Number>
- inline Vector<Number>::~Vector()
+ template <typename Number, typename MemorySpace>
+ inline Vector<Number, MemorySpace>::~Vector()
{
try
{
- template <typename Number>
- inline Vector<Number> &
- Vector<Number>::operator=(const Vector<Number> &c)
+ template <typename Number, typename MemorySpace>
+ inline Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator=(const Vector<Number, MemorySpace> &c)
{
#ifdef _MSC_VER
return this->operator=<Number>(c);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename Number2>
- inline Vector<Number> &
- Vector<Number>::operator=(const Vector<Number2> &c)
+ inline Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::
+ operator=(const Vector<Number2, MemorySpace> &c)
{
Assert(c.partitioner.get() != nullptr, ExcNotInitialized());
const size_type this_size = partitioner->local_size();
if (this_size > 0)
{
- dealii::internal::VectorOperations::Vector_copy<Number, Number2>
- copier(c.values.get(), values.get());
- dealii::internal::VectorOperations::parallel_for(
- copier, 0, this_size, thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number2, MemorySpace>::copy(
+ thread_loop_partitioner, this_size, c.data, data);
}
if (must_update_ghost_values)
- template <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename Number2>
void
- Vector<Number>::copy_locally_owned_data_from(const Vector<Number2> &src)
+ Vector<Number, MemorySpace>::copy_locally_owned_data_from(
+ const Vector<Number2, MemorySpace> &src)
{
AssertDimension(partitioner->local_size(), src.partitioner->local_size());
if (partitioner->local_size() > 0)
{
- dealii::internal::VectorOperations::Vector_copy<Number, Number2>
- copier(src.values.get(), values.get());
- dealii::internal::VectorOperations::parallel_for(
- copier, 0, partitioner->local_size(), thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number2, MemorySpace>::copy(
+ thread_loop_partitioner,
+ partitioner->local_size(),
+ src.data,
+ data);
}
}
}
} // namespace petsc_helpers
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator=(const PETScWrappers::MPI::Vector &petsc_vec)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::
+ 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
#ifdef DEAL_II_WITH_TRILINOS
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::
+ operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec)
{
# ifdef DEAL_II_WITH_MPI
IndexSet combined_set = partitioner->locally_owned_range();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::compress(::dealii::VectorOperation::values operation)
+ Vector<Number, MemorySpace>::compress(
+ ::dealii::VectorOperation::values operation)
{
compress_start(0, operation);
compress_finish(operation);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::update_ghost_values() const
+ Vector<Number, MemorySpace>::update_ghost_values() const
{
update_ghost_values_start();
update_ghost_values_finish();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::zero_out_ghosts() const
+ Vector<Number, MemorySpace>::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 <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::compress_start(const unsigned int counter,
- ::dealii::VectorOperation::values operation)
+ Vector<Number, MemorySpace>::compress_start(
+ const unsigned int counter,
+ ::dealii::VectorOperation::values operation)
{
(void)counter;
(void)operation;
import_data =
std_cxx14::make_unique<Number[]>(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<MemorySpace, ::dealii::MemorySpace::CUDA>::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<Number>(values.get() + partitioner->local_size(),
+ ArrayView<Number>(data.values.get() + partitioner->local_size(),
partitioner->n_ghost_indices()),
ArrayView<Number>(import_data.get(), partitioner->n_import_indices()),
compress_requests);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::compress_finish(::dealii::VectorOperation::values operation)
+ Vector<Number, MemorySpace>::compress_finish(
+ ::dealii::VectorOperation::values operation)
{
#ifdef DEAL_II_WITH_MPI
vector_is_ghosted = false;
operation,
ArrayView<const Number>(import_data.get(),
partitioner->n_import_indices()),
- ArrayView<Number>(values.get(), partitioner->local_size()),
- ArrayView<Number>(values.get() + partitioner->local_size(),
+ ArrayView<Number>(data.values.get(), partitioner->local_size()),
+ ArrayView<Number>(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<MemorySpace, ::dealii::MemorySpace::CUDA>::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
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::update_ghost_values_start(const unsigned int counter) const
+ Vector<Number, MemorySpace>::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.
import_data =
std_cxx14::make_unique<Number[]>(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<MemorySpace, ::dealii::MemorySpace::CUDA>::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<const Number>(values.get(), partitioner->local_size()),
+ ArrayView<const Number>(data.values.get(), partitioner->local_size()),
ArrayView<Number>(import_data.get(), partitioner->n_import_indices()),
- ArrayView<Number>(values.get() + partitioner->local_size(),
+ ArrayView<Number>(data.values.get() + partitioner->local_size(),
partitioner->n_ghost_indices()),
update_ghost_values_requests);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::update_ghost_values_finish() const
+ Vector<Number, MemorySpace>::update_ghost_values_finish() const
{
#ifdef DEAL_II_WITH_MPI
// wait for both sends and receives to complete, even though only
Threads::Mutex::ScopedLock lock(mutex);
partitioner->export_to_ghosted_array_finish(
- ArrayView<Number>(values.get() + partitioner->local_size(),
+ ArrayView<Number>(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<MemorySpace, ::dealii::MemorySpace::CUDA>::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 <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::import(
+ Vector<Number, MemorySpace>::import(
const ReadWriteVector<Number> & V,
VectorOperation::values operation,
std::shared_ptr<const CommunicationPatternBase> communication_pattern)
"Utilities::MPI::Partitioner."));
}
Vector<Number> 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
}
tmp_vector.compress(operation);
- std::copy(tmp_vector.begin(), tmp_vector.end(), begin());
+ data.copy_from(tmp_vector.begin(), local_size());
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::swap(Vector<Number> &v)
+ Vector<Number, MemorySpace>::swap(Vector<Number, MemorySpace> &v)
{
#ifdef DEAL_II_WITH_MPI
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 <typename Number>
- Vector<Number> &
- Vector<Number>::operator=(const Number s)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator=(const Number s)
{
const size_type this_size = local_size();
if (this_size > 0)
{
- dealii::internal::VectorOperations::Vector_set<Number> setter(
- s, values.get());
-
- dealii::internal::VectorOperations::parallel_for(
- setter, 0, this_size, thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::set(thread_loop_partitioner,
+ this_size,
+ s,
+ data);
}
// if we call Vector::operator=0, we want to zero out all the entries
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::reinit(const VectorSpaceVector<Number> &V,
- const bool omit_zeroing_entries)
+ Vector<Number, MemorySpace>::reinit(const VectorSpaceVector<Number> &V,
+ const bool omit_zeroing_entries)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&V) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
+ const VectorType &down_V = dynamic_cast<const VectorType &>(V);
reinit(down_V, omit_zeroing_entries);
}
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator+=(const VectorSpaceVector<Number> &vv)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator+=(const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_add_v<Number>
- 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<Number, Number, MemorySpace>::add_vector(
+ thread_loop_partitioner, partitioner->local_size(), v.data, data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator-=(const VectorSpaceVector<Number> &vv)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator-=(const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_subtract_v<Number>
- 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<Number, Number, MemorySpace>::subtract_vector(
+ thread_loop_partitioner, partitioner->local_size(), v.data, data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::add(const Number a)
+ Vector<Number, MemorySpace>::add(const Number a)
{
AssertIsFinite(a);
- dealii::internal::VectorOperations::Vectorization_add_factor<Number>
- vector_add(values.get(), a);
- dealii::internal::VectorOperations::parallel_for(
- vector_add, 0, partitioner->local_size(), thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::add_factor(
+ thread_loop_partitioner, partitioner->local_size(), a, data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::add_local(const Number a,
- const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::add_local(const Number a,
+ const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertIsFinite(a);
AssertDimension(local_size(), v.local_size());
if (a == Number(0.))
return;
- dealii::internal::VectorOperations::Vectorization_add_av<Number>
- 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<Number, Number, MemorySpace>::add_av(
+ thread_loop_partitioner, partitioner->local_size(), a, v.data, data);
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::add(const Number a,
+ const VectorSpaceVector<Number> &vv)
{
add_local(a, vv);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::add(const Number a,
- const VectorSpaceVector<Number> &vv,
- const Number b,
- const VectorSpaceVector<Number> &ww)
+ Vector<Number, MemorySpace>::add(const Number a,
+ const VectorSpaceVector<Number> &vv,
+ const Number b,
+ const VectorSpaceVector<Number> &ww)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
- Assert(dynamic_cast<const Vector<Number> *>(&ww) != nullptr,
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
+ Assert(dynamic_cast<const VectorType *>(&ww) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &w = dynamic_cast<const Vector<Number> &>(ww);
+ const VectorType &w = dynamic_cast<const VectorType &>(ww);
AssertIsFinite(a);
AssertIsFinite(b);
AssertDimension(local_size(), v.local_size());
AssertDimension(local_size(), w.local_size());
- dealii::internal::VectorOperations::Vectorization_add_avpbw<Number>
- 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<Number, Number, MemorySpace>::add_avpbw(
+ thread_loop_partitioner,
+ partitioner->local_size(),
+ a,
+ b,
+ v.data,
+ w.data,
+ data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::add(const std::vector<size_type> &indices,
- const std::vector<Number> & values)
+ Vector<Number, MemorySpace>::add(const std::vector<size_type> &indices,
+ const std::vector<Number> & values)
{
for (std::size_t i = 0; i < indices.size(); ++i)
{
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::sadd(const Number x, const Vector<Number> &v)
+ Vector<Number, MemorySpace>::sadd(const Number x,
+ const Vector<Number, MemorySpace> &v)
{
AssertIsFinite(x);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_sadd_xv<Number>
- 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<Number, Number, MemorySpace>::sadd_xv(
+ thread_loop_partitioner, partitioner->local_size(), x, v.data, data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::sadd_local(const Number x,
- const Number a,
- const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::sadd_local(const Number x,
+ const Number a,
+ const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert((dynamic_cast<const VectorType *>(&vv) != nullptr),
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertIsFinite(x);
AssertIsFinite(a);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_sadd_xav<Number>
- 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<Number, Number, MemorySpace>::sadd_xav(
+ thread_loop_partitioner,
+ partitioner->local_size(),
+ x,
+ a,
+ v.data,
+ data);
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::sadd(const Number x,
- const Number a,
- const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::sadd(const Number x,
+ const Number a,
+ const VectorSpaceVector<Number> &vv)
{
sadd_local(x, a, vv);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::sadd(const Number x,
- const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &w)
+ Vector<Number, MemorySpace>::sadd(const Number x,
+ const Number a,
+ const Vector<Number, MemorySpace> &v,
+ const Number b,
+ const Vector<Number, MemorySpace> &w)
{
AssertIsFinite(x);
AssertIsFinite(a);
AssertDimension(local_size(), v.local_size());
AssertDimension(local_size(), w.local_size());
- dealii::internal::VectorOperations::Vectorization_sadd_xavbw<Number>
- 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<Number, Number, MemorySpace>::sadd_xavbw(
+ thread_loop_partitioner,
+ partitioner->local_size(),
+ x,
+ a,
+ b,
+ v.data,
+ w.data,
+ data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator*=(const Number factor)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator*=(const Number factor)
{
AssertIsFinite(factor);
- dealii::internal::VectorOperations::Vectorization_multiply_factor<Number>
- vector_multiply(values.get(), factor);
- dealii::internal::VectorOperations::parallel_for(
- vector_multiply, 0, partitioner->local_size(), thread_loop_partitioner);
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::multiply_factor(
+ thread_loop_partitioner, partitioner->local_size(), factor, data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
- Vector<Number> &
- Vector<Number>::operator/=(const Number factor)
+ template <typename Number, typename MemorySpace>
+ Vector<Number, MemorySpace> &
+ Vector<Number, MemorySpace>::operator/=(const Number factor)
{
operator*=(static_cast<Number>(1.) / factor);
return *this;
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::scale(const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::scale(const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_scale<Number>
- 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<Number, Number, MemorySpace>::scale(thread_loop_partitioner,
+ local_size(),
+ v.data,
+ data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &vv)
+ Vector<Number, MemorySpace>::equ(const Number a,
+ const VectorSpaceVector<Number> &vv)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
AssertIsFinite(a);
AssertDimension(local_size(), v.local_size());
- dealii::internal::VectorOperations::Vectorization_equ_au<Number>
- 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<Number, Number, MemorySpace>::equ_au(
+ thread_loop_partitioner, partitioner->local_size(), a, v.data, data);
+
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::equ(const Number a,
- const Vector<Number> &v,
- const Number b,
- const Vector<Number> &w)
+ Vector<Number, MemorySpace>::equ(const Number a,
+ const Vector<Number, MemorySpace> &v,
+ const Number b,
+ const Vector<Number, MemorySpace> &w)
{
AssertIsFinite(a);
AssertIsFinite(b);
AssertDimension(local_size(), v.local_size());
AssertDimension(local_size(), w.local_size());
- dealii::internal::VectorOperations::Vectorization_equ_aubv<Number>
- 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<Number, Number, MemorySpace>::equ_aubv(
+ thread_loop_partitioner,
+ partitioner->local_size(),
+ a,
+ b,
+ v.data,
+ w.data,
+ data);
if (vector_is_ghosted)
update_ghost_values();
- template <typename Number>
+ template <typename Number, typename MemorySpace>
bool
- Vector<Number>::all_zero_local() const
+ Vector<Number, MemorySpace>::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 <typename Number>
- bool
- Vector<Number>::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<int>(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 <typename Number>
+ template <typename Number, typename MemorySpace>
template <typename Number2>
Number
- Vector<Number>::inner_product_local(const Vector<Number2> &v) const
+ Vector<Number, MemorySpace>::inner_product_local(
+ const Vector<Number2, MemorySpace> &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<Number, Number2> 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<Number, Number2, MemorySpace>::dot(thread_loop_partitioner,
+ partitioner->local_size(),
+ v.data,
+ data);
}
- template <typename Number>
- Number Vector<Number>::operator*(const VectorSpaceVector<Number> &vv) const
+ template <typename Number, typename MemorySpace>
+ Number Vector<Number, MemorySpace>::
+ operator*(const VectorSpaceVector<Number> &vv) const
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert((dynamic_cast<const VectorType *>(&vv) != nullptr),
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
Number local_result = inner_product_local(v);
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::norm_sqr_local() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::norm_sqr_local() const
{
- real_type sum;
- dealii::internal::VectorOperations::Norm2<Number, real_type> norm2(
- values.get());
- dealii::internal::VectorOperations::parallel_reduce(
- norm2, 0, partitioner->local_size(), sum, thread_loop_partitioner);
+ real_type sum;
+
+
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::norm_2(
+ thread_loop_partitioner, partitioner->local_size(), sum, data);
+
AssertIsFinite(sum);
return sum;
- template <typename Number>
+ template <typename Number, typename MemorySpace>
Number
- Vector<Number>::mean_value_local() const
+ Vector<Number, MemorySpace>::mean_value_local() const
{
Assert(size() != 0, ExcEmptyObject());
if (partitioner->local_size() == 0)
return Number();
- Number sum;
- dealii::internal::VectorOperations::MeanValue<Number> mean(values.get());
- dealii::internal::VectorOperations::parallel_reduce(
- mean, 0, partitioner->local_size(), sum, thread_loop_partitioner);
+ Number sum = ::dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::mean_value(
+ thread_loop_partitioner, partitioner->local_size(), data);
return sum / real_type(partitioner->local_size());
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
Number
- Vector<Number>::mean_value() const
+ Vector<Number, MemorySpace>::mean_value() const
{
Number local_result = mean_value_local();
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::l1_norm_local() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::l1_norm_local() const
{
- real_type sum;
- dealii::internal::VectorOperations::Norm1<Number, real_type> norm1(
- values.get());
- dealii::internal::VectorOperations::parallel_reduce(
- norm1, 0, partitioner->local_size(), sum, thread_loop_partitioner);
+ real_type sum;
+
+ dealii::internal::VectorOperations::
+ functions<Number, Number, MemorySpace>::norm_1(
+ thread_loop_partitioner, partitioner->local_size(), sum, data);
return sum;
}
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::l1_norm() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::l1_norm() const
{
real_type local_result = l1_norm_local();
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::norm_sqr() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::norm_sqr() const
{
real_type local_result = norm_sqr_local();
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::l2_norm() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::l2_norm() const
{
return std::sqrt(norm_sqr());
}
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::lp_norm_local(const real_type p) const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::lp_norm_local(const real_type p) const
{
- real_type sum;
- dealii::internal::VectorOperations::NormP<Number, real_type> 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<Number, Number, MemorySpace>::norm_p(
+ thread_loop_partitioner, partitioner->local_size(), sum, p, data);
+
return std::pow(sum, 1. / p);
}
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::lp_norm(const real_type p) const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::lp_norm(const real_type p) const
{
const real_type local_result = lp_norm_local(p);
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
- typename Vector<Number>::real_type
- Vector<Number>::linfty_norm_local() const
+ template <typename Number, typename MemorySpace>
+ typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::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<Number>::abs(values[i]), max);
+ internal::la_parallel_vector_templates_functions<Number, MemorySpace>::
+ linfty_norm_local(data, local_size, max);
return max;
}
- template <typename Number>
- inline typename Vector<Number>::real_type
- Vector<Number>::linfty_norm() const
+ template <typename Number, typename MemorySpace>
+ inline typename Vector<Number, MemorySpace>::real_type
+ Vector<Number, MemorySpace>::linfty_norm() const
{
const real_type local_result = linfty_norm_local();
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
+ template <typename Number, typename MemorySpace>
Number
- Vector<Number>::add_and_dot_local(const Number a,
- const Vector<Number> &v,
- const Vector<Number> &w)
+ Vector<Number, MemorySpace>::add_and_dot_local(
+ const Number a,
+ const Vector<Number, MemorySpace> &v,
+ const Vector<Number, MemorySpace> &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<Number> 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<Number, Number, MemorySpace>::add_and_dot(
+ thread_loop_partitioner, vec_size, a, v.data, w.data, data);
+
AssertIsFinite(sum);
+
return sum;
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
Number
- Vector<Number>::add_and_dot(const Number a,
- const VectorSpaceVector<Number> &vv,
- const VectorSpaceVector<Number> &ww)
+ Vector<Number, MemorySpace>::add_and_dot(
+ const Number a,
+ const VectorSpaceVector<Number> &vv,
+ const VectorSpaceVector<Number> &ww)
{
// Downcast. Throws an exception if invalid.
- Assert(dynamic_cast<const Vector<Number> *>(&vv) != nullptr,
+ using VectorType = Vector<Number, MemorySpace>;
+ Assert((dynamic_cast<const VectorType *>(&vv) != nullptr),
ExcVectorTypeNotCompatible());
- const Vector<Number> &v = dynamic_cast<const Vector<Number> &>(vv);
- Assert(dynamic_cast<const Vector<Number> *>(&ww) != nullptr,
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
+ Assert((dynamic_cast<const VectorType *>(&ww) != nullptr),
ExcVectorTypeNotCompatible());
- const Vector<Number> &w = dynamic_cast<const Vector<Number> &>(ww);
+ const VectorType &w = dynamic_cast<const VectorType &>(ww);
Number local_result = add_and_dot_local(a, v, w);
if (partitioner->n_mpi_processes() > 1)
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline bool
- Vector<Number>::partitioners_are_compatible(
+ Vector<Number, MemorySpace>::partitioners_are_compatible(
const Utilities::MPI::Partitioner &part) const
{
return partitioner->is_compatible(part);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline bool
- Vector<Number>::partitioners_are_globally_compatible(
+ Vector<Number, MemorySpace>::partitioners_are_globally_compatible(
const Utilities::MPI::Partitioner &part) const
{
return partitioner->is_globally_compatible(part);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
std::size_t
- Vector<Number>::memory_consumption() const
+ Vector<Number, MemorySpace>::memory_consumption() const
{
std::size_t memory = sizeof(*this);
memory += sizeof(Number) * static_cast<std::size_t>(allocated_size);
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
- Vector<Number>::print(std::ostream & out,
- const unsigned int precision,
- const bool scientific,
- const bool across) const
+ Vector<Number, MemorySpace>::print(std::ostream & out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool across) const
{
Assert(partitioner.get() != nullptr, ExcInternalError());
AssertThrow(out, ExcIO());
}
#endif
+ std::vector<Number> 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
<< "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)
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;
out.precision(old_precision);
}
- } // end of namespace distributed
-
+ } // namespace distributed
} // namespace LinearAlgebra
#include <deal.II/base/config.h>
+#include <deal.II/base/memory_space.h>
#include <deal.II/base/parallel.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/template_constraints.h>
{
namespace distributed
{
- template <typename number>
+ template <typename, typename>
class Vector;
- }
+ } // namespace distributed
} // namespace LinearAlgebra
}
// selection for diagonal matrix around parallel deal.II vector
- template <typename Number>
+ template <typename Number, typename MemorySpace>
inline void
vector_updates(
- const LinearAlgebra::distributed::Vector<Number> & src,
- const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number>> &jacobi,
- const bool start_zero,
- const double factor1,
- const double factor2,
- LinearAlgebra::distributed::Vector<Number> &update1,
- LinearAlgebra::distributed::Vector<Number> &update2,
- LinearAlgebra::distributed::Vector<Number> &,
- LinearAlgebra::distributed::Vector<Number> &dst)
+ const LinearAlgebra::distributed::Vector<Number, MemorySpace> &src,
+ const DiagonalMatrix<
+ LinearAlgebra::distributed::Vector<Number, MemorySpace>> &jacobi,
+ const bool start_zero,
+ const double factor1,
+ const double factor2,
+ LinearAlgebra::distributed::Vector<Number, MemorySpace> & update1,
+ LinearAlgebra::distributed::Vector<Number, MemorySpace> & update2,
+ LinearAlgebra::distributed::Vector<Number, MemorySpace> &,
+ LinearAlgebra::distributed::Vector<Number, MemorySpace> &dst)
{
VectorUpdater<Number> upd(src.begin(),
jacobi.get_vector().begin(),
vector.add(-mean_value);
}
- template <typename Number>
+ template <typename Number, typename MemorySpace>
void
set_initial_guess(
- ::dealii::LinearAlgebra::distributed::Vector<Number> &vector)
+ ::dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace> &vector)
{
// Choose a high-frequency mode consisting of numbers between 0 and 1
// that is cheap to compute (cheaper than random numbers) but avoids
(std::is_same<VectorType,
dealii::Vector<typename VectorType::value_type>>::value ==
false &&
- std::is_same<VectorType,
- LinearAlgebra::distributed::Vector<
- typename VectorType::value_type>>::value == false))
+ ((std::is_same<
+ VectorType,
+ LinearAlgebra::distributed::Vector<typename VectorType::value_type,
+ MemorySpace::Host>>::value ==
+ false) ||
+ (std::is_same<
+ VectorType,
+ LinearAlgebra::distributed::Vector<typename VectorType::value_type,
+ MemorySpace::CUDA>>::value ==
+ false))))
update3.reinit(src, true);
const_cast<
class CommunicationPatternBase;
namespace distributed
{
- template <typename>
+ template <typename, typename>
class Vector;
- }
+ } // namespace distributed
} // namespace LinearAlgebra
#ifdef DEAL_II_WITH_PETSC
* be used if the same communication pattern is used multiple times. This
* can be used to improve performance.
*/
+ template <typename MemorySpace>
void
- import(const distributed::Vector<Number> &vec,
- VectorOperation::values operation,
+ import(const distributed::Vector<Number, MemorySpace> &vec,
+ VectorOperation::values operation,
const std::shared_ptr<const CommunicationPatternBase>
&communication_pattern =
std::shared_ptr<const CommunicationPatternBase>());
# include <Epetra_Import.h>
#endif
-#if defined(DEAL_II_WITH_CUDA)
+#ifdef DEAL_II_WITH_CUDA
# include <deal.II/lac/cuda_vector.h>
-
-# include <cuda_runtime_api.h>
#endif
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 <typename Number>
+ Number
+ get_min(const Number a, const Number b)
+ {
+ return std::min(a, b);
+ }
+
+ template <typename Number>
+ std::complex<Number>
+ get_min(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::min not "
+ "implemented for complex numbers"));
+ return a;
+ }
+
+ template <typename Number>
+ Number
+ get_max(const Number a, const Number b)
+ {
+ return std::max(a, b);
+ }
+
+ template <typename Number>
+ std::complex<Number>
+ get_max(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::max not "
+ "implemented for complex numbers"));
+ return a;
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ struct read_write_vector_functions
+ {
+ static void
+ import(const std::shared_ptr<const ::dealii::Utilities::MPI::Partitioner>
+ & communication_pattern,
+ const Number * values,
+ const ::dealii::VectorOperation::values operation,
+ ::dealii::LinearAlgebra::ReadWriteVector<Number> &rw_vector)
+ {
+ (void)communication_pattern;
+ (void)values;
+ (void)operation;
+ (void)rw_vector;
+
+ static_assert(
+ std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value ||
+ std::is_same<MemorySpace, ::dealii::MemorySpace::CUDA>::value,
+ "MemorySpace should be Host or CUDA");
+ }
+ };
+
+
+
+ template <typename Number>
+ struct read_write_vector_functions<Number, ::dealii::MemorySpace::Host>
+ {
+ using size_type = types::global_dof_index;
+
+
+ static void
+ import(const std::shared_ptr<const ::dealii::Utilities::MPI::Partitioner>
+ & communication_pattern,
+ const Number * values,
+ const ::dealii::VectorOperation::values operation,
+ ::dealii::LinearAlgebra::ReadWriteVector<Number> &rw_vector)
+ {
+ distributed::Vector<Number, ::dealii::MemorySpace::Host> 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 <typename Number>
+ struct read_write_vector_functions<Number, ::dealii::MemorySpace::CUDA>
+ {
+ using size_type = types::global_dof_index;
+
+ static void
+ import(const std::shared_ptr<const ::dealii::Utilities::MPI::Partitioner>
+ & communication_pattern,
+ const Number * values,
+ const ::dealii::VectorOperation::values operation,
+ ::dealii::LinearAlgebra::ReadWriteVector<Number> &rw_vector)
+ {
+ distributed::Vector<Number, ::dealii::MemorySpace::Host> 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 <typename Number>
void
ReadWriteVector<Number>::resize_val(const size_type new_alloc_size)
ReadWriteVector<Number>::apply(const Functor &func)
{
FunctorTemplate<Functor> 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);
}
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 <typename Number>
- Number
- get_min(const Number a, const Number b)
- {
- return std::min(a, b);
- }
-
- template <typename Number>
- std::complex<Number>
- get_min(const std::complex<Number> a, const std::complex<Number>)
- {
- AssertThrow(false,
- ExcMessage("VectorOperation::min not "
- "implemented for complex numbers"));
- return a;
- }
-
- template <typename Number>
- Number
- get_max(const Number a, const Number b)
- {
- return std::max(a, b);
- }
- template <typename Number>
- std::complex<Number>
- get_max(const std::complex<Number> a, const std::complex<Number>)
- {
- AssertThrow(false,
- ExcMessage("VectorOperation::max not "
- "implemented for complex numbers"));
- return a;
- }
- } // namespace internal
template <typename Number>
+ template <typename MemorySpace>
void
ReadWriteVector<Number>::import(
- const distributed::Vector<Number> &vec,
- VectorOperation::values operation,
+ const distributed::Vector<Number, MemorySpace> &vec,
+ VectorOperation::values operation,
const std::shared_ptr<const CommunicationPatternBase>
&communication_pattern)
{
ExcMessage("The communication pattern is not of type "
"Utilities::MPI::Partitioner."));
}
- distributed::Vector<Number> 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<Number, MemorySpace>::import(
+ comm_pattern, vec.begin(), operation, *this);
}
-#if defined(DEAL_II_WITH_CUDA)
+#ifdef DEAL_II_COMPILER_CUDA_AWARE
template <typename Number>
void
ReadWriteVector<Number>::import(
#ifndef dealii_vector_operations_internal_h
#define dealii_vector_operations_internal_h
+#include <deal.II/base/memory_space.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/parallel.h>
#include <deal.II/base/thread_management.h>
+#include <deal.II/base/types.h>
#include <deal.II/base/vectorization.h>
+#include <deal.II/lac/cuda_kernels.h>
+#include <deal.II/lac/cuda_kernels.templates.h>
+
#include <cstdio>
#include <cstring>
template <typename Number>
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)
}
}
- Number *val;
- Number *v_val;
- Number factor;
+ Number * val;
+ const Number *v_val;
+ Number factor;
};
template <typename Number>
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)
}
}
- Number *val;
- Number *v_val;
- Number a;
- Number x;
+ Number * val;
+ const Number *v_val;
+ Number a;
+ Number x;
};
template <typename Number>
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)
{}
}
}
- Number *val;
- Number *v_val;
+ Number * val;
+ const Number *v_val;
};
template <typename Number>
template <typename Number>
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)
{}
}
}
- Number *val;
- Number *v_val;
+ Number * val;
+ const Number *v_val;
};
template <typename Number>
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)
}
}
- 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 <typename Number>
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)
}
}
- Number *val;
- Number *v_val;
- Number x;
+ Number * val;
+ const Number *v_val;
+ Number x;
};
template <typename Number>
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)
}
}
- 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 <typename Number>
struct Vectorization_scale
{
- Vectorization_scale(Number *val, Number *v_val)
+ Vectorization_scale(Number *val, const Number *v_val)
: val(val)
, v_val(v_val)
{}
}
}
- Number *val;
- Number *v_val;
+ Number * val;
+ const Number *v_val;
};
template <typename Number>
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)
}
}
- Number *val;
- Number *u_val;
- Number a;
+ Number * val;
+ const Number *u_val;
+ Number a;
};
template <typename Number>
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)
}
}
- 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 <typename Number>
(void)partitioner;
#endif
}
+
+
+ template <typename Number, typename Number2, typename MemorySpace>
+ struct functions
+ {
+ static void
+ copy(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number2, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {
+ static_assert(
+ std::is_same<MemorySpace, ::dealii::MemorySpace::CUDA>::value &&
+ std::is_same<Number, Number2>::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<Number, MemorySpace> & /*data*/)
+ {}
+
+ static void
+ add_vector(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {}
+
+ static void
+ subtract_vector(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {}
+
+ static void
+ add_factor(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ Number /*a*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*w_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*w_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace> & /*data*/)
+ {}
+
+ static void
+ scale(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*w_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {}
+
+ static Number
+ dot(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number2, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {
+ return Number();
+ }
+
+ template <typename real_type>
+ 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<Number, MemorySpace>
+ & /*v_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {}
+
+ static Number
+ mean_value(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ /*thread_loop_partitioner*/,
+ const size_type /*size*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*data*/)
+ {
+ return Number();
+ }
+
+ template <typename real_type>
+ 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 <typename real_type>
+ 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<Number, MemorySpace> & /*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<Number, MemorySpace>
+ & /*v_data*/,
+ const ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace>
+ & /*w_data*/,
+ ::dealii::MemorySpace::MemorySpaceData<Number, MemorySpace> & /*data*/)
+ {
+ return Number();
+ }
+ };
+
+
+
+ template <typename Number, typename Number2>
+ struct functions<Number, Number2, ::dealii::MemorySpace::Host>
+ {
+ static void
+ copy(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ thread_loop_partitioner,
+ const size_type size,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number2, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vector_copy<Number, Number2> 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<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vector_set<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_add_v<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vectorization_subtract_v<Number> 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<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_add_factor<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vectorization_add_av<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_add_avpbw<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_sadd_xv<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_sadd_xav<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_sadd_xavbw<Number> 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<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vectorization_multiply_factor<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vectorization_scale<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Vectorization_equ_au<Number> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Vectorization_equ_aubv<Number> 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<Number2, ::dealii::MemorySpace::Host> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Number sum;
+ dealii::internal::VectorOperations::Dot<Number, Number2> 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 <typename real_type>
+ static void
+ norm_2(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ thread_loop_partitioner,
+ const size_type size,
+ real_type & sum,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Norm2<Number, real_type> 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<Number, ::dealii::MemorySpace::Host> &data)
+ {
+ Number sum;
+ MeanValue<Number> mean(data.values.get());
+ parallel_reduce(mean, 0, size, sum, thread_loop_partitioner);
+
+ return sum;
+ }
+
+ template <typename real_type>
+ static void
+ norm_1(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ thread_loop_partitioner,
+ const size_type size,
+ real_type & sum,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Norm1<Number, real_type> norm1(data.values.get());
+ parallel_reduce(norm1, 0, size, sum, thread_loop_partitioner);
+ }
+
+ template <typename real_type>
+ 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<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ NormP<Number, real_type> 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<Number, ::dealii::MemorySpace::Host> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::Host> &w_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::Host>
+ &data)
+ {
+ Number sum;
+ AddAndDot<Number> 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 <typename Number>
+ struct functions<Number, Number, ::dealii::MemorySpace::CUDA>
+ {
+ 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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &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<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::set<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_add<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aVbW<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(
+ 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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(
+ 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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_scale<Number>
+ <<<n_blocks, block_size>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::scale<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &w_data,
+ ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &data)
+ {
+ const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+ ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &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<Number>>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
+ data.values_dev.get(),
+ v_data.values_dev.get(),
+ static_cast<unsigned int>(
+ 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 <typename real_type>
+ static void
+ norm_2(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>
+ thread_loop_partitioner,
+ const size_type size,
+ real_type & sum,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &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<Number, ::dealii::MemorySpace::CUDA> &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<Number>>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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 <typename real_type>
+ static void
+ norm_1(std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>,
+ const size_type size,
+ real_type & sum,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &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<Number>>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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 <typename real_type>
+ static void
+ norm_p(
+ std::shared_ptr<::dealii::parallel::internal::TBBPartitioner>,
+ const size_type,
+ real_type &,
+ real_type,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA> &)
+ {
+ 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<Number, ::dealii::MemorySpace::CUDA> &v_data,
+ const ::dealii::MemorySpace::
+ MemorySpaceData<Number, ::dealii::MemorySpace::CUDA> &w_data,
+ ::dealii::MemorySpace::MemorySpaceData<Number,
+ ::dealii::MemorySpace::CUDA>
+ &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<Number>
+ <<<dim3(n_blocks, 1), dim3(block_size)>>>(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
{
namespace distributed
{
- template <typename>
+ template <typename, typename>
class Vector;
}
} // namespace LinearAlgebra
#include <deal.II/base/cuda.h>
#include <deal.II/base/exceptions.h>
+#include <deal.II/base/memory_space.h>
#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/vector_memory.templates.h>
DEAL_II_NAMESPACE_OPEN
dealii::GrowingVectorMemory<
LinearAlgebra::CUDAWrappers::Vector<double>>::release_unused_memory();
+ dealii::GrowingVectorMemory<
+ LinearAlgebra::distributed::Vector<float, dealii::MemorySpace::CUDA>>::
+ release_unused_memory();
+ dealii::GrowingVectorMemory<
+ LinearAlgebra::distributed::Vector<double, dealii::MemorySpace::CUDA>>::
+ release_unused_memory();
+
cusolverStatus_t cusolver_error_code =
cusolverDnDestroy(cusolver_dn_handle);
AssertCusolver(cusolver_error_code);
cuda_solver_direct.cu
cuda_sparse_matrix.cu
cuda_vector.cu
+ la_parallel_vector.cu
+ read_write_vector.cu
)
ENDIF()
{
namespace distributed
{
-#define TEMPL_COPY_CONSTRUCTOR(S1, S2) \
- template Vector<S1> &Vector<S1>::operator=<S2>(const Vector<S2> &)
+#define TEMPL_COPY_CONSTRUCTOR(S1, S2) \
+ template Vector<S1, ::dealii::MemorySpace::Host> \
+ &Vector<S1, ::dealii::MemorySpace::Host>::operator= \
+ <S2>(const Vector<S2, ::dealii::MemorySpace::Host> &)
TEMPL_COPY_CONSTRUCTOR(double, float);
TEMPL_COPY_CONSTRUCTOR(float, double);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace LinearAlgebra
+{
+ namespace distributed
+ {
+ template class Vector<float, ::dealii::MemorySpace::CUDA>;
+ template class Vector<double, ::dealii::MemorySpace::CUDA>;
+ } // namespace distributed
+} // namespace LinearAlgebra
+
+
+DEAL_II_NAMESPACE_CLOSE
\{
namespace distributed
\{
- template class Vector<SCALAR>;
+ template class Vector<SCALAR, ::dealii::MemorySpace::Host>;
\}
\}
}
namespace distributed
\{
template void
- Vector<S1>::reinit<S2>(const Vector<S2> &, const bool);
+ Vector<S1, ::dealii::MemorySpace::Host>::reinit<S2>(
+ const Vector<S2, ::dealii::MemorySpace::Host> &,
+ const bool);
template S1
- Vector<S1>::inner_product_local<S2>(const Vector<S2> &) const;
+ Vector<S1, ::dealii::MemorySpace::Host>::inner_product_local<S2>(
+ const Vector<S2, ::dealii::MemorySpace::Host> &) const;
template void
- Vector<S1>::copy_locally_owned_data_from<S2>(const Vector<S2> &);
+ Vector<S1, ::dealii::MemorySpace::Host>::copy_locally_owned_data_from<
+ S2>(const Vector<S2, ::dealii::MemorySpace::Host> &);
\}
\}
}
namespace distributed
\{
template void
- Vector<S1>::reinit<S2>(const Vector<S2> &, const bool);
+ Vector<S1, ::dealii::MemorySpace::Host>::reinit<S2>(
+ const Vector<S2, ::dealii::MemorySpace::Host> &,
+ const bool);
template S1
- Vector<S1>::inner_product_local<S2>(const Vector<S2> &) const;
+ Vector<S1, ::dealii::MemorySpace::Host>::inner_product_local<S2>(
+ const Vector<S2, ::dealii::MemorySpace::Host> &) const;
template void
- Vector<S1>::copy_locally_owned_data_from<S2>(const Vector<S2> &);
+ Vector<S1, ::dealii::MemorySpace::Host>::copy_locally_owned_data_from<
+ S2>(const Vector<S2, ::dealii::MemorySpace::Host> &);
\}
\}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/read_write_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace LinearAlgebra
+{
+ template void
+ ReadWriteVector<float>::import(
+ const CUDAWrappers::Vector<float> &,
+ VectorOperation::values,
+ const std::shared_ptr<const CommunicationPatternBase> &);
+ template void
+ ReadWriteVector<float>::import(
+ const distributed::Vector<float, ::dealii::MemorySpace::CUDA> &,
+ VectorOperation::values,
+ const std::shared_ptr<const CommunicationPatternBase> &);
+
+ template void
+ ReadWriteVector<double>::import(
+ const CUDAWrappers::Vector<double> &,
+ VectorOperation::values,
+ const std::shared_ptr<const CommunicationPatternBase> &);
+ template void
+ ReadWriteVector<double>::import(
+ const distributed::Vector<double, ::dealii::MemorySpace::CUDA> &,
+ VectorOperation::values,
+ const std::shared_ptr<const CommunicationPatternBase> &);
+} // namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff