]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add eigensolver routines using MRRR algorithm
authorBenjamin Brands <benjamin.brands@fau.de>
Mon, 19 Mar 2018 14:10:22 +0000 (15:10 +0100)
committerBenjamin Brands <benjamin.brands@fau.de>
Mon, 19 Mar 2018 14:10:22 +0000 (15:10 +0100)
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc

index 19dca92b6e52bc74cd64e5f4e81c9610c4344d30..5a4aa99ac545fc51b1ea9abcec965eb612240fb3 100644 (file)
@@ -421,6 +421,34 @@ public:
   std::vector<NumberType> eigenpairs_symmetric_by_value(const std::pair<NumberType,NumberType> &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<NumberType> eigenpairs_symmetric_by_index_MRRR(const std::pair<unsigned int,unsigned int> &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<NumberType> eigenpairs_symmetric_by_value_MRRR(const std::pair<NumberType,NumberType> &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<NumberType,NumberType> &value_limits=
                                                  std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::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<NumberType> eigenpairs_symmetric_MRRR(const bool compute_eigenvectors,
+                                                    const std::pair<unsigned int,unsigned int> &index_limits=
+                                                      std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int),
+                                                    const std::pair<NumberType,NumberType> &value_limits=
+                                                      std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),std::numeric_limits<NumberType>::quiet_NaN()));
+
   /*
    * Stores the distributed matrix in @p filename
    * using serial routines
index 824c7f198526e772aed053f12058fe5c2a815885..97dd75a5b222de878bf3bf5de4a12318ace8c198 100644 (file)
@@ -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 <typename number>
+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
index 6ea55a6f9646d746405f7c44355db62f36765d31..6ab1b4439ad294f2e5ddb37f681b5045a567f2d1 100644 (file)
@@ -910,6 +910,185 @@ ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric(const bool compute_eigenvector
 
 
 
+template <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_index_MRRR(const std::pair<unsigned int,unsigned int> &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<unsigned int,unsigned int> 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 <typename NumberType>
+std::vector<NumberType> ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_by_value_MRRR(const std::pair<NumberType,NumberType> &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<unsigned int,unsigned int> indices = std::make_pair(numbers::invalid_unsigned_int,numbers::invalid_unsigned_int);
+
+  return eigenpairs_symmetric_MRRR(compute_eigenvectors,indices,value_limits);
+}
+
+
+
+template <typename NumberType>
+std::vector<NumberType>
+ScaLAPACKMatrix<NumberType>::eigenpairs_symmetric_MRRR(const bool compute_eigenvectors,
+                                                       const std::pair<unsigned int, unsigned int> &eigenvalue_idx,
+                                                       const std::pair<NumberType,NumberType> &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<ScaLAPACKMatrix<NumberType>> eigenvectors = compute_eigenvectors ?
+                                                              std_cxx14::make_unique<ScaLAPACKMatrix<NumberType>>(n_rows,grid,row_block_size) :
+                                                              std_cxx14::make_unique<ScaLAPACKMatrix<NumberType>>(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<NumberType> 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 <typename NumberType>
 std::vector<NumberType> ScaLAPACKMatrix<NumberType>::compute_SVD(ScaLAPACKMatrix<NumberType> *U,
     ScaLAPACKMatrix<NumberType> *VT)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.