From: Daniel Arndt Date: Fri, 17 Aug 2018 13:13:36 +0000 (+0200) Subject: SolverMinRes X-Git-Tag: v9.1.0-rc1~790^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=07de41b4598bd153a82a4ea4a0f081d6bcf8cc31;p=dealii.git SolverMinRes --- diff --git a/include/deal.II/lac/cuda_vector.h b/include/deal.II/lac/cuda_vector.h index b6dfb083f0..78ee585668 100644 --- 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 &&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 + * swap(u,v) that simply calls u.swap(v), 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 &v); + /** * Reinit functionality. The flag omit_zeroing_entries * 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 +inline void +swap(LinearAlgebra::CUDAWrappers::Vector &u, + LinearAlgebra::CUDAWrappers::Vector &v) +{ + u.swap(v); +} - // ---------------------------- Inline functions -------------------------- +namespace LinearAlgebra +{ + namespace CUDAWrappers + { template inline Number * Vector::get_values() const @@ -351,6 +390,16 @@ namespace LinearAlgebra { return complete_index_set(n_elements); } + + + + template + inline void + Vector::swap(Vector &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 cda9155757..845928b743 100644 --- a/include/deal.II/lac/solver_minres.h +++ b/include/deal.II/lac/solver_minres.h @@ -329,13 +329,6 @@ SolverMinRes::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 index 0000000000..ea8561c29a --- /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 +#include + +#include +#include +#include +#include +#include +#include + +#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 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 sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + minres_host.solve(A, sol_host, rhs_host, prec_no); + + // Solve on the device + CUDAWrappers::SparseMatrix A_dev(cuda_handle, A); + LinearAlgebra::CUDAWrappers::Vector sol_dev(size); + LinearAlgebra::CUDAWrappers::Vector rhs_dev(size); + LinearAlgebra::ReadWriteVector rw_vector(size); + for (unsigned int i = 0; i < size; ++i) + rw_vector[i] = static_cast(i); + rhs_dev.import(rw_vector, VectorOperation::insert); + SolverMinRes> 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 index 0000000000..448f7326c7 --- /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