From: willems Date: Tue, 14 Feb 2012 13:19:49 +0000 (+0000) Subject: Additional LAPACK routine for solving generalized eigenvalue problems added. The... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bd0d0a623e17a16017f7822f3473a487142a2479;p=dealii-svn.git Additional LAPACK routine for solving generalized eigenvalue problems added. The new function compute_generalized_eigenvalues_symmetric provides more functionality than then one with the same name already included in the library. In particular one can specify an interval where eigenvalues should be computed. Also one can provide the accuracy up to which eigenvalues are computed. git-svn-id: https://svn.dealii.org/trunk@25073 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/configure b/deal.II/configure index 6c2b77c6c6..3c90705ccf 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -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" diff --git a/deal.II/configure.in b/deal.II/configure.in index 1a23b9c1d8..57b8ea697e 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -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 diff --git a/deal.II/include/deal.II/base/config.h.in b/deal.II/include/deal.II/base/config.h.in index 299c749fae..458e567d5a 100644 --- a/deal.II/include/deal.II/base/config.h.in +++ b/deal.II/include/deal.II/base/config.h.in @@ -351,6 +351,9 @@ /* 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_ 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 b49ef6b772..5fa0219a3d 100644 --- a/deal.II/include/deal.II/lac/lapack_full_matrix.h +++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h @@ -326,7 +326,7 @@ class LAPACKFullMatrix : public TransposeTable /** * 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 * 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 & B, + const number lower_bound, + const number upper_bound, + const number abs_accuracy, + Vector & eigenvalues, + std::vector > & 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 * 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 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 76e78e7216..cffee7371a 100644 --- a/deal.II/include/deal.II/lac/lapack_templates.h.in +++ b/deal.II/include/deal.II/lac/lapack_templates.h.in @@ -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, diff --git a/deal.II/source/lac/lapack_full_matrix.cc b/deal.II/source/lac/lapack_full_matrix.cc index 6d17829fc4..3fad978202 100644 --- a/deal.II/source/lac/lapack_full_matrix.cc +++ b/deal.II/source/lac/lapack_full_matrix.cc @@ -434,6 +434,117 @@ LAPACKFullMatrix::compute_eigenvalues( } +template +void +LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( + LAPACKFullMatrix & B, + const number lower_bound, + const number upper_bound, + const number abs_accuracy, + Vector & eigenvalues, + std::vector > & eigenvectors, + const int itype) +{ + Assert(state == matrix, ExcState(state)); + const int nn = (this->n_cols() > 0 ? this->n_cols() : 1); + Assert(static_cast(nn) == this->n_rows(), ExcNotQuadratic()); + Assert(B.n_rows() == B.n_cols(), ExcNotQuadratic()); + Assert(static_cast(nn) == B.n_cols(), + ExcDimensionMismatch (nn, B.n_cols())); + + wr.resize(nn); + LAPACKFullMatrix matrix_eigenvectors(nn, nn); + + number* values_A = const_cast (this->data()); + number* values_B = const_cast (B.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 ?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 (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 (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 (nn); ++j) + { + eigenvectors[i](j) = values_eigenvectors[col_begin+j]; + } + } + + state = LAPACKSupport::State(unusable); +} + + template void LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric (