]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test for CUDA support of matrix-free.
authorBruno Turcksin <bruno.turcksin@gmail.com>
Sun, 14 May 2017 03:05:05 +0000 (23:05 -0400)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Mon, 15 May 2017 21:18:06 +0000 (17:18 -0400)
tests/cuda/cuda_evaluate_1d_shape.cu [new file with mode: 0644]
tests/cuda/cuda_evaluate_1d_shape.output [new file with mode: 0644]
tests/cuda/cuda_evaluate_2d_shape.cu [new file with mode: 0644]
tests/cuda/cuda_evaluate_2d_shape.output [new file with mode: 0644]
tests/cuda/matrix_free_matrix_vector_01.cu [new file with mode: 0644]
tests/cuda/matrix_free_matrix_vector_01.output [new file with mode: 0644]
tests/cuda/matrix_free_matrix_vector_04.cu [new file with mode: 0644]
tests/cuda/matrix_free_matrix_vector_04.output [new file with mode: 0644]
tests/cuda/matrix_vector_common.cuh [new file with mode: 0644]
tests/cuda/matrix_vector_mf.cuh [new file with mode: 0644]

diff --git a/tests/cuda/cuda_evaluate_1d_shape.cu b/tests/cuda/cuda_evaluate_1d_shape.cu
new file mode 100644 (file)
index 0000000..040a2fe
--- /dev/null
@@ -0,0 +1,193 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the 1d evaluation functions used in
+// CUDAWrappers::FEEvaluation. These functions are marked 'internal' but it is
+// much easier to check their correctness directly rather than from the results
+// in dependent functions
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/matrix_free/cuda_fe_evaluation.cuh>
+
+namespace CUDA = LinearAlgebra::CUDAWrappers;
+
+template <int M, int N, int type, bool add, bool dof_to_quad>
+__global__ void evaluate_tensor_product(double *dst,
+                                        double *src)
+{
+  CUDAWrappers::internal::EvaluatorTensorProduct<
+  CUDAWrappers::internal::evaluate_general,1,M-1,N,double> evaluator;
+
+  if (type == 0)
+    evaluator.template values<0,dof_to_quad,add,false> (src, dst);
+  if (type == 1)
+    evaluator.template gradients<0,dof_to_quad,add,false> (src, dst);
+
+}
+
+template <int M, int N, int type, bool add>
+void test()
+{
+  deallog << "Test " << M << " x " << N << std::endl;
+  LinearAlgebra::ReadWriteVector<double> shape_host(M*N);
+  for (unsigned int i=0; i<(M+1)/2; ++i)
+    for (unsigned int j=0; j<N; ++j)
+      {
+        shape_host[i*N+j] = -1. + 2. *
+                            static_cast<double>(Testing::rand())/RAND_MAX;
+        if (type == 1)
+          shape_host[(M-1-i)*N+N-1-j] = -shape_host[i*N+j];
+        else
+          shape_host[(M-1-i)*N+N-1-j] = shape_host[i*N+j];
+      }
+  if (type == 0 && M%2 == 1 && N%2 == 1)
+    {
+      for (unsigned int i=0; i<M; ++i)
+        shape_host[i*N+N/2] = 0.;
+      shape_host[M/2*N+N/2] = 1.;
+    }
+  if (type == 1 && M%2 == 1 && N%2 == 1)
+    shape_host[M/2*N+N/2] = 0.;
+
+  LinearAlgebra::ReadWriteVector<double> x_host(N), x_ref(N), y_host(M), y_ref(M);
+  for (unsigned int i=0; i<N; ++i)
+    x_host[i] = static_cast<double>(Testing::rand())/RAND_MAX;
+
+  // Compute reference
+  for (unsigned int i=0; i<M; ++i)
+    {
+      y_host[i] = 1.;
+      y_ref[i] = add ? y_host[i] : 0.;
+      for (unsigned int j=0; j<N; ++j)
+        y_ref[i] += shape_host[i*N+j] * x_host[j];
+    }
+
+  // Copy data to the GPU.
+  CUDA::Vector<double> x_dev(N), y_dev(M);
+  x_dev.import(x_host, VectorOperation::insert);
+  y_dev.import(y_host, VectorOperation::insert);
+
+  unsigned int size_shape_values = M*N*sizeof(double);
+
+  cudaError_t cuda_error = cudaMemcpyToSymbol(CUDAWrappers::internal::global_shape_values,
+                                              shape_host.begin(),
+                                              size_shape_values,
+                                              0,
+                                              cudaMemcpyHostToDevice);
+  AssertCuda(cuda_error);
+
+  cuda_error = cudaMemcpyToSymbol(CUDAWrappers::internal::global_shape_gradients,
+                                  shape_host.begin(),
+                                  size_shape_values,
+                                  0,
+                                  cudaMemcpyHostToDevice);
+  AssertCuda(cuda_error);
+
+  // Launch the kernel
+  evaluate_tensor_product<M,N,type,add,false><<<1,M>>>(y_dev.get_values(),
+                                                       x_dev.get_values());
+
+  // Check the results on the host
+  y_host.import(y_dev, VectorOperation::insert);
+  deallog << "Errors no transpose: ";
+  for (unsigned int i=0; i<M; ++i)
+    deallog << y_host[i] - y_ref[i] << " ";
+  deallog << std::endl;
+
+  for (unsigned int i=0; i<M; ++i)
+    y_host[i] = static_cast<double>(Testing::rand())/RAND_MAX;
+
+  // Copy y_host to the device
+  y_dev.import(y_host, VectorOperation::insert);
+
+  // Compute reference
+  for (unsigned int i=0; i<N; ++i)
+    {
+      x_host[i] = 2.;
+      x_ref[i] = add ? x_host[i] : 0.;
+      for (unsigned int j=0; j<M; ++j)
+        x_ref[i] += shape_host[j*N+i] * y_host[j];
+    }
+
+  // Copy x_host to the device
+  x_dev.import(x_host, VectorOperation::insert);
+
+  // Launch the kernel
+  evaluate_tensor_product<M,N,type,add,true><<<1,M>>>(x_dev.get_values(),
+                                                      y_dev.get_values());
+
+  // Check the results on the host
+  x_host.import(x_dev, VectorOperation::insert);
+  deallog << "Errors transpose:    ";
+  for (unsigned int i=0; i<N; ++i)
+    deallog << x_host[i] - x_ref[i] << " ";
+  deallog << std::endl;
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-14);
+
+  deallog.push("values");
+  test<4,4,0,false>();
+  test<3,3,0,false>();
+  // Not supported right now
+  // test<4,3,0,false>();
+  // test<3,4,0,false>();
+  // test<3,5,0,false>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,false>();
+  test<3,3,1,false>();
+  // Not supported right now
+  // test<4,3,1,false>();
+  // test<3,4,1,false>();
+  // test<3,5,1,false>();
+  deallog.pop();
+
+  deallog.push("add");
+
+  deallog.push("values");
+  test<4,4,0,true>();
+  test<3,3,0,true>();
+  // Not supported right now
+  // test<4,3,0,true>();
+  // test<3,4,0,true>();
+  // test<3,5,0,true>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,true>();
+  test<3,3,1,true>();
+  // Not supported right now
+  // test<4,3,1,true>();
+  // test<3,4,1,true>();
+  // test<3,5,1,true>();
+  deallog.pop();
+
+  deallog.pop();
+
+  return 0;
+}
diff --git a/tests/cuda/cuda_evaluate_1d_shape.output b/tests/cuda/cuda_evaluate_1d_shape.output
new file mode 100644 (file)
index 0000000..5aa8395
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 
diff --git a/tests/cuda/cuda_evaluate_2d_shape.cu b/tests/cuda/cuda_evaluate_2d_shape.cu
new file mode 100644 (file)
index 0000000..06863b7
--- /dev/null
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the 2d evaluation functions used in
+// CUDAWrappers::FEEvaluation. These functions are marked 'internal' but it is
+// much easier to check their correctness directly rather than from the results
+// in dependent functions
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/matrix_free/cuda_fe_evaluation.cuh>
+
+namespace CUDA = LinearAlgebra::CUDAWrappers;
+
+template <int M, int N, int type, bool add, bool dof_to_quad>
+__global__ void evaluate_tensor_product(double *dst,
+                                        double *src)
+{
+  CUDAWrappers::internal::EvaluatorTensorProduct<
+  CUDAWrappers::internal::evaluate_general,2,M-1,N,double> evaluator;
+
+  if (type == 0)
+    {
+      evaluator.template values<0,dof_to_quad,false,false> (src, src);
+      __syncthreads();
+      evaluator.template values<1,dof_to_quad,add,false> (src, dst);
+    }
+  if (type == 1)
+    {
+      evaluator.template gradients<0,dof_to_quad,false,false> (src, src);
+      __syncthreads();
+      evaluator.template gradients<1,dof_to_quad,add,false> (src, dst);
+    }
+}
+
+template <int M, int N, int type, bool add>
+void test()
+{
+  deallog << "Test " << M << " x " << N << std::endl;
+  LinearAlgebra::ReadWriteVector<double> shape_host(M*N);
+  for (unsigned int i=0; i<(M+1)/2; ++i)
+    for (unsigned int j=0; j<N; ++j)
+      {
+        shape_host[i*N+j] =  -1. + 2. *
+                             static_cast<double>(Testing::rand())/RAND_MAX;
+        if (type == 1)
+          shape_host[(M-1-i)*N+N-1-j] = -shape_host[i*N+j];
+        else
+          shape_host[(M-1-i)*N+N-1-j] = shape_host[i*N+j];
+      }
+  if (type == 0 && M%2 == 1 && N%2 == 1)
+    {
+      for (unsigned int i=0; i<M; ++i)
+        shape_host[i*N+N/2] = 0.;
+      shape_host[M/2*N+N/2] = 1.;
+    }
+  if (type == 1 && M%2 == 1 && N%2 == 1)
+    shape_host[M/2*N+N/2] = 0.;
+
+  constexpr int M_2d = M*M;
+  constexpr int N_2d = N*N;
+  LinearAlgebra::ReadWriteVector<double> x_host(N_2d), x_ref(N_2d),
+                y_host(M_2d), y_ref(M_2d);
+  for (unsigned int i=0; i<N_2d; ++i)
+    x_host[i] =  static_cast<double>(Testing::rand())/RAND_MAX;
+
+  FullMatrix<double> shape_2d(M_2d, N_2d);
+  for (unsigned int i=0; i<M; ++i)
+    {
+      for (unsigned int j=0; j<N; ++j)
+        {
+          const double shape_val = shape_host[i*N+j];
+          for (unsigned int m=0; m<M; ++m)
+            for (unsigned int n=0; n<N; ++n)
+              shape_2d(i*M+m,j*N+n) = shape_val*shape_host[m*N+n];
+        }
+    }
+
+  // Compute reference
+  for (unsigned int i=0; i<M_2d; ++i)
+    {
+      y_host[i] = 1.;
+      y_ref[i] = add ? y_host[i] : 0.;
+      for (unsigned int j=0; j<N_2d; ++j)
+        y_ref[i] += shape_2d(i,j) * x_host[j];
+    }
+
+  // Copy data to the GPU.
+  CUDA::Vector<double> x_dev(N_2d), y_dev(M_2d);
+  x_dev.import(x_host, VectorOperation::insert);
+  y_dev.import(y_host, VectorOperation::insert);
+
+  unsigned int size_shape_values = M*N*sizeof(double);
+
+  cudaError_t cuda_error = cudaMemcpyToSymbol(CUDAWrappers::internal::global_shape_values,
+                                              shape_host.begin(),
+                                              size_shape_values,
+                                              0,
+                                              cudaMemcpyHostToDevice);
+  AssertCuda(cuda_error);
+
+  cuda_error = cudaMemcpyToSymbol(CUDAWrappers::internal::global_shape_gradients,
+                                  shape_host.begin(),
+                                  size_shape_values,
+                                  0,
+                                  cudaMemcpyHostToDevice);
+  AssertCuda(cuda_error);
+
+  // Launch the kernel
+  dim3 block_dim(M,N);
+  evaluate_tensor_product<M,N,type,add,false><<<1,block_dim>>>(y_dev.get_values(),
+      x_dev.get_values());
+
+  // Check the results on the host
+  y_host.import(y_dev, VectorOperation::insert);
+  deallog << "Errors no transpose: ";
+
+  for (unsigned int i=0; i<M_2d; ++i)
+    deallog << y_host[i] - y_ref[i] << " ";
+  deallog << std::endl;
+
+  for (unsigned int i=0; i<M_2d; ++i)
+    y_host[i] = static_cast<double>(Testing::rand())/RAND_MAX;
+
+  // Copy y_host to the device
+  y_dev.import(y_host, VectorOperation::insert);
+
+  // Compute reference
+  for (unsigned int i=0; i<N_2d; ++i)
+    {
+      x_host[i] = 2.;
+      x_ref[i] = add ? x_host[i] : 0.;
+      for (unsigned int j=0; j<M_2d; ++j)
+        x_ref[i] += shape_2d(j,i) * y_host[j];
+    }
+
+  // Copy x_host to the device
+  x_dev.import(x_host, VectorOperation::insert);
+
+  // Launch the kernel
+  evaluate_tensor_product<M,N,type,add,true><<<1,block_dim>>>(x_dev.get_values(),
+                                                              y_dev.get_values());
+
+  // Check the results on the host
+  x_host.import(x_dev, VectorOperation::insert);
+  deallog << "Errors transpose:    ";
+  for (unsigned int i=0; i<N_2d; ++i)
+    deallog << x_host[i] - x_ref[i] << " ";
+  deallog << std::endl;
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-14);
+
+  deallog.push("values");
+  test<4,4,0,false>();
+  test<3,3,0,false>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,false>();
+  test<3,3,1,false>();
+  deallog.pop();
+
+  deallog.push("add");
+
+  deallog.push("values");
+  test<4,4,0,true>();
+  test<3,3,0,true>();
+  deallog.pop();
+
+  deallog.push("gradients");
+  test<4,4,1,true>();
+  test<3,3,1,true>();
+  deallog.pop();
+
+  deallog.pop();
+
+  return 0;
+}
diff --git a/tests/cuda/cuda_evaluate_2d_shape.output b/tests/cuda/cuda_evaluate_2d_shape.output
new file mode 100644 (file)
index 0000000..312a78a
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0
+DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0
diff --git a/tests/cuda/matrix_free_matrix_vector_01.cu b/tests/cuda/matrix_free_matrix_vector_01.cu
new file mode 100644 (file)
index 0000000..1a32446
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. The mesh uses a hypercube mesh with no hanging nodes and no other
+// constraints
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.cuh"
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(5-dim);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  do_test<dim, fe_degree, double, fe_degree+1> (dof, constraints);
+}
diff --git a/tests/cuda/matrix_free_matrix_vector_01.output b/tests/cuda/matrix_free_matrix_vector_01.output
new file mode 100644 (file)
index 0000000..43ef13c
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(3)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/cuda/matrix_free_matrix_vector_04.cu b/tests/cuda/matrix_free_matrix_vector_04.cu
new file mode 100644 (file)
index 0000000..3db14f8
--- /dev/null
@@ -0,0 +1,49 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. The mesh uses a hypershell mesh without hanging nodes (only cell
+// type: 2)
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.cuh"
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+
+  const SphericalManifold<dim> manifold;
+  Triangulation<dim> tria;
+  GridGenerator::hyper_shell (tria, Point<dim>(),
+                              0.5, 1., 96, true);
+  tria.set_all_manifold_ids(0);
+  tria.set_manifold (0, manifold);
+  if (dim == 2)
+    tria.refine_global (2);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  do_test<dim, fe_degree, double, fe_degree+1> (dof, constraints);
+}
diff --git a/tests/cuda/matrix_free_matrix_vector_04.output b/tests/cuda/matrix_free_matrix_vector_04.output
new file mode 100644 (file)
index 0000000..43ef13c
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(3)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/cuda/matrix_vector_common.cuh b/tests/cuda/matrix_vector_common.cuh
new file mode 100644 (file)
index 0000000..7a93226
--- /dev/null
@@ -0,0 +1,173 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// this is a template for matrix-vector products with the Helmholtz equation
+// (zero and first derivatives) on different kinds of meshes (Cartesian,
+// general, with and without hanging nodes). It also tests the multithreading
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.cuh"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+
+// forward declare this function. will be implemented in .cc files
+template <int dim, int fe_degree>
+void test();
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+void do_test(const DoFHandler<dim> &dof,
+             const ConstraintMatrix &constraints)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MappingQGeneric<dim> mapping(fe_degree);
+  CUDAWrappers::MatrixFree<dim,Number> mf_data;
+  typename CUDAWrappers::MatrixFree<dim,Number>::AdditionalData additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients |
+                                         update_JxW_values | update_quadrature_points;
+  const QGauss<1> quad (n_q_points_1d);
+  mf_data.reinit (mapping, dof, constraints, quad, additional_data);
+
+  const unsigned int n_dofs = dof.n_dofs();
+  MatrixFreeTest<dim,fe_degree,Number,n_q_points_1d> mf(mf_data);
+  Vector<Number> in_host(n_dofs), out_host(n_dofs);
+  LinearAlgebra::ReadWriteVector<Number> in(n_dofs), out(n_dofs);
+  LinearAlgebra::CUDAWrappers::Vector<Number> in_device(n_dofs);
+  LinearAlgebra::CUDAWrappers::Vector<Number> out_device(n_dofs);
+
+  for (unsigned int i=0; i<n_dofs; ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+    }
+
+  in_device.import(in, VectorOperation::insert);
+  mf.vmult(out_device, in_device);
+  cudaDeviceSynchronize();
+  out.import(out_device, VectorOperation::insert);
+
+
+  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern csp(n_dofs, n_dofs);
+    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+    sparsity.copy_from(csp);
+  }
+  SparseMatrix<double> sparse_matrix (sparsity);
+  {
+    QGauss<dim> quadrature_formula(n_q_points_1d);
+
+    FEValues<dim> fe_values (mapping, dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        fe_values.reinit (cell);
+
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+                                      fe_values.shape_grad(j,q_point)
+                                      +
+                                      10. *
+                                      fe_values.shape_value(i,q_point) *
+                                      fe_values.shape_value(j,q_point)) *
+                                     fe_values.JxW(q_point));
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global (cell_matrix,
+                                                local_dof_indices,
+                                                sparse_matrix);
+      }
+  }
+  for (unsigned i = 0; i<n_dofs; ++i)
+    in_host[i] = in[i];
+  sparse_matrix.vmult(out_host, in_host);
+
+  Number out_dist_cpu_norm = 0.;
+  Number out_norm = 0.;
+  for (unsigned i = 0; i<n_dofs; ++i)
+    {
+      out_norm += std::pow(out[i] - out_host[i], 2);
+      out_dist_cpu_norm += std::pow(out_host[i], 2);
+    }
+  const double diff_norm = out_norm / out_dist_cpu_norm;
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+int main()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog.threshold_double(5.e-11);
+    deallog.push("2d");
+    test<2,1>();
+    test<2,2>();
+    test<2,3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    test<3,3>();
+    deallog.pop();
+  }
+
+  return 0;
+}
diff --git a/tests/cuda/matrix_vector_mf.cuh b/tests/cuda/matrix_vector_mf.cuh
new file mode 100644 (file)
index 0000000..08c5ef4
--- /dev/null
@@ -0,0 +1,118 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// This is a template for matrix-vector products with the Helmholtz equation
+// (zero and first derivatives) on different kinds of meshes (Cartesian,
+// general, with and without hanging nodes).
+
+#include "../tests.h"
+
+#include <deal.II/matrix_free/cuda_matrix_free.h>
+#include <deal.II/matrix_free/cuda_fe_evaluation.cuh>
+
+#include <deal.II/lac/cuda_vector.h>
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+class HelmholtzOperatorQuad
+{
+public :
+  __device__ void operator()(
+    CUDAWrappers::FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> *fe_eval,
+    const unsigned int q_point) const;
+};
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+__device__ void
+HelmholtzOperatorQuad<dim,fe_degree,Number,n_q_points_1d>::operator()(
+  CUDAWrappers::FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> *fe_eval,
+  const unsigned int q) const
+{
+  fe_eval->submit_value(Number(10)*fe_eval->get_value(q),q);
+  fe_eval->submit_gradient(fe_eval->get_gradient(q),q);
+}
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+class HelmholtzOperator
+{
+public:
+  __device__ void operator()(
+    const unsigned int cell,
+    const typename CUDAWrappers::MatrixFree<dim,Number>::Data *gpu_data,
+    CUDAWrappers::SharedData<dim,Number>                *shared_data,
+    const Number                                        *src,
+    Number                                              *dst) const;
+
+  static const unsigned int n_dofs_1d = fe_degree+1;
+  static const unsigned int n_local_dofs =
+    dealii::Utilities::fixed_int_power<fe_degree+1,dim>::value;
+  static const unsigned int n_q_points =
+    dealii::Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+};
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+__device__ void
+HelmholtzOperator<dim,fe_degree,Number,n_q_points_1d>::operator()(
+  const unsigned int                                        cell,
+  const typename CUDAWrappers::MatrixFree<dim,Number>::Data *gpu_data,
+  CUDAWrappers::SharedData<dim,Number>                *shared_data,
+  const Number                                        *src,
+  Number                                              *dst) const
+{
+  CUDAWrappers::FEEvaluation<dim,fe_degree,n_q_points_1d,1,Number> fe_eval(
+    cell, gpu_data, shared_data);
+  fe_eval.read_dof_values(src);
+  fe_eval.evaluate(true, true);
+  fe_eval.apply_quad_point_operations(
+    HelmholtzOperatorQuad<dim,fe_degree,Number,n_q_points_1d>());
+  fe_eval.integrate(true, true);
+  fe_eval.distribute_local_to_global (dst);
+}
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d =
+          fe_degree+1>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const CUDAWrappers::MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void vmult(LinearAlgebra::CUDAWrappers::Vector<Number>       &dst,
+             const LinearAlgebra::CUDAWrappers::Vector<Number> &src);
+
+private:
+  const CUDAWrappers::MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+void MatrixFreeTest<dim,fe_degree,Number,n_q_points_1d>::vmult(
+  LinearAlgebra::CUDAWrappers::Vector<Number>       &dst,
+  const LinearAlgebra::CUDAWrappers::Vector<Number> &src)
+{
+  dst = static_cast<Number>(0.);
+  HelmholtzOperator<dim,fe_degree,Number,n_q_points_1d> helmholtz_operator;
+  data.cell_loop(helmholtz_operator, src, dst);
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.