From 41f6285b6a3ad3ff04c3168082adf30ef72f82c4 Mon Sep 17 00:00:00 2001 From: willems Date: Mon, 26 Apr 2010 14:45:03 +0000 Subject: [PATCH] implementation of computation of generalized eigenvalues for symmetric lapack matrices included git-svn-id: https://svn.dealii.org/trunk@21038 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/source/lapack_full_matrix.cc | 96 +++++++++++++++++++++++- 1 file changed, 95 insertions(+), 1 deletion(-) diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 386b85703f..5d37b9a90c 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -99,7 +99,8 @@ LAPACKFullMatrix::vmult ( const Vector &v, const bool adding) const { - Assert (state == matrix, ExcState(state)); + Assert ((state == matrix) || (state == inverse_matrix), + ExcState(state)); const int mm = this->n_rows(); const int nn = this->n_cols(); @@ -307,6 +308,99 @@ LAPACKFullMatrix::compute_eigenvalues( state = LAPACKSupport::State(eigenvalues | unusable); } + +template +void +LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( + LAPACKFullMatrix & B, + std::vector > & eigenvectors, + const int itype) +{ + Assert(state == matrix, ExcState(state)); + const int nn = this->n_cols(); + Assert(nn == this->n_rows(), ExcNotQuadratic()); + Assert(B.n_rows() == B.n_cols(), ExcNotQuadratic()); + Assert(nn == B.n_cols(), ExcDimensionMismatch (nn, B.n_cols())); + Assert(eigenvectors.size() <= nn, ExcMessage ("eigenvectors.size() > matrix.n_cols()")); + + wr.resize(nn); + wi.resize(nn); //This is set purley for consistency reasons with the + //eigenvalues() function. + + number* values_A = const_cast (this->data()); + number* values_B = const_cast (B.data()); + + int info = 0; + int lwork = 1; + const char * const jobz = (eigenvectors.size() > 0) ? (&V) : (&N); + const char * const uplo = (&U); + + // Optimal workspace query: + + // The LAPACK routine DSYGV 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); + + sygv (&itype, jobz, uplo, &nn, values_A, &nn, + values_B, &nn, + &wr[0], &work[0], &lwork, &info); + // sygv 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 = 3*nn-1 > 1 ? 3*nn-1 : 1; // no query mode +#endif + // resize workspace array + work.resize((unsigned int) lwork); + + // Finally compute the generalized + // eigenvalues. + sygv (&itype, jobz, uplo, &nn, values_A, &nn, + values_B, &nn, + &wr[0], &work[0], &lwork, &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; + + for (unsigned int i=0; i < eigenvectors.size(); ++i) + { + unsigned int col_begin(i*nn); + eigenvectors[i].reinit(nn, true); + for (unsigned int j=0; j < nn; ++j) + { + eigenvectors[i](j) = values_A[col_begin+j]; + } + } + state = LAPACKSupport::State(eigenvalues | unusable); +} + #else template -- 2.39.5