From: kronbichler Date: Tue, 22 Oct 2013 09:50:02 +0000 (+0000) Subject: Do symmetric tensor product evaluation with even-odd formula. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fa4f9f2391754d877130be15446be9a3a787e7ba;p=dealii-svn.git Do symmetric tensor product evaluation with even-odd formula. git-svn-id: https://svn.dealii.org/trunk@31380 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index edcb3095cf..aa3ff08f55 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -1304,6 +1304,10 @@ protected: void apply_hessians (const VectorizedArray in [], VectorizedArray out []); + VectorizedArray shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; + VectorizedArray shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; + VectorizedArray shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; + /** * Friend declarations. */ @@ -4390,6 +4394,204 @@ namespace internal + // This method implements a different approach to the symmetric case for + // values, gradients, and Hessians also treated with the above functions: It + // is possible to reduce the cost per dimension from N^2 to N^2/2, where N + // is the number of 1D dofs (there are only N^2/2 different entries in the + // shape matrix, so this is plausible). The approach is based on the idea of + // applying the operator on the even and odd part of the input vectors + // separately, given that the shape functions evaluated on quadrature points + // are symmetric. This method is presented e.g. in the book "Implementing + // Spectral Methods for Partial Differential Equations" by David A. Kopriva, + // Springer, 2009, section 3.5.3 (Even-Odd-Decomposition). Even though the + // experiments in the book say that the method is not efficient for N<20, it + // is more efficient in the context where the loop bounds are compile-time + // constants (templates). + template + inline + void + apply_tensor_product_evenodd (const Number shapes [][(n_q_points_1d+1)/2], + const Number in [], + Number out []) + { + AssertIndexRange (type, 3); + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + const int n_cols = nn / 2; + const int mid = mm / 2; + + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; + + // this code may look very inefficient at first sight due to the many + // different cases with if's at the innermost loop part, but all of the + // conditionals can be evaluated at compile time because they are + // templates, so the compiler should optimize everything away + for (int i2=0; i2 0) + { + if (dof_to_quad == true) + { + r0 = shapes[0][col] * xp[0]; + r1 = shapes[fe_degree][col] * xm[0]; + } + else + { + r0 = shapes[col][0] * xp[0]; + r1 = shapes[fe_degree-col][0] * xm[0]; + } + for (int ind=1; ind 0)) + r0 += shapes[col][mid] * in[stride*mid]; + + if (add == false) + { + out[stride*col] = r0 + r1; + if (type == 1 && dof_to_quad == false) + out[stride*(nn-1-col)] = r1 - r0; + else + out[stride*(nn-1-col)] = r0 - r1; + } + else + { + out[stride*col] += r0 + r1; + if (type == 1 && dof_to_quad == false) + out[stride*(nn-1-col)] += r1 - r0; + else + out[stride*(nn-1-col)] += r0 - r1; + } + } + if ( type == 0 && dof_to_quad == true && nn%2==1 && mm%2==1 ) + { + if (add==false) + out[stride*n_cols] = in[stride*mid]; + else + out[stride*n_cols] += in[stride*mid]; + } + else if (dof_to_quad == true && nn%2==1) + { + Number r0; + if (mid > 0) + { + r0 = shapes[0][n_cols] * xp[0]; + for (int ind=1; ind 0) + { + if (type == 1) + { + r0 = shapes[n_cols][0] * xm[0]; + for (int ind=1; ind j-1][0]) < zero_tol, ExcMessage(error_message)); #endif + + // Compute symmetric and skew-symmetric part of shape values for even-odd + // decomposition + for (unsigned int i=0; i<(fe_degree+1)/2; ++i) + for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) + { + shape_val_evenodd[i][q] = + 0.5 * (this->data.shape_values[i*n_q_points_1d+q] + + this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); + shape_val_evenodd[fe_degree-i][q] = + 0.5 * (this->data.shape_values[i*n_q_points_1d+q] - + this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); + + shape_gra_evenodd[i][q] = + 0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] + + this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); + shape_gra_evenodd[fe_degree-i][q] = + 0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] - + this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); + + shape_hes_evenodd[i][q] = + 0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] + + this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); + shape_hes_evenodd[fe_degree-i][q] = + 0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] - + this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); + } + if (fe_degree % 2 == 0) + for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) + { + shape_val_evenodd[fe_degree/2][q] = + this->data.shape_values[(fe_degree/2)*n_q_points_1d+q]; + shape_gra_evenodd[fe_degree/2][q] = + this->data.shape_gradients[(fe_degree/2)*n_q_points_1d+q]; + shape_hes_evenodd[fe_degree/2][q] = + this->data.shape_hessians[(fe_degree/2)*n_q_points_1d+q]; + } } @@ -5129,9 +5369,9 @@ FEEvaluation ::apply_values (const VectorizedArray in [], VectorizedArray out []) { - internal::apply_tensor_product_values, direction, dof_to_quad, add> - (this->data.shape_values.begin(), in, out); + internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 0> + (shape_val_evenodd, in, out); } @@ -5145,9 +5385,9 @@ FEEvaluation ::apply_gradients (const VectorizedArray in [], VectorizedArray out []) { - internal::apply_tensor_product_gradients, direction, dof_to_quad, add> - (this->data.shape_gradients.begin(), in, out); + internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 1> + (shape_gra_evenodd, in, out); } @@ -5164,9 +5404,9 @@ FEEvaluation ::apply_hessians (const VectorizedArray in [], VectorizedArray out []) { - internal::apply_tensor_product_hessians, direction, dof_to_quad, add> - (this->data.shape_hessians.begin(), in, out); + internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 2> + (shape_hes_evenodd, in, out); } diff --git a/tests/matrix_free/evaluate_1d_shape.cc b/tests/matrix_free/evaluate_1d_shape.cc new file mode 100644 index 0000000000..b19a94c4be --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape.cc @@ -0,0 +1,174 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 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 +// 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 + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + double shape[M][N]; + for (unsigned int i=0; i<(M+1)/2; ++i) + for (unsigned int j=0; j + (&shape[0][0],x,y); + if (type == 1) + internal::apply_tensor_product_gradients<1,M-1,N,double,0,false,add> + (&shape[0][0],x,y); + if (type == 2) + internal::apply_tensor_product_hessians<1,M-1,N,double,0,false,add> + (&shape[0][0],x,y); + + deallog << "Errors no transpose: "; + for (unsigned int i=0; i + (&shape[0][0],y,x); + if (type == 1) + internal::apply_tensor_product_gradients<1,M-1,N,double,0,true,add> + (&shape[0][0],y,x); + if (type == 2) + internal::apply_tensor_product_hessians<1,M-1,N,double,0,true,add> + (&shape[0][0],y,x); + + deallog << "Errors transpose: "; + for (unsigned int i=0; i(); + test<3,3,0,false>(); + 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>(); + test<4,3,1,false>(); + test<3,4,1,false>(); + test<3,5,1,false>(); + deallog.pop(); + + deallog.push("hessians"); + test<4,4,2,false>(); + test<3,3,2,false>(); + test<4,3,2,false>(); + test<3,4,2,false>(); + test<3,5,2,false>(); + deallog.pop(); + + deallog.push("add"); + + deallog.push("values"); + test<4,4,0,true>(); + test<3,3,0,true>(); + 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>(); + test<4,3,1,true>(); + test<3,4,1,true>(); + test<3,5,1,true>(); + deallog.pop(); + + deallog.push("hessians"); + test<4,4,2,true>(); + test<3,3,2,true>(); + test<4,3,2,true>(); + test<3,4,2,true>(); + test<3,5,2,true>(); + deallog.pop(); + + deallog.pop(); + + return 0; +} + diff --git a/tests/matrix_free/evaluate_1d_shape.output b/tests/matrix_free/evaluate_1d_shape.output new file mode 100644 index 0000000000..2bd9530d6c --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape.output @@ -0,0 +1,91 @@ + +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:values::Test 4 x 3 +DEAL:values::Errors no transpose: 0 0 0 0 +DEAL:values::Errors transpose: 0 0 0 +DEAL:values::Test 3 x 4 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 +DEAL:values::Test 3 x 5 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 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:gradients::Test 4 x 3 +DEAL:gradients::Errors no transpose: 0 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 +DEAL:gradients::Test 3 x 4 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 +DEAL:gradients::Test 3 x 5 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 +DEAL:hessians::Test 4 x 4 +DEAL:hessians::Errors no transpose: 0 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 0 +DEAL:hessians::Test 3 x 3 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 +DEAL:hessians::Test 4 x 3 +DEAL:hessians::Errors no transpose: 0 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 +DEAL:hessians::Test 3 x 4 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 0 +DEAL:hessians::Test 3 x 5 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 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:values::Test 4 x 3 +DEAL:add:values::Errors no transpose: 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 +DEAL:add:values::Test 3 x 4 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 +DEAL:add:values::Test 3 x 5 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 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 +DEAL:add:gradients::Test 4 x 3 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 +DEAL:add:gradients::Test 3 x 4 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 +DEAL:add:gradients::Test 3 x 5 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 +DEAL:add:hessians::Test 4 x 4 +DEAL:add:hessians::Errors no transpose: 0 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 +DEAL:add:hessians::Test 3 x 3 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 +DEAL:add:hessians::Test 4 x 3 +DEAL:add:hessians::Errors no transpose: 0 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 +DEAL:add:hessians::Test 3 x 4 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 +DEAL:add:hessians::Test 3 x 5 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 0 diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.cc b/tests/matrix_free/evaluate_1d_shape_evenodd.cc new file mode 100644 index 0000000000..2eab66476d --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_evenodd.cc @@ -0,0 +1,172 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 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 +// FEEvaluation. These functions are marked 'internal' but it is much easier +// to check their correctness directly rather than from the results in +// dependent functions. this function tests the even-odd path of the +// evaluation functions + +#include "../tests.h" +#include +#include + +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + double shape[M][N]; + for (unsigned int i=0; i<(M+1)/2; ++i) + for (unsigned int j=0; j(shape_sym,x,y); + + deallog << "Errors no transpose: "; + for (unsigned int i=0; i(shape_sym,y,x); + + deallog << "Errors transpose: "; + for (unsigned int i=0; i(); + test<3,3,0,false>(); + 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>(); + test<4,3,1,false>(); + test<3,4,1,false>(); + test<3,5,1,false>(); + deallog.pop(); + + deallog.push("hessians"); + test<4,4,2,false>(); + test<3,3,2,false>(); + test<4,3,2,false>(); + test<3,4,2,false>(); + test<3,5,2,false>(); + deallog.pop(); + + deallog.push("add"); + + deallog.push("values"); + test<4,4,0,true>(); + test<3,3,0,true>(); + 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>(); + test<4,3,1,true>(); + test<3,4,1,true>(); + test<3,5,1,true>(); + deallog.pop(); + + deallog.push("hessians"); + test<4,4,2,true>(); + test<3,3,2,true>(); + test<4,3,2,true>(); + test<3,4,2,true>(); + test<3,5,2,true>(); + deallog.pop(); + + deallog.pop(); + + return 0; +} + diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.output b/tests/matrix_free/evaluate_1d_shape_evenodd.output new file mode 100644 index 0000000000..2bd9530d6c --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_evenodd.output @@ -0,0 +1,91 @@ + +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:values::Test 4 x 3 +DEAL:values::Errors no transpose: 0 0 0 0 +DEAL:values::Errors transpose: 0 0 0 +DEAL:values::Test 3 x 4 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 +DEAL:values::Test 3 x 5 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 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:gradients::Test 4 x 3 +DEAL:gradients::Errors no transpose: 0 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 +DEAL:gradients::Test 3 x 4 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 +DEAL:gradients::Test 3 x 5 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 +DEAL:hessians::Test 4 x 4 +DEAL:hessians::Errors no transpose: 0 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 0 +DEAL:hessians::Test 3 x 3 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 +DEAL:hessians::Test 4 x 3 +DEAL:hessians::Errors no transpose: 0 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 +DEAL:hessians::Test 3 x 4 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 0 0 +DEAL:hessians::Test 3 x 5 +DEAL:hessians::Errors no transpose: 0 0 0 +DEAL:hessians::Errors transpose: 0 0 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:values::Test 4 x 3 +DEAL:add:values::Errors no transpose: 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 +DEAL:add:values::Test 3 x 4 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 +DEAL:add:values::Test 3 x 5 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 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 +DEAL:add:gradients::Test 4 x 3 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 +DEAL:add:gradients::Test 3 x 4 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 +DEAL:add:gradients::Test 3 x 5 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 +DEAL:add:hessians::Test 4 x 4 +DEAL:add:hessians::Errors no transpose: 0 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 +DEAL:add:hessians::Test 3 x 3 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 +DEAL:add:hessians::Test 4 x 3 +DEAL:add:hessians::Errors no transpose: 0 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 +DEAL:add:hessians::Test 3 x 4 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 +DEAL:add:hessians::Test 3 x 5 +DEAL:add:hessians::Errors no transpose: 0 0 0 +DEAL:add:hessians::Errors transpose: 0 0 0 0 0