]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Additional LAPACK routine for solving generalized eigenvalue problems added. The...
authorwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Feb 2012 13:19:49 +0000 (13:19 +0000)
committerwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Feb 2012 13:19:49 +0000 (13:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@25073 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/configure
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 6c2b77c6c65295b2a96b7bd09aae5028cb1a5b96..3c90705ccf74b461c3faf1b4c4795b50cdb6280c 100755 (executable)
@@ -1,5 +1,5 @@
 #! /bin/sh
-# From configure.in Revision: 24994 .
+# From configure.in Revision: 24995 .
 # Guess values for system-dependent variables and create Makefiles.
 # Generated by GNU Autoconf 2.68 for deal.II 7.2.pre.
 #
@@ -13502,7 +13502,7 @@ _ACEOF
 fi
 done
 
-for ac_func in dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_
+for ac_func in dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dsygvx_ ssygvx_
 do :
   as_ac_var=`$as_echo "ac_cv_func_$ac_func" | $as_tr_sh`
 ac_fn_cxx_check_func "$LINENO" "$ac_func" "$as_ac_var"
index 1a23b9c1d8c4b1e4acbc018c5f84576d13efdbb2..57b8ea697e35519c33f24d0ab24a8d8a270c4637 100644 (file)
@@ -670,7 +670,7 @@ dnl exception.
 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_])
+AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dsygvx_ ssygvx_])
 dnl Singular value decomposition
 AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgesdd_ sgesdd_ dgelsd_ sgelsd_])
 dnl Check Bessel functions in GNU libc
index 299c749fae0a29e1f4d44b20a7088b2b73184a37..458e567d5afc833192424986d34cbf6d3b7cf828 100644 (file)
 /* Define to 1 if you have the `dstev_' function. */
 #undef HAVE_DSTEV_
 
+/* Define to 1 if you have the `dsygvx_' function. */
+#undef HAVE_DSYGVX_
+
 /* Define to 1 if you have the `dsygv_' function. */
 #undef HAVE_DSYGV_
 
index b49ef6b77229fcea6bb16c9172909209a4f0e399..5fa0219a3d40de09facf71fe01d75da9a75fcfbf 100644 (file)
@@ -326,7 +326,7 @@ class LAPACKFullMatrix : public TransposeTable<number>
 
                                     /**
                                      * Compute generalized eigenvalues
-                                     * and (optionally) eigenvectors of
+                                     * and eigenvectors of
                                      * a real generalized symmetric
                                      * eigenproblem of the form
                                      * itype = 1: $Ax=\lambda B x$
@@ -335,17 +335,53 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      * where A is this matrix.
                                      * A and B are assumed to be symmetric,
                                      * and B has to be positive definite.
+                                     * 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, eigenvalues can
-                                     * be retrieved using the
-                                     * eigenvalue() function. The
-                                     * matrix itself will be
-                                     * LAPACKSupport::unusable after
-                                     * this operation. The number of
-                                     * computed eigenvectors is equal
-                                     * to eigenvectors.size()
+                                     * been called, all eigenvalues in
+                                     * (lower_bound, upper_bound] will be
+                                     * stored in eigenvalues and the
+                                     * corresponding eigenvectors will be stored
+                                     * in eigenvectors, whose dimension is set
+                                     * accordingly.
                                      *
-                                     * Note that the function does
+                                     * @note Calls the LAPACK
+                                     * function Xsygvx. For this to
+                                     * work, ./configure has to
+                                     * be told to use LAPACK.
+                                     */
+    void compute_generalized_eigenvalues_symmetric(
+                              LAPACKFullMatrix<number> & B,
+                              const number lower_bound,
+                              const number upper_bound,
+                              const number abs_accuracy,
+                              Vector<number> & eigenvalues,
+                              std::vector<Vector<number> > & eigenvectors,
+                              const int itype = 1);
+
+                                    /**
+                                     * Same as the other
+                                     * compute_generalized_eigenvalues_symmetric
+                                     * function except that all
+                                     * eigenvalues are computed
+                                     * and the tolerance is set
+                                     * automatically.
+                                     * Note that this function does
                                      * not return the computed
                                      * eigenvalues right away since
                                      * that involves copying data
@@ -362,7 +398,12 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      * eigenvalues and have a
                                      * separate function that returns
                                      * whatever eigenvalue is
-                                     * requested.
+                                     * requested. Eigenvalues can
+                                     * be retrieved using the
+                                     * eigenvalue() function.
+                                     * The number of computed 
+                                     * eigenvectors is equal
+                                     * to eigenvectors.size()
                                      *
                                      * @note Calls the LAPACK
                                      * function Xsygv. For this to
index 76e78e7216a5295526ddcebbf41e213874662faf..cffee7371ae5fcde0b935507030ab4165e86fb10 100644 (file)
@@ -70,13 +70,24 @@ 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);
-// General eigenvalues and eigenvectors of
+// 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
 void dsygv_ (const int* itype, const char* jobz, const char* uplo,
              const int* n, double* A, const int* lda, double* B,
              const int* ldb, double* w, double* work,
              const int* lwork, int* info);
+// Same functionality as dsygv_ 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 dsygvx_ (const int* itype, const char* jobz, const char* range,
+              const char* uplo, const int* n, double* A, const int* lda,
+              double* B, const int* ldb, 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);
 
 // Compute singular value decomposition using divide and conquer
 void dgesdd_ (const char* jobz,
index 6d17829fc4c26d46af519c389751b280ebb71d08..3fad9782026e210bb8ccba108564993b87793bec 100644 (file)
@@ -434,6 +434,117 @@ LAPACKFullMatrix<number>::compute_eigenvalues(
 }
 
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
+  LAPACKFullMatrix<number> & B,
+  const number lower_bound,
+  const number upper_bound,
+  const number abs_accuracy,
+  Vector<number> & eigenvalues,
+  std::vector<Vector<number> > & eigenvectors,
+  const int itype)
+{
+  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());
+  Assert(B.n_rows() == B.n_cols(), ExcNotQuadratic());
+  Assert(static_cast<unsigned int>(nn) == B.n_cols(), 
+        ExcDimensionMismatch (nn, B.n_cols()));
+
+  wr.resize(nn);
+  LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
+
+  number* values_A = const_cast<number*> (this->data());
+  number* values_B = const_cast<number*> (B.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 ?SYGVX 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);
+  
+  sygvx (&itype, jobz, range, uplo, &nn, values_A, &nn,
+         values_B, &nn, &lower_bound, &upper_bound,
+         dummy, dummy, &abs_accuracy, &n_eigenpairs,
+         &wr[0], values_eigenvectors, &nn, &work[0],
+         &lwork, &iwork[0], &ifail[0], &info);
+                                  // sygvx 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 generalized
+                                  // eigenvalues.
+  sygvx (&itype, jobz, range, uplo, &nn, values_A, &nn,
+         values_B, &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 sygv" << std::endl;
+
+  eigenvalues.reinit(n_eigenpairs);
+  eigenvectors.resize(n_eigenpairs);
+  for(unsigned int i=0; i < static_cast<unsigned int> (n_eigenpairs); ++i)
+  {
+    eigenvalues(i) = wr[i];
+    unsigned int col_begin(i*nn);
+    eigenvectors[i].reinit(nn, true);
+    for (unsigned int j=0; j < static_cast<unsigned int> (nn); ++j)
+    {
+      eigenvectors[i](j) = values_eigenvectors[col_begin+j];
+    }
+  }
+
+  state = LAPACKSupport::State(unusable);
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (

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.