From 9ef60a4455fb77927aa125ccef070b5c381d4f64 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 5 Dec 2018 22:12:30 +0100 Subject: [PATCH] make LAPACKFullMatrix::transpose() do conjugate transpose for complex data types --- include/deal.II/lac/lapack_full_matrix.h | 2 ++ source/lac/lapack_full_matrix.cc | 5 +++-- tests/lapack/full_matrix_36.cc | 23 ++++++++++++++++++----- tests/lapack/full_matrix_36.output | 3 +++ 4 files changed, 26 insertions(+), 7 deletions(-) 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::transpose(LAPACKFullMatrix &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::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 +std::string +to_string(const T &t) +{ + std::ostringstream s; + s << t; + return s.str(); +} + template 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::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(11, 27); test(15, 4); test(10, 10); + + test>(11, 27); + test>(15, 4); + test>(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 -- 2.39.5