--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}