the Arpack solver in deal.II did not output the eigenvectors, but rather the Schur...

diff --git a/include/deal.II/lac/arpack_solver.h b/include/deal.II/lac/arpack_solver.h
index 3900e1cdf9b21456b31912fc6112c0db3b29d4e6..70409a6c735d95c2495a0bcba6bcb6fbaa98da1b 100644 (file)
--- a/include/deal.II/lac/arpack_solver.h
+++ b/include/deal.II/lac/arpack_solver.h
@@ -76,7 +76,16 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
  * <code>dnaupd</code> work and also how to set the parameters appropriately
  * please take a look into the ARPACK manual.
  *
- * @author Baerbel Janssen, Agnieszka Miedlar, 2010.
+ * @note Whenever you eliminate degrees of freedom using
+ * ConstraintMatrix, you generate spurious eigenvalues and
+ * eigenvectors. If you make sure that the diagonals of eliminated
+ * matrix rows are all equal to one, you get a single additional
+ * eigenvalue. But beware that some functions in deal.II set these
+ * diagonals to rather arbitrary (from the point of view of
+ * eigenvalue problems) values. See also @ref step_36 "step-36" for an
+ * example.
+ *
+ * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015
  */
 class ArpackSolver : public Subscriptor
 {
@@ -130,7 +139,40 @@ public:
 
   /**
    * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling
-   * the <code>dneupd</code> and <code>dnaupd</code> functions of ARPACK.
+   * the <code>dneupd</code> and <code>dnaupd</code> functions of
+   * ARPACK.
+   *
+   * The function returns a vector of eigenvalues of length <i>n</i>
+   * and a vector of eigenvectors, where the latter should be twice
+   * the size of the eigenvalue vector. The first <i>n</i> vectors in
+   * <code>eigenvectors</code> will be the real parts of the
+   * eigenvectors, the second <i>n</i> the imaginary parts.
+   *
+   * @param A The operator for which we want to compute
+   * eigenvalues. Actually, this parameter is entirely unused.
+   *
+   * @param B The inner product of the underlying space, typically the
+   * mass matrix. For constrained problems, it can be a partial mass
+   * matrix, like for instance the velocity mass matrix of a Stokes
+   * problem. Only its function <code>vmult()</code> is used.
+   *
+   * @param inverse This is the possibly shifted inverse that is
+   * actually used instead of <code>A</code>. Only its function
+   * <code>vmult()</code> is used.
+   *
+   * @param eigenvalues is a vector of complex numbers in which the
+   * eigenvalues are returned.
+   *
+   * @param eigenvectors is a <b>real</b> vector of eigenvectors,
+   * containing alternatingly the real parts and the imaginary parts of the
+   * eigenvectors. Therefore, its length should be twice the number of
+   * eigenvalues. The vectors have to be initialized to match the
+   * matrices.
+   *
+   * @param n_eigenvalues The purpose of this parameter is not clear,
+   * but it is safe to set it to the size of <code>eigenvalues</code>
+   * or greater. Leave it at its default zero, which will be reset to the size
+   * of <code>eigenvalues</code> internally.
    */
   template <typename VECTOR, typename MATRIX1,
             typename MATRIX2, typename INVERSE>
@@ -140,7 +182,7 @@ public:
     const INVERSE &inverse,
     std::vector<std::complex<double> > &eigenvalues,
     std::vector<VECTOR> &eigenvectors,
-    const unsigned int n_eigenvalues);
+    const unsigned int n_eigenvalues = 0);
 
 protected:
 
@@ -233,23 +275,26 @@ void ArpackSolver::solve (
   //values change magically, so store
   //them here
 
-  const unsigned int n = system_matrix.m();
-  const unsigned int n_inside_arpack = system_matrix.m();
+  const unsigned int n = eigenvectors[0].size();
+  const unsigned int n_inside_arpack = eigenvectors[0].size();
 
+  // Number of eigenvalues for arpack
+  const unsigned int nev = (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
+  AssertIndexRange(eigenvalues.size()-1, nev);
   /*
   if(n < 0 || nev <0 || p < 0 || maxit < 0 )
        std:cout << "All input parameters have to be positive.\n";
        */
   Assert (n_eigenvalues < n,
-          ExcInvalidNumberofEigenvalues(n_eigenvalues, n));
+          ExcInvalidNumberofEigenvalues(nev, n));
 
   Assert (additional_data.number_of_arnoldi_vectors < n,
           ExcInvalidNumberofArnoldiVectors(
             additional_data.number_of_arnoldi_vectors, n));
 
-  Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
+  Assert (additional_data.number_of_arnoldi_vectors > 2*nev+1,
           ExcSmallNumberofArnoldiVectors(
-            additional_data.number_of_arnoldi_vectors, n_eigenvalues));
+            additional_data.number_of_arnoldi_vectors, nev));
   // ARPACK mode for dnaupd, here only mode 3
   int mode = 3;
 
@@ -341,7 +386,6 @@ void ArpackSolver::solve (
   //information out of the iteration
   int info = 1;
 
-  const unsigned int nev = n_eigenvalues;
   while (ido != 99)
     {
       // call of ARPACK dnaupd routine
@@ -445,9 +489,9 @@ void ArpackSolver::solve (
       int rvec = 1;
 
       // which eigenvectors
-      char howmany[4] = "All";
+      char howmany = 'A';
 
-      std::vector<int> select (ncv, 0);
+      std::vector<int> select (ncv, 1);
 
       int ldz = n;
 
@@ -459,11 +503,11 @@ void ArpackSolver::solve (
       int lworkev = 3*ncv;
       std::vector<double> workev (lworkev, 0.);
 
-      std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
-      std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
+      std::vector<double> eigenvalues_real (nev, 0.);
+      std::vector<double> eigenvalues_im (nev, 0.);
 
       // call of ARPACK dneupd routine
-      dneupd_(&rvec, howmany, &select[0], &eigenvalues_real[0],
+      dneupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0],
               &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai,
               &workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
               &resid[0], &ncv, &v[0], &ldv,
@@ -482,9 +526,11 @@ void ArpackSolver::solve (
           Assert (false, ExcArpackInfodneupd(info));
         }
 
-      for (size_type i=0; i<eigenvectors.size(); ++i)
+
+      const unsigned int n_eigenvecs = eigenvectors.size();
+      for (size_type i=0; i<n_eigenvecs; ++i)
         for (unsigned int j=0; j<n; ++j)
-          eigenvectors[i](j) = v[i*n+j];
+          eigenvectors[i](j) = z[i*n+j];
 
       delete[] workd;