From 96993d6c85a5171904b9e808785de9411ec45a0b Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Sat, 13 May 2017 23:05:05 -0400 Subject: [PATCH] Add test for CUDA support of matrix-free. --- tests/cuda/cuda_evaluate_1d_shape.cu | 193 +++++++++++++++++ tests/cuda/cuda_evaluate_1d_shape.output | 25 +++ tests/cuda/cuda_evaluate_2d_shape.cu | 201 ++++++++++++++++++ tests/cuda/cuda_evaluate_2d_shape.output | 25 +++ tests/cuda/matrix_free_matrix_vector_01.cu | 44 ++++ .../cuda/matrix_free_matrix_vector_01.output | 19 ++ tests/cuda/matrix_free_matrix_vector_04.cu | 49 +++++ .../cuda/matrix_free_matrix_vector_04.output | 19 ++ tests/cuda/matrix_vector_common.cuh | 173 +++++++++++++++ tests/cuda/matrix_vector_mf.cuh | 118 ++++++++++ 10 files changed, 866 insertions(+) create mode 100644 tests/cuda/cuda_evaluate_1d_shape.cu create mode 100644 tests/cuda/cuda_evaluate_1d_shape.output create mode 100644 tests/cuda/cuda_evaluate_2d_shape.cu create mode 100644 tests/cuda/cuda_evaluate_2d_shape.output create mode 100644 tests/cuda/matrix_free_matrix_vector_01.cu create mode 100644 tests/cuda/matrix_free_matrix_vector_01.output create mode 100644 tests/cuda/matrix_free_matrix_vector_04.cu create mode 100644 tests/cuda/matrix_free_matrix_vector_04.output create mode 100644 tests/cuda/matrix_vector_common.cuh create mode 100644 tests/cuda/matrix_vector_mf.cuh diff --git a/tests/cuda/cuda_evaluate_1d_shape.cu b/tests/cuda/cuda_evaluate_1d_shape.cu new file mode 100644 index 0000000000..040a2fedc3 --- /dev/null +++ b/tests/cuda/cuda_evaluate_1d_shape.cu @@ -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 +#include + +#include +#include +#include + +namespace CUDA = LinearAlgebra::CUDAWrappers; + +template +__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 +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + LinearAlgebra::ReadWriteVector shape_host(M*N); + for (unsigned int i=0; i<(M+1)/2; ++i) + for (unsigned int j=0; j(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 x_host(N), x_ref(N), y_host(M), y_ref(M); + for (unsigned int i=0; i(Testing::rand())/RAND_MAX; + + // Compute reference + for (unsigned int i=0; i 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<<<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(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<<<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(); + 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 index 0000000000..5aa83953a8 --- /dev/null +++ b/tests/cuda/cuda_evaluate_1d_shape.output @@ -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 index 0000000000..06863b731e --- /dev/null +++ b/tests/cuda/cuda_evaluate_2d_shape.cu @@ -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 +#include + +#include +#include +#include + +namespace CUDA = LinearAlgebra::CUDAWrappers; + +template +__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 +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + LinearAlgebra::ReadWriteVector shape_host(M*N); + for (unsigned int i=0; i<(M+1)/2; ++i) + for (unsigned int j=0; j(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 x_host(N_2d), x_ref(N_2d), + y_host(M_2d), y_ref(M_2d); + for (unsigned int i=0; i(Testing::rand())/RAND_MAX; + + FullMatrix shape_2d(M_2d, N_2d); + for (unsigned int i=0; i 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<<<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(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<<<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(); + 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 index 0000000000..312a78a60e --- /dev/null +++ b/tests/cuda/cuda_evaluate_2d_shape.output @@ -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 index 0000000000..1a32446ba3 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_01.cu @@ -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 +void test () +{ + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.refine_global(5-dim); + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + constraints.close(); + + do_test (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 index 0000000000..43ef13cd75 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_01.output @@ -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 index 0000000000..3db14f8876 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_04.cu @@ -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 +void test () +{ + + const SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_shell (tria, Point(), + 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 fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + constraints.close(); + + do_test (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 index 0000000000..43ef13cd75 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_04.output @@ -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 index 0000000000..7a93226338 --- /dev/null +++ b/tests/cuda/matrix_vector_common.cuh @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +// forward declare this function. will be implemented in .cc files +template +void test(); + + + +template +void do_test(const DoFHandler &dof, + const ConstraintMatrix &constraints) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MappingQGeneric mapping(fe_degree); + CUDAWrappers::MatrixFree mf_data; + typename CUDAWrappers::MatrixFree::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 mf(mf_data); + Vector in_host(n_dofs), out_host(n_dofs); + LinearAlgebra::ReadWriteVector in(n_dofs), out(n_dofs); + LinearAlgebra::CUDAWrappers::Vector in_device(n_dofs); + LinearAlgebra::CUDAWrappers::Vector out_device(n_dofs); + + for (unsigned int i=0; i sparse_matrix (sparsity); + { + QGauss quadrature_formula(n_q_points_1d); + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + for (unsigned i = 0; i(); + 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 index 0000000000..08c5ef40aa --- /dev/null +++ b/tests/cuda/matrix_vector_mf.cuh @@ -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 +#include + +#include + +template +class HelmholtzOperatorQuad +{ +public : + __device__ void operator()( + CUDAWrappers::FEEvaluation *fe_eval, + const unsigned int q_point) const; +}; + + + +template +__device__ void +HelmholtzOperatorQuad::operator()( + CUDAWrappers::FEEvaluation *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 +class HelmholtzOperator +{ +public: + __device__ void operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + CUDAWrappers::SharedData *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::value; + static const unsigned int n_q_points = + dealii::Utilities::fixed_int_power::value; +}; + + + +template +__device__ void +HelmholtzOperator::operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + CUDAWrappers::SharedData *shared_data, + const Number *src, + Number *dst) const +{ + CUDAWrappers::FEEvaluation 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()); + fe_eval.integrate(true, true); + fe_eval.distribute_local_to_global (dst); +} + + + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in): + data (data_in) + {}; + + void vmult(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src); + +private: + const CUDAWrappers::MatrixFree &data; +}; + + + +template +void MatrixFreeTest::vmult( + LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) +{ + dst = static_cast(0.); + HelmholtzOperator helmholtz_operator; + data.cell_loop(helmholtz_operator, src, dst); +} -- 2.39.5