From 26bf11c1f126900d32fa0103190643e8263f794a Mon Sep 17 00:00:00 2001
From: Bruno Turcksin <bruno.turcksin@gmail.com>
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 <deal.II/base/cuda.h>
+#include <deal.II/lac/cuda_solver_direct.h>
+#include <deal.II/lac/cuda_sparse_matrix.h>
+#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/vector.h>
+
+
+void test(Utilities::CUDA::Handle &cuda_handle)
+{
+  // Create the matrix on the host.
+  dealii::SparsityPattern sparsity_pattern;
+  dealii::SparseMatrix<double> matrix;
+  unsigned int const size = 30;
+  std::vector<std::vector<unsigned int>> 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<double> sol_ref(size);
+  std::default_random_engine generator;
+  std::normal_distribution<> distribution(10., 2.);
+  for (auto &val : sol_ref)
+    val = distribution(generator);
+
+  dealii::Vector<double> rhs(size);
+  matrix.vmult(rhs, sol_ref);
+
+  // Move the matrix and the rhs to the host
+  CUDAWrappers::SparseMatrix<double> matrix_dev(cuda_handle, matrix);
+
+  LinearAlgebra::CUDAWrappers::Vector<double> rhs_dev(size);
+  LinearAlgebra::ReadWriteVector<double> rhs_host(size);
+  std::copy(rhs.begin(), rhs.end(), rhs_host.begin());
+  rhs_dev.import(rhs_host, VectorOperation::insert);
+
+  LinearAlgebra::CUDAWrappers::Vector<double> solution_dev(size);
+
+  for (auto solver_type: {"Cholesky", "LU_dense", "LU_host"})
+    {
+      // Solve on the device
+      CUDAWrappers::SolverDirect<double>::AdditionalData data(solver_type);
+      SolverControl solver_control;
+
+      CUDAWrappers::SolverDirect<double> solver(cuda_handle, solver_control,
+                                                data);
+      solver.solve(matrix_dev, solution_dev, rhs_dev);
+
+      // Move the result back to the host
+      LinearAlgebra::ReadWriteVector<double> 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<<solver_type<<std::endl;
+    }
+}
+
+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_02.output b/tests/cuda/solver_02.output
new file mode 100644
index 0000000000..d0ba62aa9b
--- /dev/null
+++ b/tests/cuda/solver_02.output
@@ -0,0 +1,8 @@
+
+DEAL::Starting value 0.00000
+DEAL::Convergence step 0 value 0.00000
+DEAL::Starting value 0.00000
+DEAL::Convergence step 0 value 0.00000
+DEAL:cg::Starting value 0.000
+DEAL::Convergence step 0 value 0.00000
+DEAL:cg::Starting value 0.000
-- 
2.39.5