From: Denis Davydov Date: Wed, 5 Dec 2018 17:46:01 +0000 (+0100) Subject: Add LAPACKFullMatrix::transpose() to perform out-of-place transposition. X-Git-Tag: v9.1.0-rc1~504^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=14f2c61f17ef1b9e4b718d6fddc7372496df9e94;p=dealii.git Add LAPACKFullMatrix::transpose() to perform out-of-place transposition. --- diff --git a/doc/news/changes/minor/20181205DenisDavydov-2 b/doc/news/changes/minor/20181205DenisDavydov-2 new file mode 100644 index 0000000000..8e2a5415c8 --- /dev/null +++ b/doc/news/changes/minor/20181205DenisDavydov-2 @@ -0,0 +1,3 @@ +New: Add LAPACKFullMatrix::transpose() to perform out-of-place transposition. +
+(Denis Davydov, 2018/12/05) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index a38caaf24d..1e88525bd2 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -527,6 +527,16 @@ public: const LAPACKFullMatrix &B, const bool adding = false) const; + /** + * Performs out-place transposition. + * Matrix @p B should be appropriately sized. + * + * @note If deal.II is configured with Intel-MKL, `mkl_?omatcopy` will be used, + * otherwise transposition is done element by element. + */ + void + transpose(LAPACKFullMatrix &B) const; + /** * Scale rows of this matrix by @p V . This is equivalent to premultiplication * with a diagonal matrix $\mathbf A\leftarrow {\rm diag}(\mathbf V)\mathbf diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index fde6c6ed31..62d5ebd86c 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -25,6 +25,15 @@ # include #endif +// Intel-MKL specific functions +#ifdef DEAL_II_LAPACK_WITH_MKL +// see +// https://software.intel.com/en-us/mkl-windows-developer-guide-using-complex-types-in-c-c +# define MKL_Complex8 std::complex +# define MKL_Complex16 std::complex +# include "mkl_trans.h" +#endif + extern "C" { void @@ -1409,6 +1418,138 @@ extern "C" DEAL_II_NAMESPACE_OPEN +template +inline void +mkl_omatcopy(char, + char, + dealii::types::blas_int, + dealii::types::blas_int, + const number1, + const number2 *, + dealii::types::blas_int, + number3 *, + dealii::types::blas_int) +{ + Assert(false, ExcNotImplemented()); +} + + + +inline void +mkl_omatcopy(char ordering, + char trans, + dealii::types::blas_int rows, + dealii::types::blas_int cols, + const float alpha, + const float * A, + dealii::types::blas_int lda, + float * B, + dealii::types::blas_int ldb) +{ +#ifdef DEAL_II_LAPACK_WITH_MKL + mkl_somatcopy(ordering, trans, rows, cols, alpha, A, lda, B, ldb); +#else + (void)ordering; + (void)trans; + (void)rows; + (void)cols; + (void)alpha; + (void)A; + (void)lda; + (void)B; + (void)ldb; + Assert(false, LAPACKSupport::ExcMissing("mkl_somatcopy")); +#endif // DEAL_II_LAPACK_WITH_MKL +} + + + +inline void +mkl_omatcopy(char ordering, + char trans, + dealii::types::blas_int rows, + dealii::types::blas_int cols, + const double alpha, + const double * A, + dealii::types::blas_int lda, + double * B, + dealii::types::blas_int ldb) +{ +#ifdef DEAL_II_LAPACK_WITH_MKL + mkl_domatcopy(ordering, trans, rows, cols, alpha, A, lda, B, ldb); +#else + (void)ordering; + (void)trans; + (void)rows; + (void)cols; + (void)alpha; + (void)A; + (void)lda; + (void)B; + (void)ldb; + Assert(false, LAPACKSupport::ExcMissing("mkl_domatcopy")); +#endif // DEAL_II_LAPACK_WITH_MKL +} + + + +inline void +mkl_omatcopy(char ordering, + char trans, + dealii::types::blas_int rows, + dealii::types::blas_int cols, + const std::complex alpha, + const std::complex *A, + dealii::types::blas_int lda, + std::complex * B, + dealii::types::blas_int ldb) +{ +#ifdef DEAL_II_LAPACK_WITH_MKL + mkl_comatcopy(ordering, trans, rows, cols, alpha, A, lda, B, ldb); +#else + (void)ordering; + (void)trans; + (void)rows; + (void)cols; + (void)alpha; + (void)A; + (void)lda; + (void)B; + (void)ldb; + Assert(false, LAPACKSupport::ExcMissing("mkl_comatcopy")); +#endif // DEAL_II_LAPACK_WITH_MKL +} + + + +inline void +mkl_omatcopy(char ordering, + char trans, + dealii::types::blas_int rows, + dealii::types::blas_int cols, + const std::complex alpha, + const std::complex *A, + dealii::types::blas_int lda, + std::complex * B, + dealii::types::blas_int ldb) +{ +#ifdef DEAL_II_LAPACK_WITH_MKL + mkl_zomatcopy(ordering, trans, rows, cols, alpha, A, lda, B, ldb); +#else + (void)ordering; + (void)trans; + (void)rows; + (void)cols; + (void)alpha; + (void)A; + (void)lda; + (void)B; + (void)ldb; + Assert(false, LAPACKSupport::ExcMissing("mkl_zomatcopy")); +#endif // DEAL_II_LAPACK_WITH_MKL +} + + template inline void diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index e075424f79..904fac49c4 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -1031,6 +1031,25 @@ LAPACKFullMatrix::Tmmult(LAPACKFullMatrix & C, +template +void +LAPACKFullMatrix::transpose(LAPACKFullMatrix &B) const +{ + const LAPACKFullMatrix &A = *this; + AssertDimension(A.m(), B.n()); + AssertDimension(A.n(), B.m()); + 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); +#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); +#endif +} + + template void LAPACKFullMatrix::scale_rows(const Vector &V) diff --git a/tests/lapack/full_matrix_36.cc b/tests/lapack/full_matrix_36.cc new file mode 100644 index 0000000000..b26397d89b --- /dev/null +++ b/tests/lapack/full_matrix_36.cc @@ -0,0 +1,57 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Tests LAPACKFullMatrix::transpose(...) const + +#include + +#include +#include + +#include "../tests.h" +#include "create_matrix.h" + +template +void +test(const unsigned int n = 3, const unsigned int k = 6) +{ + LAPACKFullMatrix A(n, k); + LAPACKFullMatrix At(k, n); + + create_random(A); + A.transpose(At); + 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); + AssertThrow(at == a, + ExcMessage(std::to_string(a) + "!=" + std::to_string(at) + + " for (" + std::to_string(i) + "," + + std::to_string(j) + ")")); + } + + deallog << "OK" << std::endl; +} + +int +main() +{ + initlog(); + 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 new file mode 100644 index 0000000000..fb71de2867 --- /dev/null +++ b/tests/lapack/full_matrix_36.output @@ -0,0 +1,4 @@ + +DEAL::OK +DEAL::OK +DEAL::OK