From 05eeee5b2156d3c164281f5ef3591ae4ae1c1380 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Mon, 19 Mar 2018 15:10:22 +0100 Subject: [PATCH] add eigensolver routines using MRRR algorithm --- include/deal.II/lac/scalapack.h | 44 ++++++ include/deal.II/lac/scalapack.templates.h | 141 +++++++++++++++++ source/lac/scalapack.cc | 179 ++++++++++++++++++++++ 3 files changed, 364 insertions(+) diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 19dca92b6e..5a4aa99ac5 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -421,6 +421,34 @@ public: std::vector eigenpairs_symmetric_by_value(const std::pair &value_limits, const bool compute_eigenvectors); + /** + * Computing selected eigenvalues and, optionally, the eigenvectors of the real symmetric + * matrix $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. + * + * 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. + * + * If all eigenvalues/eigenvectors have to be computed, pass the closed interval $ \left[ 0, M-1 \right] $ in @p index_limits. + * + * Pass the closed interval $ \left[ M-r, M-1 \right] $ if the $r$ largest eigenvalues/eigenvectors are desired. + */ + std::vector eigenpairs_symmetric_by_index_MRRR(const std::pair &index_limits, + const bool compute_eigenvectors); + + /** + * Computing selected eigenvalues and, optionally, the eigenvectors 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. + * The eigenvectors are stored in the columns of the matrix, thereby + * overwriting the original content of the matrix. + */ + std::vector eigenpairs_symmetric_by_value_MRRR(const std::pair &value_limits, + const bool compute_eigenvectors); + /** * Computing the singular value decomposition (SVD) of a * matrix $A \in \mathbb{R}^{M \times N}$, optionally computing the left and/or right @@ -606,6 +634,22 @@ private: const std::pair &value_limits= std::make_pair(std::numeric_limits::quiet_NaN(),std::numeric_limits::quiet_NaN())); + /** + * Computing selected eigenvalues and, optionally, the eigenvectors 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. + * 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= + std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int), + const std::pair &value_limits= + std::make_pair(std::numeric_limits::quiet_NaN(),std::numeric_limits::quiet_NaN())); + /* * Stores the distributed matrix in @p filename * using serial routines diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index 824c7f1985..97dd75a5b2 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -743,6 +743,61 @@ extern "C" const int *IC, const int *JC, const int *DESCC); + + /** + * psyevr computes selected eigenvalues and, optionally, eigenvectors + * of a real symmetric matrix A using a parallel implementation of the MRR algorithm. + * Eigenvalues/vectors can be selected by specifying a range of values + * or a range of indices for the desired eigenvalues. + */ + void pdsyevr_(const char *jobz, + const char *range, + const char *uplo, + const int *n, + double *A, + const int *IA, + const int *JA, + const int *DESCA, + const double *VL, + const double *VU, + const int *IL, + const int *IU, + int *m, + int *nz, + double *w, + double *Z, + const int *IZ, + const int *JZ, + const int *DESCZ, + double *work, + int *lwork, + int *iwork, + int *liwork, + int *info); + void pssyevr_(const char *jobz, + const char *range, + const char *uplo, + const int *n, + float *A, + const int *IA, + const int *JA, + const int *DESCA, + const float *VL, + const float *VU, + const int *IL, + const int *IU, + int *m, + int *nz, + float *w, + float *Z, + const int *IZ, + const int *JZ, + const int *DESCZ, + float *work, + int *lwork, + int *iwork, + int *liwork, + int *info); } @@ -1652,6 +1707,92 @@ inline void ptran(const int *m, pstran_(m,n,alpha,A,IA,JA,DESCA,beta,C,IC,JC,DESCC); } + +template +inline void psyevr(const char * /*jobz*/, + const char * /*range*/, + const char * /*uplo*/, + const int * /*n*/, + number * /*A*/, + const int * /*IA*/, + const int * /*JA*/, + const int * /*DESCA*/, + const number * /*VL*/, + const number * /*VU*/, + const int * /*IL*/, + const int * /*IU*/, + int * /*m*/, + int * /*nz*/, + number * /*w*/, + number * /*Z*/, + const int * /*IZ*/, + const int * /*JZ*/, + const int * /*DESCZ*/, + number * /*work*/, + int * /*lwork*/, + int * /*iwork*/, + int * /*liwork*/, + int * /*info*/) +{ + Assert (false, dealii::ExcNotImplemented()); +} + +inline void psyevr(const char *jobz, + const char *range, + const char *uplo, + const int *n, + double *A, + const int *IA, + const int *JA, + const int *DESCA, + const double *VL, + const double *VU, + const int *IL, + const int *IU, + int *m, + int *nz, + double *w, + double *Z, + const int *IZ, + const int *JZ, + const int *DESCZ, + double *work, + int *lwork, + int *iwork, + int *liwork, + int *info) +{ + pdsyevr_(jobz,range,uplo,n,A,IA,JA,DESCA,VL,VU,IL,IU,m,nz,w,Z,IZ,JZ,DESCZ,work,lwork,iwork,liwork,info); +} + +inline void psyevr(const char *jobz, + const char *range, + const char *uplo, + const int *n, + float *A, + const int *IA, + const int *JA, + const int *DESCA, + const float *VL, + const float *VU, + const int *IL, + const int *IU, + int *m, + int *nz, + float *w, + float *Z, + const int *IZ, + const int *JZ, + const int *DESCZ, + float *work, + int *lwork, + int *iwork, + int *liwork, + int *info) +{ + pssyevr_(jobz,range,uplo,n,A,IA,JA,DESCA,VL,VU,IL,IU,m,nz,w,Z,IZ,JZ,DESCZ,work,lwork,iwork,liwork,info); +} + #endif // DEAL_II_WITH_SCALAPACK #endif // dealii_scalapack_templates_h diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 6ea55a6f96..6ab1b4439a 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -910,6 +910,185 @@ ScaLAPACKMatrix::eigenpairs_symmetric(const bool compute_eigenvector +template +std::vector ScaLAPACKMatrix::eigenpairs_symmetric_by_index_MRRR(const std::pair &index_limits, + 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)); + + 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) + return eigenpairs_symmetric_MRRR(compute_eigenvectors); + else + return eigenpairs_symmetric_MRRR(compute_eigenvectors,idx); +} + + + +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")); + + std::pair indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int); + + return eigenpairs_symmetric_MRRR(compute_eigenvectors,indices,value_limits); +} + + + +template +std::vector +ScaLAPACKMatrix::eigenpairs_symmetric_MRRR(const bool compute_eigenvectors, + const std::pair &eigenvalue_idx, + const std::pair &eigenvalue_limits) +{ + Assert (state == LAPACKSupport::matrix, + ExcMessage("Matrix has to be in Matrix state before calling this function.")); + Assert (property == LAPACKSupport::symmetric, + ExcMessage("Matrix has to be symmetric for this operation.")); + + Threads::Mutex::ScopedLock lock(mutex); + + 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")); + + // If computation of eigenvectors is not required, use a sufficiently small distributed matrix. + std::unique_ptr> eigenvectors = compute_eigenvectors ? + std_cxx14::make_unique>(n_rows,grid,row_block_size) : + std_cxx14::make_unique>(grid->n_process_rows,grid->n_process_columns,grid,1,1); + + eigenvectors->property = property; + // Number of eigenvalues to be returned from psyevr; upon successful exit ev contains the m seclected eigenvalues in ascending order. + int m = n_rows; + std::vector ev(n_rows); + + // Number of eigenvectors to be returned; + // Upon successful exit the first m=nz columns contain the selected eigenvectors (only if jobz=='V'). + int nz=0; + + if (grid->mpi_process_is_active) + { + int info = 0; + /* + * 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; + NumberType vl=NumberType(),vu=NumberType(); + int il=1,iu=1; + + // Index range for eigenvalues is not specified. + if (!use_indices) + { + // Interval for eigenvalues is not specified and consequently all eigenvalues/eigenpairs will be computed. + 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); + } + } + 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; + } + NumberType *A_loc = &this->values[0]; + + /* + * By setting lwork to -1 a workspace query for optimal length of work is performed. + */ + int lwork=-1; + int liwork=-1; + NumberType *eigenvectors_loc = (compute_eigenvectors ? &eigenvectors->values[0] : nullptr); + work.resize(1); + iwork.resize (1); + + psyevr(&jobz, &range, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, + &vl, &vu, &il, &iu, &m, &nz, ev.data(), + eigenvectors_loc, &eigenvectors->submatrix_row, &eigenvectors->submatrix_column, eigenvectors->descriptor, + work.data(), &lwork, iwork.data(), &liwork, &info); + + AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyevr", info)); + + lwork=work[0]; + work.resize(lwork); + liwork = iwork[0]; + iwork.resize(liwork); + + psyevr(&jobz, &range, &uplo, &n_rows, A_loc, &submatrix_row, &submatrix_column, descriptor, + &vl, &vu, &il, &iu, &m, &nz, ev.data(), + eigenvectors_loc, &eigenvectors->submatrix_row, &eigenvectors->submatrix_column, eigenvectors->descriptor, + work.data(), &lwork, iwork.data(), &liwork, &info); + + AssertThrow (info==0, LAPACKSupport::ExcErrorCode("psyevr", info)); + + if (compute_eigenvectors) + AssertThrow(m==nz,ExcMessage("psyevr failed to compute all eigenvectors for the selected eigenvalues")); + + // If eigenvectors are queried, copy eigenvectors to original matrix. + // As the temporary matrix eigenvectors has identical dimensions and + // block-cyclic distribution we simply swap the local array. + if (compute_eigenvectors) + this->values.swap(eigenvectors->values); + + // Adapt the size of ev to fit m upon return. + while ((int)ev.size() > m) + ev.pop_back(); + } + /* + * Send number of computed eigenvalues to inactive processes. + */ + grid->send_to_inactive(&m, 1); + + /* + * Inactive processes have to resize array of eigenvalues. + */ + if (! grid->mpi_process_is_active) + ev.resize(m); + /* + * Send the eigenvalues to processors not being part of the process grid. + */ + grid->send_to_inactive(ev.data(), ev.size()); + + /* + * If only eigenvalues are queried, the content of the matrix will be destroyed. + * If the eigenpairs are queried, matrix A on exit stores the eigenvectors in the columns. + */ + if (compute_eigenvectors) + { + property = LAPACKSupport::Property::general; + state = LAPACKSupport::eigenvalues; + } + else + state = LAPACKSupport::unusable; + + return ev; +} + + + template std::vector ScaLAPACKMatrix::compute_SVD(ScaLAPACKMatrix *U, ScaLAPACKMatrix *VT) -- 2.39.5