SolverMinRes

diff --git a/include/deal.II/lac/cuda_vector.h b/include/deal.II/lac/cuda_vector.h
index b6dfb083f024297c935856555d59e9272141d383..78ee5856680d5030bb36574ef99efa77f3212072 100644 (file)
--- a/include/deal.II/lac/cuda_vector.h
+++ b/include/deal.II/lac/cuda_vector.h
@@ -96,6 +96,24 @@ namespace LinearAlgebra
       Vector &
       operator=(Vector<Number> &&v) = default;
 
+      /**
+       * Swap the contents of this vector and the other vector @p v. One could do
+       * this operation with a temporary variable and copying over the data
+       * elements, but this function is significantly more efficient since it
+       * only swaps the pointers to the data of the two vectors and therefore
+       * does not need to allocate temporary storage and move data around.
+       *
+       * This function is analogous to the @p swap function of all C++
+       * standard containers. Also, there is a global function
+       * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+       * analogy to standard functions.
+       *
+       * This function is virtual in order to allow for derived classes to
+       * handle memory separately.
+       */
+      virtual void
+      swap(Vector<Number> &v);
+
       /**
        * Reinit functionality. The flag <tt>omit_zeroing_entries</tt>
        * determines whether the vector should be filled with zero (false) or
@@ -324,10 +342,31 @@ namespace LinearAlgebra
        */
       size_type n_elements;
     };
+  } // namespace CUDAWrappers
+} // namespace LinearAlgebra
 
+// ---------------------------- Inline functions --------------------------
 
+/**
+ * Global function @p swap which overloads the default implementation of the
+ * C++ standard library which uses a temporary object. The function simply
+ * exchanges the data of the two vectors.
+ *
+ * @relatesalso Vector
+ * @author Daniel Arndt, 2018
+ */
+template <typename Number>
+inline void
+swap(LinearAlgebra::CUDAWrappers::Vector<Number> &u,
+     LinearAlgebra::CUDAWrappers::Vector<Number> &v)
+{
+  u.swap(v);
+}
 
-    // ---------------------------- Inline functions --------------------------
+namespace LinearAlgebra
+{
+  namespace CUDAWrappers
+  {
     template <typename Number>
     inline Number *
     Vector<Number>::get_values() const
@@ -351,6 +390,16 @@ namespace LinearAlgebra
     {
       return complete_index_set(n_elements);
     }
+
+
+
+    template <typename Number>
+    inline void
+    Vector<Number>::swap(Vector<Number> &v)
+    {
+      std::swap(val, v.val);
+      std::swap(n_elements, v.n_elements);
+    }
   } // namespace CUDAWrappers
 } // namespace LinearAlgebra
 

diff --git a/include/deal.II/lac/solver_minres.h b/include/deal.II/lac/solver_minres.h
index cda91557577610f7c6a392c59a08a852126e8dec..845928b743bc2f4ab8fcc0a21120e95213f531ca 100644 (file)
--- a/include/deal.II/lac/solver_minres.h
+++ b/include/deal.II/lac/solver_minres.h
@@ -329,13 +329,6 @@ SolverMinRes<VectorType>::solve(const MatrixType &        A,
       // All vectors have to be shifted
       // one iteration step.
       // This should be changed one time.
-      //
-      // the previous code was like this:
-      //   m[2] = m[1];
-      //   m[1] = m[0];
-      // but it can be made more efficient,
-      // since the value of m[0] is no more
-      // needed in the next iteration
       swap(*m[2], *m[1]);
       swap(*m[1], *m[0]);
 

diff --git a/tests/cuda/solver_07.cu b/tests/cuda/solver_07.cu
new file mode 100644 (file)
index 0000000..ea8561c
--- /dev/null
+++ b/tests/cuda/solver_07.cu
@@ -0,0 +1,85 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Check that dealii::SolverMinRes works with CUDAWrappers::SparseMatrix
+
+#include <deal.II/base/cuda.h>
+#include <deal.II/base/exceptions.h>
+
+#include <deal.II/lac/cuda_sparse_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_minres.h>
+#include <deal.II/lac/vector.h>
+
+#include "../testmatrix.h"
+#include "../tests.h"
+
+void
+test(Utilities::CUDA::Handle &cuda_handle)
+{
+  // Build the sparse matrix on the host
+  const unsigned int   problem_size = 10;
+  unsigned int         size         = (problem_size - 1) * (problem_size - 1);
+  FDMatrix             testproblem(problem_size, problem_size);
+  SparsityPattern      structure(size, size, 5);
+  SparseMatrix<double> A;
+  testproblem.five_point_structure(structure);
+  structure.compress();
+  A.reinit(structure);
+  testproblem.five_point(A);
+
+  // Solve on the host
+  PreconditionIdentity prec_no;
+  SolverControl        control(100, 1.e-10);
+  SolverMinRes<>       minres_host(control);
+  Vector<double>       sol_host(size);
+  Vector<double>       rhs_host(size);
+  for (unsigned int i = 0; i < size; ++i)
+    rhs_host[i] = static_cast<double>(i);
+  minres_host.solve(A, sol_host, rhs_host, prec_no);
+
+  // Solve on the device
+  CUDAWrappers::SparseMatrix<double>          A_dev(cuda_handle, A);
+  LinearAlgebra::CUDAWrappers::Vector<double> sol_dev(size);
+  LinearAlgebra::CUDAWrappers::Vector<double> rhs_dev(size);
+  LinearAlgebra::ReadWriteVector<double>      rw_vector(size);
+  for (unsigned int i = 0; i < size; ++i)
+    rw_vector[i] = static_cast<double>(i);
+  rhs_dev.import(rw_vector, VectorOperation::insert);
+  SolverMinRes<LinearAlgebra::CUDAWrappers::Vector<double>> minres_dev(control);
+  minres_dev.solve(A_dev, sol_dev, rhs_dev, prec_no);
+
+  // Check the result
+  rw_vector.import(sol_dev, VectorOperation::insert);
+  for (unsigned int i = 0; i < size; ++i)
+    AssertThrow(std::fabs(rw_vector[i] - sol_host[i]) < 1e-8,
+                ExcInternalError());
+}
+
+int
+main()
+{
+  initlog();
+  deallog.depth_console(0);
+
+  Utilities::CUDA::Handle cuda_handle;
+  test(cuda_handle);
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}

diff --git a/tests/cuda/solver_07.output b/tests/cuda/solver_07.output
new file mode 100644 (file)
index 0000000..448f732
--- /dev/null
+++ b/tests/cuda/solver_07.output
@@ -0,0 +1,6 @@
+
+DEAL:minres::Starting value 416.989
+DEAL:minres::Convergence step 31 value 8.71122e-12
+DEAL:minres::Starting value 416.989
+DEAL:minres::Convergence step 31 value 8.71119e-12
+DEAL::OK