From: Denis Davydov <davydden@gmail.com>
Date: Wed, 5 Dec 2018 21:12:30 +0000 (+0100)
Subject: make LAPACKFullMatrix::transpose() do conjugate transpose for complex data types
X-Git-Tag: v9.1.0-rc1~504^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9ef60a4455fb77927aa125ccef070b5c381d4f64;p=dealii.git

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 1e88525bd2..f80b549486 100644
--- 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 904fac49c4..837191d34d 100644
--- 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 b26397d89b..62ca1372d0 100644
--- 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 fb71de2867..0aa61ff573 100644
--- 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