]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add support for CUDA memory space in LA::distributed::Vector
authorBruno Turcksin <bruno.turcksin@gmail.com>
Thu, 2 Aug 2018 14:51:24 +0000 (14:51 +0000)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 22 Aug 2018 13:18:03 +0000 (13:18 +0000)
Add support for CUDA memory space to LinearAlgebra::distributed::Vector by
adding a new template parameter with a default value of Host.

16 files changed:
include/deal.II/base/memory_space.h [new file with mode: 0644]
include/deal.II/base/types.h
include/deal.II/lac/la_parallel_block_vector.templates.h
include/deal.II/lac/la_parallel_vector.h
include/deal.II/lac/la_parallel_vector.templates.h
include/deal.II/lac/precondition.h
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/vector_operations_internal.h
include/deal.II/matrix_free/fe_evaluation.h
source/base/cuda.cu
source/lac/CMakeLists.txt
source/lac/la_parallel_vector.cc
source/lac/la_parallel_vector.cu [new file with mode: 0644]
source/lac/la_parallel_vector.inst.in
source/lac/read_write_vector.cu [new file with mode: 0644]

diff --git a/include/deal.II/base/memory_space.h b/include/deal.II/base/memory_space.h
new file mode 100644 (file)
index 0000000..f67b50c
--- /dev/null
@@ -0,0 +1,222 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_memory_space_h
+#define dealii_memory_space_h
+
+#include <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
index a6724db1331595d6a3cd1707154bbc3be47ac12f..5d94c7d8026ebacedfb0c49bb6ae0843ab259246 100644 (file)
@@ -276,7 +276,6 @@ namespace numbers
     static_cast<types::subdomain_id>(-2);
 } // namespace numbers
 
-
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index fecbee6055f25922e87b42be8cec30158f70b000..c1fb7d0f1786b1045d9f7faf834554d35f43dadc 100644 (file)
@@ -573,7 +573,7 @@ namespace LinearAlgebra
       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(
index 2467590a27f99f75af01149365e595d3c85ffa36..b7ea75208542a157c697f7c86d77beb5968f1b1d 100644 (file)
@@ -18,6 +18,7 @@
 
 #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>
@@ -175,13 +176,54 @@ namespace LinearAlgebra
      * fail in some circumstances. Therefore, it is strongly recommended to
      * not rely on this class to automatically detect the unsupported case.
      *
-     * @author Katharina Kormann, Martin Kronbichler, 2010, 2011, 2016
+     * <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 *;
@@ -192,6 +234,11 @@ namespace LinearAlgebra
       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
        */
@@ -204,7 +251,7 @@ namespace LinearAlgebra
       /**
        * 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
@@ -270,8 +317,8 @@ namespace LinearAlgebra
        */
       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
@@ -326,7 +373,7 @@ namespace LinearAlgebra
        * handle memory separately.
        */
       void
-      swap(Vector<Number> &v);
+      swap(Vector<Number, MemorySpace> &v);
 
       /**
        * Assigns the vector to the parallel partitioning of the input vector
@@ -339,8 +386,8 @@ namespace LinearAlgebra
        * at all, the vector will also update its ghost values in analogy to
        * the respective setting the Trilinos and PETSc vectors.
        */
-      Vector<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
@@ -354,8 +401,8 @@ namespace LinearAlgebra
        * 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
       /**
@@ -369,7 +416,7 @@ namespace LinearAlgebra
        * ReadWriteVector instead.
        */
       DEAL_II_DEPRECATED
-      Vector<Number> &
+      Vector<Number, MemorySpace> &
       operator=(const PETScWrappers::MPI::Vector &petsc_vec);
 #endif
 
@@ -386,7 +433,7 @@ namespace LinearAlgebra
        * ReadWriteVector instead.
        */
       DEAL_II_DEPRECATED
-      Vector<Number> &
+      Vector<Number, MemorySpace> &
       operator=(const TrilinosWrappers::MPI::Vector &trilinos_vec);
 #endif
       //@}
@@ -564,7 +611,7 @@ namespace LinearAlgebra
        */
       template <typename Number2>
       void
-      copy_locally_owned_data_from(const Vector<Number2> &src);
+      copy_locally_owned_data_from(const Vector<Number2, MemorySpace> &src);
 
       //@}
 
@@ -584,25 +631,25 @@ namespace LinearAlgebra
       /**
        * 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;
 
       /**
@@ -612,6 +659,9 @@ namespace LinearAlgebra
        * current elements. The last parameter can be used if the same
        * communication pattern is used multiple times. This can be used to
        * improve performance.
+       *
+       * @note If the MemorySpace is CUDA, the data in the ReadWriteVector will
+       * be moved to the device.
        */
       virtual void
       import(
@@ -777,7 +827,7 @@ namespace LinearAlgebra
        * zero, also ghost elements are set to zero, otherwise they remain
        * unchanged.
        */
-      virtual Vector<Number> &
+      virtual Vector<Number, MemorySpace> &
       operator=(const Number s) override;
 
       /**
@@ -804,7 +854,7 @@ namespace LinearAlgebra
        * 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.
@@ -813,11 +863,11 @@ namespace LinearAlgebra
        */
       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>.
@@ -826,10 +876,10 @@ namespace LinearAlgebra
        */
       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);
 
       //@}
 
@@ -904,6 +954,9 @@ namespace LinearAlgebra
        * 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();
@@ -911,6 +964,9 @@ namespace LinearAlgebra
       /**
        * Return constant iterator to the start of the locally owned elements
        * of the vector.
+       *
+       * @note For the CUDA memory space, the iterator points to memory on the
+       * device.
        */
       const_iterator
       begin() const;
@@ -918,6 +974,9 @@ namespace LinearAlgebra
       /**
        * Return an iterator pointing to the element past the end of the array
        * of locally owned entries.
+       *
+       * @note For the CUDA memory space, the iterator points to memory on the
+       * device.
        */
       iterator
       end();
@@ -925,6 +984,9 @@ namespace LinearAlgebra
       /**
        * Return a constant iterator pointing to the element past the end of
        * the array of the locally owned entries.
+       *
+       * @note For the CUDA memory space, the iterator points to memory on the
+       * device.
        */
       const_iterator
       end() const;
@@ -992,6 +1054,15 @@ namespace LinearAlgebra
       Number &
       local_element(const size_type local_index);
 
+      /**
+       * Return the pointer to the underlying raw array.
+       *
+       * @note For the CUDA memory space, the pointer points to memory on the
+       * device.
+       */
+      Number *
+      get_values() const;
+
       /**
        * Instead of getting individual elements of a vector via operator(),
        * this function allows getting a whole set of elements at once. The
@@ -1006,6 +1077,8 @@ namespace LinearAlgebra
        * @endcode
        *
        * @pre The sizes of the @p indices and @p values arrays must be identical.
+       *
+       * @note This function is not implemented for CUDA memory space.
        */
       template <typename OtherNumber>
       void
@@ -1125,6 +1198,13 @@ namespace LinearAlgebra
        */
       DeclException0(ExcVectorTypeNotCompatible);
 
+      /**
+       * Attempt to perform an operation not implemented on the device.
+       *
+       * @ingroup Exceptions
+       */
+      DeclException0(ExcNotAllowedForCuda);
+
       /**
        * Exception
        */
@@ -1171,18 +1251,12 @@ namespace LinearAlgebra
                  const Number                     a,
                  const VectorSpaceVector<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().
@@ -1220,9 +1294,9 @@ namespace LinearAlgebra
        * 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
@@ -1237,13 +1311,9 @@ namespace LinearAlgebra
       size_type allocated_size;
 
       /**
-       * Pointer to the array of local elements of this vector.
-       *
-       * Because we allocate these arrays via Utilities::System::posix_memalign,
-       * we need to use a custom deleter for this object that does not call
-       * <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
@@ -1309,7 +1379,7 @@ namespace LinearAlgebra
       /*
        * Make all other vector types friends.
        */
-      template <typename Number2>
+      template <typename Number2, typename MemorySpace2>
       friend class Vector;
 
       /**
@@ -1325,56 +1395,142 @@ namespace LinearAlgebra
 
 #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());
 
@@ -1386,73 +1542,79 @@ namespace LinearAlgebra
 
 
 
-    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),
@@ -1465,15 +1627,18 @@ namespace LinearAlgebra
                vector_is_ghosted == true,
              ExcMessage("You tried to read a ghost element of this vector, "
                         "but it has not imported its ghost values."));
-      return values[partitioner->global_to_local(global_index)];
+      return data.values[partitioner->global_to_local(global_index)];
     }
 
 
 
-    template <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),
@@ -1487,31 +1652,37 @@ namespace LinearAlgebra
       // (then, the compiler picks this method according to the C++ rule book
       // even if a human would pick the const method when this subsequent use
       // is just a read)
-      return values[partitioner->global_to_local(global_index)];
+      return data.values[partitioner->global_to_local(global_index)];
     }
 
 
 
-    template <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());
@@ -1519,28 +1690,44 @@ namespace LinearAlgebra
       Assert(local_index < local_size() || vector_is_ghosted == true,
              ExcMessage("You tried to read a ghost element of this vector, "
                         "but it has not imported its ghost values."));
-      return values[local_index];
+
+      return data.values[local_index];
     }
 
 
 
-    template <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]);
@@ -1548,12 +1735,13 @@ namespace LinearAlgebra
 
 
 
-    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)
         {
@@ -1565,11 +1753,12 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1584,12 +1773,12 @@ namespace LinearAlgebra
 
 
 
-    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)
         {
@@ -1603,18 +1792,18 @@ namespace LinearAlgebra
 
 
 
-    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;
     }
@@ -1633,10 +1822,10 @@ namespace LinearAlgebra
  * @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);
 }
@@ -1647,12 +1836,13 @@ swap(LinearAlgebra::distributed::Vector<Number> &u,
  *
  * @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
index 6619314f9c7d741d9c7aaa75ed77a35343a00d23..a476d61b9e6a10abd0aae02aea8dab44df462f22 100644 (file)
@@ -19,6 +19,8 @@
 
 #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>
@@ -40,7 +42,7 @@ namespace LinearAlgebra
     {
       // In the import_from_ghosted_array_finish we might need to calculate the
       // maximal and minimal value for the given number type, which is not
-      // straight forward for complex numbers. Therefore, comparison of complex
+      // straightforward for complex numbers. Therefore, comparison of complex
       // numbers is prohibited and throws an exception.
       template <typename Number>
       Number
@@ -75,13 +77,304 @@ namespace LinearAlgebra
                                "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++)
@@ -101,39 +394,23 @@ namespace LinearAlgebra
 
 
 
-    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();
 
@@ -155,11 +432,11 @@ namespace LinearAlgebra
 
 
 
-    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());
@@ -192,11 +469,11 @@ namespace LinearAlgebra
 
 
 
-    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(
@@ -208,10 +485,10 @@ namespace LinearAlgebra
 
 
 
-    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(
@@ -221,9 +498,9 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -249,22 +526,20 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -274,21 +549,19 @@ namespace LinearAlgebra
       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);
@@ -296,11 +569,10 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -308,10 +580,9 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -319,11 +590,10 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -331,8 +601,8 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    inline Vector<Number>::~Vector()
+    template <typename Number, typename MemorySpace>
+    inline Vector<Number, MemorySpace>::~Vector()
     {
       try
         {
@@ -344,9 +614,9 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -357,10 +627,11 @@ namespace LinearAlgebra
 
 
 
-    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());
 
@@ -411,10 +682,9 @@ namespace LinearAlgebra
       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)
@@ -426,18 +696,21 @@ namespace LinearAlgebra
 
 
 
-    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);
         }
     }
 
@@ -475,9 +748,10 @@ namespace LinearAlgebra
       }
     } // 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
@@ -517,9 +791,10 @@ namespace LinearAlgebra
 
 #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();
@@ -541,9 +816,10 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -551,9 +827,9 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -561,23 +837,35 @@ namespace LinearAlgebra
 
 
 
-    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;
@@ -593,10 +881,31 @@ namespace LinearAlgebra
         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);
@@ -605,9 +914,10 @@ namespace LinearAlgebra
 
 
 
-    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;
@@ -625,10 +935,26 @@ namespace LinearAlgebra
         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
@@ -636,9 +962,10 @@ namespace LinearAlgebra
 
 
 
-    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.
@@ -654,11 +981,34 @@ namespace LinearAlgebra
         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);
 
@@ -669,9 +1019,9 @@ namespace LinearAlgebra
 
 
 
-    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
@@ -685,19 +1035,34 @@ namespace LinearAlgebra
           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)
@@ -725,7 +1090,8 @@ namespace LinearAlgebra
                                  "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
@@ -774,12 +1140,12 @@ namespace LinearAlgebra
         }
       tmp_vector.compress(operation);
 
-      std::copy(tmp_vector.begin(), tmp_vector.end(), begin());
+      data.copy_from(tmp_vector.begin(), local_size());
     }
 
-    template <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
 
@@ -822,25 +1188,25 @@ namespace LinearAlgebra
       std::swap(partitioner, v.partitioner);
       std::swap(thread_loop_partitioner, v.thread_loop_partitioner);
       std::swap(allocated_size, v.allocated_size);
-      std::swap(values, v.values);
+      std::swap(data, v.data);
       std::swap(import_data, v.import_data);
       std::swap(vector_is_ghosted, v.vector_is_ghosted);
     }
 
 
 
-    template <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
@@ -853,36 +1219,37 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -892,21 +1259,21 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -916,16 +1283,15 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -933,15 +1299,16 @@ namespace LinearAlgebra
 
 
 
-    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());
@@ -950,17 +1317,17 @@ namespace LinearAlgebra
       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);
 
@@ -970,20 +1337,21 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -991,10 +1359,15 @@ namespace LinearAlgebra
       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();
@@ -1002,10 +1375,10 @@ namespace LinearAlgebra
 
 
 
-    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)
         {
@@ -1015,17 +1388,17 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -1033,34 +1406,39 @@ namespace LinearAlgebra
 
 
 
-    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);
 
@@ -1070,13 +1448,13 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -1085,10 +1463,16 @@ namespace LinearAlgebra
       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();
@@ -1096,16 +1480,15 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -1115,9 +1498,9 @@ namespace LinearAlgebra
 
 
 
-    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;
@@ -1125,21 +1508,23 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -1147,22 +1532,24 @@ namespace LinearAlgebra
 
 
 
-    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();
@@ -1170,12 +1557,12 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -1183,10 +1570,15 @@ namespace LinearAlgebra
       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();
@@ -1194,67 +1586,44 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1266,15 +1635,17 @@ namespace LinearAlgebra
 
 
 
-    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;
@@ -1282,28 +1653,27 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1317,24 +1687,24 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1346,9 +1716,9 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1360,32 +1730,33 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1399,24 +1770,24 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1428,40 +1799,43 @@ namespace LinearAlgebra
 
 
 
-    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)
@@ -1473,9 +1847,9 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -1483,9 +1857,9 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -1493,9 +1867,9 @@ namespace LinearAlgebra
 
 
 
-    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);
@@ -1514,12 +1888,12 @@ namespace LinearAlgebra
 
 
 
-    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());
@@ -1544,6 +1918,9 @@ namespace LinearAlgebra
           }
 #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
@@ -1551,10 +1928,10 @@ namespace LinearAlgebra
           << "Vector data:" << std::endl;
       if (across)
         for (size_type i = 0; i < partitioner->local_size(); ++i)
-          out << local_element(i) << ' ';
+          out << stored_elements[i] << ' ';
       else
         for (size_type i = 0; i < partitioner->local_size(); ++i)
-          out << local_element(i) << std::endl;
+          out << stored_elements[i] << std::endl;
       out << std::endl;
 
       if (vector_is_ghosted)
@@ -1563,13 +1940,13 @@ namespace LinearAlgebra
           if (across)
             for (size_type i = 0; i < partitioner->n_ghost_indices(); ++i)
               out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
-                  << '/' << local_element(partitioner->local_size() + i)
+                  << '/' << stored_elements[partitioner->local_size() + i]
                   << ") ";
           else
             for (size_type i = 0; i < partitioner->n_ghost_indices(); ++i)
               out << '(' << partitioner->ghost_indices().nth_index_in_set(i)
-                  << '/' << local_element(partitioner->local_size() + i) << ")"
-                  << std::endl;
+                  << '/' << stored_elements[partitioner->local_size() + i]
+                  << ")" << std::endl;
           out << std::endl;
         }
       out << std::flush;
@@ -1596,8 +1973,7 @@ namespace LinearAlgebra
       out.precision(old_precision);
     }
 
-  } // end of namespace distributed
-
+  } // namespace distributed
 } // namespace LinearAlgebra
 
 
index 9800cb55a695b7f8d47bca17af6663d461a710c8..96f2e06aef52296a8fd3ec44cb50a84b4918073f 100644 (file)
@@ -20,6 +20,7 @@
 
 #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>
@@ -42,9 +43,9 @@ namespace LinearAlgebra
 {
   namespace distributed
   {
-    template <typename number>
+    template <typename, typename>
     class Vector;
-  }
+  } // namespace distributed
 } // namespace LinearAlgebra
 
 
@@ -1915,18 +1916,19 @@ namespace internal
     }
 
     // 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(),
@@ -2012,10 +2014,10 @@ namespace internal
       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
@@ -2222,9 +2224,16 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
       (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<
index 8cc7b8263047c54eb692610069896c1c737a9236..9e67f981787d07044c2b59161d1b9bad600bd79a 100644 (file)
@@ -47,9 +47,9 @@ namespace LinearAlgebra
   class CommunicationPatternBase;
   namespace distributed
   {
-    template <typename>
+    template <typename, typename>
     class Vector;
-  }
+  } // namespace distributed
 } // namespace LinearAlgebra
 
 #ifdef DEAL_II_WITH_PETSC
@@ -290,9 +290,10 @@ namespace LinearAlgebra
      * be used if the same communication pattern is used multiple times. This
      * can be used to improve performance.
      */
+    template <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>());
index 664d468040cbc46f863d8d498a78a8604ef5991d..362ac1f1a5d6abfd0d001c303edfbc5d6978a187 100644 (file)
 #  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
@@ -51,6 +49,163 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace LinearAlgebra
 {
+  namespace internal
+  {
+    // In the import_from_ghosted_array_finish we need to calculate the
+    // maximal and minimal value for the given number type, which is not
+    // straight forward for complex numbers. Therefore, comparison of complex
+    // numbers is prohibited and throws an assert.
+    template <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)
@@ -178,10 +333,10 @@ namespace LinearAlgebra
   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);
   }
 
 
@@ -250,52 +405,14 @@ namespace LinearAlgebra
     return *this;
   }
 
-  namespace internal
-  {
-    // In the import_from_ghosted_array_finish we need to calculate the maximal
-    // and minimal value for the given number type, which is not straight
-    // forward for complex numbers. Therefore, comparison of complex numbers is
-    // prohibited and throws an assert.
-    template <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)
   {
@@ -318,28 +435,10 @@ namespace LinearAlgebra
                     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);
   }
 
 
@@ -570,7 +669,7 @@ namespace LinearAlgebra
 
 
 
-#if defined(DEAL_II_WITH_CUDA)
+#ifdef DEAL_II_COMPILER_CUDA_AWARE
   template <typename Number>
   void
   ReadWriteVector<Number>::import(
index 9fd318574f84d465f2d60a821b3f1f2b2d6aa7ba..9d629d327958c16cef9fb0798e1236a43a24a969 100644 (file)
 #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>
 
@@ -307,7 +312,7 @@ namespace internal
     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)
@@ -329,15 +334,18 @@ namespace internal
           }
       }
 
-      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)
@@ -360,16 +368,16 @@ namespace internal
           }
       }
 
-      Number *val;
-      Number *v_val;
-      Number  a;
-      Number  x;
+      Number *      val;
+      const Number *v_val;
+      Number        a;
+      Number        x;
     };
 
     template <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)
       {}
@@ -390,8 +398,8 @@ namespace internal
           }
       }
 
-      Number *val;
-      Number *v_val;
+      Number *      val;
+      const Number *v_val;
     };
 
     template <typename Number>
@@ -425,7 +433,7 @@ namespace internal
     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)
       {}
@@ -446,18 +454,18 @@ namespace internal
           }
       }
 
-      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)
@@ -481,17 +489,17 @@ namespace internal
           }
       }
 
-      Number *val;
-      Number *v_val;
-      Number *w_val;
-      Number  a;
-      Number  b;
+      Number *      val;
+      const Number *v_val;
+      const Number *w_val;
+      Number        a;
+      Number        b;
     };
 
     template <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)
@@ -513,20 +521,20 @@ namespace internal
           }
       }
 
-      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)
@@ -551,18 +559,18 @@ namespace internal
           }
       }
 
-      Number *val;
-      Number *v_val;
-      Number *w_val;
-      Number  x;
-      Number  a;
-      Number  b;
+      Number *      val;
+      const Number *v_val;
+      const Number *w_val;
+      Number        x;
+      Number        a;
+      Number        b;
     };
 
     template <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)
       {}
@@ -583,14 +591,14 @@ namespace internal
           }
       }
 
-      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)
@@ -612,19 +620,19 @@ namespace internal
           }
       }
 
-      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)
@@ -648,11 +656,11 @@ namespace internal
           }
       }
 
-      Number *val;
-      Number *u_val;
-      Number *v_val;
-      Number  a;
-      Number  b;
+      Number *      val;
+      const Number *u_val;
+      const Number *v_val;
+      Number        a;
+      Number        b;
     };
 
     template <typename Number>
@@ -1385,6 +1393,1060 @@ namespace internal
       (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
 
index 1459d498d88f2766f36af2b54b79cf093b81e8f0..54d3ff65004516c6a343f03acb6156cbfeccffae 100644 (file)
@@ -45,7 +45,7 @@ namespace LinearAlgebra
 {
   namespace distributed
   {
-    template <typename>
+    template <typename, typename>
     class Vector;
   }
 } // namespace LinearAlgebra
index ac746b62e210b56998eb1673ed96ad39866708ac..19a7ea64b514b3f6e5c5b536f47444bec99d9e4f 100644 (file)
 
 #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
 
@@ -47,6 +50,13 @@ namespace Utilities
       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);
index 0abdfcba5d7ad7479752680581bf0899fe0dfbac..27729b8e89f7fd0af825650a646a848e1774c4ac 100644 (file)
@@ -150,6 +150,8 @@ IF(DEAL_II_WITH_CUDA)
     cuda_solver_direct.cu
     cuda_sparse_matrix.cu
     cuda_vector.cu
+    la_parallel_vector.cu
+    read_write_vector.cu
   )
 ENDIF()
 
index 0ef8028af40687a07d6edb96b24c4f2fcba73095..b56ed27e328dfdf9a63e27619c2ab0559d978029 100644 (file)
@@ -29,8 +29,10 @@ namespace LinearAlgebra
 {
   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);
diff --git a/source/lac/la_parallel_vector.cu b/source/lac/la_parallel_vector.cu
new file mode 100644 (file)
index 0000000..352cf2c
--- /dev/null
@@ -0,0 +1,32 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <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
index 856cdcf4c402e84a2c434dec8c72e6630ead1253..12c85c8a2a0ecfd4ca3b67707a0c373547adf55e 100644 (file)
@@ -21,7 +21,7 @@ for (SCALAR : REAL_AND_COMPLEX_SCALARS)
     \{
       namespace distributed
       \{
-        template class Vector<SCALAR>;
+        template class Vector<SCALAR, ::dealii::MemorySpace::Host>;
       \}
     \}
   }
@@ -33,11 +33,15 @@ for (S1 : REAL_AND_COMPLEX_SCALARS; S2 : REAL_SCALARS)
       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> &);
       \}
     \}
   }
@@ -50,11 +54,15 @@ for (S1, S2 : COMPLEX_SCALARS)
       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> &);
       \}
     \}
   }
diff --git a/source/lac/read_write_vector.cu b/source/lac/read_write_vector.cu
new file mode 100644 (file)
index 0000000..dd39d52
--- /dev/null
@@ -0,0 +1,47 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <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

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.