From: Peter Munch Date: Sat, 3 Sep 2022 12:52:50 +0000 (+0200) Subject: Specialize EvaluatorTensorProduct for non-templated execution X-Git-Tag: v9.5.0-rc1~987^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8e6e6eb08137e98295ad874aa59028b63711e4f3;p=dealii.git Specialize EvaluatorTensorProduct for non-templated execution --- diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index a83591563d..10971a272f 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -1567,6 +1567,218 @@ namespace internal + template + inline void + even_odd_apply(const int n_rows_in, + const int n_columns_in, + const Number2 *DEAL_II_RESTRICT shapes, + const Number * in, + Number * out) + { + static_assert(type < 3, "Only three variants type=0,1,2 implemented"); + static_assert(one_line == false || direction == dim - 1, + "Single-line evaluation only works for direction=dim-1."); + + const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static; + const int n_columns = + n_columns_static == -1 ? n_columns_in : n_columns_static; + + Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || + in != out, + ExcMessage("In-place operation only supported for " + "n_rows==n_columns or single-line interpolation")); + + // We cannot statically assert that direction is less than dim, so must do + // an additional dynamic check + AssertIndexRange(direction, dim); + + const int nn = contract_over_rows ? n_columns : n_rows; + const int mm = contract_over_rows ? n_rows : n_columns; + constexpr int mm_static = + contract_over_rows ? n_rows_static : n_columns_static; + const int n_cols = nn / 2; + const int mid = mm / 2; + constexpr int mid_static = mm_static / 2; + constexpr int max_mid = 15; // for non-templated execution + + Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid, + ExcNotImplemented()); + + const int stride = Utilities::pow(n_columns, direction); + const int n_blocks1 = one_line ? 1 : stride; + const int n_blocks2 = + Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + + const int offset = (n_columns + 1) / 2; + + // 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) + { + constexpr unsigned int mid_size = + (n_rows_static == -1 || n_columns_static == -1) ? + max_mid : + (mid_static > 0 ? mid_static : 1); + Number xp[mid_size], xm[mid_size]; + for (int i = 0; i < mid; ++i) + { + if (contract_over_rows == 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)]; + } + } + Number xmid = in[stride * mid]; + for (int col = 0; col < n_cols; ++col) + { + Number r0, r1; + if (mid > 0) + { + if (contract_over_rows == true) + { + r0 = shapes[col] * xp[0]; + r1 = shapes[(n_rows - 1) * offset + col] * xm[0]; + } + else + { + r0 = shapes[col * offset] * xp[0]; + r1 = shapes[(n_rows - 1 - col) * offset] * xm[0]; + } + for (int ind = 1; ind < mid; ++ind) + { + if (contract_over_rows == true) + { + r0 += shapes[ind * offset + col] * xp[ind]; + r1 += shapes[(n_rows - 1 - ind) * offset + col] * + xm[ind]; + } + else + { + r0 += shapes[col * offset + ind] * xp[ind]; + r1 += shapes[(n_rows - 1 - col) * offset + ind] * + xm[ind]; + } + } + } + else + r0 = r1 = Number(); + if (mm % 2 == 1 && contract_over_rows == true) + { + if (type == 1) + r1 += shapes[mid * offset + col] * xmid; + else + r0 += shapes[mid * offset + col] * xmid; + } + else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3)) + r0 += shapes[col * offset + mid] * xmid; + + if (add) + { + out[stride * col] += r0 + r1; + if (type == 1 && contract_over_rows == 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 && contract_over_rows == false) + out[stride * (nn - 1 - col)] = r1 - r0; + else + out[stride * (nn - 1 - col)] = r0 - r1; + } + } + if (type == 0 && contract_over_rows == true && nn % 2 == 1 && + mm % 2 == 1 && mm > 3) + { + if (add) + out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid; + else + out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid; + } + else if (contract_over_rows == true && nn % 2 == 1) + { + Number r0; + if (mid > 0) + { + r0 = shapes[n_cols] * xp[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += shapes[ind * offset + n_cols] * xp[ind]; + } + else + r0 = Number(); + if (type != 1 && mm % 2 == 1) + r0 += shapes[mid * offset + n_cols] * xmid; + + if (add) + out[stride * n_cols] += r0; + else + out[stride * n_cols] = r0; + } + else if (contract_over_rows == false && nn % 2 == 1) + { + Number r0; + if (mid > 0) + { + if (type == 1) + { + r0 = shapes[n_cols * offset] * xm[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += shapes[n_cols * offset + ind] * xm[ind]; + } + else + { + r0 = shapes[n_cols * offset] * xp[0]; + for (int ind = 1; ind < mid; ++ind) + r0 += shapes[n_cols * offset + ind] * xp[ind]; + } + } + else + r0 = Number(); + + if ((type == 0 || type == 2) && mm % 2 == 1) + r0 += shapes[n_cols * offset + mid] * xmid; + + if (add) + out[stride * n_cols] += r0; + else + out[stride * n_cols] = r0; + } + if (one_line == false) + { + in += 1; + out += 1; + } + } + if (one_line == false) + { + in += stride * (mm - 1); + out += stride * (nn - 1); + } + } + } + + + /** * Internal evaluator for 1d-3d shape function using the tensor product form * of the basis functions. @@ -1751,7 +1963,19 @@ namespace internal static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number * in, - Number * out); + Number * out) + { + even_odd_apply(n_rows, n_columns, shape_data, in, out); + } private: const Number2 *shape_values; @@ -1760,203 +1984,133 @@ namespace internal }; - - template - template - inline void - EvaluatorTensorProduct::apply(const Number2 *DEAL_II_RESTRICT shapes, - const Number * in, - Number * out) + /** + * Internal evaluator for shape function using the tensor product form + * of the basis functions. The same as the other templated class but + * without making use of template arguments and variable loop bounds + * instead. + */ + template + struct EvaluatorTensorProduct { - static_assert(type < 3, "Only three variants type=0,1,2 implemented"); - static_assert(one_line == false || direction == dim - 1, - "Single-line evaluation only works for direction=dim-1."); - Assert(dim == direction + 1 || one_line == true || n_rows == n_columns || - in != out, - ExcMessage("In-place operation only supported for " - "n_rows==n_columns or single-line interpolation")); + EvaluatorTensorProduct() + : shape_values(nullptr) + , shape_gradients(nullptr) + , shape_hessians(nullptr) + , n_rows(numbers::invalid_unsigned_int) + , n_columns(numbers::invalid_unsigned_int) + {} - // We cannot statically assert that direction is less than dim, so must do - // an additional dynamic check - AssertIndexRange(direction, dim); + EvaluatorTensorProduct(const AlignedVector &shape_values, + const unsigned int n_rows = 0, + const unsigned int n_columns = 0) + : shape_values(shape_values.begin()) + , shape_gradients(nullptr) + , shape_hessians(nullptr) + { + AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); + } - constexpr int nn = contract_over_rows ? n_columns : n_rows; - constexpr int mm = contract_over_rows ? n_rows : n_columns; - constexpr int n_cols = nn / 2; - constexpr int mid = mm / 2; + EvaluatorTensorProduct(const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &shape_hessians, + const unsigned int n_rows = 0, + const unsigned int n_columns = 0) + : shape_values(shape_values.begin()) + , shape_gradients(shape_gradients.begin()) + , shape_hessians(shape_hessians.begin()) + , n_rows(n_rows) + , n_columns(n_columns) + { + if (!shape_values.empty()) + AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2)); + if (!shape_gradients.empty()) + AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2)); + if (!shape_hessians.empty()) + AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2)); + } - constexpr int stride = Utilities::pow(n_columns, direction); - constexpr int n_blocks1 = one_line ? 1 : stride; - constexpr int n_blocks2 = - Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + template + void + values(const Number in[], Number out[]) const + { + Assert(shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } - constexpr int offset = (n_columns + 1) / 2; + template + void + gradients(const Number in[], Number out[]) const + { + Assert(shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } - // 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 > 0 ? mid : 1], xm[mid > 0 ? mid : 1]; - for (int i = 0; i < mid; ++i) - { - if (contract_over_rows == 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)]; - } - } - Number xmid = in[stride * mid]; - for (int col = 0; col < n_cols; ++col) - { - Number r0, r1; - if (mid > 0) - { - if (contract_over_rows == true) - { - r0 = shapes[col] * xp[0]; - r1 = shapes[(n_rows - 1) * offset + col] * xm[0]; - } - else - { - r0 = shapes[col * offset] * xp[0]; - r1 = shapes[(n_rows - 1 - col) * offset] * xm[0]; - } - for (int ind = 1; ind < mid; ++ind) - { - if (contract_over_rows == true) - { - r0 += shapes[ind * offset + col] * xp[ind]; - r1 += shapes[(n_rows - 1 - ind) * offset + col] * - xm[ind]; - } - else - { - r0 += shapes[col * offset + ind] * xp[ind]; - r1 += shapes[(n_rows - 1 - col) * offset + ind] * - xm[ind]; - } - } - } - else - r0 = r1 = Number(); - if (mm % 2 == 1 && contract_over_rows == true) - { - if (type == 1) - r1 += shapes[mid * offset + col] * xmid; - else - r0 += shapes[mid * offset + col] * xmid; - } - else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3)) - r0 += shapes[col * offset + mid] * xmid; + template + void + hessians(const Number in[], Number out[]) const + { + Assert(shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, in, out); + } - if (add) - { - out[stride * col] += r0 + r1; - if (type == 1 && contract_over_rows == 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 && contract_over_rows == false) - out[stride * (nn - 1 - col)] = r1 - r0; - else - out[stride * (nn - 1 - col)] = r0 - r1; - } - } - if (type == 0 && contract_over_rows == true && nn % 2 == 1 && - mm % 2 == 1 && mm > 3) - { - if (add) - out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid; - else - out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid; - } - else if (contract_over_rows == true && nn % 2 == 1) - { - Number r0; - if (mid > 0) - { - r0 = shapes[n_cols] * xp[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[ind * offset + n_cols] * xp[ind]; - } - else - r0 = Number(); - if (type != 1 && mm % 2 == 1) - r0 += shapes[mid * offset + n_cols] * xmid; + template + void + values_one_line(const Number in[], Number out[]) const + { + Assert(shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } - if (add) - out[stride * n_cols] += r0; - else - out[stride * n_cols] = r0; - } - else if (contract_over_rows == false && nn % 2 == 1) - { - Number r0; - if (mid > 0) - { - if (type == 1) - { - r0 = shapes[n_cols * offset] * xm[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[n_cols * offset + ind] * xm[ind]; - } - else - { - r0 = shapes[n_cols * offset] * xp[0]; - for (int ind = 1; ind < mid; ++ind) - r0 += shapes[n_cols * offset + ind] * xp[ind]; - } - } - else - r0 = Number(); + template + void + gradients_one_line(const Number in[], Number out[]) const + { + Assert(shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, + in, + out); + } - if ((type == 0 || type == 2) && mm % 2 == 1) - r0 += shapes[n_cols * offset + mid] * xmid; + template + void + hessians_one_line(const Number in[], Number out[]) const + { + Assert(shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, + in, + out); + } - if (add) - out[stride * n_cols] += r0; - else - out[stride * n_cols] = r0; - } - if (one_line == false) - { - in += 1; - out += 1; - } - } - if (one_line == false) - { - in += stride * (mm - 1); - out += stride * (nn - 1); - } - } - } + template + void + apply(const Number2 *DEAL_II_RESTRICT shape_data, + const Number * in, + Number * out) const + { + even_odd_apply(n_rows, n_columns, shape_data, in, out); + } + + private: + const Number2 * shape_values; + const Number2 * shape_gradients; + const Number2 * shape_hessians; + const unsigned int n_rows; + const unsigned int n_columns; + };