void apply_hessians (const VectorizedArray<Number> in [],
VectorizedArray<Number> out []);
+ VectorizedArray<Number> shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+ VectorizedArray<Number> shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+ VectorizedArray<Number> shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2];
+
/**
* Friend declarations.
*/
+ // 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 <int dim, int fe_degree, int n_q_points_1d, typename Number,
+ int direction, bool dof_to_quad, bool add, int type>
+ 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<nn,direction>::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<n_blocks2; ++i2)
+ {
+ for (int i1=0; i1<n_blocks1; ++i1)
+ {
+ Number xp[mid], xm[mid];
+ for (int i=0; i<mid; ++i)
+ {
+ if (dof_to_quad == true && type == 1)
+ {
+ xp[i] = in[stride*i] - in[stride*(mm-1-i)];
+ xm[i] = in[stride*i] + in[stride*(mm-1-i)];
+ }
+ else
+ {
+ xp[i] = in[stride*i] + in[stride*(mm-1-i)];
+ xm[i] = in[stride*i] - in[stride*(mm-1-i)];
+ }
+ }
+ for (int col=0; col<n_cols; ++col)
+ {
+ Number r0, r1;
+ if (mid > 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<mid; ++ind)
+ {
+ if (dof_to_quad == true)
+ {
+ r0 += shapes[ind][col] * xp[ind];
+ r1 += shapes[fe_degree-ind][col] * xm[ind];
+ }
+ else
+ {
+ r0 += shapes[col][ind] * xp[ind];
+ r1 += shapes[fe_degree-col][ind] * xm[ind];
+ }
+ }
+ }
+ else
+ r0 = r1 = Number();
+ if (mm % 2 == 1 && dof_to_quad == true)
+ {
+ if (type == 1)
+ r1 += shapes[mid][col] * in[stride*mid];
+ else
+ r0 += shapes[mid][col] * in[stride*mid];
+ }
+ else if (mm % 2 == 1 && (nn % 2 == 0 || type > 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<mid; ++ind)
+ r0 += shapes[ind][n_cols] * xp[ind];
+ }
+ else
+ r0 = Number();
+ if (type != 1 && mm % 2 == 1)
+ r0 += shapes[mid][n_cols] * in[stride*mid];
+
+ if (add == false)
+ out[stride*n_cols] = r0;
+ else
+ out[stride*n_cols] += r0;
+ }
+ else if (dof_to_quad == false && nn%2 == 1)
+ {
+ Number r0;
+ if (mid > 0)
+ {
+ if (type == 1)
+ {
+ r0 = shapes[n_cols][0] * xm[0];
+ for (int ind=1; ind<mid; ++ind)
+ r0 += shapes[n_cols][ind] * xm[ind];
+ }
+ else
+ {
+ r0 = shapes[n_cols][0] * xp[0];
+ for (int ind=1; ind<mid; ++ind)
+ r0 += shapes[n_cols][ind] * xp[ind];
+ }
+ }
+ else
+ r0 = Number();
+
+ if (type == 0 && mm % 2 == 1)
+ r0 += in[stride*mid];
+ else if (type == 2 && mm % 2 == 1)
+ r0 += shapes[n_cols][mid] * in[stride*mid];
+
+ if (add == false)
+ out[stride*n_cols] = r0;
+ else
+ out[stride*n_cols] += r0;
+ }
+
+ // increment: in regular case, just go to the next point in
+ // x-direction. If we are at the end of one chunk in x-dir, need to
+ // jump over to the next layer in z-direction
+ switch (direction)
+ {
+ case 0:
+ in += mm;
+ out += nn;
+ break;
+ case 1:
+ case 2:
+ ++in;
+ ++out;
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ if (direction == 1)
+ {
+ in += nn*(mm-1);
+ out += nn*(nn-1);
+ }
+ }
+ }
+
+
+
// evaluates the given shape data in 1d-3d using the tensor product
// form assuming the symmetries of unit cell shape gradients for
// finite elements in FEEvaluationGL
fe_eval.dof_values_initialized = true;
#endif
}
-}
+
+} // end of namespace internal
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];
+ }
}
::apply_values (const VectorizedArray<Number> in [],
VectorizedArray<Number> out [])
{
- internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_values.begin(), in, out);
+ internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+ VectorizedArray<Number>, direction, dof_to_quad, add, 0>
+ (shape_val_evenodd, in, out);
}
::apply_gradients (const VectorizedArray<Number> in [],
VectorizedArray<Number> out [])
{
- internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_gradients.begin(), in, out);
+ internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+ VectorizedArray<Number>, direction, dof_to_quad, add, 1>
+ (shape_gra_evenodd, in, out);
}
::apply_hessians (const VectorizedArray<Number> in [],
VectorizedArray<Number> out [])
{
- internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_hessians.begin(), in, out);
+ internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
+ VectorizedArray<Number>, direction, dof_to_quad, add, 2>
+ (shape_hes_evenodd, in, out);
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <iostream>
+#include <fstream>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+template <int M, int N, int type, bool add>
+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<N; ++j)
+ {
+ shape[i][j] = -1. + 2. * (double)rand()/RAND_MAX;
+ if (type == 1)
+ shape[M-1-i][N-1-j] = -shape[i][j];
+ else
+ shape[M-1-i][N-1-j] = shape[i][j];
+ }
+ if (type == 0 && M%2 == 1 && N%2 == 1)
+ {
+ for (unsigned int i=0; i<M; ++i)
+ shape[i][N/2] = 0.;
+ shape[M/2][N/2] = 1;
+ }
+ if (type == 1 && M%2 == 1 && N%2 == 1)
+ shape[M/2][N/2] = 0.;
+
+ double x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = (double)rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i][j] * x[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ internal::apply_tensor_product_values<1,M-1,N,double,0,false,add>
+ (&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<M; ++i)
+ deallog << y[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ y[i] = (double)rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j][i] * y[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ internal::apply_tensor_product_values<1,M-1,N,double,0,true,add>
+ (&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<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1e-14);
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ 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;
+}
+
--- /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: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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <iostream>
+#include <fstream>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+template <int M, int N, int type, bool add>
+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<N; ++j)
+ {
+ shape[i][j] = -1. + 2. * (double)rand()/RAND_MAX;
+ if (type == 1)
+ shape[M-1-i][N-1-j] = -shape[i][j];
+ else
+ shape[M-1-i][N-1-j] = shape[i][j];
+ }
+ if (type == 0 && M%2 == 1 && N%2 == 1)
+ {
+ for (unsigned int i=0; i<M; ++i)
+ shape[i][N/2] = 0.;
+ shape[M/2][N/2] = 1;
+ }
+ if (type == 1 && M%2 == 1 && N%2 == 1)
+ shape[M/2][N/2] = 0.;
+
+ // create symmetrized shape array exactly as expected by the evenodd
+ // function
+ double shape_sym[M][(N+1)/2];
+ for (unsigned int i=0; i<M/2; ++i)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ {
+ shape_sym[i][q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
+ shape_sym[M-1-i][q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
+ }
+ if (M % 2 == 1)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ shape_sym[(M-1)/2][q] = shape[(M-1)/2][q];
+
+ double x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = (double)rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i][j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::apply_tensor_product_evenodd<1,M-1,N,double,0,false,add,type>(shape_sym,x,y);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ deallog << y[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ y[i] = (double)rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j][i] * y[j];
+ }
+
+ // apply function for tensor product
+ internal::apply_tensor_product_evenodd<1,M-1,N,double,0,true,add,type>(shape_sym,y,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1e-14);
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ 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;
+}
+
--- /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: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
projects / dealii-svn.git / commitdiff