From 06b7c1e66110e8e8d2dab4d9621bf6a81d7b2988 Mon Sep 17 00:00:00 2001 From: willems Date: Mon, 19 Mar 2012 13:55:09 +0000 Subject: [PATCH] eigenvalue solver for symmetric LAPACK matrices included. This eigensolver in particular includes the option of specifying the accuracy used for computing eigenvalues. Additionally, an interval in which eigenvalues are computed can be specified git-svn-id: https://svn.dealii.org/trunk@25300 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/configure.in | 1 + deal.II/include/deal.II/base/config.h.in | 12 +- .../include/deal.II/lac/lapack_full_matrix.h | 39 +++++++ .../include/deal.II/lac/lapack_templates.h.in | 11 ++ deal.II/source/lac/lapack_full_matrix.cc | 109 +++++++++++++++++- 5 files changed, 168 insertions(+), 4 deletions(-) diff --git a/deal.II/configure.in b/deal.II/configure.in index 57b8ea697e..6db52e46fe 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -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 diff --git a/deal.II/include/deal.II/base/config.h.in b/deal.II/include/deal.II/base/config.h.in index 66621f78e3..379ac491d5 100644 --- a/deal.II/include/deal.II/base/config.h.in +++ b/deal.II/include/deal.II/base/config.h.in @@ -244,9 +244,6 @@ 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 @@ -354,6 +351,9 @@ /* 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_ @@ -491,6 +491,9 @@ /* 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_ @@ -668,6 +671,9 @@ >= \ (major)*10000 + (minor)*100 + (subminor)) +#define DEAL_II_PETSC_VERSION_DEV() \ + (DEAL_II_USE_PETSC_DEV) + #include #include diff --git a/deal.II/include/deal.II/lac/lapack_full_matrix.h b/deal.II/include/deal.II/lac/lapack_full_matrix.h index 5fa0219a3d..9d14579fe7 100644 --- a/deal.II/include/deal.II/lac/lapack_full_matrix.h +++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h @@ -323,6 +323,45 @@ class LAPACKFullMatrix : public TransposeTable */ 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 & eigenvalues, + FullMatrix & eigenvectors); /** * Compute generalized eigenvalues diff --git a/deal.II/include/deal.II/lac/lapack_templates.h.in b/deal.II/include/deal.II/lac/lapack_templates.h.in index cffee7371a..7eb37440cc 100644 --- a/deal.II/include/deal.II/lac/lapack_templates.h.in +++ b/deal.II/include/deal.II/lac/lapack_templates.h.in @@ -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 diff --git a/deal.II/source/lac/lapack_full_matrix.cc b/deal.II/source/lac/lapack_full_matrix.cc index 3fad978202..aacd5baeda 100644 --- a/deal.II/source/lac/lapack_full_matrix.cc +++ b/deal.II/source/lac/lapack_full_matrix.cc @@ -434,6 +434,113 @@ LAPACKFullMatrix::compute_eigenvalues( } +template +void +LAPACKFullMatrix::compute_eigenvalues_symmetric( + const number lower_bound, + const number upper_bound, + const number abs_accuracy, + Vector & eigenvalues, + FullMatrix & eigenvectors) +{ + Assert(state == matrix, ExcState(state)); + const int nn = (this->n_cols() > 0 ? this->n_cols() : 1); + Assert(static_cast(nn) == this->n_rows(), ExcNotQuadratic()); + + wr.resize(nn); + LAPACKFullMatrix matrix_eigenvectors(nn, nn); + + number* values_A = const_cast (this->data()); + number* values_eigenvectors = const_cast (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 iwork(static_cast (5*nn)); + std::vector ifail(static_cast (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 (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 (n_eigenpairs); ++i) + { + eigenvalues(i) = wr[i]; + unsigned int col_begin(i*nn); + for (unsigned int j=0; j < static_cast (nn); ++j) + { + eigenvectors(j,i) = values_eigenvectors[col_begin+j]; + } + } + + state = LAPACKSupport::State(unusable); +} + + template void LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( @@ -526,7 +633,7 @@ LAPACKFullMatrix::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); -- 2.39.5