From 26bf11c1f126900d32fa0103190643e8263f794a Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Mon, 23 Apr 2018 18:38:13 -0400 Subject: [PATCH] Add test for direct solvers on the device --- tests/cuda/solver_02.cu | 109 ++++++++++++++++++++++++++++++++++++ tests/cuda/solver_02.output | 8 +++ 2 files changed, 117 insertions(+) create mode 100644 tests/cuda/solver_02.cu create mode 100644 tests/cuda/solver_02.output diff --git a/tests/cuda/solver_02.cu b/tests/cuda/solver_02.cu new file mode 100644 index 0000000000..296afe3549 --- /dev/null +++ b/tests/cuda/solver_02.cu @@ -0,0 +1,109 @@ +// --------------------------------------------------------------------- +// +// 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check that CUDA direct solvers work + +#include "../tests.h" +#include "../testmatrix.h" + +#include +#include +#include +#include +#include +#include +#include + + +void test(Utilities::CUDA::Handle &cuda_handle) +{ + // Create the matrix on the host. + dealii::SparsityPattern sparsity_pattern; + dealii::SparseMatrix matrix; + unsigned int const size = 30; + std::vector> column_indices(size); + for (unsigned int i = 0; i < size; ++i) + { + unsigned int j_max = std::min(size, i + 2); + unsigned int j_min = (i == 0) ? 0 : i - 1; + for (unsigned int j = j_min; j < j_max; ++j) + column_indices[i].emplace_back(j); + } + sparsity_pattern.copy_from(size, size, column_indices.begin(), + column_indices.end()); + matrix.reinit(sparsity_pattern); + for (unsigned int i = 0; i < size; ++i) + { + unsigned int j_max = std::min(size - 1, i + 1); + unsigned int j_min = (i == 0) ? 0 : i - 1; + matrix.set(i, j_min, -1.); + matrix.set(i, j_max, -1.); + matrix.set(i, i, 4.); + } + + // Generate a random solution and then compute the rhs + dealii::Vector sol_ref(size); + std::default_random_engine generator; + std::normal_distribution<> distribution(10., 2.); + for (auto &val : sol_ref) + val = distribution(generator); + + dealii::Vector rhs(size); + matrix.vmult(rhs, sol_ref); + + // Move the matrix and the rhs to the host + CUDAWrappers::SparseMatrix matrix_dev(cuda_handle, matrix); + + LinearAlgebra::CUDAWrappers::Vector rhs_dev(size); + LinearAlgebra::ReadWriteVector rhs_host(size); + std::copy(rhs.begin(), rhs.end(), rhs_host.begin()); + rhs_dev.import(rhs_host, VectorOperation::insert); + + LinearAlgebra::CUDAWrappers::Vector solution_dev(size); + + for (auto solver_type: {"Cholesky", "LU_dense", "LU_host"}) + { + // Solve on the device + CUDAWrappers::SolverDirect::AdditionalData data(solver_type); + SolverControl solver_control; + + CUDAWrappers::SolverDirect solver(cuda_handle, solver_control, + data); + solver.solve(matrix_dev, solution_dev, rhs_dev); + + // Move the result back to the host + LinearAlgebra::ReadWriteVector solution_host(size); + solution_host.import(solution_dev, VectorOperation::insert); + + // Check the result + for (unsigned int i = 0; i < size; ++i) + AssertThrow(std::abs(solution_host[i] - sol_ref[i]) < 1e-12, + ExcInternalError()); + std::cout<