From f8984293cd145e14c721b8809b0ff8048eda8ca5 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Mon, 19 Mar 2018 16:57:42 +0100 Subject: [PATCH] Corrections --- include/deal.II/lac/lapack_templates.h | 69 -------------------------- include/deal.II/lac/scalapack.h | 16 ++++-- source/lac/scalapack.cc | 25 ++++------ tests/scalapack/scalapack_15.cc | 2 +- tests/scalapack/scalapack_15_b.cc | 2 +- 5 files changed, 24 insertions(+), 90 deletions(-) diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index cab2707058..e3c11a310e 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -1403,36 +1403,7 @@ syevr(const char *jobz, { dsyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info); } -#else -inline void -syevr(const char * /*jobz*/, - const char * /*range*/, - const char * /*uplo*/, - const types::blas_int * /*n*/, - double * /*A*/, - const types::blas_int * /*lda*/, - const double * /*vl*/, - const double * /*vu*/, - const types::blas_int * /*il*/, - const types::blas_int * /*iu*/, - const double * /*abstol*/, - types::blas_int * /*m*/, - double * /*w*/, - double * /*z*/, - const types::blas_int * /*ldz*/, - types::blas_int * /*isuppz*/, - double * /*work*/, - types::blas_int * /*lwork*/, - types::blas_int * /*iwork*/, - types::blas_int * /*liwork*/, - types::blas_int * /*info*/) -{ - Assert (false, LAPACKSupport::ExcMissing("dsyevr")); -} -#endif - -#ifdef DEAL_II_WITH_LAPACK inline void syevr(const char *jobz, const char *range, @@ -1458,32 +1429,6 @@ syevr(const char *jobz, { ssyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info); } -#else -inline void -syevr(const char * /*jobz*/, - const char * /*range*/, - const char * /*uplo*/, - const types::blas_int * /*n*/, - float * /*A*/, - const types::blas_int * /*lda*/, - const float * /*vl*/, - const float * /*vu*/, - const types::blas_int * /*il*/, - const types::blas_int * /*iu*/, - const float * /*abstol*/, - types::blas_int * /*m*/, - float * /*w*/, - float * /*z*/, - const types::blas_int * /*ldz*/, - types::blas_int * /*isuppz*/, - float * /*work*/, - types::blas_int * /*lwork*/, - types::blas_int * /*iwork*/, - types::blas_int * /*liwork*/, - types::blas_int * /*info*/) -{ - Assert (false, LAPACKSupport::ExcMissing("ssyevr")); -} #endif @@ -1795,26 +1740,12 @@ lamch(const char *chmach, double &precision) { precision = dlamch_(chmach); } -#else -inline void -lamch(const char * /*chmach*/, double & /*precision*/) -{ - Assert(false, LAPACKSupport::ExcMissing("dlamch")); -} -#endif -#ifdef DEAL_II_WITH_LAPACK inline void lamch(const char *chmach, float &precision) { precision = slamch_(chmach); } -#else -inline void -lamch(const char * /*chmach*/, float & /*precision*/) -{ - Assert(false, LAPACKSupport::ExcMissing("slamch")); -} #endif DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 5a4aa99ac5..0f0a1bf3e7 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -423,7 +423,7 @@ public: /** * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric - * matrix $A \in \mathbb{R}^{M \times M}$ using the MRRR algorithm. + * matrix $\mathbf{A} \in \mathbb{R}^{M \times M}$ using the MRRR algorithm. * * The eigenvalues/eigenvectors are selected by prescribing a range of indices @p index_limits. * @@ -439,7 +439,8 @@ public: const bool compute_eigenvectors); /** - * Computing selected eigenvalues and, optionally, the eigenvectors using the MRRR algorithm. + * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric + * matrix $\mathbf{A} \in \mathbb{R}^{M \times M}$ using the MRRR algorithm. * The eigenvalues/eigenvectors are selected by prescribing a range of values @p value_limits for the eigenvalues. * * If successful, the computed eigenvalues are arranged in ascending order. @@ -635,14 +636,19 @@ private: std::make_pair(std::numeric_limits::quiet_NaN(),std::numeric_limits::quiet_NaN())); /** - * Computing selected eigenvalues and, optionally, the eigenvectors using the MRRR algorithm. + * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric + * matrix $\mathbf{A} \in \mathbb{R}^{M \times M}$ using the MRRR algorithm. * The eigenvalues/eigenvectors are selected by either prescribing a range of indices @p index_limits * or a range of values @p value_limits for the eigenvalues. The function will throw an exception * if both ranges are prescribed (meaning that both ranges differ from the default value) * as this ambiguity is prohibited. + * + * By calling this function the original content of the matrix will be overwritten. + * If requested, the eigenvectors are stored in the columns of the matrix. + * Also in the case that just the eigenvalues are required, + * the content of the matrix will be overwritten. + * * If successful, the computed eigenvalues are arranged in ascending order. - * The eigenvectors are stored in the columns of the matrix, thereby - * overwriting the original content of the matrix. */ std::vector eigenpairs_symmetric_MRRR(const bool compute_eigenvectors, const std::pair &index_limits= diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 6ab1b4439a..97ad95f20a 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -915,14 +915,14 @@ std::vector ScaLAPACKMatrix::eigenpairs_symmetric_by_ind const bool compute_eigenvectors) { // Check validity of index limits. - Assert (index_limits.first < (unsigned int)n_rows,ExcIndexRange(index_limits.first,0,n_rows)); - Assert (index_limits.second < (unsigned int)n_rows,ExcIndexRange(index_limits.second,0,n_rows)); + AssertIndexRange(index_limits.first,static_cast(n_rows)); + AssertIndexRange(index_limits.second,static_cast(n_rows)); - std::pair idx = std::make_pair(std::min(index_limits.first,index_limits.second), - std::max(index_limits.first,index_limits.second)); + const std::pair idx = std::make_pair(std::min(index_limits.first,index_limits.second), + std::max(index_limits.first,index_limits.second)); // Compute all eigenvalues/eigenvectors. - if (idx.first==0 && idx.second==(unsigned int)n_rows-1) + if (idx.first==0 && idx.second==static_cast(n_rows-1)) return eigenpairs_symmetric_MRRR(compute_eigenvectors); else return eigenpairs_symmetric_MRRR(compute_eigenvectors,idx); @@ -934,10 +934,10 @@ template std::vector ScaLAPACKMatrix::eigenpairs_symmetric_by_value_MRRR(const std::pair &value_limits, const bool compute_eigenvectors) { - Assert (!std::isnan(value_limits.first),ExcMessage("value_limits.first is NaN")); - Assert (!std::isnan(value_limits.second),ExcMessage("value_limits.second is NaN")); + AssertIsFinite(value_limits.first); + AssertIsFinite(value_limits.second); - std::pair indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int); + const std::pair indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int); return eigenpairs_symmetric_MRRR(compute_eigenvectors,indices,value_limits); } @@ -960,7 +960,8 @@ ScaLAPACKMatrix::eigenpairs_symmetric_MRRR(const bool compute_eigenv const bool use_values = (std::isnan(eigenvalue_limits.first) || std::isnan(eigenvalue_limits.second)) ? false : true; const bool use_indices = ((eigenvalue_idx.first==numbers::invalid_unsigned_int) || (eigenvalue_idx.second==numbers::invalid_unsigned_int)) ? false : true; - Assert(!(use_values && use_indices),ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous")); + Assert(!(use_values && use_indices), + ExcMessage("Prescribing both the index and value range for the eigenvalues is ambiguous")); // If computation of eigenvectors is not required, use a sufficiently small distributed matrix. std::unique_ptr> eigenvectors = compute_eigenvectors ? @@ -983,9 +984,8 @@ ScaLAPACKMatrix::eigenpairs_symmetric_MRRR(const bool compute_eigenv * For jobz==N only eigenvalues are computed, for jobz='V' also the eigenvectors of the matrix are computed. */ char jobz = compute_eigenvectors ? 'V' : 'N'; - char range='A'; // Default value is to compute all eigenvalues and optionally eigenvectors. - bool all_eigenpairs=true; + char range='A'; NumberType vl=NumberType(),vu=NumberType(); int il=1,iu=1; @@ -996,12 +996,10 @@ ScaLAPACKMatrix::eigenpairs_symmetric_MRRR(const bool compute_eigenv if (!use_values) { range = 'A'; - all_eigenpairs = true; } else { range = 'V'; - all_eigenpairs = false; vl = std::min(eigenvalue_limits.first,eigenvalue_limits.second); vu = std::max(eigenvalue_limits.first,eigenvalue_limits.second); } @@ -1009,7 +1007,6 @@ ScaLAPACKMatrix::eigenpairs_symmetric_MRRR(const bool compute_eigenv else { range = 'I'; - all_eigenpairs = false; // As Fortran starts counting/indexing from 1 unlike C/C++, where it starts from 0. il = std::min(eigenvalue_idx.first,eigenvalue_idx.second) + 1; iu = std::max(eigenvalue_idx.first,eigenvalue_idx.second) + 1; diff --git a/tests/scalapack/scalapack_15.cc b/tests/scalapack/scalapack_15.cc index 3b1ef5ceb5..af60c275d7 100644 --- a/tests/scalapack/scalapack_15.cc +++ b/tests/scalapack/scalapack_15.cc @@ -132,7 +132,7 @@ void test(const unsigned int size, const unsigned int block_size, const NumberTy for (unsigned int i=0; i