]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
implementation of computation of generalized eigenvalues for symmetric lapack matrice...
authorwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 26 Apr 2010 14:45:03 +0000 (14:45 +0000)
committerwillems <willems@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 26 Apr 2010 14:45:03 +0000 (14:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@21038 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/source/lapack_full_matrix.cc

index 386b85703feba12c68d1a32b52f2c79ccd1425e2..5d37b9a90cd9e6cc581469f1fa1d90f29b31f8ef 100644 (file)
@@ -99,7 +99,8 @@ LAPACKFullMatrix<number>::vmult (
   const Vector<number> &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<number>::compute_eigenvalues(
   state = LAPACKSupport::State(eigenvalues | unusable);
 }
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (
+  LAPACKFullMatrix<number> & B,
+  std::vector<Vector<number> > & 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<number*> (this->data());
+  number* values_B = const_cast<number*> (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 <typename 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.