]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
eigenvalue solver for symmetric LAPACK matrices included. This eigensolver in particu...
authorwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Mar 2012 13:55:09 +0000 (13:55 +0000)
committerwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 19 Mar 2012 13:55:09 +0000 (13:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@25300 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/configure.in
deal.II/include/deal.II/base/config.h.in
deal.II/include/deal.II/lac/lapack_full_matrix.h
deal.II/include/deal.II/lac/lapack_templates.h.in
deal.II/source/lac/lapack_full_matrix.cc

index 57b8ea697e35519c33f24d0ab24a8d8a270c4637..6db52e46fe13ab9c932bfe9f71093ac96975409d 100644 (file)
@@ -671,6 +671,7 @@ AC_CHECK_FUNCS([daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_])
 AC_CHECK_FUNCS([dgemm_ sgemm_ dgetrf_ sgetrf_ dgetri_ sgetri_])
 AC_CHECK_FUNCS([dgeqrf_ sgeqrf_ dormqr_ sormqr_ dorgqr_ sorgqr_ dtrtrs_ strtrs_])
 AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dsygvx_ ssygvx_])
+AC_CHECK_FUNCS([dsyevx_ ssyevx_])
 dnl Singular value decomposition
 AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgesdd_ sgesdd_ dgelsd_ sgelsd_])
 dnl Check Bessel functions in GNU libc
index 66621f78e37251fc263bf93d79fea6370dce2c05..379ac491d5bcd058a9e0a79df3a7d3b42df80a3b 100644 (file)
    provided by the compiler. */
 #undef DEAL_II_USE_DIRECT_ERRNO_H
 
-/* Defined if deal.II is configured with an external Boost library */
-#undef DEAL_II_USE_EXTERNAL_BOOST
-
 /* Defined if a Metis installation was found and is going to be used */
 #undef DEAL_II_USE_METIS
 
 /* Define to 1 if you have the `dstev_' function. */
 #undef HAVE_DSTEV_
 
+/* Define to 1 if you have the `dsyevx_' function. */
+#undef HAVE_DSYEVX_
+
 /* Define to 1 if you have the `dsygvx_' function. */
 #undef HAVE_DSYGVX_
 
 /* Define to 1 if you have the `sstev_' function. */
 #undef HAVE_SSTEV_
 
+/* Define to 1 if you have the `ssyevx_' function. */
+#undef HAVE_SSYEVX_
+
 /* Define to 1 if you have the `ssygvx_' function. */
 #undef HAVE_SSYGVX_
 
     >=  \
     (major)*10000 + (minor)*100 + (subminor))
 
+#define DEAL_II_PETSC_VERSION_DEV() \
+  (DEAL_II_USE_PETSC_DEV)
+
 #include <deal.II/base/numbers.h>
 #include <deal.II/base/types.h>
 
index 5fa0219a3d40de09facf71fe01d75da9a75fcfbf..9d14579fe7fd38f51d75e5411668bc3bc31e2e08 100644 (file)
@@ -323,6 +323,45 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      */
     void compute_eigenvalues (const bool right_eigenvectors = false,
                              const bool left_eigenvectors = false);
+    
+                                    /**
+                                     * Compute eigenvalues and
+                                     * eigenvectors of a real symmetric
+                                     * matrix. Only eigenvalues in the
+                                     * interval (lower_bound, upper_bound]
+                                     * are computed with the absolute
+                                     * tolerance abs_accuracy. An approximate
+                                     * eigenvalue is accepted as converged
+                                     * when it is determined to lie in an
+                                     * interval [a,b] of width less than or
+                                     * equal to abs_accuracy + eps * max( |a|,|b| ),
+                                      * where eps is the machine precision.
+                                      * If abs_accuracy is less than
+                                      * or equal to zero, then  eps*|t| will 
+                                      * be used in its place, where |t| is the 
+                                      * 1-norm of the tridiagonal matrix obtained
+                                      * by reducing A to tridiagonal form.
+                                      * Eigenvalues will be computed most accurately
+                                      * when abs_accuracy is set to twice the 
+                                      * underflow threshold, not zero.
+                                     * After this routine has
+                                     * been called, all eigenvalues in
+                                     * (lower_bound, upper_bound] will be
+                                     * stored in eigenvalues and the
+                                     * corresponding eigenvectors will be stored
+                                     * in the columns of eigenvectors, whose 
+                                     * dimension is set accordingly.
+                                     * 
+                                     * @note Calls the LAPACK function
+                                     * Xsyevx. For this to work, ./configure
+                                     * has to be told to use LAPACK.
+                                     */ 
+    void compute_eigenvalues_symmetric(
+                             const number lower_bound,
+                             const number upper_bound,
+                             const number abs_accuracy,
+                             Vector<number> & eigenvalues,
+                             FullMatrix<number> & eigenvectors);
 
                                     /**
                                      * Compute generalized eigenvalues
index cffee7371ae5fcde0b935507030ab4165e86fb10..7eb37440cc547b6fcf396d05183d2218bd065451 100644 (file)
@@ -70,6 +70,17 @@ void dgeevx_ (const char* balanc, const char* jobvl, const char* jobvr,
 void dsyev_ (const char *jobz, const char *uplo, const int *n,
             double *A, const int *lda, double *w,
             double *work, const int *lwork, int *info);
+// Same functionality as dsyev_ but with more options: E.g.
+// Compute only eigenvalues in a specific interval,
+// Compute only eigenvalues with a specific index,
+// Set tolerance for eigenvalue computation
+void dsyevx_ (const char* jobz, const char* range,
+              const char* uplo, const int* n, double* A, const int* lda,
+              const double* vl, const double* vu,
+              const int* il, const int* iu, const double* abstol,
+              int* m, double* w, double* z,
+              const int* ldz, double* work, const int* lwork, int* iwork,
+              int* ifail, int* info);
 // Generalized eigenvalues and eigenvectors of
 // 1: A*x = lambda*B*x; 2: A*B*x = lambda*x; 3: B*A*x = lambda*x
 // A and B are symmetric and B is definite
index 3fad9782026e210bb8ccba108564993b87793bec..aacd5baedaa02eb67d720bb817d035ffe29c7ddf 100644 (file)
@@ -434,6 +434,113 @@ LAPACKFullMatrix<number>::compute_eigenvalues(
 }
 
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
+  const number lower_bound,
+  const number upper_bound,
+  const number abs_accuracy,
+  Vector<number> & eigenvalues,
+  FullMatrix<number> & eigenvectors)
+{
+  Assert(state == matrix, ExcState(state));
+  const int nn = (this->n_cols() > 0 ? this->n_cols() : 1);
+  Assert(static_cast<unsigned int>(nn) == this->n_rows(), ExcNotQuadratic());
+
+  wr.resize(nn);
+  LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
+
+  number* values_A = const_cast<number*> (this->data());
+  number* values_eigenvectors = const_cast<number*> (matrix_eigenvectors.data());
+  
+  int info(0),
+      lwork(1),
+      n_eigenpairs(0);
+  const char * const jobz(&V);
+  const char * const uplo(&U);
+  const char * const range(&V);
+  const int * const  dummy(&one);
+  std::vector<int> iwork(static_cast<unsigned int> (5*nn));
+  std::vector<int> ifail(static_cast<unsigned int> (nn));
+  
+  
+                                  // Optimal workspace query:
+
+                                  // The LAPACK routine ?SYEVX requires
+                                  // a sufficient large workspace variable,
+                                  // minimum requirement is
+                                  //    work.size>=3*nn-1.
+                                  // However, to improve performance, a
+                                  // somewhat larger workspace may be needed.
+  
+                                  // SOME implementations of the LAPACK routine
+                                  // provide a workspace query call,
+                                  //   info:=0, lwork:=-1
+                                  // which returns an optimal value for the
+                                  // size of the workspace array
+                                  // (the PETSc 2.3.0 implementation does NOT
+                                  // provide this functionality!).
+
+                                  // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
+                                  // disable the workspace query.
+#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
+  lwork = -1;
+  work.resize(1);
+  
+  syevx (jobz, range,
+        uplo, &nn, values_A, &nn,
+        &lower_bound, &upper_bound,
+         dummy, dummy, &abs_accuracy,
+        &n_eigenpairs, &wr[0], values_eigenvectors,
+        &nn, &work[0], &lwork, &iwork[0],
+        &ifail[0], &info);
+                                  // syevx returns info=0 on
+                                  // success. Since we only queried
+                                  // the optimal size for work,
+                                  // everything else would not be
+                                  // acceptable.
+  Assert (info == 0, ExcInternalError());
+                                  // Allocate working array according
+                                  // to suggestion.
+  lwork = (int) (work[0]+.1);
+#else
+  lwork = 8*nn > 1 ? 8*nn : 1; // no query mode
+#endif
+                                  // resize workspace arrays
+  work.resize(static_cast<unsigned int> (lwork));
+
+                                  // Finally compute the eigenvalues.
+  syevx (jobz, range,
+        uplo, &nn, values_A, &nn,
+        &lower_bound, &upper_bound,
+         dummy, dummy, &abs_accuracy,
+        &n_eigenpairs, &wr[0], values_eigenvectors,
+        &nn, &work[0], &lwork, &iwork[0],
+        &ifail[0], &info);
+
+                                  // Negative return value implies a
+                                  // wrong argument. This should be
+                                  // internal.
+  Assert (info >=0, ExcInternalError());
+  if (info != 0)
+    std::cerr << "LAPACK error in syevx" << std::endl;
+
+  eigenvalues.reinit(n_eigenpairs);
+  eigenvectors.reinit(nn, n_eigenpairs, true);
+  for(unsigned int i=0; i < static_cast<unsigned int> (n_eigenpairs); ++i)
+  {
+    eigenvalues(i) = wr[i];
+    unsigned int col_begin(i*nn);
+    for (unsigned int j=0; j < static_cast<unsigned int> (nn); ++j)
+    {
+      eigenvectors(j,i) = values_eigenvectors[col_begin+j];
+    }
+  }
+
+  state = LAPACKSupport::State(unusable);
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
@@ -526,7 +633,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
                                   // internal.
   Assert (info >=0, ExcInternalError());
   if (info != 0)
-    std::cerr << "LAPACK error in sygv" << std::endl;
+    std::cerr << "LAPACK error in sygvx" << std::endl;
 
   eigenvalues.reinit(n_eigenpairs);
   eigenvectors.resize(n_eigenpairs);

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.