deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_product_kernels.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_tensor_product_kernels_h
17#define dealii_matrix_free_tensor_product_kernels_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/aligned_vector.h>
22#include <deal.II/base/utilities.h>
23
24#include <deal.II/matrix_free/shape_info.h>
25
26
27DEAL_II_NAMESPACE_OPEN
28
29
30
31namespace internal
32{
68
69
70
74 enum class EvaluatorQuantity
75 {
79 value,
83 gradient,
87 hessian
88 };
89
90
91
107 template <EvaluatorVariant variant,
108 EvaluatorQuantity quantity,
109 int n_rows,
110 int n_columns,
111 int stride_in,
112 int stride_out,
113 bool transpose_matrix,
114 bool add,
115 typename Number,
116 typename Number2>
117 std::enable_if_t<(variant == evaluate_general), void>
118 apply_matrix_vector_product(const Number2 *matrix,
119 const Number *in,
120 Number *out)
121 {
122 // We can only statically assert that one argument is non-zero because
123 // face evaluation might instantiate some functions, so we need to use the
124 // run-time assert to verify that we do not end up involuntarily.
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,
128 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
129 std::to_string(n_rows) + ", " +
130 std::to_string(n_columns) + " was passed!"));
131 static_assert(quantity == EvaluatorQuantity::value,
132 "This function should only use EvaluatorQuantity::value");
133
134 constexpr int mm = transpose_matrix ? n_rows : n_columns,
135 nn = transpose_matrix ? n_columns : n_rows;
136
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)
141 {
142 Number res0;
143 if (transpose_matrix == true)
144 {
145 res0 = matrix[col] * x[0];
146 for (int i = 1; i < mm; ++i)
147 {
148 const Number2 mji = matrix[i * n_columns + col];
149 if constexpr (numbers::NumberTraits<Number>::is_complex &&
150 numbers::NumberTraits<Number2>::is_complex)
151 {
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());
156 }
157 else
158 res0 += mji * x[i];
159 }
160 }
161 else
162 {
163 res0 = matrix[col * n_columns] * x[0];
164 for (int i = 1; i < mm; ++i)
165 {
166 const Number2 mij = matrix[col * n_columns + i];
167 if constexpr (numbers::NumberTraits<Number>::is_complex &&
168 numbers::NumberTraits<Number2>::is_complex)
169 {
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());
174 }
175 else
176 res0 += mij * x[i];
177 }
178 }
179 if (add)
180 out[stride_out * col] += res0;
181 else
182 out[stride_out * col] = res0;
183 }
184 }
185
186
187
192 template <EvaluatorVariant variant,
193 EvaluatorQuantity quantity,
194 bool transpose_matrix,
195 bool add,
196 bool consider_strides,
197 typename Number,
198 typename Number2,
199 int n_components = 1>
200 std::enable_if_t<(variant == evaluate_general), void>
201 apply_matrix_vector_product(const Number2 *matrix,
202 const Number *in,
203 Number *out,
204 const int n_rows,
205 const int n_columns,
206 const int stride_in_given,
207 const int stride_out_given)
208 {
209 const int mm = transpose_matrix ? n_rows : n_columns,
210 nn = transpose_matrix ? n_columns : n_rows;
211 Assert(n_rows > 0 && n_columns > 0,
212 ExcInternalError("Empty evaluation task!"));
213 Assert(n_rows > 0 && n_columns > 0,
214 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
215 std::to_string(n_rows) + ", " +
216 std::to_string(n_columns) + " was passed!"));
217
218 static_assert(quantity == EvaluatorQuantity::value,
219 "This function should only use EvaluatorQuantity::value");
220
221 Assert(consider_strides || (stride_in_given == 1 && stride_out_given == 1),
222 ExcInternalError());
223 const int stride_in = consider_strides ? stride_in_given : 1;
224 const int stride_out = consider_strides ? stride_out_given : 1;
225
226 static_assert(n_components > 0 && n_components < 4,
227 "Invalid number of components");
228
229 // specialization for n_rows = 2 that manually unrolls the innermost loop
230 // to make the operation perform better (not completely as good as the
231 // templated one, but much better than the generic version down below,
232 // because the loop over col can be more effectively unrolled by the
233 // compiler)
234 if (transpose_matrix && n_rows == 2 && n_components == 1)
235 {
236 const Number2 *matrix_1 = matrix + n_columns;
237 const Number x0 = in[0], x1 = in[stride_in];
238 for (int col = 0; col < nn; ++col)
239 {
240 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
241 if (add)
242 out[stride_out * col] += result;
243 else
244 out[stride_out * col] = result;
245 }
246 }
247 else if (transpose_matrix && n_rows == 3 && n_components == 1)
248 {
249 const Number2 *matrix_1 = matrix + n_columns;
250 const Number2 *matrix_2 = matrix_1 + n_columns;
251 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
252 for (int col = 0; col < nn; ++col)
253 {
254 const Number result =
255 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
256 if (add)
257 out[stride_out * col] += result;
258 else
259 out[stride_out * col] = result;
260 }
261 }
262 else if (std::abs(in - out) < std::min(stride_out * nn, stride_in * mm) &&
263 n_components == 1)
264 {
265 Assert(mm <= 128,
266 ExcNotImplemented("For large sizes, arrays may not overlap"));
267 std::array<Number, 129> x;
268 for (int i = 0; i < mm; ++i)
269 x[i] = in[stride_in * i];
270
271 for (int col = 0; col < nn; ++col)
272 {
273 Number res0;
274 if (transpose_matrix == true)
275 {
276 res0 = matrix[col] * x[0];
277 for (int i = 1; i < mm; ++i)
278 res0 += matrix[i * n_columns + col] * x[i];
279 }
280 else
281 {
282 res0 = matrix[col * n_columns] * x[0];
283 for (int i = 1; i < mm; ++i)
284 res0 += matrix[col * n_columns + i] * x[i];
285 }
286 if (add)
287 out[stride_out * col] += res0;
288 else
289 out[stride_out * col] = res0;
290 }
291 }
292 else
293 {
294 const Number *in0 = in;
295 const Number *in1 = n_components > 1 ? in + mm : nullptr;
296 const Number *in2 = n_components > 2 ? in + 2 * mm : nullptr;
297
298 Number *out0 = out;
299 Number *out1 = n_components > 1 ? out + nn : nullptr;
300 Number *out2 = n_components > 2 ? out + 2 * nn : nullptr;
301
302 int nn_regular = (nn / 4) * 4;
303 for (int col = 0; col < nn_regular; col += 4)
304 {
305 Number res[12];
306 if (transpose_matrix == true)
307 {
308 const Number2 *matrix_ptr = matrix + col;
309 const Number a = in0[0];
310 res[0] = matrix_ptr[0] * a;
311 res[1] = matrix_ptr[1] * a;
312 res[2] = matrix_ptr[2] * a;
313 res[3] = matrix_ptr[3] * a;
314
315 if (n_components > 1)
316 {
317 const Number b = in1[0];
318 res[4] = matrix_ptr[0] * b;
319 res[5] = matrix_ptr[1] * b;
320 res[6] = matrix_ptr[2] * b;
321 res[7] = matrix_ptr[3] * b;
322 }
323
324 if (n_components > 2)
325 {
326 const Number c = in2[0];
327 res[8] = matrix_ptr[0] * c;
328 res[9] = matrix_ptr[1] * c;
329 res[10] = matrix_ptr[2] * c;
330 res[11] = matrix_ptr[3] * c;
331 }
332
333 matrix_ptr += n_columns;
334 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
335 {
336 const Number a = in0[stride_in * i];
337 res[0] += matrix_ptr[0] * a;
338 res[1] += matrix_ptr[1] * a;
339 res[2] += matrix_ptr[2] * a;
340 res[3] += matrix_ptr[3] * a;
341
342 if (n_components > 1)
343 {
344 const Number b = in1[stride_in * i];
345 res[4] += matrix_ptr[0] * b;
346 res[5] += matrix_ptr[1] * b;
347 res[6] += matrix_ptr[2] * b;
348 res[7] += matrix_ptr[3] * b;
349 }
350 if (n_components > 2)
351 {
352 const Number c = in2[stride_in * i];
353 res[8] += matrix_ptr[0] * c;
354 res[9] += matrix_ptr[1] * c;
355 res[10] += matrix_ptr[2] * c;
356 res[11] += matrix_ptr[3] * c;
357 }
358 }
359 }
360 else
361 {
362 const Number2 *matrix_0 = matrix + col * n_columns;
363 const Number2 *matrix_1 = matrix + (col + 1) * n_columns;
364 const Number2 *matrix_2 = matrix + (col + 2) * n_columns;
365 const Number2 *matrix_3 = matrix + (col + 3) * n_columns;
366
367 const Number a = in0[0];
368 res[0] = matrix_0[0] * a;
369 res[1] = matrix_1[0] * a;
370 res[2] = matrix_2[0] * a;
371 res[3] = matrix_3[0] * a;
372
373 if (n_components > 1)
374 {
375 const Number b = in1[0];
376 res[4] = matrix_0[0] * b;
377 res[5] = matrix_1[0] * b;
378 res[6] = matrix_2[0] * b;
379 res[7] = matrix_3[0] * b;
380 }
381
382 if (n_components > 2)
383 {
384 const Number c = in2[0];
385 res[8] = matrix_0[0] * c;
386 res[9] = matrix_1[0] * c;
387 res[10] = matrix_2[0] * c;
388 res[11] = matrix_3[0] * c;
389 }
390
391 for (int i = 1; i < mm; ++i)
392 {
393 const Number a = in0[stride_in * i];
394 res[0] += matrix_0[i] * a;
395 res[1] += matrix_1[i] * a;
396 res[2] += matrix_2[i] * a;
397 res[3] += matrix_3[i] * a;
398
399 if (n_components > 1)
400 {
401 const Number b = in1[stride_in * i];
402 res[4] += matrix_0[i] * b;
403 res[5] += matrix_1[i] * b;
404 res[6] += matrix_2[i] * b;
405 res[7] += matrix_3[i] * b;
406 }
407
408 if (n_components > 2)
409 {
410 const Number c = in2[stride_in * i];
411 res[8] += matrix_0[i] * c;
412 res[9] += matrix_1[i] * c;
413 res[10] += matrix_2[i] * c;
414 res[11] += matrix_3[i] * c;
415 }
416 }
417 }
418 if (add)
419 {
420 out0[0] += res[0];
421 out0[stride_out] += res[1];
422 out0[2 * stride_out] += res[2];
423 out0[3 * stride_out] += res[3];
424 if (n_components > 1)
425 {
426 out1[0] += res[4];
427 out1[stride_out] += res[5];
428 out1[2 * stride_out] += res[6];
429 out1[3 * stride_out] += res[7];
430 }
431 if (n_components > 2)
432 {
433 out2[0] += res[8];
434 out2[stride_out] += res[9];
435 out2[2 * stride_out] += res[10];
436 out2[3 * stride_out] += res[11];
437 }
438 }
439 else
440 {
441 out0[0] = res[0];
442 out0[stride_out] = res[1];
443 out0[2 * stride_out] = res[2];
444 out0[3 * stride_out] = res[3];
445 if (n_components > 1)
446 {
447 out1[0] = res[4];
448 out1[stride_out] = res[5];
449 out1[2 * stride_out] = res[6];
450 out1[3 * stride_out] = res[7];
451 }
452 if (n_components > 2)
453 {
454 out2[0] = res[8];
455 out2[stride_out] = res[9];
456 out2[2 * stride_out] = res[10];
457 out2[3 * stride_out] = res[11];
458 }
459 }
460 out0 += 4 * stride_out;
461 if (n_components > 1)
462 out1 += 4 * stride_out;
463 if (n_components > 2)
464 out2 += 4 * stride_out;
465 }
466 if (nn - nn_regular == 3)
467 {
468 Number res0, res1, res2, res3, res4, res5, res6, res7, res8;
469 if (transpose_matrix == true)
470 {
471 const Number2 *matrix_ptr = matrix + nn_regular;
472 res0 = matrix_ptr[0] * in0[0];
473 res1 = matrix_ptr[1] * in0[0];
474 res2 = matrix_ptr[2] * in0[0];
475 if (n_components > 1)
476 {
477 res3 = matrix_ptr[0] * in1[0];
478 res4 = matrix_ptr[1] * in1[0];
479 res5 = matrix_ptr[2] * in1[0];
480 }
481 if (n_components > 2)
482 {
483 res6 = matrix_ptr[0] * in2[0];
484 res7 = matrix_ptr[1] * in2[0];
485 res8 = matrix_ptr[2] * in2[0];
486 }
487 matrix_ptr += n_columns;
488 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
489 {
490 res0 += matrix_ptr[0] * in0[stride_in * i];
491 res1 += matrix_ptr[1] * in0[stride_in * i];
492 res2 += matrix_ptr[2] * in0[stride_in * i];
493 if (n_components > 1)
494 {
495 res3 += matrix_ptr[0] * in1[stride_in * i];
496 res4 += matrix_ptr[1] * in1[stride_in * i];
497 res5 += matrix_ptr[2] * in1[stride_in * i];
498 }
499 if (n_components > 2)
500 {
501 res6 += matrix_ptr[0] * in2[stride_in * i];
502 res7 += matrix_ptr[1] * in2[stride_in * i];
503 res8 += matrix_ptr[2] * in2[stride_in * i];
504 }
505 }
506 }
507 else
508 {
509 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
510 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
511 const Number2 *matrix_2 = matrix + (nn_regular + 2) * n_columns;
512
513 res0 = matrix_0[0] * in0[0];
514 res1 = matrix_1[0] * in0[0];
515 res2 = matrix_2[0] * in0[0];
516 if (n_components > 1)
517 {
518 res3 = matrix_0[0] * in1[0];
519 res4 = matrix_1[0] * in1[0];
520 res5 = matrix_2[0] * in1[0];
521 }
522 if (n_components > 2)
523 {
524 res6 = matrix_0[0] * in2[0];
525 res7 = matrix_1[0] * in2[0];
526 res8 = matrix_2[0] * in2[0];
527 }
528 for (int i = 1; i < mm; ++i)
529 {
530 res0 += matrix_0[i] * in0[stride_in * i];
531 res1 += matrix_1[i] * in0[stride_in * i];
532 res2 += matrix_2[i] * in0[stride_in * i];
533 if (n_components > 1)
534 {
535 res3 += matrix_0[i] * in1[stride_in * i];
536 res4 += matrix_1[i] * in1[stride_in * i];
537 res5 += matrix_2[i] * in1[stride_in * i];
538 }
539 if (n_components > 2)
540 {
541 res6 += matrix_0[i] * in2[stride_in * i];
542 res7 += matrix_1[i] * in2[stride_in * i];
543 res8 += matrix_2[i] * in2[stride_in * i];
544 }
545 }
546 }
547 if (add)
548 {
549 out0[0] += res0;
550 out0[stride_out] += res1;
551 out0[2 * stride_out] += res2;
552 if (n_components > 1)
553 {
554 out1[0] += res3;
555 out1[stride_out] += res4;
556 out1[2 * stride_out] += res5;
557 }
558 if (n_components > 2)
559 {
560 out2[0] += res6;
561 out2[stride_out] += res7;
562 out2[2 * stride_out] += res8;
563 }
564 }
565 else
566 {
567 out0[0] = res0;
568 out0[stride_out] = res1;
569 out0[2 * stride_out] = res2;
570 if (n_components > 1)
571 {
572 out1[0] = res3;
573 out1[stride_out] = res4;
574 out1[2 * stride_out] = res5;
575 }
576 if (n_components > 2)
577 {
578 out2[0] = res6;
579 out2[stride_out] = res7;
580 out2[2 * stride_out] = res8;
581 }
582 }
583 }
584 else if (nn - nn_regular == 2)
585 {
586 Number res0, res1, res2, res3, res4, res5;
587 if (transpose_matrix == true)
588 {
589 const Number2 *matrix_ptr = matrix + nn_regular;
590 res0 = matrix_ptr[0] * in0[0];
591 res1 = matrix_ptr[1] * in0[0];
592 if (n_components > 1)
593 {
594 res2 = matrix_ptr[0] * in1[0];
595 res3 = matrix_ptr[1] * in1[0];
596 }
597 if (n_components > 2)
598 {
599 res4 = matrix_ptr[0] * in2[0];
600 res5 = matrix_ptr[1] * in2[0];
601 }
602 matrix_ptr += n_columns;
603 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
604 {
605 res0 += matrix_ptr[0] * in0[stride_in * i];
606 res1 += matrix_ptr[1] * in0[stride_in * i];
607 if (n_components > 1)
608 {
609 res2 += matrix_ptr[0] * in1[stride_in * i];
610 res3 += matrix_ptr[1] * in1[stride_in * i];
611 }
612 if (n_components > 2)
613 {
614 res4 += matrix_ptr[0] * in2[stride_in * i];
615 res5 += matrix_ptr[1] * in2[stride_in * i];
616 }
617 }
618 }
619 else
620 {
621 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
622 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
623
624 res0 = matrix_0[0] * in0[0];
625 res1 = matrix_1[0] * in0[0];
626 if (n_components > 1)
627 {
628 res2 = matrix_0[0] * in1[0];
629 res3 = matrix_1[0] * in1[0];
630 }
631 if (n_components > 2)
632 {
633 res4 = matrix_0[0] * in2[0];
634 res5 = matrix_1[0] * in2[0];
635 }
636 for (int i = 1; i < mm; ++i)
637 {
638 res0 += matrix_0[i] * in0[stride_in * i];
639 res1 += matrix_1[i] * in0[stride_in * i];
640 if (n_components > 1)
641 {
642 res2 += matrix_0[i] * in1[stride_in * i];
643 res3 += matrix_1[i] * in1[stride_in * i];
644 }
645 if (n_components > 2)
646 {
647 res4 += matrix_0[i] * in2[stride_in * i];
648 res5 += matrix_1[i] * in2[stride_in * i];
649 }
650 }
651 }
652 if (add)
653 {
654 out0[0] += res0;
655 out0[stride_out] += res1;
656 if (n_components > 1)
657 {
658 out1[0] += res2;
659 out1[stride_out] += res3;
660 }
661 if (n_components > 2)
662 {
663 out2[0] += res4;
664 out2[stride_out] += res5;
665 }
666 }
667 else
668 {
669 out0[0] = res0;
670 out0[stride_out] = res1;
671 if (n_components > 1)
672 {
673 out1[0] = res2;
674 out1[stride_out] = res3;
675 }
676 if (n_components > 2)
677 {
678 out2[0] = res4;
679 out2[stride_out] = res5;
680 }
681 }
682 }
683 else if (nn - nn_regular == 1)
684 {
685 Number res0, res1, res2;
686 if (transpose_matrix == true)
687 {
688 const Number2 *matrix_ptr = matrix + nn_regular;
689 res0 = matrix_ptr[0] * in0[0];
690 if (n_components > 1)
691 res1 = matrix_ptr[0] * in1[0];
692 if (n_components > 2)
693 res2 = matrix_ptr[0] * in2[0];
694 matrix_ptr += n_columns;
695 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
696 {
697 res0 += matrix_ptr[0] * in0[stride_in * i];
698 if (n_components > 1)
699 res1 += matrix_ptr[0] * in1[stride_in * i];
700 if (n_components > 2)
701 res2 += matrix_ptr[0] * in2[stride_in * i];
702 }
703 }
704 else
705 {
706 const Number2 *matrix_ptr = matrix + nn_regular * n_columns;
707 res0 = matrix_ptr[0] * in0[0];
708 if (n_components > 1)
709 res1 = matrix_ptr[0] * in1[0];
710 if (n_components > 2)
711 res2 = matrix_ptr[0] * in2[0];
712 for (int i = 1; i < mm; ++i)
713 {
714 res0 += matrix_ptr[i] * in0[stride_in * i];
715 if (n_components > 1)
716 res1 += matrix_ptr[i] * in1[stride_in * i];
717 if (n_components > 2)
718 res2 += matrix_ptr[i] * in2[stride_in * i];
719 }
720 }
721 if (add)
722 {
723 out0[0] += res0;
724 if (n_components > 1)
725 out1[0] += res1;
726 if (n_components > 2)
727 out2[0] += res2;
728 }
729 else
730 {
731 out0[0] = res0;
732 if (n_components > 1)
733 out1[0] = res1;
734 if (n_components > 2)
735 out2[0] = res2;
736 }
737 }
738 }
739 }
740
741
742
749 template <EvaluatorVariant variant,
750 EvaluatorQuantity quantity,
751 int n_rows,
752 int n_columns,
753 int stride_in,
754 int stride_out,
755 bool transpose_matrix,
756 bool add,
757 typename Number,
758 typename Number2>
759 std::enable_if_t<(variant == evaluate_symmetric), void>
760 apply_matrix_vector_product(const Number2 *matrix,
761 const Number *in,
762 Number *out)
763 {
764 // We can only statically assert that one argument is non-zero because
765 // face evaluation might instantiate some functions, so we need to use the
766 // run-time assert to verify that we do not end up involuntarily.
767 static_assert(n_rows > 0 || n_columns > 0,
768 "Specialization only for n_rows, n_columns > 0");
769 Assert(n_rows > 0 && n_columns > 0,
770 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
771 std::to_string(n_rows) + ", " +
772 std::to_string(n_columns) + " was passed!"));
773
774 constexpr int mm = transpose_matrix ? n_rows : n_columns,
775 nn = transpose_matrix ? n_columns : n_rows;
776 constexpr int n_cols = nn / 2;
777 constexpr int mid = mm / 2;
778
779 std::array<Number, mm> x;
780 for (int i = 0; i < mm; ++i)
781 x[i] = in[stride_in * i];
782
783 if (quantity == EvaluatorQuantity::value)
784 {
785 // In this case, the 1d shape values read (sorted lexicographically,
786 // rows run over 1d dofs, columns over quadrature points):
787 // Q2 --> [ 0.687 0 -0.087 ]
788 // [ 0.4 1 0.4 ]
789 // [-0.087 0 0.687 ]
790 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
791 // [ 0.521 1.005 -0.01 -0.230 ]
792 // [-0.230 -0.01 1.005 0.521 ]
793 // [ 0.049 0.002 0.003 0.66 ]
794 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
795 // [ 0.608 1.059 0 0.039 0.176 ]
796 // [-0.409 -0.113 1 -0.113 -0.409 ]
797 // [ 0.176 0.039 0 1.059 0.608 ]
798 // [-0.032 -0.007 0 0.022 0.658 ]
799 //
800 // In these matrices, we want to use avoid computations involving
801 // zeros and ones and use the symmetry in entries starting from (1,1)
802 // forward and (N,N) backward, respectively to reduce the number of
803 // read operations.
804 for (int col = 0; col < n_cols; ++col)
805 {
806 Number2 val0, val1;
807 Number res0, res1;
808 if (transpose_matrix == true)
809 {
810 val0 = matrix[col];
811 val1 = matrix[nn - 1 - col];
812 }
813 else
814 {
815 val0 = matrix[col * n_columns];
816 val1 = matrix[(col + 1) * n_columns - 1];
817 }
818 if (mid > 0)
819 {
820 res0 = val0 * x[0];
821 res1 = val1 * x[0];
822 res0 += val1 * x[mm - 1];
823 res1 += val0 * x[mm - 1];
824 for (int ind = 1; ind < mid; ++ind)
825 {
826 if (transpose_matrix == true)
827 {
828 val0 = matrix[ind * n_columns + col];
829 val1 = matrix[ind * n_columns + nn - 1 - col];
830 }
831 else
832 {
833 val0 = matrix[col * n_columns + ind];
834 val1 = matrix[(col + 1) * n_columns - 1 - ind];
835 }
836 res0 += val0 * x[ind];
837 res1 += val1 * x[ind];
838 res0 += val1 * x[mm - 1 - ind];
839 res1 += val0 * x[mm - 1 - ind];
840 }
841 }
842 else
843 res0 = res1 = Number();
844 if (transpose_matrix == true)
845 {
846 if (mm % 2 == 1)
847 {
848 const Number tmp = matrix[mid * n_columns + col] * x[mid];
849 res0 += tmp;
850 res1 += tmp;
851 }
852 }
853 else
854 {
855 if (mm % 2 == 1 && nn % 2 == 0)
856 {
857 const Number tmp = matrix[col * n_columns + mid] * x[mid];
858 res0 += tmp;
859 res1 += tmp;
860 }
861 }
862 if (add)
863 {
864 out[stride_out * col] += res0;
865 out[stride_out * (nn - 1 - col)] += res1;
866 }
867 else
868 {
869 out[stride_out * col] = res0;
870 out[stride_out * (nn - 1 - col)] = res1;
871 }
872 }
873 if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
874 {
875 if (add)
876 out[stride_out * n_cols] += x[mid];
877 else
878 out[stride_out * n_cols] = x[mid];
879 }
880 else if (transpose_matrix == true && nn % 2 == 1)
881 {
882 Number res0;
883 if (mid > 0)
884 {
885 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
886 for (int ind = 1; ind < mid; ++ind)
887 {
888 const Number2 val0 = matrix[ind * n_columns + n_cols];
889 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
890 }
891 }
892 else
893 res0 = Number();
894 if (add)
895 out[stride_out * n_cols] += res0;
896 else
897 out[stride_out * n_cols] = res0;
898 }
899 else if (transpose_matrix == false && nn % 2 == 1)
900 {
901 Number res0;
902 if (mid > 0)
903 {
904 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
905 for (int ind = 1; ind < mid; ++ind)
906 {
907 const Number2 val0 = matrix[n_cols * n_columns + ind];
908 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
909 ;
910 }
911 if (mm % 2)
912 res0 += x[mid];
913 }
914 else
915 res0 = in[0];
916 if (add)
917 out[stride_out * n_cols] += res0;
918 else
919 out[stride_out * n_cols] = res0;
920 }
921 }
922 else if (quantity == EvaluatorQuantity::gradient)
923 {
924 // For the specialized loop used for gradient computations we again
925 // exploit symmetries according to the following entries (sorted
926 // lexicographically, rows run over 1d dofs, columns over quadrature
927 // points):
928 // Q2 --> [-2.549 -1 0.549 ]
929 // [ 3.098 0 -3.098 ]
930 // [-0.549 1 2.549 ]
931 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
932 // [ 6.07 -1.44 -2.97 2.196 ]
933 // [-2.196 2.97 1.44 -6.07 ]
934 // [ 0.44 -0.5 1.03 4.315 ]
935 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
936 // [10.111 -2.76 -2.667 2.066 -2.306 ]
937 // [-5.688 5.773 0 -5.773 5.688 ]
938 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
939 // [-0.413 0.353 -0.333 -0.353 0.413 ]
940 for (int col = 0; col < n_cols; ++col)
941 {
942 Number2 val0, val1;
943 Number res0, res1;
944 if (transpose_matrix == true)
945 {
946 val0 = matrix[col];
947 val1 = matrix[nn - 1 - col];
948 }
949 else
950 {
951 val0 = matrix[col * n_columns];
952 val1 = matrix[(nn - col - 1) * n_columns];
953 }
954 if (mid > 0)
955 {
956 res0 = val0 * x[0];
957 res1 = val1 * x[0];
958 res0 -= val1 * x[mm - 1];
959 res1 -= val0 * x[mm - 1];
960 for (int ind = 1; ind < mid; ++ind)
961 {
962 if (transpose_matrix == true)
963 {
964 val0 = matrix[ind * n_columns + col];
965 val1 = matrix[ind * n_columns + nn - 1 - col];
966 }
967 else
968 {
969 val0 = matrix[col * n_columns + ind];
970 val1 = matrix[(nn - col - 1) * n_columns + ind];
971 }
972 res0 += val0 * x[ind];
973 res1 += val1 * x[ind];
974 res0 -= val1 * x[mm - 1 - ind];
975 res1 -= val0 * x[mm - 1 - ind];
976 }
977 }
978 else
979 res0 = res1 = Number();
980 if (mm % 2 == 1)
981 {
982 if (transpose_matrix == true)
983 val0 = matrix[mid * n_columns + col];
984 else
985 val0 = matrix[col * n_columns + mid];
986 const Number tmp = val0 * x[mid];
987 res0 += tmp;
988 res1 -= tmp;
989 }
990 if (add)
991 {
992 out[stride_out * col] += res0;
993 out[stride_out * (nn - 1 - col)] += res1;
994 }
995 else
996 {
997 out[stride_out * col] = res0;
998 out[stride_out * (nn - 1 - col)] = res1;
999 }
1000 }
1001 if (nn % 2 == 1)
1002 {
1003 Number2 val0;
1004 Number res0;
1005 if (transpose_matrix == true)
1006 val0 = matrix[n_cols];
1007 else
1008 val0 = matrix[n_cols * n_columns];
1009 res0 = val0 * (x[0] - x[mm - 1]);
1010 for (int ind = 1; ind < mid; ++ind)
1011 {
1012 if (transpose_matrix == true)
1013 val0 = matrix[ind * n_columns + n_cols];
1014 else
1015 val0 = matrix[n_cols * n_columns + ind];
1016 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
1017 res0 += in1;
1018 }
1019 if (add)
1020 out[stride_out * n_cols] += res0;
1021 else
1022 out[stride_out * n_cols] = res0;
1023 }
1024 }
1025 else
1026 {
1027 // Hessians are almost the same as values, apart from some missing '1'
1028 // entries
1029 for (int col = 0; col < n_cols; ++col)
1030 {
1031 Number2 val0, val1;
1032 Number res0, res1;
1033 if (transpose_matrix == true)
1034 {
1035 val0 = matrix[col];
1036 val1 = matrix[nn - 1 - col];
1037 }
1038 else
1039 {
1040 val0 = matrix[col * n_columns];
1041 val1 = matrix[(col + 1) * n_columns - 1];
1042 }
1043 if (mid > 0)
1044 {
1045 res0 = val0 * x[0];
1046 res1 = val1 * x[0];
1047 res0 += val1 * x[mm - 1];
1048 res1 += val0 * x[mm - 1];
1049 for (int ind = 1; ind < mid; ++ind)
1050 {
1051 if (transpose_matrix == true)
1052 {
1053 val0 = matrix[ind * n_columns + col];
1054 val1 = matrix[ind * n_columns + nn - 1 - col];
1055 }
1056 else
1057 {
1058 val0 = matrix[col * n_columns + ind];
1059 val1 = matrix[(col + 1) * n_columns - 1 - ind];
1060 }
1061 res0 += val0 * x[ind];
1062 res1 += val1 * x[ind];
1063 res0 += val1 * x[mm - 1 - ind];
1064 res1 += val0 * x[mm - 1 - ind];
1065 }
1066 }
1067 else
1068 res0 = res1 = Number();
1069 if (mm % 2 == 1)
1070 {
1071 if (transpose_matrix == true)
1072 val0 = matrix[mid * n_columns + col];
1073 else
1074 val0 = matrix[col * n_columns + mid];
1075 const Number tmp = val0 * x[mid];
1076 res0 += tmp;
1077 res1 += tmp;
1078 }
1079 if (add)
1080 {
1081 out[stride_out * col] += res0;
1082 out[stride_out * (nn - 1 - col)] += res1;
1083 }
1084 else
1085 {
1086 out[stride_out * col] = res0;
1087 out[stride_out * (nn - 1 - col)] = res1;
1088 }
1089 }
1090 if (nn % 2 == 1)
1091 {
1092 Number2 val0;
1093 Number res0;
1094 if (transpose_matrix == true)
1095 val0 = matrix[n_cols];
1096 else
1097 val0 = matrix[n_cols * n_columns];
1098 if (mid > 0)
1099 {
1100 res0 = val0 * (x[0] + x[mm - 1]);
1101 for (int ind = 1; ind < mid; ++ind)
1102 {
1103 if (transpose_matrix == true)
1104 val0 = matrix[ind * n_columns + n_cols];
1105 else
1106 val0 = matrix[n_cols * n_columns + ind];
1107 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
1108 res0 += in1;
1109 }
1110 }
1111 else
1112 res0 = Number();
1113 if (mm % 2 == 1)
1114 {
1115 if (transpose_matrix == true)
1116 val0 = matrix[mid * n_columns + n_cols];
1117 else
1118 val0 = matrix[n_cols * n_columns + mid];
1119 res0 += val0 * x[mid];
1120 }
1121 if (add)
1122 out[stride_out * n_cols] += res0;
1123 else
1124 out[stride_out * n_cols] = res0;
1125 }
1126 }
1127 }
1128
1129
1130
1149 template <EvaluatorVariant variant,
1150 EvaluatorQuantity quantity,
1151 int n_rows_static,
1152 int n_columns_static,
1153 int stride_in_static,
1154 int stride_out_static,
1155 bool transpose_matrix,
1156 bool add,
1157 typename Number,
1158 typename Number2>
1159#ifndef DEBUG
1160 inline DEAL_II_ALWAYS_INLINE
1161#endif
1162 std::enable_if_t<(variant == evaluate_evenodd), void>
1163 apply_matrix_vector_product(const Number2 *DEAL_II_RESTRICT matrix,
1164 const Number *in,
1165 Number *out,
1166 int n_rows_runtime = 0,
1167 int n_columns_runtime = 0,
1168 int stride_in_runtime = 0,
1169 int stride_out_runtime = 0)
1170 {
1171 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
1172 "Negative loop ranges are not allowed!");
1173
1174 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
1175 const int n_columns =
1176 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
1177 const int stride_in =
1178 stride_in_static == 0 ? stride_in_runtime : stride_in_static;
1179 const int stride_out =
1180 stride_out_static == 0 ? stride_out_runtime : stride_out_static;
1181
1182 Assert(n_rows > 0 && n_columns > 0,
1183 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
1184 std::to_string(n_rows) + ", " +
1185 std::to_string(n_columns) + " was passed!"));
1186
1187 const int mm = transpose_matrix ? n_rows : n_columns,
1188 nn = transpose_matrix ? n_columns : n_rows;
1189 const int n_half = nn / 2;
1190 const int m_half = mm / 2;
1191
1192 constexpr int array_length =
1193 (n_rows_static == 0) ?
1194 16 // for non-templated execution
1195 :
1196 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
1197 const int offset = (n_columns + 1) / 2;
1198
1199 Assert(m_half <= array_length, ExcNotImplemented());
1200
1201 std::array<Number, array_length> xp, xm;
1202 for (int i = 0; i < m_half; ++i)
1203 {
1204 if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
1205 {
1206 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1207 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1208 }
1209 else
1210 {
1211 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1212 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1213 }
1214 }
1215 Number xmid = in[stride_in * m_half];
1216 for (int col = 0; col < n_half; ++col)
1217 {
1218 Number r0, r1;
1219 if (m_half > 0)
1220 {
1221 if (transpose_matrix == true)
1222 {
1223 r0 = matrix[col] * xp[0];
1224 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
1225 }
1226 else
1227 {
1228 r0 = matrix[col * offset] * xp[0];
1229 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
1230 }
1231 for (int ind = 1; ind < m_half; ++ind)
1232 {
1233 if (transpose_matrix == true)
1234 {
1235 r0 += matrix[ind * offset + col] * xp[ind];
1236 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
1237 }
1238 else
1239 {
1240 r0 += matrix[col * offset + ind] * xp[ind];
1241 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
1242 }
1243 }
1244 }
1245 else
1246 r0 = r1 = Number();
1247 if (mm % 2 == 1 && transpose_matrix == true)
1248 {
1249 if (quantity == EvaluatorQuantity::gradient)
1250 r1 += matrix[m_half * offset + col] * xmid;
1251 else
1252 r0 += matrix[m_half * offset + col] * xmid;
1253 }
1254 else if (mm % 2 == 1 &&
1255 (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
1256 mm == 3))
1257 r0 += matrix[col * offset + m_half] * xmid;
1258
1259 if (add)
1260 {
1261 out[stride_out * col] += r0 + r1;
1262 if (quantity == EvaluatorQuantity::gradient &&
1263 transpose_matrix == false)
1264 out[stride_out * (nn - 1 - col)] += r1 - r0;
1265 else
1266 out[stride_out * (nn - 1 - col)] += r0 - r1;
1267 }
1268 else
1269 {
1270 out[stride_out * col] = r0 + r1;
1271 if (quantity == EvaluatorQuantity::gradient &&
1272 transpose_matrix == false)
1273 out[stride_out * (nn - 1 - col)] = r1 - r0;
1274 else
1275 out[stride_out * (nn - 1 - col)] = r0 - r1;
1276 }
1277 }
1278 if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
1279 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
1280 {
1281 if (add)
1282 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
1283 else
1284 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
1285 }
1286 else if (transpose_matrix == true && nn % 2 == 1)
1287 {
1288 Number r0;
1289 if (m_half > 0)
1290 {
1291 r0 = matrix[n_half] * xp[0];
1292 for (int ind = 1; ind < m_half; ++ind)
1293 r0 += matrix[ind * offset + n_half] * xp[ind];
1294 }
1295 else
1296 r0 = Number();
1297 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1298 r0 += matrix[m_half * offset + n_half] * xmid;
1299
1300 if (add)
1301 out[stride_out * n_half] += r0;
1302 else
1303 out[stride_out * n_half] = r0;
1304 }
1305 else if (transpose_matrix == false && nn % 2 == 1)
1306 {
1307 Number r0;
1308 if (m_half > 0)
1309 {
1310 if (quantity == EvaluatorQuantity::gradient)
1311 {
1312 r0 = matrix[n_half * offset] * xm[0];
1313 for (int ind = 1; ind < m_half; ++ind)
1314 r0 += matrix[n_half * offset + ind] * xm[ind];
1315 }
1316 else
1317 {
1318 r0 = matrix[n_half * offset] * xp[0];
1319 for (int ind = 1; ind < m_half; ++ind)
1320 r0 += matrix[n_half * offset + ind] * xp[ind];
1321 }
1322 }
1323 else
1324 r0 = Number();
1325
1326 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1327 r0 += matrix[n_half * offset + m_half] * xmid;
1328
1329 if (add)
1330 out[stride_out * n_half] += r0;
1331 else
1332 out[stride_out * n_half] = r0;
1333 }
1334 }
1335
1336
1337
1342 template <EvaluatorVariant variant,
1343 EvaluatorQuantity quantity,
1344 bool transpose_matrix,
1345 bool add,
1346 bool consider_strides,
1347 typename Number,
1348 typename Number2>
1349 std::enable_if_t<(variant == evaluate_evenodd), void>
1350 apply_matrix_vector_product(const Number2 *matrix,
1351 const Number *in,
1352 Number *out,
1353 int n_rows,
1354 int n_columns,
1355 int stride_in,
1356 int stride_out)
1357 {
1358 apply_matrix_vector_product<evaluate_evenodd,
1359 quantity,
1360 0,
1361 0,
1362 consider_strides ? 0 : 1,
1363 consider_strides ? 0 : 1,
1364 transpose_matrix,
1365 add>(
1366 matrix, in, out, n_rows, n_columns, stride_in, stride_out);
1367 }
1368
1369
1370
1386 template <EvaluatorVariant variant,
1387 EvaluatorQuantity quantity,
1388 int n_rows,
1389 int n_columns,
1390 int stride_in,
1391 int stride_out,
1392 bool transpose_matrix,
1393 bool add,
1394 typename Number,
1395 typename Number2>
1396 std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
1397 apply_matrix_vector_product(const Number2 *matrix,
1398 const Number *in,
1399 Number *out)
1400 {
1401 static_assert(n_rows > 0 && n_columns > 0,
1402 "Specialization requires n_rows, n_columns > 0");
1403
1404 constexpr bool evaluate_antisymmetric =
1405 (quantity == EvaluatorQuantity::gradient);
1406
1407 constexpr int mm = transpose_matrix ? n_rows : n_columns,
1408 nn = transpose_matrix ? n_columns : n_rows;
1409 constexpr int n_half = nn / 2;
1410 constexpr int m_half = mm / 2;
1411
1412 if (transpose_matrix)
1413 {
1414 std::array<Number, mm> x;
1415 for (unsigned int i = 0; i < mm; ++i)
1416 x[i] = in[stride_in * i];
1417 for (unsigned int col = 0; col < n_half; ++col)
1418 {
1419 Number r0, r1;
1420 if (m_half > 0)
1421 {
1422 r0 = matrix[col] * x[0];
1423 r1 = matrix[col + n_columns] * x[1];
1424 for (unsigned int ind = 1; ind < m_half; ++ind)
1425 {
1426 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
1427 r1 +=
1428 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
1429 }
1430 }
1431 else
1432 r0 = r1 = Number();
1433 if (mm % 2 == 1)
1434 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
1435 if (add)
1436 {
1437 out[stride_out * col] += r0 + r1;
1438 if (evaluate_antisymmetric)
1439 out[stride_out * (nn - 1 - col)] += r1 - r0;
1440 else
1441 out[stride_out * (nn - 1 - col)] += r0 - r1;
1442 }
1443 else
1444 {
1445 out[stride_out * col] = r0 + r1;
1446 if (evaluate_antisymmetric)
1447 out[stride_out * (nn - 1 - col)] = r1 - r0;
1448 else
1449 out[stride_out * (nn - 1 - col)] = r0 - r1;
1450 }
1451 }
1452 if (nn % 2 == 1)
1453 {
1454 Number r0;
1455 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1456 if (m_half > 0)
1457 {
1458 r0 = matrix[n_half + shift * n_columns] * x[shift];
1459 for (unsigned int ind = 1; ind < m_half; ++ind)
1460 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1461 x[2 * ind + shift];
1462 }
1463 else
1464 r0 = 0;
1465 if (!evaluate_antisymmetric && mm % 2 == 1)
1466 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1467 if (add)
1468 out[stride_out * n_half] += r0;
1469 else
1470 out[stride_out * n_half] = r0;
1471 }
1472 }
1473 else
1474 {
1475 std::array<Number, m_half + 1> xp, xm;
1476 for (int i = 0; i < m_half; ++i)
1477 if (!evaluate_antisymmetric)
1478 {
1479 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1480 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1481 }
1482 else
1483 {
1484 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1485 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1486 }
1487 if (mm % 2 == 1)
1488 xp[m_half] = in[stride_in * m_half];
1489 for (unsigned int col = 0; col < n_half; ++col)
1490 {
1491 Number r0, r1;
1492 if (m_half > 0)
1493 {
1494 r0 = matrix[2 * col * n_columns] * xp[0];
1495 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1496 for (unsigned int ind = 1; ind < m_half; ++ind)
1497 {
1498 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1499 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1500 }
1501 }
1502 else
1503 r0 = r1 = Number();
1504 if (mm % 2 == 1)
1505 {
1506 if (evaluate_antisymmetric)
1507 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1508 else
1509 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1510 }
1511 if (add)
1512 {
1513 out[stride_out * (2 * col)] += r0;
1514 out[stride_out * (2 * col + 1)] += r1;
1515 }
1516 else
1517 {
1518 out[stride_out * (2 * col)] = r0;
1519 out[stride_out * (2 * col + 1)] = r1;
1520 }
1521 }
1522 if (nn % 2 == 1)
1523 {
1524 Number r0;
1525 if (m_half > 0)
1526 {
1527 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1528 for (unsigned int ind = 1; ind < m_half; ++ind)
1529 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1530 }
1531 else
1532 r0 = Number();
1533 if (mm % 2 == 1 && !evaluate_antisymmetric)
1534 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1535 if (add)
1536 out[stride_out * (nn - 1)] += r0;
1537 else
1538 out[stride_out * (nn - 1)] = r0;
1539 }
1540 }
1541 }
1542
1543
1544
1567 template <EvaluatorVariant variant,
1568 int dim,
1569 int n_rows,
1570 int n_columns,
1571 typename Number,
1572 typename Number2 = Number>
1573 struct EvaluatorTensorProduct
1574 {
1575 static constexpr unsigned int n_rows_of_product =
1576 Utilities::pow(n_rows, dim);
1577 static constexpr unsigned int n_columns_of_product =
1578 Utilities::pow(n_columns, dim);
1579
1584 EvaluatorTensorProduct()
1585 : shape_values(nullptr)
1586 , shape_gradients(nullptr)
1587 , shape_hessians(nullptr)
1588 {}
1589
1593 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1594 const AlignedVector<Number2> &shape_gradients,
1595 const AlignedVector<Number2> &shape_hessians,
1596 const unsigned int = 0,
1597 const unsigned int = 0)
1598 : shape_values(shape_values.begin())
1599 , shape_gradients(shape_gradients.begin())
1600 , shape_hessians(shape_hessians.begin())
1601 {
1602 if (variant == evaluate_evenodd)
1603 {
1604 if (!shape_values.empty())
1605 AssertDimension(shape_values.size(),
1606 n_rows * ((n_columns + 1) / 2));
1607 if (!shape_gradients.empty())
1608 AssertDimension(shape_gradients.size(),
1609 n_rows * ((n_columns + 1) / 2));
1610 if (!shape_hessians.empty())
1611 AssertDimension(shape_hessians.size(),
1612 n_rows * ((n_columns + 1) / 2));
1613 }
1614 else
1615 {
1616 Assert(shape_values.empty() ||
1617 shape_values.size() == n_rows * n_columns,
1618 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1619 Assert(shape_gradients.empty() ||
1620 shape_gradients.size() == n_rows * n_columns,
1621 ExcDimensionMismatch(shape_gradients.size(),
1622 n_rows * n_columns));
1623 Assert(shape_hessians.empty() ||
1624 shape_hessians.size() == n_rows * n_columns,
1625 ExcDimensionMismatch(shape_hessians.size(),
1626 n_rows * n_columns));
1627 }
1628 }
1629
1633 EvaluatorTensorProduct(const Number2 *shape_values,
1634 const Number2 *shape_gradients,
1635 const Number2 *shape_hessians,
1636 const unsigned int dummy1 = 0,
1637 const unsigned int dummy2 = 0)
1638 : shape_values(shape_values)
1639 , shape_gradients(shape_gradients)
1640 , shape_hessians(shape_hessians)
1641 {
1642 (void)dummy1;
1643 (void)dummy2;
1644 }
1645
1671 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1672 void
1673 values(const Number in[], Number out[]) const
1674 {
1675 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1676 apply<direction, contract_over_rows, add, false, value_type, stride>(
1677 shape_values, in, out);
1678 }
1679
1685 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1686 void
1687 gradients(const Number in[], Number out[]) const
1688 {
1689 constexpr EvaluatorQuantity gradient_type =
1690 (variant == evaluate_general ? EvaluatorQuantity::value :
1691 EvaluatorQuantity::gradient);
1692 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1693 shape_gradients, in, out);
1694 }
1695
1701 template <int direction, bool contract_over_rows, bool add>
1702 void
1703 hessians(const Number in[], Number out[]) const
1704 {
1705 constexpr EvaluatorQuantity hessian_type =
1706 (((variant == evaluate_general) |
1707 (variant == evaluate_symmetric_hierarchical)) ?
1708 EvaluatorQuantity::value :
1709 EvaluatorQuantity::hessian);
1710 apply<direction, contract_over_rows, add, false, hessian_type>(
1711 shape_hessians, in, out);
1712 }
1713
1721 template <int direction, bool contract_over_rows, bool add>
1722 void
1723 values_one_line(const Number in[], Number out[]) const
1724 {
1725 Assert(shape_values != nullptr, ExcNotInitialized());
1726 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1727 shape_values, in, out);
1728 }
1729
1737 template <int direction, bool contract_over_rows, bool add>
1738 void
1739 gradients_one_line(const Number in[], Number out[]) const
1740 {
1741 Assert(shape_gradients != nullptr, ExcNotInitialized());
1742 constexpr EvaluatorQuantity gradient_type =
1743 (variant == evaluate_general ? EvaluatorQuantity::value :
1744 EvaluatorQuantity::gradient);
1745 apply<direction, contract_over_rows, add, true, gradient_type>(
1746 shape_gradients, in, out);
1747 }
1748
1756 template <int direction, bool contract_over_rows, bool add>
1757 void
1758 hessians_one_line(const Number in[], Number out[]) const
1759 {
1760 Assert(shape_hessians != nullptr, ExcNotInitialized());
1761 constexpr EvaluatorQuantity hessian_type =
1762 (((variant == evaluate_general) |
1763 (variant == evaluate_symmetric_hierarchical)) ?
1764 EvaluatorQuantity::value :
1765 EvaluatorQuantity::hessian);
1766 apply<direction, contract_over_rows, add, true, hessian_type>(
1767 shape_hessians, in, out);
1768 }
1769
1806 template <int direction,
1807 bool contract_over_rows,
1808 bool add,
1809 bool one_line = false,
1810 EvaluatorQuantity quantity = EvaluatorQuantity::value,
1811 int stride = 1>
1812 static void
1813 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1814 const Number *in,
1815 Number *out);
1816
1817 private:
1818 const Number2 *shape_values;
1819 const Number2 *shape_gradients;
1820 const Number2 *shape_hessians;
1821 };
1822
1823
1824
1825 template <EvaluatorVariant variant,
1826 int dim,
1827 int n_rows,
1828 int n_columns,
1829 typename Number,
1830 typename Number2>
1831 template <int direction,
1832 bool contract_over_rows,
1833 bool add,
1834 bool one_line,
1835 EvaluatorQuantity quantity,
1836 int stride>
1837 inline void
1838 EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
1839 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1840 const Number *in,
1841 Number *out)
1842 {
1843 static_assert(one_line == false || direction == dim - 1,
1844 "Single-line evaluation only works for direction=dim-1.");
1845 Assert(shape_data != nullptr,
1846 ExcMessage(
1847 "The given array shape_data must not be the null pointer!"));
1848 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1849 in != out,
1850 ExcMessage("In-place operation only supported for "
1851 "n_rows==n_columns or single-line interpolation"));
1852 AssertIndexRange(direction, dim);
1853 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1854 nn = contract_over_rows ? n_columns : n_rows;
1855
1856 constexpr int stride_operation = Utilities::pow(n_columns, direction);
1857 constexpr int n_blocks1 = one_line ? 1 : stride_operation;
1858 constexpr int n_blocks2 =
1859 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1860
1861 constexpr int stride_in = !contract_over_rows ? stride : 1;
1862 constexpr int stride_out = contract_over_rows ? stride : 1;
1863 for (int i2 = 0; i2 < n_blocks2; ++i2)
1864 {
1865 for (int i1 = 0; i1 < n_blocks1; ++i1)
1866 {
1867 apply_matrix_vector_product<variant,
1868 quantity,
1869 n_rows,
1870 n_columns,
1871 stride_operation * stride_in,
1872 stride_operation * stride_out,
1873 contract_over_rows,
1874 add>(shape_data, in, out);
1875
1876 if (one_line == false)
1877 {
1878 in += stride_in;
1879 out += stride_out;
1880 }
1881 }
1882 if (one_line == false)
1883 {
1884 in += stride_operation * (mm - 1) * stride_in;
1885 out += stride_operation * (nn - 1) * stride_out;
1886 }
1887 }
1888 }
1889
1890
1891
1905 template <EvaluatorVariant variant,
1906 int dim,
1907 typename Number,
1908 typename Number2>
1909 struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
1910 {
1911 static constexpr unsigned int n_rows_of_product =
1912 numbers::invalid_unsigned_int;
1913 static constexpr unsigned int n_columns_of_product =
1914 numbers::invalid_unsigned_int;
1915
1920 EvaluatorTensorProduct()
1921 : shape_values(nullptr)
1922 , shape_gradients(nullptr)
1923 , shape_hessians(nullptr)
1924 , n_rows(numbers::invalid_unsigned_int)
1925 , n_columns(numbers::invalid_unsigned_int)
1926 {}
1927
1931 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1932 const AlignedVector<Number2> &shape_gradients,
1933 const AlignedVector<Number2> &shape_hessians,
1934 const unsigned int n_rows = 0,
1935 const unsigned int n_columns = 0)
1936 : shape_values(shape_values.begin())
1937 , shape_gradients(shape_gradients.begin())
1938 , shape_hessians(shape_hessians.begin())
1939 , n_rows(n_rows)
1940 , n_columns(n_columns)
1941 {
1942 if (variant == evaluate_evenodd)
1943 {
1944 if (!shape_values.empty())
1945 AssertDimension(shape_values.size(),
1946 n_rows * ((n_columns + 1) / 2));
1947 if (!shape_gradients.empty())
1948 AssertDimension(shape_gradients.size(),
1949 n_rows * ((n_columns + 1) / 2));
1950 if (!shape_hessians.empty())
1951 AssertDimension(shape_hessians.size(),
1952 n_rows * ((n_columns + 1) / 2));
1953 }
1954 else
1955 {
1956 Assert(shape_values.empty() ||
1957 shape_values.size() == n_rows * n_columns,
1958 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1959 Assert(shape_gradients.empty() ||
1960 shape_gradients.size() == n_rows * n_columns,
1961 ExcDimensionMismatch(shape_gradients.size(),
1962 n_rows * n_columns));
1963 Assert(shape_hessians.empty() ||
1964 shape_hessians.size() == n_rows * n_columns,
1965 ExcDimensionMismatch(shape_hessians.size(),
1966 n_rows * n_columns));
1967 }
1968 }
1969
1973 EvaluatorTensorProduct(const Number2 *shape_values,
1974 const Number2 *shape_gradients,
1975 const Number2 *shape_hessians,
1976 const unsigned int n_rows = 0,
1977 const unsigned int n_columns = 0)
1978 : shape_values(shape_values)
1979 , shape_gradients(shape_gradients)
1980 , shape_hessians(shape_hessians)
1981 , n_rows(n_rows)
1982 , n_columns(n_columns)
1983 {}
1984
1985 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1986 void
1987 values(const Number *in, Number *out) const
1988 {
1989 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1990 apply<direction, contract_over_rows, add, false, value_type, stride>(
1991 shape_values, in, out);
1992 }
1993
1994 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1995 void
1996 gradients(const Number *in, Number *out) const
1997 {
1998 constexpr EvaluatorQuantity gradient_type =
1999 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2000 EvaluatorQuantity::gradient);
2001 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
2002 shape_gradients, in, out);
2003 }
2004
2005 template <int direction, bool contract_over_rows, bool add>
2006 void
2007 hessians(const Number *in, Number *out) const
2008 {
2009 constexpr EvaluatorQuantity hessian_type =
2010 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2011 EvaluatorQuantity::hessian);
2012 apply<direction, contract_over_rows, add, false, hessian_type>(
2013 shape_hessians, in, out);
2014 }
2015
2016 template <int direction, bool contract_over_rows, bool add>
2017 void
2018 values_one_line(const Number in[], Number out[]) const
2019 {
2020 Assert(shape_values != nullptr, ExcNotInitialized());
2021 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
2022 shape_values, in, out);
2023 }
2024
2025 template <int direction, bool contract_over_rows, bool add>
2026 void
2027 gradients_one_line(const Number in[], Number out[]) const
2028 {
2029 Assert(shape_gradients != nullptr, ExcNotInitialized());
2030 constexpr EvaluatorQuantity gradient_type =
2031 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2032 EvaluatorQuantity::gradient);
2033 apply<direction, contract_over_rows, add, true, gradient_type>(
2034 shape_gradients, in, out);
2035 }
2036
2037 template <int direction, bool contract_over_rows, bool add>
2038 void
2039 hessians_one_line(const Number in[], Number out[]) const
2040 {
2041 Assert(shape_hessians != nullptr, ExcNotInitialized());
2042 constexpr EvaluatorQuantity hessian_type =
2043 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2044 EvaluatorQuantity::hessian);
2045 apply<direction, contract_over_rows, add, true, hessian_type>(
2046 shape_hessians, in, out);
2047 }
2048
2049 template <int direction,
2050 bool contract_over_rows,
2051 bool add,
2052 bool one_line = false,
2053 EvaluatorQuantity quantity = EvaluatorQuantity::value,
2054 int stride = 1>
2055 void
2056 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2057 const Number *in,
2058 Number *out) const;
2059
2060 const Number2 *shape_values;
2061 const Number2 *shape_gradients;
2062 const Number2 *shape_hessians;
2063 const unsigned int n_rows;
2064 const unsigned int n_columns;
2065 };
2066
2067
2068
2069 template <EvaluatorVariant variant,
2070 int dim,
2071 typename Number,
2072 typename Number2>
2073 template <int direction,
2074 bool contract_over_rows,
2075 bool add,
2076 bool one_line,
2077 EvaluatorQuantity quantity,
2078 int stride>
2079 inline void
2080 EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply(
2081 const Number2 *DEAL_II_RESTRICT shape_data,
2082 const Number *in,
2083 Number *out) const
2084 {
2085 static_assert(one_line == false || direction == dim - 1,
2086 "Single-line evaluation only works for direction=dim-1.");
2087 Assert(shape_data != nullptr,
2088 ExcMessage(
2089 "The given array shape_data must not be the null pointer!"));
2090 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2091 in != out,
2092 ExcMessage("In-place operation only supported for "
2093 "n_rows==n_columns or single-line interpolation"));
2094 AssertIndexRange(direction, dim);
2095 const int mm = contract_over_rows ? n_rows : n_columns,
2096 nn = contract_over_rows ? n_columns : n_rows;
2097
2098 const int stride_operation =
2099 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
2100 const int n_blocks1 = one_line ? 1 : stride_operation;
2101 const int n_blocks2 = direction >= dim - 1 ?
2102 1 :
2103 Utilities::fixed_power<dim - direction - 1>(n_rows);
2104 Assert(n_rows <= 128, ExcNotImplemented());
2105
2106 constexpr int stride_in = !contract_over_rows ? stride : 1;
2107 constexpr int stride_out = contract_over_rows ? stride : 1;
2108 for (int i2 = 0; i2 < n_blocks2; ++i2)
2109 {
2110 for (int i1 = 0; i1 < n_blocks1; ++i1)
2111 {
2112 // the empty template case can only run the general evaluator or
2113 // evenodd
2114 constexpr EvaluatorVariant restricted_variant =
2115 variant == evaluate_evenodd ? evaluate_evenodd : evaluate_general;
2116 apply_matrix_vector_product<restricted_variant,
2117 quantity,
2118 contract_over_rows,
2119 add,
2120 (direction != 0 || stride != 1)>(
2121 shape_data,
2122 in,
2123 out,
2124 n_rows,
2125 n_columns,
2126 stride_operation * stride_in,
2127 stride_operation * stride_out);
2128
2129 if (one_line == false)
2130 {
2131 in += stride_in;
2132 out += stride_out;
2133 }
2134 }
2135 if (one_line == false)
2136 {
2137 in += stride_operation * (mm - 1) * stride_in;
2138 out += stride_operation * (nn - 1) * stride_out;
2139 }
2140 }
2141 }
2142
2143
2144
2145 template <int dim,
2146 int fe_degree,
2147 int n_q_points_1d,
2148 bool contract_over_rows,
2149 bool symmetric_evaluate = true>
2150 struct EvaluatorTensorProductAnisotropic
2151 {
2152 template <int direction,
2153 int stride = 1,
2154 typename Number = double,
2155 typename Number2 = double>
2156 static void
2157 normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
2158 const Number *in,
2159 Number *out,
2160 const bool add_into_result = false,
2161 const int subface_index_1d = 0)
2162 {
2163 AssertIndexRange(direction, dim);
2164 AssertDimension(fe_degree, data.fe_degree);
2165 AssertDimension(n_q_points_1d, data.n_q_points_1d);
2166 constexpr int n_rows = fe_degree + 1;
2167 constexpr int n_columns = n_q_points_1d;
2168 constexpr int mm = contract_over_rows ? n_rows : n_columns;
2169 constexpr int nn = contract_over_rows ? n_columns : n_rows;
2170 const Number2 *shape_data =
2171 symmetric_evaluate ?
2172 data.shape_values_eo.data() :
2173 data.values_within_subface[subface_index_1d].data();
2174 Assert(shape_data != nullptr, ExcNotInitialized());
2175 Assert(contract_over_rows == false || !add_into_result,
2176 ExcMessage("Cannot add into result if contract_over_rows = true"));
2177
2178 constexpr int n_blocks1 = Utilities::pow(fe_degree, direction);
2179 constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1);
2180 constexpr int stride_in = contract_over_rows ? 1 : stride;
2181 constexpr int stride_out = contract_over_rows ? stride : 1;
2182 constexpr EvaluatorVariant variant =
2183 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
2184
2185 for (int i2 = 0; i2 < n_blocks2; ++i2)
2186 {
2187 for (int i1 = 0; i1 < n_blocks1; ++i1)
2188 {
2189 if (contract_over_rows == false && add_into_result)
2190 apply_matrix_vector_product<variant,
2191 EvaluatorQuantity::value,
2192 n_rows,
2193 n_columns,
2194 n_blocks1 * stride_in,
2195 n_blocks1 * stride_out,
2196 contract_over_rows,
2197 true>(shape_data, in, out);
2198 else
2199 apply_matrix_vector_product<variant,
2200 EvaluatorQuantity::value,
2201 n_rows,
2202 n_columns,
2203 n_blocks1 * stride_in,
2204 n_blocks1 * stride_out,
2205 contract_over_rows,
2206 false>(shape_data, in, out);
2207
2208 in += stride_in;
2209 out += stride_out;
2210 }
2211 in += n_blocks1 * (mm - 1) * stride_in;
2212 out += n_blocks1 * (nn - 1) * stride_out;
2213 }
2214 }
2215
2216 template <int direction,
2217 int normal_direction,
2218 int stride = 1,
2219 typename Number = double,
2220 typename Number2 = double>
2221 static void
2222 tangential(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
2223 const Number *in,
2224 Number *out,
2225 const int subface_index_1d = 0)
2226 {
2227 AssertIndexRange(direction, dim);
2228 AssertDimension(fe_degree - 1, data.fe_degree);
2229 AssertDimension(n_q_points_1d, data.n_q_points_1d);
2230 static_assert(direction != normal_direction,
2231 "Cannot interpolate tangentially in normal direction");
2232
2233 constexpr int n_rows = std::max(fe_degree, 0);
2234 constexpr int n_columns = n_q_points_1d;
2235 const Number2 *shape_data =
2236 symmetric_evaluate ?
2237 data.shape_values_eo.data() :
2238 data.values_within_subface[subface_index_1d].data();
2239 Assert(shape_data != nullptr, ExcNotInitialized());
2240
2241 constexpr int n_blocks1 =
2242 (direction > normal_direction) ?
2243 Utilities::pow(n_q_points_1d, direction) :
2244 (direction > 0 ?
2245 (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
2246 1);
2247 constexpr int n_blocks2 =
2248 (direction > normal_direction) ?
2249 Utilities::pow(fe_degree, dim - 1 - direction) :
2250 ((direction + 1 < dim) ?
2251 (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
2252 1);
2253
2254 constexpr EvaluatorVariant variant =
2255 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
2256
2257 // Since we may perform an in-place interpolation, we must run the step
2258 // expanding the size of the basis backward ('contract_over_rows' aka
2259 // 'evaluate' case), so shift the pointers and decrement during the loop
2260 if (contract_over_rows)
2261 {
2262 in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
2263 out +=
2264 stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
2265 for (int i2 = 0; i2 < n_blocks2; ++i2)
2266 {
2267 for (int i1 = 0; i1 < n_blocks1; ++i1)
2268 {
2269 apply_matrix_vector_product<variant,
2270 EvaluatorQuantity::value,
2271 n_rows,
2272 n_columns,
2273 n_blocks1,
2274 n_blocks1 * stride,
2275 true,
2276 false>(shape_data, in, out);
2277
2278 --in;
2279 out -= stride;
2280 }
2281 in -= n_blocks1 * (n_rows - 1);
2282 out -= n_blocks1 * (n_columns - 1) * stride;
2283 }
2284 }
2285 else
2286 {
2287 for (int i2 = 0; i2 < n_blocks2; ++i2)
2288 {
2289 for (int i1 = 0; i1 < n_blocks1; ++i1)
2290 {
2291 apply_matrix_vector_product<variant,
2292 EvaluatorQuantity::value,
2293 n_rows,
2294 n_columns,
2295 n_blocks1 * stride,
2296 n_blocks1,
2297 false,
2298 false>(shape_data, in, out);
2299
2300 in += stride;
2301 ++out;
2302 }
2303 in += n_blocks1 * (n_columns - 1) * stride;
2304 out += n_blocks1 * (n_rows - 1);
2305 }
2306 }
2307 }
2308 };
2309
2310
2311
2357 template <int n_rows_template,
2358 int stride_template,
2359 bool contract_onto_face,
2360 bool add,
2361 int max_derivative,
2362 typename Number,
2363 typename Number2>
2364 inline std::enable_if_t<contract_onto_face, void>
2365 interpolate_to_face(const Number2 *shape_values,
2366 const std::array<int, 2> &n_blocks,
2367 const std::array<int, 2> &steps,
2368 const Number *input,
2369 Number *DEAL_II_RESTRICT output,
2370 const int n_rows_runtime = 0,
2371 const int stride_runtime = 1)
2372 {
2373 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2374 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2375
2376 Number *output1 = output + n_blocks[0] * n_blocks[1];
2377 Number *output2 = output1 + n_blocks[0] * n_blocks[1];
2378 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2379 {
2380 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2381 {
2382 Number res0 = shape_values[0] * input[0];
2383 Number res1, res2;
2384 if (max_derivative > 0)
2385 res1 = shape_values[n_rows] * input[0];
2386 if (max_derivative > 1)
2387 res2 = shape_values[2 * n_rows] * input[0];
2388 for (int ind = 1; ind < n_rows; ++ind)
2389 {
2390 res0 += shape_values[ind] * input[stride * ind];
2391 if (max_derivative > 0)
2392 res1 += shape_values[ind + n_rows] * input[stride * ind];
2393 if (max_derivative > 1)
2394 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
2395 }
2396 if (add)
2397 {
2398 output[i1] += res0;
2399 if (max_derivative > 0)
2400 output1[i1] += res1;
2401 if (max_derivative > 1)
2402 output2[i2] += res2;
2403 }
2404 else
2405 {
2406 output[i1] = res0;
2407 if (max_derivative > 0)
2408 output1[i1] = res1;
2409 if (max_derivative > 1)
2410 output2[i1] = res2;
2411 }
2412 input += steps[0];
2413 }
2414 output += n_blocks[0];
2415 if (max_derivative > 0)
2416 output1 += n_blocks[0];
2417 if (max_derivative > 1)
2418 output2 += n_blocks[0];
2419 input += steps[1];
2420 }
2421 }
2422
2423
2424
2432 constexpr bool
2433 use_collocation_evaluation(const unsigned int fe_degree,
2434 const unsigned int n_q_points_1d)
2435 {
2436 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2437 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2438 }
2439
2440
2441
2447 template <int n_rows_template,
2448 int stride_template,
2449 bool contract_onto_face,
2450 bool add,
2451 int max_derivative,
2452 typename Number,
2453 typename Number2>
2454 inline std::enable_if_t<!contract_onto_face, void>
2455 interpolate_to_face(const Number2 *shape_values,
2456 const std::array<int, 2> &n_blocks,
2457 const std::array<int, 2> &steps,
2458 const Number *input,
2459 Number *DEAL_II_RESTRICT output,
2460 const int n_rows_runtime = 0,
2461 const int stride_runtime = 1)
2462 {
2463 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2464 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2465
2466 const Number *input1 = input + n_blocks[0] * n_blocks[1];
2467 const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
2468 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2469 {
2470 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2471 {
2472 const Number in = input[i1];
2473 Number in1, in2;
2474 if (max_derivative > 0)
2475 in1 = input1[i1];
2476 if (max_derivative > 1)
2477 in2 = input2[i1];
2478 for (int col = 0; col < n_rows; ++col)
2479 {
2480 Number result =
2481 add ? (output[col * stride] + shape_values[col] * in) :
2482 (shape_values[col] * in);
2483 if (max_derivative > 0)
2484 result += shape_values[col + n_rows] * in1;
2485 if (max_derivative > 1)
2486 result += shape_values[col + 2 * n_rows] * in2;
2487
2488 output[col * stride] = result;
2489 }
2490 output += steps[0];
2491 }
2492 input += n_blocks[0];
2493 if (max_derivative > 0)
2494 input1 += n_blocks[0];
2495 if (max_derivative > 1)
2496 input2 += n_blocks[0];
2497 output += steps[1];
2498 }
2499 }
2500
2501 template <int dim, int n_points_1d_template, typename Number>
2502 inline void
2503 weight_fe_q_dofs_by_entity(const Number *weights,
2504 const unsigned int n_components,
2505 const int n_points_1d_non_template,
2506 Number *data)
2507 {
2508 const int n_points_1d = n_points_1d_template != -1 ?
2509 n_points_1d_template :
2510 n_points_1d_non_template;
2511
2512 Assert(n_points_1d > 0, ExcNotImplemented());
2513 Assert(n_points_1d < 100, ExcNotImplemented());
2514
2515 unsigned int compressed_index[100];
2516 compressed_index[0] = 0;
2517 for (int i = 1; i < n_points_1d - 1; ++i)
2518 compressed_index[i] = 1;
2519 compressed_index[n_points_1d - 1] = 2;
2520
2521 for (unsigned int c = 0; c < n_components; ++c)
2522 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2523 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2524 {
2525 const unsigned int shift =
2526 9 * compressed_index[k] + 3 * compressed_index[j];
2527 data[0] *= weights[shift];
2528 // loop bound as int avoids compiler warnings in case n_points_1d
2529 // == 1 (polynomial degree 0)
2530 const Number weight = weights[shift + 1];
2531 for (int i = 1; i < n_points_1d - 1; ++i)
2532 data[i] *= weight;
2533 data[n_points_1d - 1] *= weights[shift + 2];
2534 data += n_points_1d;
2535 }
2536 }
2537
2538
2539 template <int dim, int n_points_1d_template, typename Number>
2540 inline void
2541 weight_fe_q_dofs_by_entity_shifted(const Number *weights,
2542 const unsigned int n_components,
2543 const int n_points_1d_non_template,
2544 Number *data)
2545 {
2546 const int n_points_1d = n_points_1d_template != -1 ?
2547 n_points_1d_template :
2548 n_points_1d_non_template;
2549
2550 Assert((n_points_1d % 2) == 1,
2551 ExcMessage("The function can only with add number of points"));
2552 Assert(n_points_1d > 0, ExcNotImplemented());
2553 Assert(n_points_1d < 100, ExcNotImplemented());
2554
2555 const unsigned int n_inside_1d = n_points_1d / 2;
2556
2557 unsigned int compressed_index[100];
2558
2559 unsigned int c = 0;
2560 for (int i = 0; i < n_inside_1d; ++i)
2561 compressed_index[c++] = 0;
2562 compressed_index[c++] = 1;
2563 for (int i = 0; i < n_inside_1d; ++i)
2564 compressed_index[c++] = 2;
2565
2566 for (unsigned int c = 0; c < n_components; ++c)
2567 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2568 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2569 {
2570 const unsigned int shift =
2571 9 * compressed_index[k] + 3 * compressed_index[j];
2572
2573 unsigned int c = 0;
2574 const Number weight1 = weights[shift];
2575 for (int i = 0; i < n_inside_1d; ++i)
2576 data[c++] *= weight1;
2577 data[c++] *= weights[shift + 1];
2578 const Number weight2 = weights[shift + 2];
2579 for (int i = 0; i < n_inside_1d; ++i)
2580 data[c++] *= weight2;
2581 data += n_points_1d;
2582 }
2583 }
2584
2585
2586 template <int dim, int n_points_1d_template, typename Number>
2587 inline bool
2588 compute_weights_fe_q_dofs_by_entity(const Number *data,
2589 const unsigned int n_components,
2590 const int n_points_1d_non_template,
2591 Number *weights)
2592 {
2593 const int n_points_1d = n_points_1d_template != -1 ?
2594 n_points_1d_template :
2595 n_points_1d_non_template;
2596
2597 Assert(n_points_1d > 0, ExcNotImplemented());
2598 Assert(n_points_1d < 100, ExcNotImplemented());
2599
2600 unsigned int compressed_index[100];
2601 compressed_index[0] = 0;
2602 for (int i = 1; i < n_points_1d - 1; ++i)
2603 compressed_index[i] = 1;
2604 compressed_index[n_points_1d - 1] = 2;
2605
2606 // Insert the number data into a storage position for weight,
2607 // ensuring that the array has either not been touched before
2608 // or the previous content is the same. In case the previous
2609 // content has a different value, we exit this function and
2610 // signal to outer functions that the compression was not possible.
2611 const auto check_and_set = [](Number &weight, const Number &data) {
2612 if (weight == Number(-1.0) || weight == data)
2613 {
2614 weight = data;
2615 return true; // success for the entry
2616 }
2617
2618 return false; // failure for the entry
2619 };
2620
2621 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2622 weights[c] = Number(-1.0);
2623
2624 for (unsigned int c = 0; c < n_components; ++c)
2625 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2626 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2627 ++j, data += n_points_1d)
2628 {
2629 const unsigned int shift =
2630 9 * compressed_index[k] + 3 * compressed_index[j];
2631
2632 if (!check_and_set(weights[shift], data[0]))
2633 return false; // failure
2634
2635 for (int i = 1; i < n_points_1d - 1; ++i)
2636 if (!check_and_set(weights[shift + 1], data[i]))
2637 return false; // failure
2638
2639 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
2640 return false; // failure
2641 }
2642
2643 return true; // success
2644 }
2645
2646
2647 template <int dim, int n_points_1d_template, typename Number>
2648 inline bool
2649 compute_weights_fe_q_dofs_by_entity_shifted(
2650 const Number *data,
2651 const unsigned int n_components,
2652 const int n_points_1d_non_template,
2653 Number *weights)
2654 {
2655 const int n_points_1d = n_points_1d_template != -1 ?
2656 n_points_1d_template :
2657 n_points_1d_non_template;
2658
2659 Assert((n_points_1d % 2) == 1,
2660 ExcMessage("The function can only with add number of points"));
2661 Assert(n_points_1d > 0, ExcNotImplemented());
2662 Assert(n_points_1d < 100, ExcNotImplemented());
2663
2664 const unsigned int n_inside_1d = n_points_1d / 2;
2665
2666 unsigned int compressed_index[100];
2667
2668 unsigned int c = 0;
2669 for (int i = 0; i < n_inside_1d; ++i)
2670 compressed_index[c++] = 0;
2671 compressed_index[c++] = 1;
2672 for (int i = 0; i < n_inside_1d; ++i)
2673 compressed_index[c++] = 2;
2674
2675 // Insert the number data into a storage position for weight,
2676 // ensuring that the array has either not been touched before
2677 // or the previous content is the same. In case the previous
2678 // content has a different value, we exit this function and
2679 // signal to outer functions that the compression was not possible.
2680 const auto check_and_set = [](Number &weight, const Number &data) {
2681 if (weight == Number(-1.0) || weight == data)
2682 {
2683 weight = data;
2684 return true; // success for the entry
2685 }
2686
2687 return false; // failure for the entry
2688 };
2689
2690 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2691 weights[c] = Number(-1.0);
2692
2693 for (unsigned int comp = 0; comp < n_components; ++comp)
2694 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2695 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2696 ++j, data += n_points_1d)
2697 {
2698 const unsigned int shift =
2699 9 * compressed_index[k] + 3 * compressed_index[j];
2700
2701 unsigned int c = 0;
2702
2703 for (int i = 0; i < n_inside_1d; ++i)
2704 if (!check_and_set(weights[shift], data[c++]))
2705 return false; // failure
2706
2707 if (!check_and_set(weights[shift + 1], data[c++]))
2708 return false; // failure
2709
2710 for (int i = 0; i < n_inside_1d; ++i)
2711 if (!check_and_set(weights[shift + 2], data[c++]))
2712 return false; // failure
2713 }
2714
2715 return true; // success
2716 }
2717
2718
2719} // end of namespace internal
2720
2721
2722DEAL_II_NAMESPACE_CLOSE
2723
2724#endif
aligned_vector.h
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition config.h:110
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
data
std::vector< index_type > data
Definition mpi.cc:735
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2433
internal::weight_fe_q_dofs_by_entity
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
Definition tensor_product_kernels.h:2503
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:118
internal::weight_fe_q_dofs_by_entity_shifted
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
Definition tensor_product_kernels.h:2541
internal::interpolate_to_face
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
Definition tensor_product_kernels.h:2365
internal::compute_weights_fe_q_dofs_by_entity
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
Definition tensor_product_kernels.h:2588
internal::compute_weights_fe_q_dofs_by_entity_shifted
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
Definition tensor_product_kernels.h:2649
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
shape_info.h
internal::EvaluatorTensorProductAnisotropic::normal
static void normal(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const bool add_into_result=false, const int subface_index_1d=0)
Definition tensor_product_kernels.h:2157
internal::EvaluatorTensorProductAnisotropic::tangential
static void tangential(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const int subface_index_1d=0)
Definition tensor_product_kernels.h:2222
internal::EvaluatorTensorProduct< variant, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
Definition tensor_product_kernels.h:1931
internal::EvaluatorTensorProduct< variant, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
Definition tensor_product_kernels.h:1973
internal::EvaluatorTensorProduct::values
void values(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1673
internal::EvaluatorTensorProduct::EvaluatorTensorProduct
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition tensor_product_kernels.h:1633
internal::EvaluatorTensorProduct::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1723
internal::EvaluatorTensorProduct::gradients
void gradients(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1687
internal::EvaluatorTensorProduct::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int=0, const unsigned int=0)
Definition tensor_product_kernels.h:1593
internal::EvaluatorTensorProduct::n_rows_of_product
static constexpr unsigned int n_rows_of_product
Definition tensor_product_kernels.h:1575
internal::EvaluatorTensorProduct::hessians
void hessians(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1703
internal::EvaluatorTensorProduct::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition tensor_product_kernels.h:1839
internal::EvaluatorTensorProduct::n_columns_of_product
static constexpr unsigned int n_columns_of_product
Definition tensor_product_kernels.h:1577
internal::EvaluatorTensorProduct::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1739
internal::EvaluatorTensorProduct::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1758
Generated by doxygen 1.9.8