]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add eigenvalue support to TridiagonalMatrix
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 5 Sep 2006 23:12:47 +0000 (23:12 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 5 Sep 2006 23:12:47 +0000 (23:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@13835 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/lapack_support.h
deal.II/lac/include/lac/tridiagonal_matrix.h
deal.II/lac/source/tridiagonal_matrix.cc

index 6b4455781bb5ee25e005cdcd78243da4c0654fcc..992c958377e9950ea0b789ea6cf5b17f2b08878d 100644 (file)
@@ -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));
+  
 }
 
 
index 4fdcffecfa7d3b9c18479dcd80bf9cdd0b6b793c..135d48c97ef24f583956fee6b219a5e3afd7903f 100644 (file)
@@ -16,6 +16,7 @@
 #include <base/config.h>
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
+#include <lac/lapack_support.h>
 
 #include <vector>
 #include <iomanip>
@@ -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<number>& 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;
 };
 
 /**@}*/
index e69d43023e8bcc9d805906880ceebfd187a16292..45ab985d2a83e7d77a5d9944bdb36cc320882aff 100644 (file)
@@ -13,6 +13,9 @@
 
 #include <lac/tridiagonal_matrix.h>
 #include <lac/vector.h>
+#include <lac/lapack_templates.h>
+
+using namespace LAPACKSupport;
 
 template<typename number>
 TridiagonalMatrix<number>::TridiagonalMatrix(
@@ -22,7 +25,8 @@ TridiagonalMatrix<number>::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<number>::reinit(
   diagonal.resize(size);
   right.resize(size);
   left.resize(symmetric ? 0 : size);
+  state = matrix;
 }
 
 
@@ -43,6 +48,8 @@ template<typename number>
 bool
 TridiagonalMatrix<number>::all_zero() const
 {
+  Assert(state == matrix, ExcState(state));
+  
   typename std::vector<number>::const_iterator i;
   typename std::vector<number>::const_iterator e;
 
@@ -68,6 +75,8 @@ TridiagonalMatrix<number>::vmult (
   const Vector<number> &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<number>::Tvmult (
   const Vector<number> &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<number>::matrix_scalar_product(
   const Vector<number>& w,
   const Vector<number>& v) const
 {
+  Assert(state == matrix, ExcState(state));
+  
   const unsigned int e=n()-1;
   typename std::vector<number>::const_iterator d = diagonal.begin();
   typename std::vector<number>::const_iterator r = right.begin();
@@ -206,11 +219,34 @@ TridiagonalMatrix<number>::matrix_norm_square(
 
 template<typename number>
 void
-TridiagonalMatrix<number>::eigenvalues(Vector<number>&) const
+TridiagonalMatrix<number>::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<typename number>
+number
+TridiagonalMatrix<number>::eigenvalue(const unsigned int i) const
 {
-  Assert(false, ExcNotImplemented());
+  Assert(state == eigenvalues, ExcState(state));
+  Assert(i<n(), ExcIndexRange(i,0,n()));
+  return diagonal[i];
 }
 
+
 /*
 template<typename number>
 TridiagonalMatrix<number>::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.