]> https://gitweb.dealii.org/ - dealii.git/commitdiff
eigenvalues
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 29 Mar 2005 03:08:00 +0000 (03:08 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 29 Mar 2005 03:08:00 +0000 (03:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@10285 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/lapack_full_matrix.h
deal.II/lac/include/lac/lapack_support.h
deal.II/lac/include/lac/lapack_templates.h
deal.II/lac/source/lapack_full_matrix.cc

index c979b233772a845841f627af2553cacf7bca1ff9..6e59194e9a5ce25aa1dfe24071acd8398f386ff3 100644 (file)
@@ -16,6 +16,9 @@
 #include <base/table.h>
 #include <lac/lapack_support.h>
 
+#include <vector>
+#include <complex>
+
 // forward declarations
 template<typename number> class Vector;
 template<typename number> class FullMatrix;
@@ -203,6 +206,65 @@ class LAPACKFullMatrix : public TransposeTable<number>
     void Tvmult_add (Vector<number>       &w,
                     const Vector<number> &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<number>
+    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<number> work;
+                                    /**
+                                     * Real parts of
+                                     * eigenvalues. Filled by
+                                     * compute_eigenvalues.
+                                     */
+    std::vector<number> wr;
+    
+                                    /**
+                                     * Imaginary parts of
+                                     * eigenvalues. Filled by
+                                     * compute_eigenvalues.
+                                     */
+    std::vector<number> wi;
+    
 };
 
 
@@ -216,6 +278,8 @@ LAPACKFullMatrix<number>::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<number>::fill (
                   dst_offset_j+entry->column()-src_offset_j)
            = entry->value();
     }
+  
+  state = LAPACKSupport::matrix;
 }
 
 
+template <typename number>
+std::complex<number>
+LAPACKFullMatrix<number>::eigenvalue (unsigned int i) const
+{
+  Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
+  Assert (wr.size() == n_rows(), ExcInternalError());
+  Assert (wi.size() == n_rows(), ExcInternalError());
+  Assert (i<n_rows(), ExcIndexRange(i,0,n_rows()));
+  
+  return std::complex<number>(wr[i], wi[i]);
+}
+
 
 
 #endif
index 8525644fdb64e7f5dd64ec0ec8ad9efacb215917..8105cf3e027684ef6ed7c49b85147e3b9dab9ddf 100644 (file)
@@ -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.
                                    */
index fda099c411c4b3b399c567d6153a13d1b9e1be51..9fad4a7386f9b81be714ed58a8c1a66c92655f32 100644 (file)
@@ -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);
 }
index e2ff5d087107070b22ee6624f97c66f5ae563b55..42d60bc6d0285dd5685f11ba03e26b8231972ff2 100644 (file)
@@ -17,6 +17,8 @@
 #include <lac/full_matrix.h>
 #include <lac/vector.h>
 
+#include <iostream>
+
 using namespace LAPACKSupport;
 
 template <typename number>
@@ -39,7 +41,9 @@ template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
                :
                TransposeTable<number> (M)
-{}
+{
+  state = LAPACKSupport::matrix;
+}
 
 
 template <typename number>
@@ -47,6 +51,7 @@ LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number>& M)
 {
   TransposeTable<number>::operator=(M);
+  state = LAPACKSupport::matrix;
   return *this;
 }
 
@@ -61,6 +66,8 @@ LAPACKFullMatrix<number>::operator = (const FullMatrix<number2>& M)
   for (unsigned int i=0;i<this->n_rows();++i)
     for (unsigned int j=0;j<this->n_cols();++j)
       (*this)(i,j) = M(i,j);
+  
+  state = LAPACKSupport::matrix;
   return *this;
 }
 
@@ -72,13 +79,13 @@ LAPACKFullMatrix<number>::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 <typename number>
 void
@@ -87,6 +94,8 @@ LAPACKFullMatrix<number>::vmult (
   const Vector<number> &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<number>::Tvmult (
   const Vector<number> &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<number>::Tvmult (
 
 #endif
 
-// template <typename number>
-// LAPACKFullMatrix<number>::()
-// {}
+#ifdef HAVE_LIBLAPACK
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_eigenvalues()
+{
+  const int nn = this->n_cols();
+  wr.resize(nn);
+  wi.resize(nn);
+  number* values = const_cast<number*> (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 <typename number>
 // LAPACKFullMatrix<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.