From 15ae836802926694b5dc47b056c47d43a1b73c4b Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 29 May 2018 16:47:05 -0400 Subject: [PATCH] Fix LAPACKFullMatrix>. --- include/deal.II/lac/lapack_full_matrix.h | 10 +- source/lac/lapack_full_matrix.cc | 364 ++++++++++++++++++----- source/lac/lapack_full_matrix.inst.in | 2 +- 3 files changed, 299 insertions(+), 77 deletions(-) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index da79619726..f1c93e3416 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -912,10 +912,11 @@ private: * Real parts of eigenvalues or the singular values. Filled by * compute_eigenvalues() or compute_svd(). */ - std::vector wr; + std::vector::real_type> wr; /** - * Imaginary parts of eigenvalues. Filled by compute_eigenvalues. + * Imaginary parts of eigenvalues, or, in the complex scalar case, the + * eigenvalues themselves. Filled by compute_eigenvalues. */ std::vector wi; @@ -1121,7 +1122,10 @@ LAPACKFullMatrix::eigenvalue(const size_type i) const Assert(wi.size() == this->n_rows(), ExcInternalError()); AssertIndexRange(i, this->n_rows()); - return std::complex(wr[i], wi[i]); + if (numbers::NumberTraits::is_complex) + return std::complex(wi[i]); + else + return std::complex(wr[i], wi[i]); } diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 74ec0da5de..552fc261c4 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -13,6 +13,7 @@ // // --------------------------------------------------------------------- +#include #include #include @@ -30,6 +31,214 @@ DEAL_II_NAMESPACE_OPEN using namespace LAPACKSupport; +namespace internal +{ + namespace LAPACKFullMatrixImplementation + { + // ZGEEV/CGEEV and DGEEV/SGEEV need different work arrays and different + // output arrays for eigenvalues. This makes working with generic scalar + // types a bit difficult. To get around this, geev_helper has the same + // signature for real and complex arguments, but it ignores some + // parameters when called with a real type and ignores different + // parameters when called with a complex type. + template + void + geev_helper(const char vl, + const char vr, + AlignedVector & matrix, + const types::blas_int n_rows, + std::vector & real_part_eigenvalues, + std::vector & imag_part_eigenvalues, + std::vector & left_eigenvectors, + std::vector & right_eigenvectors, + std::vector & real_work, + std::vector & /*complex_work*/, + const types::blas_int work_flag, + types::blas_int & info) + { + static_assert(std::is_same::value || + std::is_same::value, + "Only implemented for double and float"); + Assert(matrix.size() == static_cast(n_rows * n_rows), + ExcInternalError()); + Assert(static_cast(n_rows) <= real_part_eigenvalues.size(), + ExcInternalError()); + Assert(static_cast(n_rows) <= imag_part_eigenvalues.size(), + ExcInternalError()); + if (vl == 'V') + Assert(static_cast(n_rows * n_rows) <= + left_eigenvectors.size(), + ExcInternalError()); + if (vr == 'V') + Assert(static_cast(n_rows * n_rows) <= + right_eigenvectors.size(), + ExcInternalError()); + Assert(work_flag == -1 || + static_cast(2 * n_rows) <= real_work.size(), + ExcInternalError()); + Assert(work_flag == -1 || std::max(1, 3 * n_rows) <= work_flag, + ExcInternalError()); + geev(&vl, + &vr, + &n_rows, + &matrix[0], + &n_rows, + real_part_eigenvalues.data(), + imag_part_eigenvalues.data(), + left_eigenvectors.data(), + &n_rows, + right_eigenvectors.data(), + &n_rows, + real_work.data(), + &work_flag, + &info); + } + + template + void + geev_helper(const char vl, + const char vr, + AlignedVector> &matrix, + const types::blas_int n_rows, + std::vector & /*real_part_eigenvalues*/, + std::vector> &eigenvalues, + std::vector> &left_eigenvectors, + std::vector> &right_eigenvectors, + std::vector> &complex_work, + std::vector & real_work, + const types::blas_int work_flag, + types::blas_int & info) + { + static_assert( + std::is_same::value || std::is_same::value, + "Only implemented for std::complex and std::complex"); + Assert(matrix.size() == static_cast(n_rows * n_rows), + ExcInternalError()); + Assert(static_cast(n_rows) <= eigenvalues.size(), + ExcInternalError()); + if (vl == 'V') + Assert(static_cast(n_rows * n_rows) <= + left_eigenvectors.size(), + ExcInternalError()); + if (vr == 'V') + Assert(static_cast(n_rows * n_rows) <= + right_eigenvectors.size(), + ExcInternalError()); + Assert(std::max(1, work_flag) <= real_work.size(), + ExcInternalError()); + Assert(work_flag == -1 || + std::max(1, 2 * n_rows) <= (work_flag), + ExcInternalError()); + + geev(&vl, + &vr, + &n_rows, + &matrix[0], + &n_rows, + eigenvalues.data(), + left_eigenvectors.data(), + &n_rows, + right_eigenvectors.data(), + &n_rows, + complex_work.data(), + &work_flag, + real_work.data(), + &info); + } + + + + template + void + gesdd_helper(const char job, + const types::blas_int n_rows, + const types::blas_int n_cols, + AlignedVector & matrix, + std::vector & singular_values, + AlignedVector & left_vectors, + AlignedVector & right_vectors, + std::vector & real_work, + std::vector & /*complex work*/, + std::vector &integer_work, + const types::blas_int & work_flag, + types::blas_int & info) + { + Assert(job == 'A' || job == 'S' || job == 'O' || job == 'N', + ExcInternalError()); + Assert(static_cast(n_rows * n_cols) == matrix.size(), + ExcInternalError()); + Assert(std::min(n_rows, n_cols) <= singular_values.size(), + ExcInternalError()); + Assert(8 * std::min(n_rows, n_cols) <= integer_work.size(), + ExcInternalError()); + Assert(work_flag == -1 || + static_cast(work_flag) <= real_work.size(), + ExcInternalError()); + gesdd(&job, + &n_rows, + &n_cols, + &matrix[0], + &n_rows, + singular_values.data(), + &left_vectors[0], + &n_rows, + &right_vectors[0], + &n_cols, + real_work.data(), + &work_flag, + integer_work.data(), + &info); + } + + + + template + void + gesdd_helper(const char job, + const types::blas_int n_rows, + const types::blas_int n_cols, + AlignedVector> &matrix, + std::vector & singular_values, + AlignedVector> &left_vectors, + AlignedVector> &right_vectors, + std::vector> & work, + std::vector & real_work, + std::vector & integer_work, + const types::blas_int & work_flag, + types::blas_int & info) + { + Assert(job == 'A' || job == 'S' || job == 'O' || job == 'N', + ExcInternalError()); + Assert(static_cast(n_rows * n_cols) == matrix.size(), + ExcInternalError()); + Assert(static_cast(std::min(n_rows, n_cols)) <= + singular_values.size(), + ExcInternalError()); + Assert(8 * std::min(n_rows, n_cols) <= integer_work.size(), + ExcInternalError()); + Assert(work_flag == -1 || + static_cast(work_flag) <= real_work.size(), + ExcInternalError()); + + gesdd(&job, + &n_rows, + &n_cols, + &matrix[0], + &n_rows, + singular_values.data(), + &left_vectors[0], + &n_rows, + &right_vectors[0], + &n_cols, + work.data(), + &work_flag, + real_work.data(), + integer_work.data(), + &info); + } + } // namespace LAPACKFullMatrixImplementation +} // namespace internal + template LAPACKFullMatrix::LAPACKFullMatrix(const size_type n) : TransposeTable(n, n), @@ -1313,36 +1522,47 @@ LAPACKFullMatrix::compute_svd() Assert(state == matrix, ExcState(state)); state = LAPACKSupport::unusable; - const types::blas_int mm = this->m(); - const types::blas_int nn = this->n(); - number *const values = &this->values[0]; + const types::blas_int mm = this->m(); + const types::blas_int nn = this->n(); wr.resize(std::max(mm, nn)); std::fill(wr.begin(), wr.end(), 0.); ipiv.resize(8 * mm); - svd_u = std_cxx14::make_unique>(mm, mm); - svd_vt = std_cxx14::make_unique>(nn, nn); - number *const mu = &svd_u->values[0]; - number *const mvt = &svd_vt->values[0]; + svd_u = std_cxx14::make_unique>(mm, mm); + svd_vt = std_cxx14::make_unique>(nn, nn); types::blas_int info = 0; // First determine optimal workspace size work.resize(1); types::blas_int lwork = -1; - gesdd(&LAPACKSupport::A, - &mm, - &nn, - values, - &mm, - wr.data(), - mu, - &mm, - mvt, - &nn, - work.data(), - &lwork, - ipiv.data(), - &info); + + // TODO double check size + std::vector::real_type> real_work; + if (numbers::NumberTraits::is_complex) + { + // This array is only used by the complex versions. + std::size_t min = std::min(this->m(), this->n()); + std::size_t max = std::max(this->m(), this->n()); + real_work.resize( + std::max(5 * min * min + 5 * min, 2 * max * min + 2 * min * min + min)); + } + + // make sure that the first entry in the work array is clear, in case the + // routine does not completely overwrite the memory: + work[0] = number(); + internal::LAPACKFullMatrixImplementation::gesdd_helper(LAPACKSupport::A, + mm, + nn, + this->values, + wr, + svd_u->values, + svd_vt->values, + work, + real_work, + ipiv, + lwork, + info); + 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. @@ -1350,20 +1570,18 @@ LAPACKFullMatrix::compute_svd() work.resize(lwork); // Do the actual SVD. - gesdd(&LAPACKSupport::A, - &mm, - &nn, - values, - &mm, - wr.data(), - mu, - &mm, - mvt, - &nn, - work.data(), - &lwork, - ipiv.data(), - &info); + internal::LAPACKFullMatrixImplementation::gesdd_helper(LAPACKSupport::A, + mm, + nn, + this->values, + wr, + svd_u->values, + svd_vt->values, + work, + real_work, + ipiv, + lwork, + info); AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("gesdd", info)); work.resize(0); @@ -1382,11 +1600,12 @@ LAPACKFullMatrix::compute_inverse_svd(const double threshold) Assert(state == LAPACKSupport::svd, ExcState(state)); + const typename numbers::NumberTraits::real_type one(1.0); const double lim = std::abs(wr[0]) * threshold; for (size_type i = 0; i < wr.size(); ++i) { if (std::abs(wr[i]) > lim) - wr[i] = number(1.) / wr[i]; + wr[i] = one / wr[i]; else wr[i] = 0.; } @@ -1405,9 +1624,10 @@ LAPACKFullMatrix::compute_inverse_svd_with_kernel( Assert(state == LAPACKSupport::svd, ExcState(state)); - const unsigned int n_wr = wr.size(); + const typename numbers::NumberTraits::real_type one(1.0); + const unsigned int n_wr = wr.size(); for (size_type i = 0; i < n_wr - kernel_size; ++i) - wr[i] = number(1.) / wr[i]; + wr[i] = one / wr[i]; for (size_type i = n_wr - kernel_size; i < n_wr; ++i) wr[i] = 0.; state = LAPACKSupport::inverse_svd; @@ -1620,12 +1840,10 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) if (left) vl.resize(nn * nn); - number *const values = &this->values[0]; - - types::blas_int info = 0; - types::blas_int lwork = 1; - const char *const jobvr = (right) ? (&V) : (&N); - const char *const jobvl = (left) ? (&V) : (&N); + types::blas_int info = 0; + types::blas_int lwork = 1; + const char jobvr = (right) ? V : N; + const char jobvl = (left) ? V : N; /* * The LAPACK routine xGEEV requires a sufficiently large work array; the @@ -1639,20 +1857,23 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) lwork = -1; work.resize(1); - geev(jobvl, - jobvr, - &nn, - values, - &nn, - wr.data(), - wi.data(), - vl.data(), - &nn, - vr.data(), - &nn, - work.data(), - &lwork, - &info); + std::vector::real_type> real_work; + if (numbers::NumberTraits::is_complex) + // This array is only used by the complex versions. + real_work.resize(2 * this->m()); + internal::LAPACKFullMatrixImplementation::geev_helper(jobvl, + jobvr, + this->values, + this->m(), + wr, + wi, + vl, + vr, + work, + real_work, + lwork, + info); + // geev returns info=0 on success. Since we only queried the optimal size // for work, everything else would not be acceptable. Assert(info == 0, ExcInternalError()); @@ -1664,21 +1885,18 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, const bool left) work.resize((size_type)lwork); // Finally compute the eigenvalues. - geev(jobvl, - jobvr, - &nn, - values, - &nn, - wr.data(), - wi.data(), - vl.data(), - &nn, - vr.data(), - &nn, - work.data(), - &lwork, - &info); - // Negative return value implies a wrong argument. This should be internal. + internal::LAPACKFullMatrixImplementation::geev_helper(jobvl, + jobvr, + this->values, + this->m(), + wr, + wi, + vl, + vr, + work, + real_work, + lwork, + info); Assert(info >= 0, ExcInternalError()); // TODO:[GK] What if the QR method fails? diff --git a/source/lac/lapack_full_matrix.inst.in b/source/lac/lapack_full_matrix.inst.in index dea62ebbaf..71264b6dfb 100644 --- a/source/lac/lapack_full_matrix.inst.in +++ b/source/lac/lapack_full_matrix.inst.in @@ -15,7 +15,7 @@ -for (S : REAL_SCALARS) +for (S : REAL_AND_COMPLEX_SCALARS) { template class LAPACKFullMatrix; template class PreconditionLU; -- 2.39.5