]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix (P)ARPACK interface for non-symmetric matrices 2871/head
authorJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Thu, 21 Jul 2016 13:19:58 +0000 (15:19 +0200)
committerJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Mon, 1 Aug 2016 07:49:09 +0000 (09:49 +0200)
-> test for length n+1 of the eigenvectors in the non-symmetric case, which is related to the storage scheme of complex components
-> allow ARPACK to store n+1 eigenvalues in the work arrays in case the eigenvalues n and n+1 form a complex conjugate pair
-> copy nev eigenvectors from the work arrays as nev is increased by one by dneupd in some cases
-> add set_initial_vector function
-> add tests for advection-diffusion FE matrices with complex eigenvalues
-> improve description of the tests

12 files changed:
doc/news/changes.h
include/deal.II/lac/arpack_solver.h
include/deal.II/lac/parpack_solver.h
tests/arpack/arpack_advection_diffusion.cc [new file with mode: 0644]
tests/arpack/arpack_advection_diffusion.output [new file with mode: 0644]
tests/arpack/parpack_advection_diffusion_petsc.cc [new file with mode: 0644]
tests/arpack/parpack_advection_diffusion_petsc.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output [new file with mode: 0644]
tests/arpack/parpack_advection_diffusion_trilinos.cc [new file with mode: 0644]
tests/arpack/parpack_advection_diffusion_trilinos.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output [new file with mode: 0644]
tests/arpack/step-36_ar.cc
tests/arpack/step-36_parpack.cc
tests/arpack/step-36_parpack_trilinos.cc

index d917083a7d3d5feb9ced9c8305301af9172c6e50..761e56ba88936e4748b3380c1586c5c92a79c65c 100644 (file)
@@ -395,6 +395,10 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> Fixed: (P)ARPACK interface for non-symmetric matrices.
+ <br>
+ (Joscha Gedicke, 2016/08/01)
+
  <li> Fixed: The TrilinosWrappers::SparsityPattern::print() and
  TrilinosWrappers::SparsityPattern::print_gnuplot() methods did not produce
  correct output on distributed computations. This is now fixed.
index 63751a8250365eede877d1120cded9427d9e57ee..dee70d45c75fe03308d8e9f10add1ec6ee808ee2 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
-extern "C" void dnaupd_(int *ido, char *bmat, const unsigned int *n, char *which,
-                        const unsigned int *nev, const double *tol, double *resid, int *ncv,
+extern "C" void dnaupd_(int *ido, char *bmat, unsigned int *n, char *which,
+                        unsigned int *nev, const double *tol, double *resid, int *ncv,
+                        double *v, int *ldv, int *iparam, int *ipntr,
+                        double *workd, double *workl, int *lworkl,
+                        int *info);
+
+extern "C" void dsaupd_(int *ido, char *bmat, unsigned int *n, char *which,
+                        unsigned int *nev, double *tol, double *resid, int *ncv,
                         double *v, int *ldv, int *iparam, int *ipntr,
                         double *workd, double *workl, int *lworkl,
                         int *info);
 
 extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
                         double *di, double *z, int *ldz, double *sigmar,
-                        double *sigmai, double *workev, char *bmat,const unsigned int *n, char *which,
-                        const unsigned int *nev, const double *tol, double *resid, int *ncv,
+                        double *sigmai, double *workev, char *bmat, unsigned int *n, char *which,
+                        unsigned int *nev, double *tol, double *resid, int *ncv,
+                        double *v, int *ldv, int *iparam, int *ipntr,
+                        double *workd, double *workl, int *lworkl, int *info);
+
+extern "C" void dseupd_(int *rvec, char *howmany, int *select, double *d,
+                        double *z, int *ldz, double *sigmar,
+                        char *bmat, unsigned int *n, char *which,
+                        unsigned int *nev, double *tol, double *resid, int *ncv,
                         double *v, int *ldv, int *iparam, int *ipntr,
                         double *workd, double *workl, int *lworkl, int *info);
 
 /**
  * Interface for using ARPACK. ARPACK is a collection of Fortran77 subroutines
  * designed to solve large scale eigenvalue problems.  Here we interface to
- * the routines <code>dneupd</code> and <code>dnaupd</code> of ARPACK.  The
+ * the routines <code>dnaupd</code> and <code>dneupd</code> of ARPACK.
+ * If the operator is specified to be symmetric we use the symmetric interface
+ * <code>dsaupd</code> and <code>dseupd</code> of ARPACK instead.  The
  * package is designed to compute a few eigenvalues and corresponding
  * eigenvectors of a general n by n matrix A. It is most appropriate for large
  * sparse matrices A.
@@ -72,8 +87,9 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
  * Through the AdditionalData the user can specify some of the parameters to
  * be set.
  *
- * For further information on how the ARPACK routines <code>dneupd</code> and
- * <code>dnaupd</code> work and also how to set the parameters appropriately
+ * For further information on how the ARPACK routines <code>dsaupd</code>,
+ * <code>dseupd</code>, <code>dnaupd</code> and <code>dneupd</code> work
+ * and also how to set the parameters appropriately
  * please take a look into the ARPACK manual.
  *
  * @note Whenever you eliminate degrees of freedom using ConstraintMatrix, you
@@ -85,7 +101,7 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
  * @ref step_36 "step-36"
  * for an example.
  *
- * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015
+ * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015, Joscha Gedicke 2016
  */
 class ArpackSolver : public Subscriptor
 {
@@ -121,9 +137,11 @@ public:
   {
     const unsigned int number_of_arnoldi_vectors;
     const WhichEigenvalues eigenvalue_of_interest;
+    const bool symmetric;
     AdditionalData(
       const unsigned int number_of_arnoldi_vectors = 15,
-      const WhichEigenvalues eigenvalue_of_interest = largest_magnitude);
+      const WhichEigenvalues eigenvalue_of_interest = largest_magnitude,
+      const bool symmetric = false);
   };
 
   /**
@@ -137,15 +155,28 @@ public:
   ArpackSolver(SolverControl &control,
                const AdditionalData &data = AdditionalData());
 
+  /**
+   * Set initial vector for building Krylov space.
+   */
+  template <typename VectorType>
+  void set_initial_vector(const VectorType &vec);
+
   /**
    * 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>dsaupd</code> and <code>dseupd</code> or
+   * <code>dnaupd</code> and <code>dneupd</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.
+   * vector of eigenvectors of length <i>n</i> in the symmetric case
+   * and of length <i>n+1</i> in the non-symmetric case. In the symmetric case
+   * all eigenvectors are real. In the non-symmetric case complex eigenvalues
+   * always occur as complex conjugate pairs. Therefore the eigenvector for an
+   * eigenvalue with nonzero complex part is stored by putting the real and
+   * the imaginary parts in consecutive real-valued vectors. The eigenvector
+   * of the complex conjugate eigenvalue does not need to be stored, since it
+   * is just the complex conjugate of the stored eigenvector. Thus, if the last
+   * n-th eigenvalue has a nonzero imaginary part, Arpack needs in total n+1
+   * real-valued vectors to store real and imaginary parts of the eigenvectors.
    *
    * @param A The operator for which we want to compute eigenvalues. Actually,
    * this parameter is entirely unused.
@@ -163,9 +194,16 @@ public:
    * 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.
+   * the real parts of all eigenvectors and the imaginary parts of the eigenvectors
+   * corresponding to complex conjugate eigenvalue pairs.
+   * Therefore, its length should be <i>n</i> in the symmetric case and <i>n+1</i>
+   * in the non-symmetric case. In the non-symmetric case the storage scheme
+   * leads for example to the following pattern. Suppose that the first two
+   * eigenvalues are real and the third and fourth are a complex conjugate
+   * pair. Asking for three eigenpairs results in <i>[real(v1),real(v2),
+   * real(v3),imag(v3)]</i>. Note that we get the same pattern if we ask for four
+   * eigenpairs in this example, since the fourth eigenvector is simply the
+   * complex conjugate of the third one.
    *
    * @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.
@@ -194,43 +232,70 @@ protected:
    */
   const AdditionalData additional_data;
 
+  /**
+   * Store an initial vector
+   */
+  bool initial_vector_provided;
+  std::vector<double> resid;
+
 private:
 
   /**
    * Exceptions.
    */
-  DeclException2 (ExcInvalidNumberofEigenvalues, int, int,
+  DeclException2 (ArpackExcInvalidNumberofEigenvalues, int, int,
                   << "Number of wanted eigenvalues " << arg1
                   << " is larger that the size of the matrix " << arg2);
 
-  DeclException2 (ExcInvalidNumberofArnoldiVectors, int, int,
+  DeclException2 (ArpackExcInvalidEigenvectorSize, int, int,
+                  << "Number of wanted eigenvalues " << arg1
+                  << " is larger that the size of eigenvectors " << arg2);
+
+  DeclException2 (ArpackExcInvalidEigenvectorSizeNonsymmetric, int, int,
+                  << "To store the real and complex parts of " << arg1
+                  << " eigenvectors in real-valued vectors, their size (currently set to " << arg2
+                  << ") should be greater than or equal to " << arg1+1);
+
+  DeclException2 (ArpackExcInvalidEigenvalueSize, int, int,
+                  << "Number of wanted eigenvalues " << arg1
+                  << " is larger that the size of eigenvalues " << arg2);
+
+  DeclException2 (ArpackExcInvalidNumberofArnoldiVectors, int, int,
                   << "Number of Arnoldi vectors " << arg1
                   << " is larger that the size of the matrix " << arg2);
 
-  DeclException2 (ExcSmallNumberofArnoldiVectors, int, int,
+  DeclException2 (ArpackExcSmallNumberofArnoldiVectors, int, int,
                   << "Number of Arnoldi vectors " << arg1
                   << " is too small to obtain " << arg2
                   << " eigenvalues");
 
-  DeclException1 (ExcArpackIdo, int, << "This ido " << arg1
+  DeclException1 (ArpackExcArpackIdo, int, << "This ido " << arg1
                   << " is not supported. Check documentation of ARPACK");
 
-  DeclException1 (ExcArpackMode, int, << "This mode " << arg1
+  DeclException1 (ArpackExcArpackMode, int, << "This mode " << arg1
                   << " is not supported. Check documentation of ARPACK");
 
-  DeclException1 (ExcArpackInfodsaupd, int,
+  DeclException1 (ArpackExcArpackInfodsaupd, int,
                   << "Error with dsaupd, info " << arg1
                   << ". Check documentation of ARPACK");
 
-  DeclException1 (ExcArpackInfodneupd, int,
+  DeclException1 (ArpackExcArpackInfodnaupd, int,
+                  << "Error with dnaupd, info " << arg1
+                  << ". Check documentation of ARPACK");
+
+  DeclException1 (ArpackExcArpackInfodseupd, int,
+                  << "Error with dseupd, info " << arg1
+                  << ". Check documentation of ARPACK");
+
+  DeclException1 (ArpackExcArpackInfodneupd, int,
                   << "Error with dneupd, info " << arg1
                   << ". Check documentation of ARPACK");
 
-  DeclException1 (ExcArpackInfoMaxIt, int,
+  DeclException1 (ArpackExcArpackInfoMaxIt, int,
                   << "Maximum number " << arg1
                   << " of iterations reached.");
 
-  DeclExceptionMsg (ExcArpackNoShifts,
+  DeclExceptionMsg (ArpackExcArpackNoShifts,
                     "No shifts could be applied during implicit"
                     " Arnoldi update, try increasing the number of"
                     " Arnoldi vectors.");
@@ -240,10 +305,12 @@ private:
 inline
 ArpackSolver::AdditionalData::
 AdditionalData (const unsigned int number_of_arnoldi_vectors,
-                const WhichEigenvalues eigenvalue_of_interest)
+                const WhichEigenvalues eigenvalue_of_interest,
+                const bool symmetric)
   :
   number_of_arnoldi_vectors(number_of_arnoldi_vectors),
-  eigenvalue_of_interest(eigenvalue_of_interest)
+  eigenvalue_of_interest(eigenvalue_of_interest),
+  symmetric(symmetric)
 {}
 
 
@@ -252,10 +319,21 @@ ArpackSolver::ArpackSolver (SolverControl &control,
                             const AdditionalData &data)
   :
   solver_control (control),
-  additional_data (data)
-
+  additional_data (data),
+  initial_vector_provided(false)
 {}
 
+template <typename VectorType>
+inline
+void ArpackSolver::
+set_initial_vector(const VectorType &vec)
+{
+  initial_vector_provided = true;
+  resid.resize(vec.size());
+  for (size_type i = 0; i < vec.size(); ++i)
+    resid[i] = vec[i];
+}
+
 
 template <typename VectorType, typename MatrixType1,
           typename MatrixType2, typename INVERSE>
@@ -267,31 +345,43 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
                           std::vector<VectorType>            &eigenvectors,
                           const unsigned int                  n_eigenvalues)
 {
-  //inside the routines of ARPACK the
-  //values change magically, so store
-  //them here
-
-  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(nev, n));
+  // Problem size
+  unsigned int n = eigenvectors[0].size();
+
+  // Number of eigenvalues
+  const unsigned int nev_const = (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
+  // nev for arpack, which might change by plus one during dneupd
+  unsigned int nev = nev_const;
+
+  // check input sizes
+  if (additional_data.symmetric)
+    {
+      Assert (nev <= eigenvectors.size(),
+              ArpackExcInvalidEigenvectorSize(nev, eigenvectors.size()));
+    }
+  else
+    Assert (nev+1 <= eigenvectors.size(),
+            ArpackExcInvalidEigenvectorSizeNonsymmetric(nev, eigenvectors.size()));
+
+  Assert (nev <= eigenvalues.size(),
+          ArpackExcInvalidEigenvalueSize(nev, eigenvalues.size()));
+
+  // check large enough problem size
+  Assert (nev < n,
+          ArpackExcInvalidNumberofEigenvalues(nev, n));
 
   Assert (additional_data.number_of_arnoldi_vectors < n,
-          ExcInvalidNumberofArnoldiVectors(
+          ArpackExcInvalidNumberofArnoldiVectors(
             additional_data.number_of_arnoldi_vectors, n));
 
+  // check whether we have enough Arnoldi vectors
   Assert (additional_data.number_of_arnoldi_vectors > 2*nev+1,
-          ExcSmallNumberofArnoldiVectors(
+          ArpackExcSmallNumberofArnoldiVectors(
             additional_data.number_of_arnoldi_vectors, nev));
-  // ARPACK mode for dnaupd, here only mode 3
+
+  /* ARPACK mode for dsaupd/dnaupd, here only mode 3,
+   * i.e. shift-invert mode
+   */
   int mode = 3;
 
   // reverse communication parameter
@@ -342,10 +432,11 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
     }
 
   // tolerance for ARPACK
-  const double tol = control().tolerance();
+  double tol = control().tolerance();
 
   // if the starting vector is used it has to be in resid
-  std::vector<double> resid(n, 1.);
+  if (!initial_vector_provided || resid.size() != n)
+    resid.resize(n, 1.);
 
   // number of Arnoldi basis vectors specified
   // in additional_data
@@ -371,23 +462,24 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
   std::vector<int> ipntr (14, 0);
 
   // work arrays for ARPACK
-  double *workd;
-  workd = new double[3*n];
-
-  for (unsigned int i=0; i<3*n; ++i)
-    workd[i] = 0.0;
-
-  int lworkl = 3*ncv*(ncv + 6);
+  std::vector<double> workd (3*n, 0.);
+  int lworkl = additional_data.symmetric ? ncv*ncv + 8*ncv : 3*ncv*ncv+6*ncv;
   std::vector<double> workl (lworkl, 0.);
+
   //information out of the iteration
   int info = 1;
 
   while (ido != 99)
     {
-      // call of ARPACK dnaupd routine
-      dnaupd_(&ido, bmat, &n_inside_arpack, which, &nev, &tol,
-              &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
-              workd, &workl[0], &lworkl, &info);
+      // call of ARPACK dsaupd/dnaupd routine
+      if (additional_data.symmetric)
+        dsaupd_(&ido, bmat, &n, which, &nev, &tol,
+                &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
+                &workd[0], &workl[0], &lworkl, &info);
+      else
+        dnaupd_(&ido, bmat, &n, which, &nev, &tol,
+                &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
+                &workd[0], &workl[0], &lworkl, &info);
 
       if (ido == 99)
         break;
@@ -408,7 +500,7 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
 
 
               for (size_type i=0; i<src.size(); ++i)
-                src(i) = *(workd+ipntr[0]-1+i);
+                src(i) = workd[ipntr[0]-1+i];
 
               // multiplication with mass matrix M
               mass_matrix.vmult(tmp, src);
@@ -416,7 +508,7 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
               inverse.vmult(dst,tmp);
 
               for (size_type i=0; i<dst.size(); ++i)
-                *(workd+ipntr[1]-1+i) = dst(i);
+                workd[ipntr[1]-1+i] = dst(i);
             }
             break;
 
@@ -431,14 +523,14 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
 
               for (size_type i=0; i<src.size(); ++i)
                 {
-                  src(i) = *(workd+ipntr[2]-1+i);
-                  tmp(i) = *(workd+ipntr[0]-1+i);
+                  src(i) = workd[ipntr[2]-1+i];
+                  tmp(i) = workd[ipntr[0]-1+i];
                 }
               // solving linear system
               inverse.vmult(dst,src);
 
               for (size_type i=0; i<dst.size(); ++i)
-                *(workd+ipntr[1]-1+i) = dst(i);
+                workd[ipntr[1]-1+i] = dst(i);
             }
             break;
 
@@ -450,32 +542,37 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
               dst.reinit(src);
 
               for (size_type i=0; i<src.size(); ++i)
-                src(i) = *(workd+ipntr[0]-1+i);
+                src(i) = workd[ipntr[0]-1+i];
 
               // Multiplication with mass matrix M
               mass_matrix.vmult(dst, src);
 
               for (size_type i=0; i<dst.size(); ++i)
-                *(workd+ipntr[1]-1+i) = dst(i);
+                workd[ipntr[1]-1+i] = dst(i);
 
             }
             break;
 
             default:
-              Assert (false, ExcArpackIdo(ido));
+              Assert (false, ArpackExcArpackIdo(ido));
               break;
             }
         }
         break;
         default:
-          Assert (false, ExcArpackMode(mode));
+          Assert (false, ArpackExcArpackMode(mode));
           break;
         }
     }
 
   if (info<0)
     {
-      Assert (false, ExcArpackInfodsaupd(info));
+      if (additional_data.symmetric)
+        {
+          Assert (false, ArpackExcArpackInfodsaupd(info));
+        }
+      else
+        Assert (false, ArpackExcArpackInfodnaupd(info));
     }
   else
     {
@@ -491,49 +588,54 @@ void ArpackSolver::solve (const MatrixType1                  &/*system_matrix*/,
 
       int ldz = n;
 
-      std::vector<double> z (ldz*ncv, 0.);
-
       double sigmar = 0.0; // real part of the shift
       double sigmai = 0.0; // imaginary part of the shift
 
-      int lworkev = 3*ncv;
-      std::vector<double> workev (lworkev, 0.);
-
-      std::vector<double> eigenvalues_real (nev, 0.);
-      std::vector<double> eigenvalues_im (nev, 0.);
+      std::vector<double> eigenvalues_real (nev+1, 0.);
+      std::vector<double> eigenvalues_im (nev+1, 0.);
 
-      // call of ARPACK dneupd routine
-      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,
-              &iparam[0], &ipntr[0], workd, &workl[0], &lworkl, &info);
+      // call of ARPACK dseupd/dneupd routine
+      if (additional_data.symmetric)
+        {
+          std::vector<double> z (ldz*nev, 0.);
+          dseupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0],
+                  &z[0], &ldz, &sigmar, bmat, &n, which, &nev, &tol,
+                  &resid[0], &ncv, &v[0], &ldv,
+                  &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
+        }
+      else
+        {
+          std::vector<double> workev (3*ncv, 0.);
+          dneupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0],
+                  &eigenvalues_im[0], &v[0], &ldz, &sigmar, &sigmai,
+                  &workev[0], bmat, &n, which, &nev, &tol,
+                  &resid[0], &ncv, &v[0], &ldv,
+                  &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
+        }
 
       if (info == 1)
         {
-          Assert (false, ExcArpackInfoMaxIt(control().max_steps()));
+          Assert (false, ArpackExcArpackInfoMaxIt(control().max_steps()));
         }
       else if (info == 3)
         {
-          Assert (false, ExcArpackNoShifts());
+          Assert (false, ArpackExcArpackNoShifts());
         }
       else if (info!=0)
         {
-          Assert (false, ExcArpackInfodneupd(info));
+          if (additional_data.symmetric)
+            {
+              Assert (false, ArpackExcArpackInfodseupd(info));
+            }
+          else
+            Assert (false, ArpackExcArpackInfodneupd(info));
         }
 
-
-      const unsigned int n_eigenvecs = eigenvectors.size();
-      for (size_type i=0; i<n_eigenvecs; ++i)
+      for (unsigned int i=0; i<nev; ++i)
         for (unsigned int j=0; j<n; ++j)
           eigenvectors[i](j) = v[i*n+j];
 
-      delete[] workd;
-
-      AssertDimension (eigenvalues.size(), eigenvalues_real.size());
-      AssertDimension (eigenvalues.size(), eigenvalues_im.size());
-
-      for (size_type i=0; i<eigenvalues.size(); ++i)
+      for (unsigned int i=0; i<nev_const; ++i)
         eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
                                                eigenvalues_im[i]);
     }
index de67f19d17db134977a578fc8125c97e611da575..ba1a80db232f1d519f27741b6ae9e9e38ad066b3 100644 (file)
@@ -428,6 +428,11 @@ private:
                   << "Number of wanted eigenvalues " << arg1
                   << " is larger that the size of eigenvectors " << arg2);
 
+  DeclException2 (PArpackExcInvalidEigenvectorSizeNonsymmetric, int, int,
+                  << "To store the real and complex parts of " << arg1
+                  << " eigenvectors in real-valued vectors, their size (currently set to " << arg2
+                  << ") should be greater than or equal to " << arg1+1);
+
   DeclException2 (PArpackExcInvalidEigenvalueSize, int, int,
                   << "Number of wanted eigenvalues " << arg1
                   << " is larger that the size of eigenvalues " << arg2);
@@ -603,8 +608,14 @@ void PArpackSolver<VectorType>::solve
  const unsigned int                  n_eigenvalues)
 {
 
-  Assert (n_eigenvalues <= eigenvectors.size(),
-          PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size()));
+  if (additional_data.symmetric)
+    {
+      Assert (n_eigenvalues <= eigenvectors.size(),
+              PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size()));
+    }
+  else
+    Assert (n_eigenvalues+1 <= eigenvectors.size(),
+            PArpackExcInvalidEigenvectorSizeNonsymmetric(n_eigenvalues, eigenvectors.size()));
 
   Assert (n_eigenvalues <= eigenvalues.size(),
           PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
@@ -872,8 +883,8 @@ void PArpackSolver<VectorType>::solve
       double sigmar = shift_value; // real part of the shift
       double sigmai = 0.0; // imaginary part of the shift
 
-      std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
-      std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
+      std::vector<double> eigenvalues_real (n_eigenvalues+1, 0.);
+      std::vector<double> eigenvalues_im (n_eigenvalues+1, 0.);
 
       // call of ARPACK pdneupd routine
       if (additional_data.symmetric)
@@ -884,7 +895,7 @@ void PArpackSolver<VectorType>::solve
                  &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
       else
         pdneupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
-                 &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai,
+                 &eigenvalues_im[0], &v[0], &ldz, &sigmar, &sigmai,
                  &workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
                  &resid[0], &ncv, &v[0], &ldv,
                  &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info);
@@ -902,7 +913,7 @@ void PArpackSolver<VectorType>::solve
           AssertThrow (false, PArpackExcInfoPdneupd(info));
         }
 
-      for (size_type i=0; i<n_eigenvalues; ++i)
+      for (size_type i=0; i<nev; ++i)
         {
           eigenvectors[i] = 0.0;
           Assert (i*nloc + nloc <= v.size(), dealii::ExcInternalError() );
@@ -919,7 +930,7 @@ void PArpackSolver<VectorType>::solve
     }
 
   // Throw an error if the solver did not converge.
-  AssertThrow (iparam[4] == (int)n_eigenvalues,
+  AssertThrow (iparam[4] >= (int)n_eigenvalues,
                PArpackExcConvergedEigenvectors(n_eigenvalues,iparam[4]));
 
   // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
diff --git a/tests/arpack/arpack_advection_diffusion.cc b/tests/arpack/arpack_advection_diffusion.cc
new file mode 100644 (file)
index 0000000..c08966e
--- /dev/null
@@ -0,0 +1,356 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2016 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Authors: Toby D. Young, Polish Academy of Sciences,
+ *          Wolfgang Bangerth, Texas A&M University
+ *          Joscha Gedicke, U Heidelberg
+ *
+ * This file tests the non-symmetric interface to ARPACK for an advection-diffussion
+ * operator. The advection is chosen in such a way that we compute complex eigenvalues.
+ * The most critical case when we cut a complex conjugate pair is tested, i.e. the
+ * last eigenvalue is complex but the conjugate pair is not included in the range
+ * of eigenvalues we asked for. This is a particular case in the nonsymmetric arpack
+ * interface that needs to be taken care of.
+ *
+ * We test that the computed vectors are eigenvectors and mass-normal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_i*B*x_i = 1
+ *
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/lac/arpack_solver.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include <complex>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace dealii;
+
+
+template <int dim>
+class EigenvalueProblem
+{
+public:
+  EigenvalueProblem (unsigned int n_eigenvalues);
+  void run ();
+
+private:
+  void make_grid_and_dofs ();
+  void assemble_system ();
+  void solve ();
+
+  Triangulation<dim> triangulation;
+  FE_Q<dim>          fe;
+  DoFHandler<dim>    dof_handler;
+
+  SparsityPattern                     sparsity_pattern;
+  SparseMatrix<double>                stiffness_matrix, mass_matrix;
+  std::vector<Vector<double> >        arpack_vectors;
+  std::vector<Vector<double> >        eigenvectors;
+  std::vector<std::complex<double> >  eigenvalues;
+
+  ConstraintMatrix constraints;
+
+  unsigned int n_eigenvalues;
+};
+
+
+
+template <int dim>
+EigenvalueProblem<dim>::EigenvalueProblem (unsigned int n_eigenvalues)
+  :
+  fe (1),
+  dof_handler (triangulation),
+  n_eigenvalues(n_eigenvalues)
+{
+}
+
+
+
+template <int dim>
+void EigenvalueProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (5);
+  dof_handler.distribute_dofs (fe);
+
+  DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
+  constraints.close ();
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                           dof_handler.n_dofs(),
+                           dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+  stiffness_matrix.reinit (sparsity_pattern);
+  mass_matrix.reinit (sparsity_pattern);
+
+  eigenvalues.resize (n_eigenvalues);
+
+  arpack_vectors.resize (n_eigenvalues + 1);
+  for (unsigned int i=0; i<arpack_vectors.size (); ++i)
+    arpack_vectors[i].reinit (dof_handler.n_dofs ());
+
+  eigenvectors.resize(2*n_eigenvalues);
+  for (unsigned int i = 0; i < eigenvectors.size(); ++i)
+    eigenvectors[i].reinit (dof_handler.n_dofs ());
+
+}
+
+
+
+template <int dim>
+void EigenvalueProblem<dim>::assemble_system ()
+{
+  QGauss<dim>   quadrature_formula(2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_stiffness_matrix = 0;
+      cell_mass_matrix      = 0;
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        {
+          const Point<dim> cur_point = fe_values.quadrature_point(q_point);
+          Tensor<1,dim> advection;
+          advection[0] = 10.;
+          advection[1] = 10.*cur_point[0];
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              {
+                cell_stiffness_matrix (i, j)
+                += (fe_values.shape_grad (i, q_point) *
+                    fe_values.shape_grad (j, q_point)
+                    +
+                    (advection * fe_values.shape_grad (i, q_point)) *
+                    fe_values.shape_value (j, q_point)
+                   ) * fe_values.JxW (q_point);
+
+                cell_mass_matrix (i, j)
+                += (fe_values.shape_value (i, q_point) *
+                    fe_values.shape_value (j, q_point)
+                   ) * fe_values.JxW (q_point);
+              }
+        }
+
+      cell->get_dof_indices (local_dof_indices);
+
+      constraints
+      .distribute_local_to_global (cell_stiffness_matrix,
+                                   local_dof_indices,
+                                   stiffness_matrix);
+      constraints
+      .distribute_local_to_global (cell_mass_matrix,
+                                   local_dof_indices,
+                                   mass_matrix);
+    }
+
+  stiffness_matrix.compress (VectorOperation::add);
+  mass_matrix.compress (VectorOperation::add);
+
+}
+
+
+
+template <int dim>
+void EigenvalueProblem<dim>::solve ()
+{
+  SolverControl solver_control (dof_handler.n_dofs(), 1e-10);
+  SparseDirectUMFPACK inverse;
+  inverse.initialize (stiffness_matrix);
+  const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2;
+  ArpackSolver::AdditionalData additional_data(num_arnoldi_vectors,
+                                               ArpackSolver::largest_real_part);
+  ArpackSolver eigensolver (solver_control, additional_data);
+  arpack_vectors[0] = 1.;
+  eigensolver.set_initial_vector(arpack_vectors[0]);
+  eigensolver.solve (stiffness_matrix,
+                     mass_matrix,
+                     inverse,
+                     eigenvalues,
+                     arpack_vectors,
+                     eigenvalues.size());
+
+  // extract real and complex components of eigenvectors
+  for (unsigned int i = 0; i < n_eigenvalues; ++i)
+    {
+      eigenvectors[i] = arpack_vectors[i];
+      if (eigenvalues[i].imag() != 0.)
+        {
+          eigenvectors[i + eigenvalues.size()] = arpack_vectors[i+1];
+          if ( i+1 < eigenvalues.size())
+            {
+              eigenvectors[i+1] = arpack_vectors[i];
+              eigenvectors[i+1 + eigenvalues.size()] = arpack_vectors[i+1];
+              eigenvectors[i+1 + eigenvalues.size()] *= -1;
+              ++i;
+            }
+        }
+    }
+
+  // make sure that we have eigenvectors and they are mass-normal:
+  // a) (A*x_i-\lambda*B*x_i).L2() == 0
+  // b) x_i*B*x_i = 1
+  {
+    Vector<double> Ax(eigenvectors[0]), Bx(eigenvectors[0]);
+    Vector<double> Ay(eigenvectors[0]), By(eigenvectors[0]);
+    for (unsigned int i=0; i < n_eigenvalues; ++i)
+      {
+        stiffness_matrix.vmult(Ax,eigenvectors[i]);
+        stiffness_matrix.vmult(Ay,eigenvectors[i + n_eigenvalues]);
+        mass_matrix.vmult(Bx,eigenvectors[i]);
+        mass_matrix.vmult(By,eigenvectors[i + n_eigenvalues]);
+
+        Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
+        Ax.add(std::imag(eigenvalues[i]),By);
+        Ay.add(-1.0*std::real(eigenvalues[i]),By);
+        Ay.add(-1.0*std::imag(eigenvalues[i]),Bx);
+        Vector<double> tmpx(Ax), tmpy(Ay);
+        tmpx.scale(Ax);
+        tmpy.scale(Ay);
+        tmpx+=tmpy;
+        if (std::sqrt(tmpx.l1_norm()) > 1e-8)
+          deallog << "Returned vector " << i << " is not an eigenvector!"
+                  << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm())
+                  << std::endl;
+
+        const double tmp =
+          std::abs( eigenvectors[i] * Bx +  eigenvectors[i+n_eigenvalues] * By - 1.) +
+          std::abs( eigenvectors[i+n_eigenvalues] * Bx -  eigenvectors[i] * By);
+        if ( tmp > 1e-8)
+          deallog << "Eigenvector " << i << " is not normal! failing norm is"
+                  << tmp << std::endl;
+
+
+      }
+  }
+}
+
+bool my_compare(std::complex<double> a, std::complex<double> b)
+{
+  if (a.imag()==0.)
+    return a.imag() < a.imag();
+  else
+    return a.real() < b.real();
+}
+
+template <int dim>
+void EigenvalueProblem<dim>::run ()
+{
+  make_grid_and_dofs ();
+
+  assemble_system ();
+
+  solve ();
+
+  std::sort(eigenvalues.begin(), eigenvalues.end(), my_compare);
+
+  for (unsigned int i = 0; i < n_eigenvalues; ++i)
+    deallog << "      Eigenvalue " << i
+            << " : " << eigenvalues[i]
+            << std::endl;
+}
+
+
+int main (int argc, char **argv)
+{
+
+  try
+    {
+      using namespace dealii;
+
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog.threshold_double(1.e-10);
+
+
+      EigenvalueProblem<2> problem(4);
+      problem.run ();
+    }
+
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/arpack/arpack_advection_diffusion.output b/tests/arpack/arpack_advection_diffusion.output
new file mode 100644 (file)
index 0000000..492d63b
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::      Eigenvalue 0 : (31.3555,0.00000)
+DEAL::      Eigenvalue 1 : (39.6418,-4.53224)
+DEAL::      Eigenvalue 2 : (39.6418,4.53224)
+DEAL::      Eigenvalue 3 : (49.9700,-9.19446)
diff --git a/tests/arpack/parpack_advection_diffusion_petsc.cc b/tests/arpack/parpack_advection_diffusion_petsc.cc
new file mode 100644 (file)
index 0000000..1a9a79b
--- /dev/null
@@ -0,0 +1,434 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2016 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ * This file tests the non-symmetric interface to PARPACK for an advection-diffussion
+ * operator with PETSc mpi vectors.
+ *
+ * We test that the computed vectors are eigenvectors and mass-normal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_i*B*x_i = 1
+ *
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/parpack_solver.h>
+
+#include <fstream>
+#include <iostream>
+
+const unsigned int dim = 2;//run in 2d to save time
+
+const double eps = 1e-10;
+
+template <typename DoFHandlerType>
+std::vector<dealii::IndexSet>
+locally_owned_dofs_per_subdomain (const DoFHandlerType  &dof_handler)
+{
+  std::vector< dealii::types::subdomain_id > subdomain_association (dof_handler.n_dofs ());
+  dealii::DoFTools::get_subdomain_association (dof_handler, subdomain_association);
+
+  const unsigned int n_subdomains = 1 + (*max_element (subdomain_association.begin (),
+                                                       subdomain_association.end ()   ));
+
+  std::vector<dealii::IndexSet> index_sets (n_subdomains,dealii::IndexSet(dof_handler.n_dofs()));
+
+  // loop over subdomain_association and populate IndexSet when a
+  // change in subdomain ID is found
+  dealii::types::global_dof_index i_min          = 0;
+  dealii::types::global_dof_index this_subdomain = subdomain_association[0];
+
+  for (dealii::types::global_dof_index index = 1;
+       index < subdomain_association.size (); ++index)
+    {
+      //found index different from the current one
+      if (subdomain_association[index] != this_subdomain)
+        {
+          index_sets[this_subdomain].add_range (i_min, index);
+          i_min = index;
+          this_subdomain = subdomain_association[index];
+        }
+    }
+
+  // the very last element is of different index
+  if (i_min == subdomain_association.size () - 1)
+    {
+      index_sets[this_subdomain].add_index (i_min);
+    }
+
+  // otherwise there are at least two different indices
+  else
+    {
+      index_sets[this_subdomain].add_range (
+        i_min, subdomain_association.size ());
+    }
+
+  for (unsigned int i = 0; i < n_subdomains; i++)
+    index_sets[i].compress ();
+
+  return index_sets;
+} //locally_owned_dofs_per_subdomain
+
+class PETScInverse
+{
+public:
+  PETScInverse(const dealii::PETScWrappers::MatrixBase &A, dealii::SolverControl &cn,const MPI_Comm &mpi_communicator = PETSC_COMM_SELF):
+    solver(cn,mpi_communicator),
+    matrix(A),
+    preconditioner(matrix)
+  {
+
+  }
+
+  void vmult  ( dealii::PETScWrappers::MPI::Vector         &dst,
+                const dealii::PETScWrappers::MPI::Vector   &src) const
+  {
+    ;
+    solver.solve(matrix, dst, src,preconditioner);
+  }
+
+
+private:
+  mutable dealii::PETScWrappers::SolverGMRES solver;
+  const dealii::PETScWrappers::MatrixBase &matrix;
+  PETScWrappers::PreconditionBlockJacobi preconditioner;
+
+};
+
+void test ()
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 4;
+
+  MPI_Comm mpi_communicator = MPI_COMM_WORLD;
+  const unsigned int n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
+  const unsigned int this_mpi_process = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
+
+
+  dealii::Triangulation<dim> triangulation;
+  dealii::DoFHandler<dim> dof_handler(triangulation);
+  dealii::FE_Q<dim> fe(1);
+  dealii::ConstraintMatrix constraints;
+  dealii::IndexSet locally_owned_dofs;
+  dealii::IndexSet locally_relevant_dofs;
+
+  std::vector<dealii::PETScWrappers::MPI::Vector> eigenfunctions;
+  std::vector<dealii::PETScWrappers::MPI::Vector>  arpack_vectors;
+  std::vector<std::complex<PetscScalar>>          eigenvalues;
+  dealii::PETScWrappers::MPI::SparseMatrix        stiffness_matrix, mass_matrix;
+
+  dealii::GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (global_mesh_refinement_steps);
+
+  // we do not use metis but rather partition by hand below.
+  //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation);
+  {
+    const double x0 = -1.0;
+    const double x1 =  1.0;
+    const double dL = (x1-x0) / n_mpi_processes;
+
+    dealii::Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+    for (; cell!=endc; ++cell)
+      {
+        const dealii::Point<dim> &center = cell->center();
+        const double x = center[0];
+
+        const unsigned int id = std::floor ( (x-x0)/dL);
+        cell->set_subdomain_id (id);
+      }
+  }
+
+  dof_handler.distribute_dofs (fe);
+  dealii::DoFRenumbering::subdomain_wise (dof_handler);
+  std::vector<dealii::IndexSet> locally_owned_dofs_per_processor
+    = locally_owned_dofs_per_subdomain (dof_handler);
+  locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process];
+  locally_relevant_dofs.clear();
+  dealii::DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                                   locally_relevant_dofs);
+
+  constraints.clear();
+  constraints.reinit (locally_relevant_dofs);
+  dealii::DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  dealii::VectorTools::interpolate_boundary_values (dof_handler,
+                                                    0,
+                                                    dealii::ZeroFunction<dim> (),
+                                                    constraints);
+  constraints.close ();
+
+  dealii::CompressedSimpleSparsityPattern csp (locally_relevant_dofs);
+  // Fill in ignoring all cells that are not locally owned
+  dealii::DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                           constraints,
+                                           /* keep constrained dofs */ true);
+  std::vector<dealii::types::global_dof_index> n_locally_owned_dofs(n_mpi_processes);
+  for (unsigned int i = 0; i < n_mpi_processes; i++)
+    n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements();
+
+  dealii::SparsityTools::distribute_sparsity_pattern
+  (csp,
+   n_locally_owned_dofs,
+   mpi_communicator,
+   locally_relevant_dofs);
+
+  // Initialise the stiffness and mass matrices
+  stiffness_matrix.reinit (locally_owned_dofs,
+                           locally_owned_dofs,
+                           csp,
+                           mpi_communicator);
+
+  mass_matrix.reinit (locally_owned_dofs,
+                      locally_owned_dofs,
+                      csp,
+                      mpi_communicator);
+
+  eigenvalues.resize (number_of_eigenvalues);
+
+  arpack_vectors.resize (number_of_eigenvalues+1);
+  for (unsigned int i=0; i<arpack_vectors.size (); ++i)
+    arpack_vectors[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+  eigenfunctions.resize (2*number_of_eigenvalues);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+
+  // ready for assembly
+  stiffness_matrix = 0;
+  mass_matrix = 0;
+
+  dealii::QGauss<dim>   quadrature_formula(2);
+  dealii::FEValues<dim> fe_values (fe, quadrature_formula,
+                                   dealii::update_values |
+                                   dealii::update_gradients |
+                                   dealii::update_quadrature_points |
+                                   dealii::update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  dealii::FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+  dealii::FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<dealii::types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  typename dealii::DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    if (cell->subdomain_id() == this_mpi_process)
+      {
+        fe_values.reinit (cell);
+        cell_stiffness_matrix = 0;
+        cell_mass_matrix      = 0;
+
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            const Point<dim> cur_point = fe_values.quadrature_point(q_point);
+            Tensor<1,dim> advection;
+            advection[0] = 10.;
+            advection[1] = 10.*cur_point[0];
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                {
+                  cell_stiffness_matrix (i, j)
+                  += (fe_values.shape_grad (i, q_point) *
+                      fe_values.shape_grad (j, q_point)
+                      +
+                      (advection * fe_values.shape_grad (i, q_point)) *
+                      fe_values.shape_value (j, q_point)
+                     ) * fe_values.JxW (q_point);
+
+                  cell_mass_matrix (i, j)
+                  += (fe_values.shape_value (i, q_point) *
+                      fe_values.shape_value (j, q_point)
+                     ) * fe_values.JxW (q_point);
+                }
+
+            cell->get_dof_indices (local_dof_indices);
+
+            constraints
+            .distribute_local_to_global (cell_stiffness_matrix,
+                                         local_dof_indices,
+                                         stiffness_matrix);
+            constraints
+            .distribute_local_to_global (cell_mass_matrix,
+                                         local_dof_indices,
+                                         mass_matrix);
+          }
+      }
+
+  stiffness_matrix.compress (dealii::VectorOperation::add);
+  mass_matrix.compress (dealii::VectorOperation::add);
+
+  // test Arpack
+  {
+    std::vector<std::complex<double> > lambda(eigenfunctions.size());
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      eigenfunctions[i] = PetscScalar();
+
+    dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-10,/*log_history*/false,/*log_results*/false);
+    PETScInverse inverse(stiffness_matrix,solver_control,mpi_communicator);
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2;
+
+    dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector>::AdditionalData
+    additional_data(num_arnoldi_vectors,
+                    dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector>::largest_real_part,
+                    false);
+
+    dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector> eigensolver (solver_control,
+        mpi_communicator,
+        additional_data);
+    eigensolver.reinit(locally_owned_dofs);
+    arpack_vectors[0] = 1.;
+    eigensolver.set_initial_vector(arpack_vectors[0]);
+    eigensolver.solve (stiffness_matrix,
+                       mass_matrix,
+                       inverse,
+                       eigenvalues,
+                       arpack_vectors,
+                       eigenvalues.size());
+
+    // extract real and complex components of eigenvectors
+    for (unsigned int i = 0; i < eigenvalues.size(); ++i)
+      {
+        eigenfunctions[i] = arpack_vectors[i];
+        if (eigenvalues[i].imag() != 0.)
+          {
+            eigenfunctions[i + eigenvalues.size()] = arpack_vectors[i+1];
+            if ( i+1 < eigenvalues.size())
+              {
+                eigenfunctions[i+1] = arpack_vectors[i];
+                eigenfunctions[i+1 + eigenvalues.size()] = arpack_vectors[i+1];
+                eigenfunctions[i+1 + eigenvalues.size()] *= -1;
+                ++i;
+              }
+          }
+      }
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      dealii::deallog << eigenvalues[i] << std::endl;
+
+    // make sure that we have eigenvectors and they are mass-normal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_i*B*x_i=1
+    {
+      const double precision = 1e-7;
+      PETScWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      PETScWrappers::MPI::Vector Ay(eigenfunctions[0]), By(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenvalues.size(); ++i)
+        {
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          stiffness_matrix.vmult(Ay,eigenfunctions[i + eigenvalues.size()]);
+          mass_matrix.vmult(Bx,eigenfunctions[i]);
+          mass_matrix.vmult(By,eigenfunctions[i + eigenvalues.size()]);
+
+          Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
+          Ax.add(std::imag(eigenvalues[i]),By);
+          Ay.add(-1.0*std::real(eigenvalues[i]),By);
+          Ay.add(-1.0*std::imag(eigenvalues[i]),Bx);
+          PETScWrappers::MPI::Vector tmpx(Ax), tmpy(Ay);
+          tmpx.scale(Ax);
+          tmpy.scale(Ay);
+          tmpx+=tmpy;
+          if (std::sqrt(tmpx.l1_norm()) > precision)
+            deallog << "Returned vector " << i << " is not an eigenvector!"
+                    << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm())
+                    << std::endl;
+
+          const double tmp =
+            std::abs( eigenfunctions[i] * Bx +  eigenfunctions[i+eigenvalues.size()] * By - 1.) +
+            std::abs( eigenfunctions[i+eigenvalues.size()] * Bx -  eigenfunctions[i] * By);
+          if ( tmp > precision)
+            deallog << "Eigenvector " << i << " is not normal! failing norm is "
+                    << tmp << std::endl;
+        }
+    }
+  }
+
+
+  dof_handler.clear ();
+  dealii::deallog << "Ok"<<std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  dealii::deallog.attach(logfile,/*do not print job id*/false);
+  dealii::deallog.threshold_double(eps);
+
+  try
+    {
+      dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        test ();
+      }
+
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/arpack/parpack_advection_diffusion_petsc.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output b/tests/arpack/parpack_advection_diffusion_petsc.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..0744b5c
--- /dev/null
@@ -0,0 +1,5 @@
+DEAL::(32.2359,0.00000)
+DEAL::(40.3938,-4.48216)
+DEAL::(40.3938,4.48216)
+DEAL::(50.5706,-9.07560)
+DEAL::Ok
diff --git a/tests/arpack/parpack_advection_diffusion_trilinos.cc b/tests/arpack/parpack_advection_diffusion_trilinos.cc
new file mode 100644 (file)
index 0000000..df8ef1f
--- /dev/null
@@ -0,0 +1,427 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2016 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ * This file tests the non-symmetric interface to PARPACK for an advection-diffussion
+ * operator with Trilinos mpi vectors.
+ *
+ * We test that the computed vectors are eigenvectors and mass-normal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_i*B*x_i = 1
+ *
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <deal.II/lac/iterative_inverse.h>
+
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/parpack_solver.h>
+
+#include <fstream>
+#include <iostream>
+
+const unsigned int dim = 2;//run in 2d to save time
+
+using namespace dealii;
+
+const double eps = 1e-10;
+
+template <typename DoFHandlerType>
+std::vector<IndexSet>
+locally_owned_dofs_per_subdomain (const DoFHandlerType  &dof_handler)
+{
+  std::vector< types::subdomain_id > subdomain_association (dof_handler.n_dofs ());
+  DoFTools::get_subdomain_association (dof_handler, subdomain_association);
+
+  const unsigned int n_subdomains = 1 + (*max_element (subdomain_association.begin (),
+                                                       subdomain_association.end ()   ));
+
+  std::vector<IndexSet> index_sets (n_subdomains,IndexSet(dof_handler.n_dofs()));
+
+  // loop over subdomain_association and populate IndexSet when a
+  // change in subdomain ID is found
+  types::global_dof_index i_min          = 0;
+  types::global_dof_index this_subdomain = subdomain_association[0];
+
+  for (types::global_dof_index index = 1;
+       index < subdomain_association.size (); ++index)
+    {
+      //found index different from the current one
+      if (subdomain_association[index] != this_subdomain)
+        {
+          index_sets[this_subdomain].add_range (i_min, index);
+          i_min = index;
+          this_subdomain = subdomain_association[index];
+        }
+    }
+
+  // the very last element is of different index
+  if (i_min == subdomain_association.size () - 1)
+    {
+      index_sets[this_subdomain].add_index (i_min);
+    }
+
+  // otherwise there are at least two different indices
+  else
+    {
+      index_sets[this_subdomain].add_range (
+        i_min, subdomain_association.size ());
+    }
+
+  for (unsigned int i = 0; i < n_subdomains; i++)
+    index_sets[i].compress ();
+
+  return index_sets;
+} //locally_owned_dofs_per_subdomain
+
+
+void test ()
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 4;
+
+  MPI_Comm mpi_communicator = MPI_COMM_WORLD;
+  const unsigned int n_mpi_processes = Utilities::MPI::n_mpi_processes(mpi_communicator);
+  const unsigned int this_mpi_process = Utilities::MPI::this_mpi_process(mpi_communicator);
+
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+  FE_Q<dim> fe(1);
+  ConstraintMatrix constraints;
+  IndexSet locally_owned_dofs;
+  IndexSet locally_relevant_dofs;
+
+  std::vector<TrilinosWrappers::MPI::Vector> eigenfunctions;
+  std::vector<TrilinosWrappers::MPI::Vector>  arpack_vectors;
+  std::vector<std::complex<double>>              eigenvalues;
+  TrilinosWrappers::SparseMatrix             stiffness_matrix, mass_matrix;
+
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (global_mesh_refinement_steps);
+
+  // we do not use metis but rather partition by hand below.
+  //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation);
+  {
+    const double x0 = -1.0;
+    const double x1 =  1.0;
+    const double dL = (x1-x0) / n_mpi_processes;
+
+    Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+    for (; cell!=endc; ++cell)
+      {
+        const Point<dim> &center = cell->center();
+        const double x = center[0];
+
+        const unsigned int id = std::floor ( (x-x0)/dL);
+        cell->set_subdomain_id (id);
+      }
+  }
+
+
+  dof_handler.distribute_dofs (fe);
+  DoFRenumbering::subdomain_wise (dof_handler);
+  std::vector<IndexSet> locally_owned_dofs_per_processor
+    = locally_owned_dofs_per_subdomain (dof_handler);
+  locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process];
+  locally_relevant_dofs.clear();
+  DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                           locally_relevant_dofs);
+
+  constraints.clear();
+  constraints.reinit (locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            ZeroFunction<dim> (),
+                                            constraints);
+  constraints.close ();
+
+  CompressedSimpleSparsityPattern csp (locally_relevant_dofs);
+  // Fill in ignoring all cells that are not locally owned
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                   constraints,
+                                   /* keep constrained dofs */ true);
+  std::vector<types::global_dof_index> n_locally_owned_dofs(n_mpi_processes);
+  for (unsigned int i = 0; i < n_mpi_processes; i++)
+    n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements();
+
+  SparsityTools::distribute_sparsity_pattern
+  (csp,
+   n_locally_owned_dofs,
+   mpi_communicator,
+   locally_relevant_dofs);
+
+  // Initialise the stiffness and mass matrices
+  stiffness_matrix.reinit (locally_owned_dofs,
+                           locally_owned_dofs,
+                           csp,
+                           mpi_communicator);
+
+  mass_matrix.reinit (locally_owned_dofs,
+                      locally_owned_dofs,
+                      csp,
+                      mpi_communicator);
+
+  eigenvalues.resize (number_of_eigenvalues);
+
+  arpack_vectors.resize (number_of_eigenvalues+1);
+  for (unsigned int i=0; i<arpack_vectors.size (); ++i)
+    arpack_vectors[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+  eigenfunctions.resize (2*number_of_eigenvalues);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+  // ready for assembly
+  stiffness_matrix = 0;
+  mass_matrix = 0;
+
+  QGauss<dim>   quadrature_formula(2);
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values |
+                           update_gradients |
+                           update_quadrature_points |
+                           update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active (),
+  endc = dof_handler.end ();
+  for (; cell!=endc; ++cell)
+    if (cell->subdomain_id() == this_mpi_process)
+      {
+        fe_values.reinit (cell);
+        cell_stiffness_matrix = 0;
+        cell_mass_matrix      = 0;
+
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            const Point<dim> cur_point = fe_values.quadrature_point(q_point);
+            Tensor<1,dim> advection;
+            advection[0] = 10.;
+            advection[1] = 10.*cur_point[0];
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                {
+                  cell_stiffness_matrix (i, j)
+                  += (fe_values.shape_grad (i, q_point) *
+                      fe_values.shape_grad (j, q_point)
+                      +
+                      (advection * fe_values.shape_grad (i, q_point)) *
+                      fe_values.shape_value (j, q_point)
+                     ) * fe_values.JxW (q_point);
+
+                  cell_mass_matrix (i, j)
+                  += (fe_values.shape_value (i, q_point) *
+                      fe_values.shape_value (j, q_point)
+                     ) * fe_values.JxW (q_point);
+                }
+
+            cell->get_dof_indices (local_dof_indices);
+
+            constraints
+            .distribute_local_to_global (cell_stiffness_matrix,
+                                         local_dof_indices,
+                                         stiffness_matrix);
+            constraints
+            .distribute_local_to_global (cell_mass_matrix,
+                                         local_dof_indices,
+                                         mass_matrix);
+          }
+      }
+
+  stiffness_matrix.compress (VectorOperation::add);
+  mass_matrix.compress (VectorOperation::add);
+
+  // test Arpack
+  {
+    const double shift = 4.0;
+    std::vector<std::complex<double> > lambda(eigenfunctions.size());
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      eigenfunctions[i] = 0.;
+
+    SolverControl solver_control     (dof_handler.n_dofs(), 1e-9,/*log_history*/false,/*log_results*/false);
+    SolverControl solver_control_lin (dof_handler.n_dofs(), 1e-10,/*log_history*/false,/*log_results*/false);
+
+    PArpackSolver<TrilinosWrappers::MPI::Vector>::Shift<TrilinosWrappers::SparseMatrix> shifted_matrix(stiffness_matrix,mass_matrix,shift);
+    TrilinosWrappers::PreconditionIdentity preconditioner;
+    IterativeInverse<TrilinosWrappers::MPI::Vector > shift_and_invert;
+    shift_and_invert.initialize(shifted_matrix,preconditioner);
+    shift_and_invert.solver.select("gmres");
+    static ReductionControl inner_control_c(/*maxiter*/stiffness_matrix.m(),
+                                                       /*tolerance (global)*/ 0.0,
+                                                       /*reduce (w.r.t. initial)*/ 1.e-13);
+    shift_and_invert.solver.set_control(inner_control_c);
+
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2;
+
+    PArpackSolver<TrilinosWrappers::MPI::Vector>::AdditionalData
+    additional_data(num_arnoldi_vectors,
+                    PArpackSolver<TrilinosWrappers::MPI::Vector>::largest_real_part,
+                    false);
+
+    PArpackSolver<TrilinosWrappers::MPI::Vector> eigensolver (solver_control,
+                                                              mpi_communicator,
+                                                              additional_data);
+    eigensolver.reinit(locally_owned_dofs);
+    eigensolver.set_shift(shift);
+    arpack_vectors[0] = 1.;
+    eigensolver.set_initial_vector(arpack_vectors[0]);
+    // avoid output of iterative solver:
+    const unsigned int previous_depth = deallog.depth_file(0);
+    eigensolver.solve (stiffness_matrix,
+                       mass_matrix,
+                       shift_and_invert,
+                       eigenvalues,
+                       arpack_vectors,
+                       eigenvalues.size());
+    deallog.depth_file(previous_depth);
+
+    // extract real and complex components of eigenvectors
+    for (unsigned int i = 0; i < eigenvalues.size(); ++i)
+      {
+        eigenfunctions[i] = arpack_vectors[i];
+        if (eigenvalues[i].imag() != 0.)
+          {
+            eigenfunctions[i + eigenvalues.size()] = arpack_vectors[i+1];
+            if ( i+1 < eigenvalues.size())
+              {
+                eigenfunctions[i+1] = arpack_vectors[i];
+                eigenfunctions[i+1 + eigenvalues.size()] = arpack_vectors[i+1];
+                eigenfunctions[i+1 + eigenvalues.size()] *= -1;
+                ++i;
+              }
+          }
+      }
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      deallog << eigenvalues[i] << std::endl;
+
+    // make sure that we have eigenvectors and they are mass-normal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_i*B*x_i=1
+    {
+      const double precision = 1e-7;
+      TrilinosWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      TrilinosWrappers::MPI::Vector Ay(eigenfunctions[0]), By(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenvalues.size(); ++i)
+        {
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          stiffness_matrix.vmult(Ay,eigenfunctions[i + eigenvalues.size()]);
+          mass_matrix.vmult(Bx,eigenfunctions[i]);
+          mass_matrix.vmult(By,eigenfunctions[i + eigenvalues.size()]);
+
+          Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
+          Ax.add(std::imag(eigenvalues[i]),By);
+          Ay.add(-1.0*std::real(eigenvalues[i]),By);
+          Ay.add(-1.0*std::imag(eigenvalues[i]),Bx);
+          TrilinosWrappers::MPI::Vector tmpx(Ax), tmpy(Ay);
+          tmpx.scale(Ax);
+          tmpy.scale(Ay);
+          tmpx+=tmpy;
+          if (std::sqrt(tmpx.l1_norm()) > precision)
+            deallog << "Returned vector " << i << " is not an eigenvector!"
+                    << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm())
+                    << std::endl;
+
+          const double tmp =
+            std::abs( eigenfunctions[i] * Bx +  eigenfunctions[i+eigenvalues.size()] * By - 1.) +
+            std::abs( eigenfunctions[i+eigenvalues.size()] * Bx -  eigenfunctions[i] * By);
+          if ( tmp > precision)
+            deallog << "Eigenvector " << i << " is not normal! failing norm is "
+                    << tmp << std::endl;
+        }
+
+    }
+  }
+
+
+  dof_handler.clear ();
+  deallog << "Ok"<<std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile,/*do not print job id*/false);
+  deallog.threshold_double(eps);
+
+  try
+    {
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        test ();
+      }
+
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/arpack/parpack_advection_diffusion_trilinos.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output b/tests/arpack/parpack_advection_diffusion_trilinos.with_mpi=true.with_trilinos=true.with_arpack_with_parpack=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..0744b5c
--- /dev/null
@@ -0,0 +1,5 @@
+DEAL::(32.2359,0.00000)
+DEAL::(40.3938,-4.48216)
+DEAL::(40.3938,4.48216)
+DEAL::(50.5706,-9.07560)
+DEAL::Ok
index 498793a7421647a8d8e93b4f862a4a2d4c4ec92c..45be84b73931e13a9b55ab1ec0bcc99ab2853535 100644 (file)
@@ -1,6 +1,6 @@
 /* ---------------------------------------------------------------------
  *
- * Copyright (C) 2009 - 2015 by the deal.II authors
+ * Copyright (C) 2009 - 2016 by the deal.II authors
  *
  * This file is part of the deal.II library.
  *
  *
  * Authors: Toby D. Young, Polish Academy of Sciences,
  *          Wolfgang Bangerth, Texas A&M University
+ *
+ * This file tests the ARPACK interface for a symmetric operator taken from step-36.
+ *
+ * We test that the computed vectors are eigenvectors and mass-orthonormal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_j*B*x_i = \delta_{i,j}
+ *
  */
 
 #include "../tests.h"
@@ -217,7 +224,8 @@ namespace Step36
     inverse.initialize (stiffness_matrix);
     const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2;
     ArpackSolver::AdditionalData additional_data(num_arnoldi_vectors,
-                                                 ArpackSolver::largest_real_part);
+                                                 ArpackSolver::largest_magnitude,
+                                                 true);
     ArpackSolver eigensolver (solver_control, additional_data);
     eigensolver.solve (stiffness_matrix,
                        mass_matrix,
@@ -236,27 +244,25 @@ namespace Step36
           mass_matrix.vmult(Bx,eigenfunctions[i]);
 
           for (unsigned int j=0; j < eigenfunctions.size(); j++)
-            Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< 1e-8,
-                    ExcMessage("Eigenvectors " +
-                               Utilities::int_to_string(i) +
-                               " and " +
-                               Utilities::int_to_string(j) +
-                               " are not orthonormal!"
-                               " failing norm is " +
-                               Utilities::to_string(
-                                 std::abs( eigenfunctions[j] * Bx - (i==j) )
-                               )
-                              ));
+            if ( std::abs( eigenfunctions[j] * Bx - (i==j))> 1e-8)
+              deallog << "Eigenvectors " +
+                      Utilities::int_to_string(i) +
+                      " and " +
+                      Utilities::int_to_string(j) +
+                      " are not orthonormal!"
+                      " failing norm is " +
+                      Utilities::to_string(std::abs( eigenfunctions[j] * Bx - (i==j)))
+                      << std::endl;
 
           stiffness_matrix.vmult(Ax,eigenfunctions[i]);
           Ax.add(-1.0*std::real(eigenvalues[i]),Bx);
-          Assert (Ax.l2_norm() < 1e-8,
-                  ExcMessage("Returned vector " +
-                             Utilities::int_to_string(i) +
-                             " is not an eigenvector!"
-                             " L2 norm of the residual is " +
-                             Utilities::to_string(Ax.l2_norm())
-                            ));
+          if (Ax.l2_norm() > 1e-8)
+            deallog << "Returned vector " +
+                    Utilities::int_to_string(i) +
+                    " is not an eigenvector!"
+                    " L2 norm of the residual is " +
+                    Utilities::to_string(Ax.l2_norm())
+                    << std::endl;
         }
     }
     for (unsigned int i=0; i<eigenfunctions.size(); ++i)
index 92d86848149434b44f24645f12800b9cc1d99028..51131b531ff3737dfe90fa02aa82a0a2aa1b68d7 100644 (file)
@@ -1,3 +1,28 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2016 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * This file tests the PARPACK interface for a symmetric operator taken from step-36
+ * using PETSc mpi vectors.
+ *
+ * We test that the computed vectors are eigenvectors and mass-orthonormal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_j*B*x_i = \delta_{i,j}
+ *
+ */
+
 #include "../tests.h"
 
 #include <deal.II/base/logstream.h>
index b063152c534153d88ba367912c2c2211d5db9b44..c02481b66acb438e71e9e56eaccb3893d80b8c4c 100644 (file)
@@ -1,3 +1,28 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2016 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * This file tests the PARPACK interface for a symmetric operator taken from step-36
+ * using Trilinos mpi vectors.
+ *
+ * We test that the computed vectors are eigenvectors and mass-orthonormal, i.e.
+ * a) (A*x_i-\lambda*B*x_i).L2() == 0
+ * b) x_j*B*x_i = \delta_{i,j}
+ *
+ */
+
 #include "../tests.h"
 
 #include <deal.II/base/logstream.h>

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.