From: Martin Kronbichler <martin.kronbichler@rub.de>
Date: Fri, 15 Mar 2024 09:42:46 +0000 (+0100)
Subject: Tensor product kernels: specialize a function for complex numbers
X-Git-Tag: v9.6.0-rc1~475^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2ad3412e12ceb22aad05b792afa2e7fe82b2cfed;p=dealii.git

Tensor product kernels: specialize a function for complex numbers
---

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<Number>::is_complex &&
+                              numbers::NumberTraits<Number2>::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<Number>::is_complex &&
+                              numbers::NumberTraits<Number2>::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;