make LAPACKFullMatrix::transpose() do conjugate transpose for complex data types

diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h
index 1e88525bd2532fe5fa17e8d3d2c3c6c6aab77d73..f80b549486b4dbc92943b305efb5942a4b5afdcb 100644 (file)
--- a/include/deal.II/lac/lapack_full_matrix.h
+++ b/include/deal.II/lac/lapack_full_matrix.h
@@ -531,6 +531,8 @@ public:
    * Performs out-place transposition.
    * Matrix @p B should be appropriately sized.
    *
+   * @note for complex number types, conjugate transpose will be performed.
+   *
    * @note If deal.II is configured with Intel-MKL, `mkl_?omatcopy` will be used,
    * otherwise transposition is done element by element.
    */

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 904fac49c4a07eb70a5917a0bf9b21cb733542a8..837191d34d99f19ab84d04157478d575e5884a84 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -1041,11 +1041,12 @@ LAPACKFullMatrix<number>::transpose(LAPACKFullMatrix<number> &B) const
   const types::blas_int m = B.m();
   const types::blas_int n = B.n();
 #ifdef DEAL_II_LAPACK_WITH_MKL
-  mkl_omatcopy('C', 'T', n, m, 1., &A.values[0], n, &B.values[0], m);
+  const number one = 1.;
+  mkl_omatcopy('C', 'C', n, m, one, &A.values[0], n, &B.values[0], m);
 #else
   for (types::blas_int i = 0; i < m; ++i)
     for (types::blas_int j = 0; j < n; ++j)
-      B(i, j) = A(j, i);
+      B(i, j) = numbers::NumberTraits<number>::conjugate(A(j, i));
 #endif
 }
 

diff --git a/tests/lapack/full_matrix_36.cc b/tests/lapack/full_matrix_36.cc
index b26397d89b3f2b341ec51b35932d7eec97ccadd5..62ca1372d033427d70e1a466c8506b8bce633189 100644 (file)
--- a/tests/lapack/full_matrix_36.cc
+++ b/tests/lapack/full_matrix_36.cc
@@ -24,6 +24,15 @@
 #include "../tests.h"
 #include "create_matrix.h"
 
+template <typename T>
+std::string
+to_string(const T &t)
+{
+  std::ostringstream s;
+  s << t;
+  return s.str();
+}
+
 template <typename NumberType>
 void
 test(const unsigned int n = 3, const unsigned int k = 6)
@@ -36,12 +45,12 @@ test(const unsigned int n = 3, const unsigned int k = 6)
   for (unsigned int i = 0; i < A.m(); ++i)
     for (unsigned int j = 0; j < A.n(); ++j)
       {
-        const auto &at = At(j, i);
-        const auto &a  = A(i, j);
+        const auto at = At(j, i);
+        const auto a  = numbers::NumberTraits<NumberType>::conjugate(A(i, j));
         AssertThrow(at == a,
-                    ExcMessage(std::to_string(a) + "!=" + std::to_string(at) +
-                               " for (" + std::to_string(i) + "," +
-                               std::to_string(j) + ")"));
+                    ExcMessage(to_string(a) + "!=" + to_string(at) + " for (" +
+                               std::to_string(i) + "," + std::to_string(j) +
+                               ")"));
       }
 
   deallog << "OK" << std::endl;
@@ -54,4 +63,8 @@ main()
   test<double>(11, 27);
   test<double>(15, 4);
   test<double>(10, 10);
+
+  test<std::complex<double>>(11, 27);
+  test<std::complex<double>>(15, 4);
+  test<std::complex<double>>(10, 10);
 }

diff --git a/tests/lapack/full_matrix_36.output b/tests/lapack/full_matrix_36.output
index fb71de2867cdd56077fa628bd1498d2d19382019..0aa61ff573b91e6ae8d4a06474a0678fa9ec448f 100644 (file)
--- a/tests/lapack/full_matrix_36.output
+++ b/tests/lapack/full_matrix_36.output
@@ -2,3 +2,6 @@
 DEAL::OK
 DEAL::OK
 DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK