206 const int stride_in_given,
207 const int stride_out_given)
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,
213 Assert(n_rows > 0 && n_columns > 0,
215 std::to_string(n_rows) +
", " +
216 std::to_string(n_columns) +
" was passed!"));
219 "This function should only use EvaluatorQuantity::value");
221 Assert(consider_strides || (stride_in_given == 1 && stride_out_given == 1),
223 const int stride_in = consider_strides ? stride_in_given : 1;
224 const int stride_out = consider_strides ? stride_out_given : 1;
226 static_assert(n_components > 0 && n_components < 4,
227 "Invalid number of components");
234 if (transpose_matrix && n_rows == 2 && n_components == 1)
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)
240 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
242 out[stride_out * col] += result;
244 out[stride_out * col] = result;
247 else if (transpose_matrix && n_rows == 3 && n_components == 1)
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)
254 const Number result =
255 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
257 out[stride_out * col] += result;
259 out[stride_out * col] = result;
267 std::array<Number, 129> x;
268 for (
int i = 0; i < mm; ++i)
269 x[i] = in[stride_in * i];
271 for (
int col = 0; col < nn; ++col)
274 if (transpose_matrix ==
true)
276 res0 = matrix[col] * x[0];
277 for (
int i = 1; i < mm; ++i)
278 res0 += matrix[i * n_columns + col] * x[i];
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];
287 out[stride_out * col] += res0;
289 out[stride_out * col] = res0;
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;
299 Number *out1 = n_components > 1 ? out + nn :
nullptr;
300 Number *out2 = n_components > 2 ? out + 2 * nn :
nullptr;
302 int nn_regular = (nn / 4) * 4;
303 for (
int col = 0; col < nn_regular; col += 4)
306 if (transpose_matrix ==
true)
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;
315 if (n_components > 1)
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;
324 if (n_components > 2)
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;
333 matrix_ptr += n_columns;
334 for (
int i = 1; i < mm; ++i, matrix_ptr += n_columns)
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;
342 if (n_components > 1)
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;
350 if (n_components > 2)
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;
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;
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;
373 if (n_components > 1)
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;
382 if (n_components > 2)
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;
391 for (
int i = 1; i < mm; ++i)
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;
399 if (n_components > 1)
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;
408 if (n_components > 2)
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;
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)
427 out1[stride_out] += res[5];
428 out1[2 * stride_out] += res[6];
429 out1[3 * stride_out] += res[7];
431 if (n_components > 2)
434 out2[stride_out] += res[9];
435 out2[2 * stride_out] += res[10];
436 out2[3 * stride_out] += res[11];
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)
448 out1[stride_out] = res[5];
449 out1[2 * stride_out] = res[6];
450 out1[3 * stride_out] = res[7];
452 if (n_components > 2)
455 out2[stride_out] = res[9];
456 out2[2 * stride_out] = res[10];
457 out2[3 * stride_out] = res[11];
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;
466 if (nn - nn_regular == 3)
468 Number res0, res1, res2, res3, res4, res5, res6, res7, res8;
469 if (transpose_matrix ==
true)
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)
477 res3 = matrix_ptr[0] * in1[0];
478 res4 = matrix_ptr[1] * in1[0];
479 res5 = matrix_ptr[2] * in1[0];
481 if (n_components > 2)
483 res6 = matrix_ptr[0] * in2[0];
484 res7 = matrix_ptr[1] * in2[0];
485 res8 = matrix_ptr[2] * in2[0];
487 matrix_ptr += n_columns;
488 for (
int i = 1; i < mm; ++i, matrix_ptr += n_columns)
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)
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];
499 if (n_components > 2)
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];
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;
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)
518 res3 = matrix_0[0] * in1[0];
519 res4 = matrix_1[0] * in1[0];
520 res5 = matrix_2[0] * in1[0];
522 if (n_components > 2)
524 res6 = matrix_0[0] * in2[0];
525 res7 = matrix_1[0] * in2[0];
526 res8 = matrix_2[0] * in2[0];
528 for (
int i = 1; i < mm; ++i)
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)
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];
539 if (n_components > 2)
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];
550 out0[stride_out] += res1;
551 out0[2 * stride_out] += res2;
552 if (n_components > 1)
555 out1[stride_out] += res4;
556 out1[2 * stride_out] += res5;
558 if (n_components > 2)
561 out2[stride_out] += res7;
562 out2[2 * stride_out] += res8;
568 out0[stride_out] = res1;
569 out0[2 * stride_out] = res2;
570 if (n_components > 1)
573 out1[stride_out] = res4;
574 out1[2 * stride_out] = res5;
576 if (n_components > 2)
579 out2[stride_out] = res7;
580 out2[2 * stride_out] = res8;
584 else if (nn - nn_regular == 2)
586 Number res0, res1, res2, res3, res4, res5;
587 if (transpose_matrix ==
true)
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)
594 res2 = matrix_ptr[0] * in1[0];
595 res3 = matrix_ptr[1] * in1[0];
597 if (n_components > 2)
599 res4 = matrix_ptr[0] * in2[0];
600 res5 = matrix_ptr[1] * in2[0];
602 matrix_ptr += n_columns;
603 for (
int i = 1; i < mm; ++i, matrix_ptr += n_columns)
605 res0 += matrix_ptr[0] * in0[stride_in * i];
606 res1 += matrix_ptr[1] * in0[stride_in * i];
607 if (n_components > 1)
609 res2 += matrix_ptr[0] * in1[stride_in * i];
610 res3 += matrix_ptr[1] * in1[stride_in * i];
612 if (n_components > 2)
614 res4 += matrix_ptr[0] * in2[stride_in * i];
615 res5 += matrix_ptr[1] * in2[stride_in * i];
621 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
622 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
624 res0 = matrix_0[0] * in0[0];
625 res1 = matrix_1[0] * in0[0];
626 if (n_components > 1)
628 res2 = matrix_0[0] * in1[0];
629 res3 = matrix_1[0] * in1[0];
631 if (n_components > 2)
633 res4 = matrix_0[0] * in2[0];
634 res5 = matrix_1[0] * in2[0];
636 for (
int i = 1; i < mm; ++i)
638 res0 += matrix_0[i] * in0[stride_in * i];
639 res1 += matrix_1[i] * in0[stride_in * i];
640 if (n_components > 1)
642 res2 += matrix_0[i] * in1[stride_in * i];
643 res3 += matrix_1[i] * in1[stride_in * i];
645 if (n_components > 2)
647 res4 += matrix_0[i] * in2[stride_in * i];
648 res5 += matrix_1[i] * in2[stride_in * i];
655 out0[stride_out] += res1;
656 if (n_components > 1)
659 out1[stride_out] += res3;
661 if (n_components > 2)
664 out2[stride_out] += res5;
670 out0[stride_out] = res1;
671 if (n_components > 1)
674 out1[stride_out] = res3;
676 if (n_components > 2)
679 out2[stride_out] = res5;
683 else if (nn - nn_regular == 1)
685 Number res0, res1, res2;
686 if (transpose_matrix ==
true)
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)
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];
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)
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];
724 if (n_components > 1)
726 if (n_components > 2)
732 if (n_components > 1)
734 if (n_components > 2)
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,
771 std::to_string(n_rows) +
", " +
772 std::to_string(n_columns) +
" was passed!"));
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;
779 std::array<Number, mm> x;
780 for (
int i = 0; i < mm; ++i)
781 x[i] = in[stride_in * i];
804 for (
int col = 0; col < n_cols; ++col)
808 if (transpose_matrix ==
true)
811 val1 = matrix[nn - 1 - col];
815 val0 = matrix[col * n_columns];
816 val1 = matrix[(col + 1) * n_columns - 1];
822 res0 += val1 * x[mm - 1];
823 res1 += val0 * x[mm - 1];
824 for (
int ind = 1; ind < mid; ++ind)
826 if (transpose_matrix ==
true)
828 val0 = matrix[ind * n_columns + col];
829 val1 = matrix[ind * n_columns + nn - 1 - col];
833 val0 = matrix[col * n_columns + ind];
834 val1 = matrix[(col + 1) * n_columns - 1 - ind];
836 res0 += val0 * x[ind];
837 res1 += val1 * x[ind];
838 res0 += val1 * x[mm - 1 - ind];
839 res1 += val0 * x[mm - 1 - ind];
843 res0 = res1 = Number();
844 if (transpose_matrix ==
true)
848 const Number tmp = matrix[mid * n_columns + col] * x[mid];
855 if (mm % 2 == 1 && nn % 2 == 0)
857 const Number tmp = matrix[col * n_columns + mid] * x[mid];
864 out[stride_out * col] += res0;
865 out[stride_out * (nn - 1 - col)] += res1;
869 out[stride_out * col] = res0;
870 out[stride_out * (nn - 1 - col)] = res1;
873 if (transpose_matrix ==
true && nn % 2 == 1 && mm % 2 == 1)
876 out[stride_out * n_cols] += x[mid];
878 out[stride_out * n_cols] = x[mid];
880 else if (transpose_matrix ==
true && nn % 2 == 1)
885 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
886 for (
int ind = 1; ind < mid; ++ind)
888 const Number2 val0 = matrix[ind * n_columns + n_cols];
889 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
895 out[stride_out * n_cols] += res0;
897 out[stride_out * n_cols] = res0;
899 else if (transpose_matrix ==
false && nn % 2 == 1)
904 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
905 for (
int ind = 1; ind < mid; ++ind)
907 const Number2 val0 = matrix[n_cols * n_columns + ind];
908 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
916 out[stride_out * n_cols] += res0;
918 out[stride_out * n_cols] = res0;
939 for (
int col = 0; col < n_cols; ++col)
943 if (transpose_matrix ==
true)
946 val1 = matrix[nn - 1 - col];
950 val0 = matrix[col * n_columns];
951 val1 = matrix[(nn - col - 1) * n_columns];
957 res0 -= val1 * x[mm - 1];
958 res1 -= val0 * x[mm - 1];
959 for (
int ind = 1; ind < mid; ++ind)
961 if (transpose_matrix ==
true)
963 val0 = matrix[ind * n_columns + col];
964 val1 = matrix[ind * n_columns + nn - 1 - col];
968 val0 = matrix[col * n_columns + ind];
969 val1 = matrix[(nn - col - 1) * n_columns + ind];
971 res0 += val0 * x[ind];
972 res1 += val1 * x[ind];
973 res0 -= val1 * x[mm - 1 - ind];
974 res1 -= val0 * x[mm - 1 - ind];
978 res0 = res1 = Number();
981 if (transpose_matrix ==
true)
982 val0 = matrix[mid * n_columns + col];
984 val0 = matrix[col * n_columns + mid];
985 const Number tmp = val0 * x[mid];
991 out[stride_out * col] += res0;
992 out[stride_out * (nn - 1 - col)] += res1;
996 out[stride_out * col] = res0;
997 out[stride_out * (nn - 1 - col)] = res1;
1004 if (transpose_matrix ==
true)
1005 val0 = matrix[n_cols];
1007 val0 = matrix[n_cols * n_columns];
1008 res0 = val0 * (x[0] - x[mm - 1]);
1009 for (
int ind = 1; ind < mid; ++ind)
1011 if (transpose_matrix ==
true)
1012 val0 = matrix[ind * n_columns + n_cols];
1014 val0 = matrix[n_cols * n_columns + ind];
1015 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
1019 out[stride_out * n_cols] += res0;
1021 out[stride_out * n_cols] = res0;
1028 for (
int col = 0; col < n_cols; ++col)
1032 if (transpose_matrix ==
true)
1035 val1 = matrix[nn - 1 - col];
1039 val0 = matrix[col * n_columns];
1040 val1 = matrix[(col + 1) * n_columns - 1];
1046 res0 += val1 * x[mm - 1];
1047 res1 += val0 * x[mm - 1];
1048 for (
int ind = 1; ind < mid; ++ind)
1050 if (transpose_matrix ==
true)
1052 val0 = matrix[ind * n_columns + col];
1053 val1 = matrix[ind * n_columns + nn - 1 - col];
1057 val0 = matrix[col * n_columns + ind];
1058 val1 = matrix[(col + 1) * n_columns - 1 - ind];
1060 res0 += val0 * x[ind];
1061 res1 += val1 * x[ind];
1062 res0 += val1 * x[mm - 1 - ind];
1063 res1 += val0 * x[mm - 1 - ind];
1067 res0 = res1 = Number();
1070 if (transpose_matrix ==
true)
1071 val0 = matrix[mid * n_columns + col];
1073 val0 = matrix[col * n_columns + mid];
1074 const Number tmp = val0 * x[mid];
1080 out[stride_out * col] += res0;
1081 out[stride_out * (nn - 1 - col)] += res1;
1085 out[stride_out * col] = res0;
1086 out[stride_out * (nn - 1 - col)] = res1;
1093 if (transpose_matrix ==
true)
1094 val0 = matrix[n_cols];
1096 val0 = matrix[n_cols * n_columns];
1099 res0 = val0 * (x[0] + x[mm - 1]);
1100 for (
int ind = 1; ind < mid; ++ind)
1102 if (transpose_matrix ==
true)
1103 val0 = matrix[ind * n_columns + n_cols];
1105 val0 = matrix[n_cols * n_columns + ind];
1106 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
1114 if (transpose_matrix ==
true)
1115 val0 = matrix[mid * n_columns + n_cols];
1117 val0 = matrix[n_cols * n_columns + mid];
1118 res0 += val0 * x[mid];
1121 out[stride_out * n_cols] += res0;
1123 out[stride_out * n_cols] = res0;
1165 int n_rows_runtime = 0,
1166 int n_columns_runtime = 0,
1167 int stride_in_runtime = 0,
1168 int stride_out_runtime = 0)
1170 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
1171 "Negative loop ranges are not allowed!");
1173 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
1174 const int n_columns =
1175 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
1176 const int stride_in =
1177 stride_in_static == 0 ? stride_in_runtime : stride_in_static;
1178 const int stride_out =
1179 stride_out_static == 0 ? stride_out_runtime : stride_out_static;
1181 Assert(n_rows > 0 && n_columns > 0,
1183 std::to_string(n_rows) +
", " +
1184 std::to_string(n_columns) +
" was passed!"));
1186 const int mm = transpose_matrix ? n_rows : n_columns,
1187 nn = transpose_matrix ? n_columns : n_rows;
1188 const int n_half = nn / 2;
1189 const int m_half = mm / 2;
1191 constexpr int array_length =
1192 (n_rows_static == 0) ?
1195 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
1196 const int offset = (n_columns + 1) / 2;
1200 std::array<Number, array_length> xp, xm;
1201 for (
int i = 0; i < m_half; ++i)
1205 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1206 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1210 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1211 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1214 Number xmid = in[stride_in * m_half];
1215 for (
int col = 0; col < n_half; ++col)
1220 if (transpose_matrix ==
true)
1222 r0 = matrix[col] * xp[0];
1223 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
1227 r0 = matrix[col * offset] * xp[0];
1228 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
1230 for (
int ind = 1; ind < m_half; ++ind)
1232 if (transpose_matrix ==
true)
1234 r0 += matrix[ind * offset + col] * xp[ind];
1235 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
1239 r0 += matrix[col * offset + ind] * xp[ind];
1240 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
1246 if (mm % 2 == 1 && transpose_matrix ==
true)
1249 r1 += matrix[m_half * offset + col] * xmid;
1251 r0 += matrix[m_half * offset + col] * xmid;
1253 else if (mm % 2 == 1 &&
1256 r0 += matrix[col * offset + m_half] * xmid;
1260 out[stride_out * col] += r0 + r1;
1262 transpose_matrix ==
false)
1263 out[stride_out * (nn - 1 - col)] += r1 - r0;
1265 out[stride_out * (nn - 1 - col)] += r0 - r1;
1269 out[stride_out * col] = r0 + r1;
1271 transpose_matrix ==
false)
1272 out[stride_out * (nn - 1 - col)] = r1 - r0;
1274 out[stride_out * (nn - 1 - col)] = r0 - r1;
1278 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
1281 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
1283 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
1285 else if (transpose_matrix ==
true && nn % 2 == 1)
1290 r0 = matrix[n_half] * xp[0];
1291 for (
int ind = 1; ind < m_half; ++ind)
1292 r0 += matrix[ind * offset + n_half] * xp[ind];
1297 r0 += matrix[m_half * offset + n_half] * xmid;
1300 out[stride_out * n_half] += r0;
1302 out[stride_out * n_half] = r0;
1304 else if (transpose_matrix ==
false && nn % 2 == 1)
1311 r0 = matrix[n_half * offset] * xm[0];
1312 for (
int ind = 1; ind < m_half; ++ind)
1313 r0 += matrix[n_half * offset + ind] * xm[ind];
1317 r0 = matrix[n_half * offset] * xp[0];
1318 for (
int ind = 1; ind < m_half; ++ind)
1319 r0 += matrix[n_half * offset + ind] * xp[ind];
1326 r0 += matrix[n_half * offset + m_half] * xmid;
1329 out[stride_out * n_half] += r0;
1331 out[stride_out * n_half] = r0;
1400 static_assert(n_rows > 0 && n_columns > 0,
1401 "Specialization requires n_rows, n_columns > 0");
1403 constexpr bool evaluate_antisymmetric =
1406 constexpr int mm = transpose_matrix ? n_rows : n_columns,
1407 nn = transpose_matrix ? n_columns : n_rows;
1408 constexpr int n_half = nn / 2;
1409 constexpr int m_half = mm / 2;
1411 if (transpose_matrix)
1413 std::array<Number, mm> x;
1414 for (
unsigned int i = 0; i < mm; ++i)
1415 x[i] = in[stride_in * i];
1416 for (
unsigned int col = 0; col < n_half; ++col)
1421 r0 = matrix[col] * x[0];
1422 r1 = matrix[col + n_columns] * x[1];
1423 for (
unsigned int ind = 1; ind < m_half; ++ind)
1425 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
1427 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
1433 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
1436 out[stride_out * col] += r0 + r1;
1437 if (evaluate_antisymmetric)
1438 out[stride_out * (nn - 1 - col)] += r1 - r0;
1440 out[stride_out * (nn - 1 - col)] += r0 - r1;
1444 out[stride_out * col] = r0 + r1;
1445 if (evaluate_antisymmetric)
1446 out[stride_out * (nn - 1 - col)] = r1 - r0;
1448 out[stride_out * (nn - 1 - col)] = r0 - r1;
1454 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1457 r0 = matrix[n_half + shift * n_columns] * x[shift];
1458 for (
unsigned int ind = 1; ind < m_half; ++ind)
1459 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1464 if (!evaluate_antisymmetric && mm % 2 == 1)
1465 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1467 out[stride_out * n_half] += r0;
1469 out[stride_out * n_half] = r0;
1474 std::array<Number, m_half + 1> xp, xm;
1475 for (
int i = 0; i < m_half; ++i)
1476 if (!evaluate_antisymmetric)
1478 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1479 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1483 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1484 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1487 xp[m_half] = in[stride_in * m_half];
1488 for (
unsigned int col = 0; col < n_half; ++col)
1493 r0 = matrix[2 * col * n_columns] * xp[0];
1494 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1495 for (
unsigned int ind = 1; ind < m_half; ++ind)
1497 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1498 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1505 if (evaluate_antisymmetric)
1506 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1508 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1512 out[stride_out * (2 * col)] += r0;
1513 out[stride_out * (2 * col + 1)] += r1;
1517 out[stride_out * (2 * col)] = r0;
1518 out[stride_out * (2 * col + 1)] = r1;
1526 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1527 for (
unsigned int ind = 1; ind < m_half; ++ind)
1528 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1532 if (mm % 2 == 1 && !evaluate_antisymmetric)
1533 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1535 out[stride_out * (nn - 1)] += r0;
1537 out[stride_out * (nn - 1)] = r0;