125 static_assert(n_rows > 0 || n_columns > 0,
126 "Specialization only for n_rows, n_columns > 0");
127 Assert(n_rows > 0 && n_columns > 0,
129 std::to_string(n_rows) +
", " +
130 std::to_string(n_columns) +
" was passed!"));
132 "This function should only use EvaluatorQuantity::value");
134 constexpr int mm = transpose_matrix ? n_rows : n_columns,
135 nn = transpose_matrix ? n_columns : n_rows;
137 std::array<Number, mm> x;
138 for (
int i = 0; i < mm; ++i)
139 x[i] = in[stride_in * i];
140 for (
int col = 0; col < nn; ++col)
143 if (transpose_matrix ==
true)
145 res0 = matrix[col] * x[0];
146 for (
int i = 1; i < mm; ++i)
148 const Number2 mji = matrix[i * n_columns + col];
152 res0.real(res0.real() + mji.real() * x[i].real() -
153 mji.imag() * x[i].imag());
154 res0.imag(res0.imag() + mji.imag() * x[i].real() +
155 mji.real() * x[i].imag());
163 res0 = matrix[col * n_columns] * x[0];
164 for (
int i = 1; i < mm; ++i)
166 const Number2 mij = matrix[col * n_columns + i];
170 res0.real(res0.real() + mij.real() * x[i].real() -
171 mij.imag() * x[i].imag());
172 res0.imag(res0.imag() + mij.imag() * x[i].real() +
173 mij.real() * x[i].imag());
180 out[stride_out * col] += res0;
182 out[stride_out * col] = res0;
205 const int stride_out)
207 const int mm = transpose_matrix ? n_rows : n_columns,
208 nn = transpose_matrix ? n_columns : n_rows;
209 Assert(n_rows > 0 && n_columns > 0,
211 Assert(n_rows > 0 && n_columns > 0,
213 std::to_string(n_rows) +
", " +
214 std::to_string(n_columns) +
" was passed!"));
217 "This function should only use EvaluatorQuantity::value");
224 if (transpose_matrix && n_rows == 2)
226 const Number2 *matrix_1 = matrix + n_columns;
227 const Number x0 = in[0], x1 = in[stride_in];
228 for (
int col = 0; col < nn; ++col)
230 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
232 out[stride_out * col] += result;
234 out[stride_out * col] = result;
237 else if (transpose_matrix && n_rows == 3)
239 const Number2 *matrix_1 = matrix + n_columns;
240 const Number2 *matrix_2 = matrix_1 + n_columns;
241 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
242 for (
int col = 0; col < nn; ++col)
244 const Number result =
245 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
247 out[stride_out * col] += result;
249 out[stride_out * col] = result;
254 std::array<Number, 129> x;
255 for (
int i = 0; i < mm; ++i)
256 x[i] = in[stride_in * i];
258 for (
int col = 0; col < nn; ++col)
261 if (transpose_matrix ==
true)
263 res0 = matrix[col] * x[0];
264 for (
int i = 1; i < mm; ++i)
265 res0 += matrix[i * n_columns + col] * x[i];
269 res0 = matrix[col * n_columns] * x[0];
270 for (
int i = 1; i < mm; ++i)
271 res0 += matrix[col * n_columns + i] * x[i];
274 out[stride_out * col] += res0;
276 out[stride_out * col] = res0;
283 for (
int col = 0; col < nn; ++col)
286 if (transpose_matrix ==
true)
288 res0 = matrix[col] * in[0];
289 for (
int i = 1; i < mm; ++i)
290 res0 += matrix[i * n_columns + col] * in[stride_in * i];
294 res0 = matrix[col * n_columns] * in[0];
295 for (
int i = 1; i < mm; ++i)
296 res0 += matrix[col * n_columns + i] * in[stride_in * i];
299 out[stride_out * col] += res0;
301 out[stride_out * col] = res0;
332 static_assert(n_rows > 0 || n_columns > 0,
333 "Specialization only for n_rows, n_columns > 0");
334 Assert(n_rows > 0 && n_columns > 0,
336 std::to_string(n_rows) +
", " +
337 std::to_string(n_columns) +
" was passed!"));
339 constexpr int mm = transpose_matrix ? n_rows : n_columns,
340 nn = transpose_matrix ? n_columns : n_rows;
341 constexpr int n_cols = nn / 2;
342 constexpr int mid = mm / 2;
344 std::array<Number, mm> x;
345 for (
int i = 0; i < mm; ++i)
346 x[i] = in[stride_in * i];
369 for (
int col = 0; col < n_cols; ++col)
373 if (transpose_matrix ==
true)
376 val1 = matrix[nn - 1 - col];
380 val0 = matrix[col * n_columns];
381 val1 = matrix[(col + 1) * n_columns - 1];
387 res0 += val1 * x[mm - 1];
388 res1 += val0 * x[mm - 1];
389 for (
int ind = 1; ind < mid; ++ind)
391 if (transpose_matrix ==
true)
393 val0 = matrix[ind * n_columns + col];
394 val1 = matrix[ind * n_columns + nn - 1 - col];
398 val0 = matrix[col * n_columns + ind];
399 val1 = matrix[(col + 1) * n_columns - 1 - ind];
401 res0 += val0 * x[ind];
402 res1 += val1 * x[ind];
403 res0 += val1 * x[mm - 1 - ind];
404 res1 += val0 * x[mm - 1 - ind];
408 res0 = res1 = Number();
409 if (transpose_matrix ==
true)
413 const Number tmp = matrix[mid * n_columns + col] * x[mid];
420 if (mm % 2 == 1 && nn % 2 == 0)
422 const Number tmp = matrix[col * n_columns + mid] * x[mid];
429 out[stride_out * col] += res0;
430 out[stride_out * (nn - 1 - col)] += res1;
434 out[stride_out * col] = res0;
435 out[stride_out * (nn - 1 - col)] = res1;
438 if (transpose_matrix ==
true && nn % 2 == 1 && mm % 2 == 1)
441 out[stride_out * n_cols] += x[mid];
443 out[stride_out * n_cols] = x[mid];
445 else if (transpose_matrix ==
true && nn % 2 == 1)
450 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
451 for (
int ind = 1; ind < mid; ++ind)
453 const Number2 val0 = matrix[ind * n_columns + n_cols];
454 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
460 out[stride_out * n_cols] += res0;
462 out[stride_out * n_cols] = res0;
464 else if (transpose_matrix ==
false && nn % 2 == 1)
469 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
470 for (
int ind = 1; ind < mid; ++ind)
472 const Number2 val0 = matrix[n_cols * n_columns + ind];
473 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
482 out[stride_out * n_cols] += res0;
484 out[stride_out * n_cols] = res0;
505 for (
int col = 0; col < n_cols; ++col)
509 if (transpose_matrix ==
true)
512 val1 = matrix[nn - 1 - col];
516 val0 = matrix[col * n_columns];
517 val1 = matrix[(nn - col - 1) * n_columns];
523 res0 -= val1 * x[mm - 1];
524 res1 -= val0 * x[mm - 1];
525 for (
int ind = 1; ind < mid; ++ind)
527 if (transpose_matrix ==
true)
529 val0 = matrix[ind * n_columns + col];
530 val1 = matrix[ind * n_columns + nn - 1 - col];
534 val0 = matrix[col * n_columns + ind];
535 val1 = matrix[(nn - col - 1) * n_columns + ind];
537 res0 += val0 * x[ind];
538 res1 += val1 * x[ind];
539 res0 -= val1 * x[mm - 1 - ind];
540 res1 -= val0 * x[mm - 1 - ind];
544 res0 = res1 = Number();
547 if (transpose_matrix ==
true)
548 val0 = matrix[mid * n_columns + col];
550 val0 = matrix[col * n_columns + mid];
551 const Number tmp = val0 * x[mid];
557 out[stride_out * col] += res0;
558 out[stride_out * (nn - 1 - col)] += res1;
562 out[stride_out * col] = res0;
563 out[stride_out * (nn - 1 - col)] = res1;
570 if (transpose_matrix ==
true)
571 val0 = matrix[n_cols];
573 val0 = matrix[n_cols * n_columns];
574 res0 = val0 * (x[0] - x[mm - 1]);
575 for (
int ind = 1; ind < mid; ++ind)
577 if (transpose_matrix ==
true)
578 val0 = matrix[ind * n_columns + n_cols];
580 val0 = matrix[n_cols * n_columns + ind];
581 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
585 out[stride_out * n_cols] += res0;
587 out[stride_out * n_cols] = res0;
594 for (
int col = 0; col < n_cols; ++col)
598 if (transpose_matrix ==
true)
601 val1 = matrix[nn - 1 - col];
605 val0 = matrix[col * n_columns];
606 val1 = matrix[(col + 1) * n_columns - 1];
612 res0 += val1 * x[mm - 1];
613 res1 += val0 * x[mm - 1];
614 for (
int ind = 1; ind < mid; ++ind)
616 if (transpose_matrix ==
true)
618 val0 = matrix[ind * n_columns + col];
619 val1 = matrix[ind * n_columns + nn - 1 - col];
623 val0 = matrix[col * n_columns + ind];
624 val1 = matrix[(col + 1) * n_columns - 1 - ind];
626 res0 += val0 * x[ind];
627 res1 += val1 * x[ind];
628 res0 += val1 * x[mm - 1 - ind];
629 res1 += val0 * x[mm - 1 - ind];
633 res0 = res1 = Number();
636 if (transpose_matrix ==
true)
637 val0 = matrix[mid * n_columns + col];
639 val0 = matrix[col * n_columns + mid];
640 const Number tmp = val0 * x[mid];
646 out[stride_out * col] += res0;
647 out[stride_out * (nn - 1 - col)] += res1;
651 out[stride_out * col] = res0;
652 out[stride_out * (nn - 1 - col)] = res1;
659 if (transpose_matrix ==
true)
660 val0 = matrix[n_cols];
662 val0 = matrix[n_cols * n_columns];
665 res0 = val0 * (x[0] + x[mm - 1]);
666 for (
int ind = 1; ind < mid; ++ind)
668 if (transpose_matrix ==
true)
669 val0 = matrix[ind * n_columns + n_cols];
671 val0 = matrix[n_cols * n_columns + ind];
672 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
680 if (transpose_matrix ==
true)
681 val0 = matrix[mid * n_columns + n_cols];
683 val0 = matrix[n_cols * n_columns + mid];
684 res0 += val0 * x[mid];
687 out[stride_out * n_cols] += res0;
689 out[stride_out * n_cols] = res0;
731 int n_rows_runtime = 0,
732 int n_columns_runtime = 0,
733 int stride_in_runtime = 0,
734 int stride_out_runtime = 0)
736 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
737 "Negative loop ranges are not allowed!");
739 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
740 const int n_columns =
741 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
742 const int stride_in =
743 n_rows_static == 0 ? stride_in_runtime : stride_in_static;
744 const int stride_out =
745 n_rows_static == 0 ? stride_out_runtime : stride_out_static;
747 Assert(n_rows > 0 && n_columns > 0,
749 std::to_string(n_rows) +
", " +
750 std::to_string(n_columns) +
" was passed!"));
752 const int mm = transpose_matrix ? n_rows : n_columns,
753 nn = transpose_matrix ? n_columns : n_rows;
754 const int n_half = nn / 2;
755 const int m_half = mm / 2;
757 constexpr int array_length =
758 (n_rows_static == 0) ?
761 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
762 const int offset = (n_columns + 1) / 2;
766 std::array<Number, array_length> xp, xm;
767 for (
int i = 0; i < m_half; ++i)
771 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
772 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
776 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
777 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
780 Number xmid = in[stride_in * m_half];
781 for (
int col = 0; col < n_half; ++col)
786 if (transpose_matrix ==
true)
788 r0 = matrix[col] * xp[0];
789 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
793 r0 = matrix[col * offset] * xp[0];
794 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
796 for (
int ind = 1; ind < m_half; ++ind)
798 if (transpose_matrix ==
true)
800 r0 += matrix[ind * offset + col] * xp[ind];
801 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
805 r0 += matrix[col * offset + ind] * xp[ind];
806 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
812 if (mm % 2 == 1 && transpose_matrix ==
true)
815 r1 += matrix[m_half * offset + col] * xmid;
817 r0 += matrix[m_half * offset + col] * xmid;
819 else if (mm % 2 == 1 &&
822 r0 += matrix[col * offset + m_half] * xmid;
826 out[stride_out * col] += r0 + r1;
828 transpose_matrix ==
false)
829 out[stride_out * (nn - 1 - col)] += r1 - r0;
831 out[stride_out * (nn - 1 - col)] += r0 - r1;
835 out[stride_out * col] = r0 + r1;
837 transpose_matrix ==
false)
838 out[stride_out * (nn - 1 - col)] = r1 - r0;
840 out[stride_out * (nn - 1 - col)] = r0 - r1;
844 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
847 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
849 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
851 else if (transpose_matrix ==
true && nn % 2 == 1)
856 r0 = matrix[n_half] * xp[0];
857 for (
int ind = 1; ind < m_half; ++ind)
858 r0 += matrix[ind * offset + n_half] * xp[ind];
863 r0 += matrix[m_half * offset + n_half] * xmid;
866 out[stride_out * n_half] += r0;
868 out[stride_out * n_half] = r0;
870 else if (transpose_matrix ==
false && nn % 2 == 1)
877 r0 = matrix[n_half * offset] * xm[0];
878 for (
int ind = 1; ind < m_half; ++ind)
879 r0 += matrix[n_half * offset + ind] * xm[ind];
883 r0 = matrix[n_half * offset] * xp[0];
884 for (
int ind = 1; ind < m_half; ++ind)
885 r0 += matrix[n_half * offset + ind] * xp[ind];
892 r0 += matrix[n_half * offset + m_half] * xmid;
895 out[stride_out * n_half] += r0;
897 out[stride_out * n_half] = r0;
965 static_assert(n_rows > 0 && n_columns > 0,
966 "Specialization requires n_rows, n_columns > 0");
968 constexpr bool evaluate_antisymmetric =
971 constexpr int mm = transpose_matrix ? n_rows : n_columns,
972 nn = transpose_matrix ? n_columns : n_rows;
973 constexpr int n_half = nn / 2;
974 constexpr int m_half = mm / 2;
976 if (transpose_matrix)
978 std::array<Number, mm> x;
979 for (
unsigned int i = 0; i < mm; ++i)
980 x[i] = in[stride_in * i];
981 for (
unsigned int col = 0; col < n_half; ++col)
986 r0 = matrix[col] * x[0];
987 r1 = matrix[col + n_columns] * x[1];
988 for (
unsigned int ind = 1; ind < m_half; ++ind)
990 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
992 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
998 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
1001 out[stride_out * col] += r0 + r1;
1002 if (evaluate_antisymmetric)
1003 out[stride_out * (nn - 1 - col)] += r1 - r0;
1005 out[stride_out * (nn - 1 - col)] += r0 - r1;
1009 out[stride_out * col] = r0 + r1;
1010 if (evaluate_antisymmetric)
1011 out[stride_out * (nn - 1 - col)] = r1 - r0;
1013 out[stride_out * (nn - 1 - col)] = r0 - r1;
1019 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1022 r0 = matrix[n_half + shift * n_columns] * x[shift];
1023 for (
unsigned int ind = 1; ind < m_half; ++ind)
1024 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1029 if (!evaluate_antisymmetric && mm % 2 == 1)
1030 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1032 out[stride_out * n_half] += r0;
1034 out[stride_out * n_half] = r0;
1039 std::array<Number, m_half + 1> xp, xm;
1040 for (
int i = 0; i < m_half; ++i)
1041 if (!evaluate_antisymmetric)
1043 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1044 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1048 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1049 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1052 xp[m_half] = in[stride_in * m_half];
1053 for (
unsigned int col = 0; col < n_half; ++col)
1058 r0 = matrix[2 * col * n_columns] * xp[0];
1059 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1060 for (
unsigned int ind = 1; ind < m_half; ++ind)
1062 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1063 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1070 if (evaluate_antisymmetric)
1071 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1073 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1077 out[stride_out * (2 * col)] += r0;
1078 out[stride_out * (2 * col + 1)] += r1;
1082 out[stride_out * (2 * col)] = r0;
1083 out[stride_out * (2 * col + 1)] = r1;
1091 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1092 for (
unsigned int ind = 1; ind < m_half; ++ind)
1093 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1097 if (mm % 2 == 1 && !evaluate_antisymmetric)
1098 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1100 out[stride_out * (nn - 1)] += r0;
1102 out[stride_out * (nn - 1)] = r0;
1928 const std::array<int, 2> &n_blocks,
1929 const std::array<int, 2> &steps,
1930 const Number *input,
1932 const int n_rows_runtime = 0,
1933 const int stride_runtime = 1)
1935 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1936 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1938 Number *output1 = output + n_blocks[0] * n_blocks[1];
1939 Number *output2 = output1 + n_blocks[0] * n_blocks[1];
1940 for (
int i2 = 0; i2 < n_blocks[1]; ++i2)
1942 for (
int i1 = 0; i1 < n_blocks[0]; ++i1)
1944 Number res0 = shape_values[0] * input[0];
1946 if (max_derivative > 0)
1947 res1 = shape_values[n_rows] * input[0];
1948 if (max_derivative > 1)
1949 res2 = shape_values[2 * n_rows] * input[0];
1950 for (
int ind = 1; ind < n_rows; ++ind)
1952 res0 += shape_values[ind] * input[stride * ind];
1953 if (max_derivative > 0)
1954 res1 += shape_values[ind + n_rows] * input[stride * ind];
1955 if (max_derivative > 1)
1956 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
1961 if (max_derivative > 0)
1962 output1[i1] += res1;
1963 if (max_derivative > 1)
1964 output2[i2] += res2;
1969 if (max_derivative > 0)
1971 if (max_derivative > 1)
1976 output += n_blocks[0];
1977 if (max_derivative > 0)
1978 output1 += n_blocks[0];
1979 if (max_derivative > 1)
1980 output2 += n_blocks[0];
2151 const unsigned int n_components,
2152 const int n_points_1d_non_template,
2155 const int n_points_1d = n_points_1d_template != -1 ?
2156 n_points_1d_template :
2157 n_points_1d_non_template;
2162 unsigned int compressed_index[100];
2163 compressed_index[0] = 0;
2164 for (
int i = 1; i < n_points_1d - 1; ++i)
2165 compressed_index[i] = 1;
2166 compressed_index[n_points_1d - 1] = 2;
2173 const auto check_and_set = [](Number &weight,
const Number &data) {
2174 if (weight == Number(-1.0) || weight == data)
2183 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2184 weights[c] = Number(-1.0);
2186 for (
unsigned int c = 0; c < n_components; ++c)
2187 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2188 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
2189 ++j, data += n_points_1d)
2191 const unsigned int shift =
2192 9 * compressed_index[k] + 3 * compressed_index[j];
2194 if (!check_and_set(weights[shift], data[0]))
2197 for (
int i = 1; i < n_points_1d - 1; ++i)
2198 if (!check_and_set(weights[shift + 1], data[i]))
2201 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
2213 const unsigned int n_components,
2214 const int n_points_1d_non_template,
2217 const int n_points_1d = n_points_1d_template != -1 ?
2218 n_points_1d_template :
2219 n_points_1d_non_template;
2221 Assert((n_points_1d % 2) == 1,
2222 ExcMessage(
"The function can only with add number of points"));
2226 const unsigned int n_inside_1d = n_points_1d / 2;
2228 unsigned int compressed_index[100];
2231 for (
int i = 0; i < n_inside_1d; ++i)
2232 compressed_index[c++] = 0;
2233 compressed_index[c++] = 1;
2234 for (
int i = 0; i < n_inside_1d; ++i)
2235 compressed_index[c++] = 2;
2242 const auto check_and_set = [](Number &weight,
const Number &data) {
2243 if (weight == Number(-1.0) || weight == data)
2252 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2253 weights[c] = Number(-1.0);
2255 for (
unsigned int comp = 0; comp < n_components; ++comp)
2256 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2257 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
2258 ++j, data += n_points_1d)
2260 const unsigned int shift =
2261 9 * compressed_index[k] + 3 * compressed_index[j];
2265 for (
int i = 0; i < n_inside_1d; ++i)
2266 if (!check_and_set(weights[shift], data[c++]))
2269 if (!check_and_set(weights[shift + 1], data[c++]))
2272 for (
int i = 0; i < n_inside_1d; ++i)
2273 if (!check_and_set(weights[shift + 2], data[c++]))