From f0520be89e409a7c14b2b80623e7ec4415f59507 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 29 Mar 2005 03:08:00 +0000 Subject: [PATCH] eigenvalues git-svn-id: https://svn.dealii.org/trunk@10285 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/lapack_full_matrix.h | 78 ++++++++++++++++++++ deal.II/lac/include/lac/lapack_support.h | 20 ++++- deal.II/lac/include/lac/lapack_templates.h | 66 +++++++++-------- deal.II/lac/source/lapack_full_matrix.cc | 67 +++++++++++++++-- 4 files changed, 189 insertions(+), 42 deletions(-) diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h index c979b23377..6e59194e9a 100644 --- a/deal.II/lac/include/lac/lapack_full_matrix.h +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -16,6 +16,9 @@ #include #include +#include +#include + // forward declarations template class Vector; template class FullMatrix; @@ -203,6 +206,65 @@ class LAPACKFullMatrix : public TransposeTable void Tvmult_add (Vector &w, const Vector &v) const; + /** + * Compute eigenvalues of the + * matrix. After this routine has + * been called, eigenvalues can + * be retrieved using the + * eigenvalue() function. The + * matrix itself will be + * LAPACKSupport::unusable after + * this operation. + * + * @note Calls the LAPACK + * function Xgeev. + */ + void compute_eigenvalues (); + + /** + * Retrieve eigenvalue after + * compute_eigenvalues() was + * called. + */ + std::complex + eigenvalue (unsigned int i) const; + + private: + /** + * Since LAPACK operations + * notoriously change the meaning + * of the matrix entries, we + * record the current state after + * the last operation here. + */ + LAPACKSupport::State state; + /** + * Additional properties of the + * matrix which may help to + * select more efficient LAPACK + * functions. + */ + LAPACKSupport::Properties properties; + + /** + * The working array used for + * some LAPACK functions. + */ + mutable std::vector work; + /** + * Real parts of + * eigenvalues. Filled by + * compute_eigenvalues. + */ + std::vector wr; + + /** + * Imaginary parts of + * eigenvalues. Filled by + * compute_eigenvalues. + */ + std::vector wi; + }; @@ -216,6 +278,8 @@ LAPACKFullMatrix::copy_from (const MATRIX& M) for (typename MATRIX::const_iterator entry = M.begin(); entry != end; ++entry) this->el(entry->row(), entry->column()) = entry->value(); + + state = LAPACKSupport::matrix; } @@ -243,9 +307,23 @@ LAPACKFullMatrix::fill ( dst_offset_j+entry->column()-src_offset_j) = entry->value(); } + + state = LAPACKSupport::matrix; } +template +std::complex +LAPACKFullMatrix::eigenvalue (unsigned int i) const +{ + Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState()); + Assert (wr.size() == n_rows(), ExcInternalError()); + Assert (wi.size() == n_rows(), ExcInternalError()); + Assert (i(wr[i], wi[i]); +} + #endif diff --git a/deal.II/lac/include/lac/lapack_support.h b/deal.II/lac/include/lac/lapack_support.h index 8525644fdb..8105cf3e02 100644 --- a/deal.II/lac/include/lac/lapack_support.h +++ b/deal.II/lac/include/lac/lapack_support.h @@ -25,12 +25,14 @@ namespace LAPACKSupport */ enum State { - /// Contents is something useless. - unusable, /// Contents is actually a matrix. matrix, /// Contents is an LU decomposition. - lu + lu, + /// Eigenvalue vector is filled + eigenvalues, + /// Contents is something useless. + unusable = 0x8000 }; @@ -54,6 +56,18 @@ namespace LAPACKSupport hessenberg = 8 }; + /** + * Character constant. + */ + extern const char V = 'V'; + /** + * Character constant. + */ + extern const char T = 'T'; + /** + * Character constant. + */ + extern const char N = 'N'; /** * Integer constant. */ diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index fda099c411..9fad4a7386 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -30,36 +30,38 @@ void sgemv_ (const char* trans, const int* m, const int* n, const float* x, const int* incx, const float* b, float* y, const int* incy); // Compute eigenvalues and vectors -void dgeev_ (char* jobvl, char* jobvr, - int* n, double* A, int* lda, +void dgeev_ (const char* jobvl, const char* jobvr, + const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, - double* vl, int* ldvl, - double* vr, int* ldva, - double* work, int* lwork, + double* vl, const int* ldvl, + double* vr, const int* ldva, + double* work, const int* lwork, int* info); -void sgeev_ (char* jobvl, char* jobvr, - int* n, float* A, int* lda, +void sgeev_ (const char* jobvl, const char* jobvr, + const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, - float* vl, int* ldvl, - float* vr, int* ldva, - float* work, int* lwork, + float* vl, const int* ldvl, + float* vr, const int* ldva, + float* work, const int* lwork, int* info); // Compute eigenvalues and vectors (expert) -void dgeevx_ (char* balanc, char* jobvl, char* jobvr, char* sense, - int* n, double* A, int* lda, +void dgeevx_ (const char* balanc, const char* jobvl, const char* jobvr, + const char* sense, + const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, - double* vl, int* ldvl, - double* vr, int* ldvr, + double* vl, const int* ldvl, + double* vr, const int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, int* lwork, int* iwork, int* info); -void sgeevx_ (char* balanc, char* jobvl, char* jobvr, char* sense, - int* n, float* A, int* lda, +void sgeevx_ (const char* balanc, const char* jobvl, const char* jobvr, + const char* sense, + const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, - float* vl, int* ldvl, - float* vr, int* ldvr, + float* vl, const int* ldvl, + float* vr, const int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, @@ -67,18 +69,18 @@ void sgeevx_ (char* balanc, char* jobvl, char* jobvr, char* sense, int* iwork, int* info); // Compute singular value decomposition void dgesvd_ (int* jobu, int* jobvt, - int* n, int* m, double* A, int* lda, + const int* n, const int* m, double* A, const int* lda, double* s, - double* u, int* ldu, - double* vt, int* ldvt, - double* work, int* lwork, + double* u, const int* ldu, + double* vt, const int* ldvt, + double* work, const int* lwork, int* info); void sgesvd_ (int* jobu, int* jobvt, - int* n, int* m, float* A, int* lda, + const int* n, const int* m, float* A, const int* lda, float* s, - float* u, int* ldu, - float* vt, int* ldvt, - float* work, int* lwork, + float* u, const int* ldu, + float* vt, const int* ldvt, + float* work, const int* lwork, int* info); } @@ -100,42 +102,42 @@ gemv (const char* trans, const int* m, const int* n, const float* alpha, const f inline void -geev (char* jobvl, char* jobvr, int* n, double* A, int* lda, double* lambda_re, double* lambda_im, double* vl, int* ldvl, double* vr, int* ldva, double* work, int* lwork, int* info) +geev (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldva, double* work, const int* lwork, int* info) { dgeev_ (jobvl,jobvr,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldva,work,lwork,info); } inline void -geev (char* jobvl, char* jobvr, int* n, float* A, int* lda, float* lambda_re, float* lambda_im, float* vl, int* ldvl, float* vr, int* ldva, float* work, int* lwork, int* info) +geev (const char* jobvl, const char* jobvr, const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, float* vl, const int* ldvl, float* vr, const int* ldva, float* work, const int* lwork, int* info) { sgeev_ (jobvl,jobvr,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldva,work,lwork,info); } inline void -geevx (char* balanc, char* jobvl, char* jobvr, char* sense, int* n, double* A, int* lda, double* lambda_re, double* lambda_im, double* vl, int* ldvl, double* vr, int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, int* lwork, int* iwork, int* info) +geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldvr, int* ilo, int* ihi, double* scale, double* abnrm, double* rconde, double* rcondv, double* work, int* lwork, int* iwork, int* info) { dgeevx_ (balanc,jobvl,jobvr,sense,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldvr,ilo,ihi,scale,abnrm,rconde,rcondv,work,lwork,iwork,info); } inline void -geevx (char* balanc, char* jobvl, char* jobvr, char* sense, int* n, float* A, int* lda, float* lambda_re, float* lambda_im, float* vl, int* ldvl, float* vr, int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, float* work, int* lwork, int* iwork, int* info) +geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sense, const int* n, float* A, const int* lda, float* lambda_re, float* lambda_im, float* vl, const int* ldvl, float* vr, const int* ldvr, int* ilo, int* ihi, float* scale, float* abnrm, float* rconde, float* rcondv, float* work, int* lwork, int* iwork, int* info) { sgeevx_ (balanc,jobvl,jobvr,sense,n,A,lda,lambda_re,lambda_im,vl,ldvl,vr,ldvr,ilo,ihi,scale,abnrm,rconde,rcondv,work,lwork,iwork,info); } inline void -gesvd (int* jobu, int* jobvt, int* n, int* m, double* A, int* lda, double* s, double* u, int* ldu, double* vt, int* ldvt, double* work, int* lwork, int* info) +gesvd (int* jobu, int* jobvt, const int* n, const int* m, double* A, const int* lda, double* s, double* u, const int* ldu, double* vt, const int* ldvt, double* work, const int* lwork, int* info) { dgesvd_ (jobu,jobvt,n,m,A,lda,s,u,ldu,vt,ldvt,work,lwork,info); } inline void -gesvd (int* jobu, int* jobvt, int* n, int* m, float* A, int* lda, float* s, float* u, int* ldu, float* vt, int* ldvt, float* work, int* lwork, int* info) +gesvd (int* jobu, int* jobvt, const int* n, const int* m, float* A, const int* lda, float* s, float* u, const int* ldu, float* vt, const int* ldvt, float* work, const int* lwork, int* info) { sgesvd_ (jobu,jobvt,n,m,A,lda,s,u,ldu,vt,ldvt,work,lwork,info); } diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index e2ff5d0871..42d60bc6d0 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -17,6 +17,8 @@ #include #include +#include + using namespace LAPACKSupport; template @@ -39,7 +41,9 @@ template LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) : TransposeTable (M) -{} +{ + state = LAPACKSupport::matrix; +} template @@ -47,6 +51,7 @@ LAPACKFullMatrix & LAPACKFullMatrix::operator = (const LAPACKFullMatrix& M) { TransposeTable::operator=(M); + state = LAPACKSupport::matrix; return *this; } @@ -61,6 +66,8 @@ LAPACKFullMatrix::operator = (const FullMatrix& M) for (unsigned int i=0;in_rows();++i) for (unsigned int j=0;jn_cols();++j) (*this)(i,j) = M(i,j); + + state = LAPACKSupport::matrix; return *this; } @@ -72,13 +79,13 @@ LAPACKFullMatrix::operator = (const double d) Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); if (this->n_elements() != 0) - std::fill_n (this->val, this->n_elements(), number()); + std::fill_n (this->val, this->n_elements(), number()); + state = LAPACKSupport::matrix; return *this; } -#ifdef HAVE_LIBLAPACK - +#ifdef HAVE_LIBBLAS template void @@ -87,6 +94,8 @@ LAPACKFullMatrix::vmult ( const Vector &v, const bool adding) const { + Assert (state == matrix, ExcInvalidState()); + const int mm = this->n_rows(); const int nn = this->n_cols(); const number alpha = 1.; @@ -103,6 +112,8 @@ LAPACKFullMatrix::Tvmult ( const Vector &v, const bool adding) const { + Assert (state == matrix, ExcInvalidState()); + const int mm = this->n_rows(); const int nn = this->n_cols(); const number alpha = 1.; @@ -138,10 +149,52 @@ LAPACKFullMatrix::Tvmult ( #endif -// template -// LAPACKFullMatrix::() -// {} +#ifdef HAVE_LIBLAPACK + +template +void +LAPACKFullMatrix::compute_eigenvalues() +{ + const int nn = this->n_cols(); + wr.resize(nn); + wi.resize(nn); + number* values = const_cast (this->data()); + + int info; + int lwork = -1; + work.resize(1); + geev(&N, &N, &nn, values, &nn, + &wr[0], &wi[0], + 0, &one, 0, &one, + &work[0], &lwork, &info); + // geev 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); + work.resize((unsigned int) lwork); + // Finally compute the eigenvalues. + geev(&N, &N, &nn, values, &nn, + &wr[0], &wi[0], + 0, &one, 0, &one, + &work[0], &lwork, &info); + // Negative return value implies a + // wrong argument. This should be + // internal. + Assert (info >=0, ExcInternalError()); +//TODO:[GK] What if the QR method fails? + if (info != 0) + std::cerr << "LAPACK error in geev" << std::endl; + + state = LAPACKSupport::State(eigenvalues | unusable); +} +#endif // template // LAPACKFullMatrix::() -- 2.39.5