From 17013e99e16975e3c6d87391dfde0843dbf0df47 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 20 May 2020 15:01:22 -0400 Subject: [PATCH] Remove lapack_templates.h from more public headers. As mentioned previously, this file can conflict with other packages' LAPACK wrapper headers. Fortunately we only currently use three LAPACK functions in headers - its easy enough to special-case these so that we can get rid of the header inclusion. Fixes #9135. --- include/deal.II/lac/qr.h | 66 ++++++++++---- include/deal.II/lac/utilities.h | 39 +++++++-- source/lac/CMakeLists.txt | 2 + source/lac/qr.cc | 150 ++++++++++++++++++++++++++++++++ source/lac/utilities.cc | 89 +++++++++++++++++++ 5 files changed, 322 insertions(+), 24 deletions(-) create mode 100644 source/lac/qr.cc create mode 100644 source/lac/utilities.cc diff --git a/include/deal.II/lac/qr.h b/include/deal.II/lac/qr.h index 69ceeeda32..343038c482 100644 --- a/include/deal.II/lac/qr.h +++ b/include/deal.II/lac/qr.h @@ -22,7 +22,6 @@ #include #include -#include #include #include @@ -433,6 +432,41 @@ private: // ------------------- inline and template functions ---------------- #ifndef DOXYGEN +namespace internal +{ + namespace QRImplementation + { + // We want to avoid including our own LAPACK wrapper header in any external + // headers to avoid possible conflicts with other packages that may define + // their own such header. At the same time we want to be able to call some + // LAPACK functions from the template functions below. To resolve both + // problems define some extra wrappers here that can be in the header: + template + void + call_trmv(const char uplo, + const char trans, + const char diag, + const types::blas_int n, + const Number * a, + const types::blas_int lda, + Number * x, + const types::blas_int incx); + + template + void + call_trtrs(const char uplo, + const char trans, + const char diag, + const types::blas_int n, + const types::blas_int nrhs, + const Number * a, + const types::blas_int lda, + Number * b, + const types::blas_int ldb, + types::blas_int * info); + } // namespace QRImplementation +} // namespace internal + template @@ -482,16 +516,16 @@ BaseQR::solve(Vector & x, const int N = this->current_size; const int n_rhs = 1; int info = 0; - trtrs("U", - transpose ? "T" : "N", - "N", - &N, - &n_rhs, - &this->R(0, 0), - &lda, - &x(0), - &ldb, - &info); + internal::QRImplementation::call_trtrs('U', + transpose ? 'T' : 'N', + 'N', + N, + n_rhs, + &this->R(0, 0), + lda, + &x(0), + ldb, + &info); } @@ -579,8 +613,8 @@ ImplicitQR::append_column(const VectorType &column) const int N = this->current_size; const int n_rhs = 1; int info = 0; - trtrs( - "U", "T", "N", &N, &n_rhs, &this->R(0, 0), &lda, &u(0), &ldb, &info); + internal::QRImplementation::call_trtrs( + 'U', 'T', 'N', N, n_rhs, &this->R(0, 0), lda, &u(0), ldb, &info); // finally get the diagonal element: // rho2 = |column|^2 - |u|^2 @@ -852,7 +886,8 @@ QR::multiply_with_A(VectorType &y, const Vector &x) const const int N = this->current_size; const int lda = N; const int incx = 1; - trmv("U", "N", "N", &N, &this->R(0, 0), &lda, &x1[0], &incx); + internal::QRImplementation::call_trmv( + 'U', 'N', 'N', N, &this->R(0, 0), lda, &x1[0], incx); multiply_with_Q(y, x1); } @@ -868,7 +903,8 @@ QR::multiply_with_AT(Vector &y, const VectorType &x) const const int N = this->current_size; const int lda = N; const int incx = 1; - trmv("U", "T", "N", &N, &this->R(0, 0), &lda, &y[0], &incx); + internal::QRImplementation::call_trmv( + 'U', 'T', 'N', N, &this->R(0, 0), lda, &y[0], incx); } #endif // no DOXYGEN diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h index 24f7e8b05b..34c84f103a 100644 --- a/include/deal.II/lac/utilities.h +++ b/include/deal.II/lac/utilities.h @@ -22,7 +22,6 @@ #include #include -#include #include #include @@ -211,6 +210,28 @@ namespace Utilities #ifndef DOXYGEN +namespace internal +{ + namespace UtilitiesImplementation + { + // We want to avoid including our own LAPACK wrapper header in any external + // headers to avoid possible conflicts with other packages that may define + // their own such header. At the same time we want to be able to call some + // LAPACK functions from the template functions below. To resolve both + // problems define some extra wrappers here that can be in the header: + template + void + call_stev(const char jobz, + const types::blas_int n, + Number * d, + Number * e, + Number * z, + const types::blas_int ldz, + Number * work, + types::blas_int * info); + } // namespace UtilitiesImplementation +} // namespace internal + namespace Utilities { namespace LinearAlgebra @@ -380,14 +401,14 @@ namespace Utilities std::vector work; // ^^ types::blas_int info; // call lapack_templates.h wrapper: - stev("N", - &n, - diagonal.data(), - subdiagonal.data(), - Z.data(), - &ldz, - work.data(), - &info); + internal::UtilitiesImplementation::call_stev('N', + n, + diagonal.data(), + subdiagonal.data(), + Z.data(), + ldz, + work.data(), + &info); Assert(info == 0, LAPACKSupport::ExcErrorCode("dstev", info)); diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 07e08b9851..4a90cc4f47 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -53,9 +53,11 @@ SET(_unity_include_src SET(_separate_src full_matrix.cc lapack_full_matrix.cc + qr.cc sparse_matrix.cc sparse_matrix_inst2.cc tridiagonal_matrix.cc + utilities.cc ) # Add CUDA wrapper files diff --git a/source/lac/qr.cc b/source/lac/qr.cc new file mode 100644 index 0000000000..10612fddcf --- /dev/null +++ b/source/lac/qr.cc @@ -0,0 +1,150 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2020 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. +// +// --------------------------------------------------------------------- + +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace QRImplementation + { + // see the corresponding note in the header + template + void + call_trmv(const char uplo, + const char trans, + const char diag, + const types::blas_int n, + const Number * a, + const types::blas_int lda, + Number * x, + const types::blas_int incx) + { + trmv(&uplo, &trans, &diag, &n, a, &lda, x, &incx); + } + + template + void + call_trtrs(const char uplo, + const char trans, + const char diag, + const types::blas_int n, + const types::blas_int nrhs, + const Number * a, + const types::blas_int lda, + Number * b, + const types::blas_int ldb, + types::blas_int * info) + { + trtrs(&uplo, &trans, &diag, &n, &nrhs, a, &lda, b, &ldb, info); + } + + template void + call_trmv(const char, + const char, + const char, + const types::blas_int, + const float *, + const types::blas_int, + float *, + const types::blas_int); + + template void + call_trmv(const char, + const char, + const char, + const types::blas_int, + const double *, + const types::blas_int, + double *, + const types::blas_int); + + template void + call_trmv(const char, + const char, + const char, + const types::blas_int, + const std::complex *, + const types::blas_int, + std::complex *, + const types::blas_int); + + template void + call_trmv(const char, + const char, + const char, + const types::blas_int, + const std::complex *, + const types::blas_int, + std::complex *, + const types::blas_int); + + template void + call_trtrs(const char, + const char, + const char, + const types::blas_int, + const types::blas_int, + const float *, + const types::blas_int, + float *, + const types::blas_int, + types::blas_int *); + + template void + call_trtrs(const char, + const char, + const char, + const types::blas_int, + const types::blas_int, + const double *, + const types::blas_int, + double *, + const types::blas_int, + types::blas_int *); + + template void + call_trtrs(const char, + const char, + const char, + const types::blas_int, + const types::blas_int, + const std::complex *, + const types::blas_int, + std::complex *, + const types::blas_int, + types::blas_int *); + + template void + call_trtrs(const char, + const char, + const char, + const types::blas_int, + const types::blas_int, + const std::complex *, + const types::blas_int, + std::complex *, + const types::blas_int, + types::blas_int *); + } // namespace QRImplementation +} // namespace internal + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/utilities.cc b/source/lac/utilities.cc new file mode 100644 index 0000000000..78318693de --- /dev/null +++ b/source/lac/utilities.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2020 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. +// +// --------------------------------------------------------------------- + +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace UtilitiesImplementation + { + // see the corresponding note in the header + template + void + call_stev(const char jobz, + const types::blas_int n, + Number * d, + Number * e, + Number * z, + const types::blas_int ldz, + Number * work, + types::blas_int * info) + { + stev(&jobz, &n, d, e, z, &ldz, work, info); + } + + + template void + call_stev(const char, + const types::blas_int, + float *, + float *, + float *, + const types::blas_int, + float *, + types::blas_int *); + + template void + call_stev(const char, + const types::blas_int, + double *, + double *, + double *, + const types::blas_int, + double *, + types::blas_int *); + + template void + call_stev(const char, + const types::blas_int, + std::complex *, + std::complex *, + std::complex *, + const types::blas_int, + std::complex *, + types::blas_int *); + + template void + call_stev(const char, + const types::blas_int, + std::complex *, + std::complex *, + std::complex *, + const types::blas_int, + std::complex *, + types::blas_int *); + } // namespace UtilitiesImplementation +} // namespace internal + + + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5