From 51132fc04a75c2c8b8467aca2fbf0a8e5082a35b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 31 Oct 2017 15:48:10 -0600 Subject: [PATCH] Make LapackFullMatrix also work with std::complex data types. This is slightly complicated by the fact that the LAPACK functions sometimes have a different signature for complex data types. This patch does not address this -- these functions will simply fail if the signature is wrong. However, at least we can compile things with this patch. --- include/deal.II/lac/lapack_full_matrix.h | 5 +++-- source/lac/lapack_full_matrix.cc | 20 ++++++++++---------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 8228c23426..7d75ba796b 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -117,10 +117,11 @@ public: /** * This operator assigns a scalar to a matrix. To avoid confusion with - * constructors, zero is the only value allowed for d + * constructors, zero (when cast to the @p number type) is the only + * value allowed for d */ LAPACKFullMatrix & - operator = (const double d); + operator = (const number d); /** * This operator multiplies all entries by a fixed factor. diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index c6918d2a3b..8d00a95018 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -123,10 +123,10 @@ LAPACKFullMatrix::operator = (const SparseMatrix &M) template LAPACKFullMatrix & -LAPACKFullMatrix::operator = (const double d) +LAPACKFullMatrix::operator = (const number d) { (void)d; - Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); + Assert (d==number(0), ExcScalarAssignmentOnlyForZeroValue()); if (this->n_elements() != 0) this->reset_values(); @@ -718,7 +718,7 @@ LAPACKFullMatrix::compute_svd() AssertThrow (info==0, LAPACKSupport::ExcErrorCode("gesdd", info)); // Resize the work array. Add one to the size computed by LAPACK to be on // the safe side. - lwork = static_cast(work[0] + 1); + lwork = static_cast(std::fabs(work[0]) + 1); work.resize(lwork); // Do the actual SVD. @@ -743,11 +743,11 @@ LAPACKFullMatrix::compute_inverse_svd(const double threshold) Assert (state==LAPACKSupport::svd, ExcState(state)); - const double lim = wr[0]*threshold; + const double lim = std::fabs(wr[0])*threshold; for (size_type i=0; i lim) - wr[i] = 1./wr[i]; + if (std::fabs(wr[i]) > lim) + wr[i] = number(1.)/wr[i]; else wr[i] = 0.; } @@ -911,7 +911,7 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, Assert (info == 0, ExcInternalError()); // Allocate working array according to suggestion (same strategy as was // noted in compute_svd). - lwork = static_cast(work[0] + 1); + lwork = static_cast(std::fabs(work[0]) + 1); // resize workspace array work.resize((size_type ) lwork); @@ -984,7 +984,7 @@ LAPACKFullMatrix::compute_eigenvalues_symmetric(const number lowe Assert (info == 0, ExcInternalError()); // Allocate working array according to suggestion (same strategy as was noted in // compute_svd). - lwork = static_cast(work[0] + 1); + lwork = static_cast(std::fabs(work[0]) + 1); work.resize(static_cast (lwork)); // Finally compute the eigenvalues. @@ -1075,7 +1075,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( Assert (info == 0, ExcInternalError()); // Allocate working array according to suggestion (same strategy as was // noted in compute_svd). - lwork = static_cast(work[0] + 1); + lwork = static_cast(std::fabs(work[0]) + 1); // resize workspace arrays work.resize(static_cast (lwork)); @@ -1157,7 +1157,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( Assert (info == 0, ExcInternalError()); // Allocate working array according to suggestion (same strategy as was // noted in compute_svd). - lwork = static_cast(work[0] + 1); + lwork = static_cast(std::fabs(work[0]) + 1); // resize workspace array work.resize(static_cast(lwork)); -- 2.39.5