From c602fedffa63d3beb9abedc3b785dc6a6f970084 Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 5 Sep 2006 23:12:47 +0000 Subject: [PATCH] add eigenvalue support to TridiagonalMatrix git-svn-id: https://svn.dealii.org/trunk@13835 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/lapack_support.h | 36 ++++++++++++++++- deal.II/lac/include/lac/tridiagonal_matrix.h | 26 ++++++++++-- deal.II/lac/source/tridiagonal_matrix.cc | 42 ++++++++++++++++++-- 3 files changed, 97 insertions(+), 7 deletions(-) diff --git a/deal.II/lac/include/lac/lapack_support.h b/deal.II/lac/include/lac/lapack_support.h index 6b4455781b..992c958377 100644 --- a/deal.II/lac/include/lac/lapack_support.h +++ b/deal.II/lac/include/lac/lapack_support.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005 by the deal.II authors +// Copyright (C) 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -35,6 +35,26 @@ namespace LAPACKSupport unusable = 0x8000 }; + /** + * Function printing the name of a State. + */ + inline const char* state_name(State s) + { + switch (s) + { + case matrix: + return "matrix"; + case lu: + return "lu decomposition"; + case eigenvalues: + return "eigenvalues"; + case unusable: + return "unusable"; + default: + return "unknown"; + } + return "internal error"; + } /** * A matrix can have certain features allowing for optimization, but @@ -76,6 +96,20 @@ namespace LAPACKSupport * Integer constant. */ static const int one = 1; + + /** + * Exception thrown when a matrix + * is not in a suitable state for + * an operation. For instance, a + * LAPACK routine may have left the + * matrix in an unusable state, + * then vmult does not make sense + * anymore. + */ + DeclException1(ExcState, State, + << "The function cannot be called while the matrix is in state " + << state_name(arg1)); + } diff --git a/deal.II/lac/include/lac/tridiagonal_matrix.h b/deal.II/lac/include/lac/tridiagonal_matrix.h index 4fdcffecfa..135d48c97e 100644 --- a/deal.II/lac/include/lac/tridiagonal_matrix.h +++ b/deal.II/lac/include/lac/tridiagonal_matrix.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -278,8 +279,8 @@ class TridiagonalMatrix ///@name LAPACK operations //@{ /** - * Return the eigenvalues of the - * matrix in the vector provided. + * Compute the eigenvalues of the + * symmetric tridiagonal matrix. * * @note This function requires * configuration of deal.II with @@ -287,7 +288,13 @@ class TridiagonalMatrix * the matrix must use symmetric * storage technique. */ - void eigenvalues(Vector& eigenvalues) const; + void compute_eigenvalues(); + /** + * After calling + * compute_eigenvalues(), you can + * access each eigenvalue here. + */ + number eigenvalue(const unsigned i) const; //@} ///@name Miscellanea //@{ @@ -347,6 +354,19 @@ class TridiagonalMatrix * matrix is assumed symmetric. */ bool is_symmetric; + + /** + * The state of the + * matrix. Normally, the state is + * LAPACKSupport::matrix, + * indicating that the object can + * be used for regular matrix + * operations. + * + * See explanation of this data + * type for details. + */ + LAPACKSupport::State state; }; /**@}*/ diff --git a/deal.II/lac/source/tridiagonal_matrix.cc b/deal.II/lac/source/tridiagonal_matrix.cc index e69d43023e..45ab985d2a 100644 --- a/deal.II/lac/source/tridiagonal_matrix.cc +++ b/deal.II/lac/source/tridiagonal_matrix.cc @@ -13,6 +13,9 @@ #include #include +#include + +using namespace LAPACKSupport; template TridiagonalMatrix::TridiagonalMatrix( @@ -22,7 +25,8 @@ TridiagonalMatrix::TridiagonalMatrix( diagonal(size, 0.), left((symmetric ? 0 : size), 0.), right(size, 0.), - is_symmetric(symmetric) + is_symmetric(symmetric), + state(matrix) {} @@ -36,6 +40,7 @@ TridiagonalMatrix::reinit( diagonal.resize(size); right.resize(size); left.resize(symmetric ? 0 : size); + state = matrix; } @@ -43,6 +48,8 @@ template bool TridiagonalMatrix::all_zero() const { + Assert(state == matrix, ExcState(state)); + typename std::vector::const_iterator i; typename std::vector::const_iterator e; @@ -68,6 +75,8 @@ TridiagonalMatrix::vmult ( const Vector &v, const bool adding) const { + Assert(state == matrix, ExcState(state)); + Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n())); Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n())); @@ -128,6 +137,8 @@ TridiagonalMatrix::Tvmult ( const Vector &v, const bool adding) const { + Assert(state == matrix, ExcState(state)); + Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n())); Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n())); @@ -177,6 +188,8 @@ TridiagonalMatrix::matrix_scalar_product( const Vector& w, const Vector& v) const { + Assert(state == matrix, ExcState(state)); + const unsigned int e=n()-1; typename std::vector::const_iterator d = diagonal.begin(); typename std::vector::const_iterator r = right.begin(); @@ -206,11 +219,34 @@ TridiagonalMatrix::matrix_norm_square( template void -TridiagonalMatrix::eigenvalues(Vector&) const +TridiagonalMatrix::compute_eigenvalues() +{ +#ifdef HAVE_LIBLAPACK + Assert(state == matrix, ExcState(state)); + Assert(is_symmetric, ExcNotImplemented()); + + const int nn = n(); + int info; + stev (&N, &nn, &*diagonal.begin(), &*right.begin(), 0, &one, 0, &info); + Assert(info == 0, ExcInternalError()); + + state = eigenvalues; +#else + Assert(false, ExcNeedsLAPACK()); +#endif +} + + +template +number +TridiagonalMatrix::eigenvalue(const unsigned int i) const { - Assert(false, ExcNotImplemented()); + Assert(state == eigenvalues, ExcState(state)); + Assert(i TridiagonalMatrix:: -- 2.39.5