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 (file)
index 0000000..8e2a541
--- /dev/null
+++ b/doc/news/changes/minor/20181205DenisDavydov-2
@@ -0,0 +1,3 @@
+New: Add LAPACKFullMatrix::transpose() to perform out-of-place transposition.
+<br>
+(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 a38caaf24d62cd3249b3172565210a728f58b80b..1e88525bd2532fe5fa17e8d3d2c3c6c6aab77d73 100644 (file)
--- 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<number> &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<number> &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 fde6c6ed31ede478bec6df2a7e889c588a76e89c..62d5ebd86cf8d6257a27ec57729dc06ec1c38798 100644 (file)
--- a/include/deal.II/lac/lapack_templates.h
+++ b/include/deal.II/lac/lapack_templates.h
@@ -25,6 +25,15 @@
 #  include <cfenv>
 #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<float>
+#  define MKL_Complex16 std::complex<double>
+#  include "mkl_trans.h"
+#endif
+
 extern "C"
 {
   void
@@ -1409,6 +1418,138 @@ extern "C"
 DEAL_II_NAMESPACE_OPEN
 
 
+template <typename number1, typename number2, typename number3>
+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<float>  alpha,
+             const std::complex<float> *A,
+             dealii::types::blas_int    lda,
+             std::complex<float> *      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<double>  alpha,
+             const std::complex<double> *A,
+             dealii::types::blas_int     lda,
+             std::complex<double> *      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 <typename number1, typename number2, typename number3>
 inline void

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index e075424f7975f6b562f4bff2f6d8503dfaee485e..904fac49c4a07eb70a5917a0bf9b21cb733542a8 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -1031,6 +1031,25 @@ LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &      C,
 
 
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::transpose(LAPACKFullMatrix<number> &B) const
+{
+  const LAPACKFullMatrix<number> &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 <typename number>
 void
 LAPACKFullMatrix<number>::scale_rows(const Vector<number> &V)

diff --git a/tests/lapack/full_matrix_36.cc b/tests/lapack/full_matrix_36.cc
new file mode 100644 (file)
index 0000000..b26397d
--- /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 <deal.II/lac/lapack_full_matrix.h>
+
+#include <iostream>
+#include <tuple>
+
+#include "../tests.h"
+#include "create_matrix.h"
+
+template <typename NumberType>
+void
+test(const unsigned int n = 3, const unsigned int k = 6)
+{
+  LAPACKFullMatrix<NumberType> A(n, k);
+  LAPACKFullMatrix<NumberType> 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<double>(11, 27);
+  test<double>(15, 4);
+  test<double>(10, 10);
+}

diff --git a/tests/lapack/full_matrix_36.output b/tests/lapack/full_matrix_36.output
new file mode 100644 (file)
index 0000000..fb71de2
--- /dev/null
+++ b/tests/lapack/full_matrix_36.output
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK