From 7de4f6b5471ebd82cad778124aa7dc9bbf956b67 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 17 Aug 2018 15:18:23 +0200 Subject: [PATCH] SolverFGMRES, SolverBiCGStab, SolverRichardson, SolverRelaxation --- tests/cuda/solver_04.cu | 86 +++++++++++++++++++++++ tests/cuda/solver_04.output | 6 ++ tests/cuda/solver_05.cu | 85 +++++++++++++++++++++++ tests/cuda/solver_05.output | 6 ++ tests/cuda/solver_09.cu | 114 ++++++++++++++++++++++++++++++ tests/cuda/solver_09.output | 6 ++ tests/cuda/solver_10.cu | 134 ++++++++++++++++++++++++++++++++++++ tests/cuda/solver_10.output | 6 ++ 8 files changed, 443 insertions(+) create mode 100644 tests/cuda/solver_04.cu create mode 100644 tests/cuda/solver_04.output create mode 100644 tests/cuda/solver_05.cu create mode 100644 tests/cuda/solver_05.output create mode 100644 tests/cuda/solver_09.cu create mode 100644 tests/cuda/solver_09.output create mode 100644 tests/cuda/solver_10.cu create mode 100644 tests/cuda/solver_10.output diff --git a/tests/cuda/solver_04.cu b/tests/cuda/solver_04.cu new file mode 100644 index 0000000000..0aba76d39b --- /dev/null +++ b/tests/cuda/solver_04.cu @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// 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::SolverBicgstab 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); + SolverBicgstab<> bicgstab_host(control); + Vector sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + bicgstab_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); + SolverBicgstab> bicgstab_dev( + control); + bicgstab_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_04.output b/tests/cuda/solver_04.output new file mode 100644 index 0000000000..df448144a1 --- /dev/null +++ b/tests/cuda/solver_04.output @@ -0,0 +1,6 @@ + +DEAL:Bicgstab::Starting value 416.989 +DEAL:Bicgstab::Convergence step 23 value 1.39578e-11 +DEAL:Bicgstab::Starting value 416.989 +DEAL:Bicgstab::Convergence step 23 value 1.39578e-11 +DEAL::OK diff --git a/tests/cuda/solver_05.cu b/tests/cuda/solver_05.cu new file mode 100644 index 0000000000..1a4bb7f48a --- /dev/null +++ b/tests/cuda/solver_05.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::SolverFGMRES 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); + SolverFGMRES<> fgmres_host(control); + Vector sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + fgmres_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); + SolverFGMRES> fgmres_dev(control); + fgmres_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_05.output b/tests/cuda/solver_05.output new file mode 100644 index 0000000000..792aee14c9 --- /dev/null +++ b/tests/cuda/solver_05.output @@ -0,0 +1,6 @@ + +DEAL:FGMRES::Starting value 416.989 +DEAL:FGMRES::Convergence step 31 value 3.31663e-11 +DEAL:FGMRES::Starting value 416.989 +DEAL:FGMRES::Convergence step 31 value 3.31763e-11 +DEAL::OK diff --git a/tests/cuda/solver_09.cu b/tests/cuda/solver_09.cu new file mode 100644 index 0000000000..742924bcbf --- /dev/null +++ b/tests/cuda/solver_09.cu @@ -0,0 +1,114 @@ +// --------------------------------------------------------------------- +// +// 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::SolverRichardson works with CUDAWrappers::SparseMatrix + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../testmatrix.h" +#include "../tests.h" + +template +class PreconditionOperator +{ +public: + PreconditionOperator(const MatrixType &inverse_diagonal_matrix_) + : inverse_diagonal_matrix(inverse_diagonal_matrix_) + {} + + template + void + vmult(VectorType &u, const VectorType &v) const + { + inverse_diagonal_matrix.vmult(u, v); + } + +private: + const MatrixType &inverse_diagonal_matrix; +}; + + +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); + SparseMatrix A_diagonal_inverse; + A_diagonal_inverse.reinit(structure); + for (unsigned int i = 0; i < size; ++i) + A_diagonal_inverse(i, i) = 1. / A(i, i); + + // Solve on the host + PreconditionOperator> preconditioner(A_diagonal_inverse); + SolverControl control(1000, 1.e-2); + SolverRichardson<> richardson_host(control); + Vector sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + richardson_host.solve(A, sol_host, rhs_host, preconditioner); + + // Solve on the device + CUDAWrappers::SparseMatrix A_dev(cuda_handle, A); + CUDAWrappers::SparseMatrix A_diagonal_inverse_dev(cuda_handle, + A_diagonal_inverse); + PreconditionOperator> preconditioner_dev( + A_diagonal_inverse_dev); + 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); + SolverRichardson> richardson_dev( + control); + richardson_dev.solve(A_dev, sol_dev, rhs_dev, preconditioner_dev); + + // 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_09.output b/tests/cuda/solver_09.output new file mode 100644 index 0000000000..c4357dd61c --- /dev/null +++ b/tests/cuda/solver_09.output @@ -0,0 +1,6 @@ + +DEAL:Richardson::Starting value 416.989 +DEAL:Richardson::Convergence step 207 value 0.00982595 +DEAL:Richardson::Starting value 416.989 +DEAL:Richardson::Convergence step 207 value 0.00982595 +DEAL::OK diff --git a/tests/cuda/solver_10.cu b/tests/cuda/solver_10.cu new file mode 100644 index 0000000000..f38f500402 --- /dev/null +++ b/tests/cuda/solver_10.cu @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// +// 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::SolverRelaxation works with CUDAWrappers::SparseMatrix + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../testmatrix.h" +#include "../tests.h" + + +template +class RelaxationOperator +{ +public: + RelaxationOperator(const MatrixType &system_matrix_, + const MatrixType &inverse_diagonal_matrix_) + : system_matrix(system_matrix_) + , inverse_diagonal_matrix(inverse_diagonal_matrix_) + {} + + template + void + step(VectorType &u, const VectorType &v) const + { + // u = u - omega*inverse_diagonal_matrix*(system_matrix*u-v) + const double omega = 1.; + VectorType tmp_1(v.size()); + system_matrix.vmult(tmp_1, u); + tmp_1 -= v; + VectorType tmp_2(u.size()); + inverse_diagonal_matrix.vmult(tmp_2, tmp_1); + tmp_2 *= omega; + u -= tmp_2; + } + + template + void + Tstep(VectorType &u, const VectorType &v) const + { + AssertThrow(false, ExcNotImplemented()); + } + +private: + const MatrixType &system_matrix; + const MatrixType &inverse_diagonal_matrix; +}; + + +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); + SparseMatrix A_diagonal_inverse; + A_diagonal_inverse.reinit(structure); + for (unsigned int i = 0; i < size; ++i) + A_diagonal_inverse(i, i) = 1. / A(i, i); + + // Solve on the host + RelaxationOperator> relaxation_operator( + A, A_diagonal_inverse); + SolverControl control(1000, 1.e-3); + SolverRelaxation<> relaxation_host(control); + Vector sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + relaxation_host.solve(A, sol_host, rhs_host, relaxation_operator); + + // Solve on the device + CUDAWrappers::SparseMatrix A_dev(cuda_handle, A); + CUDAWrappers::SparseMatrix A_diagonal_inverse_dev(cuda_handle, + A_diagonal_inverse); + RelaxationOperator> + relaxation_operator_dev(A_dev, A_diagonal_inverse_dev); + 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); + SolverRelaxation> relaxation_dev( + control); + relaxation_dev.solve(A_dev, sol_dev, rhs_dev, relaxation_operator_dev); + + // 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(10); + + Utilities::CUDA::Handle cuda_handle; + test(cuda_handle); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/cuda/solver_10.output b/tests/cuda/solver_10.output new file mode 100644 index 0000000000..d46239fd5e --- /dev/null +++ b/tests/cuda/solver_10.output @@ -0,0 +1,6 @@ + +DEAL:Relaxation::Starting value 416.989 +DEAL:Relaxation::Convergence step 253 value 0.000976934 +DEAL:Relaxation::Starting value 416.989 +DEAL:Relaxation::Convergence step 253 value 0.000976934 +DEAL::OK -- 2.39.5