From 2ad3412e12ceb22aad05b792afa2e7fe82b2cfed Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 15 Mar 2024 10:42:46 +0100 Subject: [PATCH] Tensor product kernels: specialize a function for complex numbers --- .../matrix_free/tensor_product_kernels.h | 28 +++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 4f2319420b..237511c5c7 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -146,13 +146,37 @@ namespace internal { res0 = matrix[col] * x[0]; for (int i = 1; i < mm; ++i) - res0 += matrix[i * n_columns + col] * x[i]; + { + const Number2 mji = matrix[i * n_columns + col]; + if constexpr (numbers::NumberTraits::is_complex && + numbers::NumberTraits::is_complex) + { + res0.real(res0.real() + mji.real() * x[i].real() - + mji.imag() * x[i].imag()); + res0.imag(res0.imag() + mji.imag() * x[i].real() + + mji.real() * x[i].imag()); + } + else + res0 += mji * x[i]; + } } else { res0 = matrix[col * n_columns] * x[0]; for (int i = 1; i < mm; ++i) - res0 += matrix[col * n_columns + i] * x[i]; + { + const Number2 mij = matrix[col * n_columns + i]; + if constexpr (numbers::NumberTraits::is_complex && + numbers::NumberTraits::is_complex) + { + res0.real(res0.real() + mij.real() * x[i].real() - + mij.imag() * x[i].imag()); + res0.imag(res0.imag() + mij.imag() * x[i].real() + + mij.real() * x[i].imag()); + } + else + res0 += mij * x[i]; + } } if (add) out[stride_out * col] += res0; -- 2.39.5