]> https://gitweb.dealii.org/ - dealii.git/commitdiff
extend pArpack interface to mode 1 and 2 4678/head
authorDenis Davydov <davydden@gmail.com>
Tue, 1 Aug 2017 18:48:59 +0000 (20:48 +0200)
committerDenis Davydov <davydden@gmail.com>
Fri, 4 Aug 2017 14:06:02 +0000 (16:06 +0200)
doc/news/changes/minor/20170801DenisDavydov [new file with mode: 0644]
include/deal.II/lac/parpack_solver.h
tests/arpack/step-36_parpack_mf_02.cc [new file with mode: 0644]
tests/arpack/step-36_parpack_mf_02.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output [new file with mode: 0644]
tests/arpack/step-36_parpack_mf_03.cc [new file with mode: 0644]
tests/arpack/step-36_parpack_mf_03.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output [new file with mode: 0644]
tests/slepc/step-36_parallel.cc
tests/slepc/step-36_parallel_02.cc [new file with mode: 0644]
tests/slepc/step-36_parallel_02.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output [new file with mode: 0644]
tests/slepc/step-36_parallel_03.cc [new file with mode: 0644]
tests/slepc/step-36_parallel_03.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170801DenisDavydov b/doc/news/changes/minor/20170801DenisDavydov
new file mode 100644 (file)
index 0000000..9751b27
--- /dev/null
@@ -0,0 +1,4 @@
+New: Extend pArpack solver to cover mode 1 (standard eigenvalue problem) and
+mode 2 (generalized eigenvalue problem without spectral transformation).
+<br>
+(Denis Davydov, 2017/08/01)
index 5548deb3f5e6631d6efc848d565294e2af7626d0..47fc83538cbf4c98b5842a8e311d16241c64af6f 100644 (file)
@@ -108,10 +108,17 @@ extern "C" {
  * eigenvector/eigenvalue pairs to solve for. Here, <code>lambda</code> is a
  * vector that will contain the eigenvalues computed, <code>x</code> a vector
  * of objects of type <code>V</code> that will contain the eigenvectors
- * computed. <code>OP</code> is an inverse operation for the matrix <code>A -
+ * computed.
+ *
+ * Currently, only three modes of (P)Arpack are implemented. In mode 3 (default),
+ * <code>OP</code> is an inverse operation for the matrix <code>A -
  * sigma * B</code>, where <code> sigma </code> is a shift value, set to zero
- * by default. Note that (P)Arpack supports other transformations, but currently
- * this class implements only shift-and-invert mode.
+ * by default. Whereas in mode 2, <code>OP</code> is an inverse of <code>M</code>.
+ * Finally, mode 1 corresponds to standard eigenvalue problem without
+ * spectral transformation $Ax=\lambda x$.
+ * The mode can be specified via AdditionalData object. Note that for
+ * shift-and-invert (mode=3), the sought eigenpairs are those after the
+ * spectral transformation is applied.
  *
  * The <code>OP</code> can be specified by using a LinearOperator:
  * @code
@@ -125,9 +132,6 @@ extern "C" {
  *   const auto op_shift_invert = inverse_operator(op_shift, cg, PreconditionIdentity ());
  * @endcode
  *
- * Through the AdditionalData the user can specify some of the parameters to
- * be set.
- *
  * The class is intended to be used with MPI and can work on arbitrary vector
  * and matrix distributed classes.  Both symmetric and non-symmetric
  * <code>A</code> are supported.
@@ -137,7 +141,7 @@ extern "C" {
  * also how to set the parameters appropriately please take a look into the
  * PARPACK manual.
  *
- * @author Denis Davydov, 2015.
+ * @author Denis Davydov, 2015, 2017
  */
 template <typename VectorType>
 class PArpackSolver : public Subscriptor
@@ -261,10 +265,12 @@ public:
     const unsigned int number_of_arnoldi_vectors;
     const WhichEigenvalues eigenvalue_of_interest;
     const bool symmetric;
+    const int mode;
     AdditionalData(
       const unsigned int number_of_arnoldi_vectors = 15,
       const WhichEigenvalues eigenvalue_of_interest = largest_magnitude,
-      const bool symmetric = false);
+      const bool symmetric = false,
+      const int mode = 3);
   };
 
   /**
@@ -307,6 +313,9 @@ public:
    * Set shift @p sigma for shift-and-invert spectral transformation.
    *
    * If this function is not called, the shift is assumed to be zero.
+   *
+   * @note only relevant for <code>mode=3</code> (see the general documentation of this
+   * class for a definition of what the different modes are).
    */
   void set_shift(const std::complex<double> sigma);
 
@@ -314,6 +323,10 @@ public:
    * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling
    * the <code>pd(n/s)eupd</code> and <code>pd(n/s)aupd</code> functions of
    * PARPACK.
+   *
+   * In <code>mode=3</code>, @p inverse should correspond to $[A-\sigma B]^{-1}$,
+   * whereas in <code>mode=2</code> it should represent $B^{-1}$. For
+   * <code>mode=1</code> both @p B and @p inverse are ignored.
    */
   template <typename MatrixType1,
             typename MatrixType2, typename INVERSE>
@@ -537,11 +550,13 @@ template <typename VectorType>
 PArpackSolver<VectorType>::AdditionalData::
 AdditionalData (const unsigned int     number_of_arnoldi_vectors,
                 const WhichEigenvalues eigenvalue_of_interest,
-                const bool             symmetric)
+                const bool             symmetric,
+                const int              mode)
   :
   number_of_arnoldi_vectors(number_of_arnoldi_vectors),
   eigenvalue_of_interest(eigenvalue_of_interest),
-  symmetric(symmetric)
+  symmetric(symmetric),
+  mode(mode)
 {
   //Check for possible options for symmetric problems
   if (symmetric)
@@ -555,6 +570,8 @@ AdditionalData (const unsigned int     number_of_arnoldi_vectors,
       Assert(eigenvalue_of_interest!=smallest_imaginary_part,
              ExcMessage("'smallest imaginary part' can only be used for non-symmetric problems!"));
     }
+  Assert (mode >= 1 && mode <= 3,
+          ExcMessage("Currently, only modes 1, 2 and 3 are supported."));
 }
 
 template <typename VectorType>
@@ -676,7 +693,7 @@ void PArpackSolver<VectorType>::reinit(const IndexSet &locally_owned_dofs,
 template <typename VectorType>
 template <typename MatrixType1,typename MatrixType2, typename INVERSE>
 void PArpackSolver<VectorType>::solve
-(const MatrixType1                  &/*system_matrix*/,
+(const MatrixType1                  &system_matrix,
  const MatrixType2                  &mass_matrix,
  const INVERSE                      &inverse,
  std::vector<std::complex<double> > &eigenvalues,
@@ -709,11 +726,8 @@ void PArpackSolver<VectorType>::solve
   Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
           PArpackExcSmallNumberofArnoldiVectors(
             additional_data.number_of_arnoldi_vectors, n_eigenvalues));
-  // ARPACK mode for dnaupd, here only
-  //  Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
-  //c           ===> OP = (inv[K - sigma*M])*M  and  B = M.
-  //c           ===> Shift-and-Invert mode
-  int mode = 3;
+
+  int mode = additional_data.mode;
 
   // reverse communication parameter
   // must be zero on the first call to pdnaupd
@@ -721,7 +735,9 @@ void PArpackSolver<VectorType>::solve
 
   // 'G' generalized eigenvalue problem
   // 'I' standard eigenvalue problem
-  char bmat[2] = "G";
+  char bmat[2];
+  bmat[0] = (mode == 1 ) ? 'I' : 'G';
+  bmat[1] = '\0';
 
   // Specify the eigenvalues of interest, possible parameters:
   // "LA" algebraically largest
@@ -781,9 +797,9 @@ void PArpackSolver<VectorType>::solve
   iparam[3] = 1;
 
   // Sets the mode of dsaupd:
-  // 1 is exact shifting,
-  // 2 is user-supplied shifts,
-  // 3 is shift-invert mode,
+  // 1 is A*x=lambda*x, OP = A, B = I
+  // 2 is A*x = lambda*M*x, OP = inv[M]*A, B = M
+  // 3 is shift-invert mode, OP = inv[A-sigma*M]*M, B = M
   // 4 is buckling mode,
   // 5 is Cayley mode.
 
@@ -815,192 +831,157 @@ void PArpackSolver<VectorType>::solve
                  &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
                  &workd[0], &workl[0], &lworkl, &info);
 
+      AssertThrow (info == 0, PArpackExcInfoPdnaupd(info));
+
+      // if we converge, we shall not modify anything in work arrays!
       if (ido == 99)
         break;
 
-      switch (mode)
-        {
-//        OP = (inv[K - sigma*M])*M
-        case 3:
+      // IPNTR(1) is the pointer into WORKD for X,
+      // IPNTR(2) is the pointer into WORKD for Y.
+      const int shift_x = ipntr[0]-1;
+      const int shift_y = ipntr[1]-1;
+      Assert (shift_x>=0, dealii::ExcInternalError() );
+      Assert (shift_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
+      Assert (shift_y>=0, dealii::ExcInternalError() );
+      Assert (shift_y+nloc <= (int)workd.size(), dealii::ExcInternalError() );
+
+      src = 0.;
+
+      // switch based on both ido and mode
+      if ((ido == -1) ||
+          (ido == 1 && mode<3))
+        // compute  Y = OP * X
         {
-          switch (ido)
-            {
-            case -1:
-              // compute  Y = OP * X  where
-              // IPNTR(1) is the pointer into WORKD for X,
-              // IPNTR(2) is the pointer into WORKD for Y.
+          src.add (nloc,
+                   &local_indices[0],
+                   &workd[0]+shift_x );
+          src.compress (VectorOperation::add);
+
+          if (mode == 3)
+            // OP = inv[K - sigma*M]*M
             {
-              const int shift_x = ipntr[0]-1;
-              const int shift_y = ipntr[1]-1;
-              Assert (shift_x>=0, dealii::ExcInternalError() );
-              Assert (shift_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-              Assert (shift_y>=0, dealii::ExcInternalError() );
-              Assert (shift_y+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-
-              src = 0.0;
-              src.add (nloc,
-                       &local_indices[0],
-                       &workd[0]+shift_x );
-              src.compress (VectorOperation::add);
-
-              // multiplication with mass matrix M
               mass_matrix.vmult(tmp, src);
-              // solving linear system
               inverse.vmult(dst,tmp);
-
-              // store the result
-              dst.extract_subvector_to (local_indices.begin(),
-                                        local_indices.end(),
-                                        &workd[0]+shift_y  );
             }
-            break;
-
-            case  1:
-              // compute  Y = OP * X where
-              // IPNTR(1) is the pointer into WORKD for X,
-              // IPNTR(2) is the pointer into WORKD for Y.
-              // In mode 3,4 and 5, the vector B * X is already
-              // available in WORKD(ipntr(3)).  It does not
-              // need to be recomputed in forming OP * X.
+          else if (mode == 2)
+            // OP = inv[M]*K
             {
-              const int shift_x   = ipntr[0]-1;
-              const int shift_y   = ipntr[1]-1;
-              const int shift_b_x = ipntr[2]-1;
-
-              Assert (shift_x>=0, dealii::ExcInternalError() );
-              Assert (shift_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-              Assert (shift_y>=0, dealii::ExcInternalError() );
-              Assert (shift_y+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-              Assert (shift_b_x>=0, dealii::ExcInternalError() );
-              Assert (shift_b_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-              Assert (shift_y>=0, dealii::ExcInternalError() );
-              Assert (shift_y+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-
-              src = 0.0; // B*X
-              src.add (nloc,
-                       &local_indices[0],
-                       &workd[0]+shift_b_x );
-
-              tmp = 0.0; // X
-              tmp.add (nloc,
-                       &local_indices[0],
-                       &workd[0]+shift_x);
-
-              src.compress (VectorOperation::add);
-              tmp.compress (VectorOperation::add);
-
-              // solving linear system
-              inverse.vmult(dst,src);
-
-              // store the result
-              dst.extract_subvector_to (local_indices.begin(),
+              system_matrix.vmult(tmp, src);
+              // store M*X in X
+              tmp.extract_subvector_to (local_indices.begin(),
                                         local_indices.end(),
-                                        &workd[0]+shift_y  );
-
+                                        &workd[0]+shift_x);
+              inverse.vmult(dst,tmp);
             }
-            break;
-
-            case  2:
-              // compute  Y = B * X  where
-              // IPNTR(1) is the pointer into WORKD for X,
-              // IPNTR(2) is the pointer into WORKD for Y.
+          else if (mode == 1)
             {
+              system_matrix.vmult(dst, src);
+            }
+          else
+            AssertThrow (false, PArpackExcMode(mode));
 
-              const int shift_x = ipntr[0]-1;
-              const int shift_y = ipntr[1]-1;
-              Assert (shift_x>=0, dealii::ExcInternalError() );
-              Assert (shift_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-              Assert (shift_y>=0, dealii::ExcInternalError() );
-              Assert (shift_y+nloc <= (int)workd.size(), dealii::ExcInternalError() );
-
-              src = 0.0;
-              src.add (nloc,
-                       &local_indices[0],
-                       &workd[0]+shift_x );
-              src.compress (VectorOperation::add);
-
-              // Multiplication with mass matrix M
-              mass_matrix.vmult(dst, src);
-
-              // store the result
-              dst.extract_subvector_to (local_indices.begin(),
-                                        local_indices.end(),
-                                        &workd[0]+shift_y);
+        }
+      else if (ido == 1 && mode >= 3)
+        // compute  Y = OP * X for mode 3, 4 and 5, where
+        // the vector B * X is already available in WORKD(ipntr(3)).
+        {
+          const int shift_b_x = ipntr[2]-1;
+          Assert (shift_b_x>=0, dealii::ExcInternalError() );
+          Assert (shift_b_x+nloc <= (int)workd.size(), dealii::ExcInternalError() );
+
+          // B*X
+          src.add (nloc,
+                   &local_indices[0],
+                   &workd[0]+shift_b_x );
+          src.compress (VectorOperation::add);
+
+          // solving linear system
+          Assert (mode == 3, ExcNotImplemented());
+          inverse.vmult(dst,src);
+        }
+      else if (ido == 2)
+        // compute  Y = B * X
+        {
+          src.add (nloc,
+                   &local_indices[0],
+                   &workd[0]+shift_x );
+          src.compress (VectorOperation::add);
 
+          // Multiplication with mass matrix M
+          if (mode == 1)
+            {
+              dst = src;
             }
-            break;
-
-            default:
-              AssertThrow (false, PArpackExcIdo(ido));
-              break;
+          else
+            // mode 2,3 and 5 have B=M
+            {
+              mass_matrix.vmult(dst, src);
             }
         }
-        break;
-        default:
-          AssertThrow (false, PArpackExcMode(mode));
-          break;
-        }
-    }
+      else
+        AssertThrow (false, PArpackExcIdo(ido));
+      // Note: IDO = 3 does not appear to be required for currently
+      // implemented modes
 
-  if (info<0)
-    {
-      AssertThrow (false, PArpackExcInfoPdnaupd(info));
-    }
-  else
-    {
-      // 1 - compute eigenvectors,
-      // 0 - only eigenvalues
-      int rvec = 1;
+      // store the result
+      dst.extract_subvector_to (local_indices.begin(),
+                                local_indices.end(),
+                                &workd[0]+shift_y);
+    } // end of pd*aupd_ loop
 
-      // which eigenvectors
-      char howmany[4] = "All";
+  // 1 - compute eigenvectors,
+  // 0 - only eigenvalues
+  int rvec = 1;
 
-      std::vector<double> eigenvalues_real (n_eigenvalues+1, 0.);
-      std::vector<double> eigenvalues_im (n_eigenvalues+1, 0.);
+  // which eigenvectors
+  char howmany[4] = "All";
 
-      // call of ARPACK pdneupd routine
-      if (additional_data.symmetric)
-        pdseupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
-                 &z[0], &ldz, &sigmar,
-                 bmat, &n_inside_arpack, which, &nev, &tol,
-                 &resid[0], &ncv, &v[0], &ldv,
-                 &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], &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);
+  std::vector<double> eigenvalues_real (n_eigenvalues+1, 0.);
+  std::vector<double> eigenvalues_im (n_eigenvalues+1, 0.);
 
-      if (info == 1)
-        {
-          AssertThrow (false, PArpackExcInfoMaxIt(control().max_steps()));
-        }
-      else if (info == 3)
-        {
-          AssertThrow (false, PArpackExcNoShifts(1));
-        }
-      else if (info!=0)
-        {
-          AssertThrow (false, PArpackExcInfoPdneupd(info));
-        }
+  // call of ARPACK pdneupd routine
+  if (additional_data.symmetric)
+    pdseupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
+             &z[0], &ldz, &sigmar,
+             bmat, &n_inside_arpack, which, &nev, &tol,
+             &resid[0], &ncv, &v[0], &ldv,
+             &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], &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);
 
-      for (int i=0; i<nev; ++i)
-        {
-          eigenvectors[i] = 0.0;
-          Assert (i*nloc + nloc <= (int)v.size(), dealii::ExcInternalError() );
+  if (info == 1)
+    {
+      AssertThrow (false, PArpackExcInfoMaxIt(control().max_steps()));
+    }
+  else if (info == 3)
+    {
+      AssertThrow (false, PArpackExcNoShifts(1));
+    }
+  else if (info!=0)
+    {
+      AssertThrow (false, PArpackExcInfoPdneupd(info));
+    }
 
-          eigenvectors[i].add (nloc,
-                               &local_indices[0],
-                               &v[i*nloc] );
-          eigenvectors[i].compress (VectorOperation::add);
-        }
+  for (int i=0; i<nev; ++i)
+    {
+      eigenvectors[i] = 0.0;
+      Assert (i*nloc + nloc <= (int)v.size(), dealii::ExcInternalError() );
 
-      for (size_type i=0; i<n_eigenvalues; ++i)
-        eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
-                                               eigenvalues_im[i]);
+      eigenvectors[i].add (nloc,
+                           &local_indices[0],
+                           &v[i*nloc] );
+      eigenvectors[i].compress (VectorOperation::add);
     }
 
+  for (size_type i=0; i<n_eigenvalues; ++i)
+    eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
+                                           eigenvalues_im[i]);
+
   // Throw an error if the solver did not converge.
   AssertThrow (iparam[4] >= (int)n_eigenvalues,
                PArpackExcConvergedEigenvectors(n_eigenvalues,iparam[4]));
diff --git a/tests/arpack/step-36_parpack_mf_02.cc b/tests/arpack/step-36_parpack_mf_02.cc
new file mode 100644 (file)
index 0000000..7b425e6
--- /dev/null
@@ -0,0 +1,249 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2017 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.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Same as step-36_parpack_mf but solve for largest eigenvalues of Laplace.
+ * For the same problem slepc/step-36_parallel_02.cc produces:
+ * 6099.84
+ * 6035.05
+ * 6035.05
+ * 5970.26
+ * 5931.69
+ */
+
+#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/distributed/tria.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/matrix_free/operators.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/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;
+
+using namespace dealii;
+
+const double eps = 1e-10;
+
+const unsigned int fe_degree = 1;
+
+void test ()
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 5;
+
+  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);
+
+  parallel::distributed::Triangulation<dim> triangulation (mpi_communicator);
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (global_mesh_refinement_steps);
+
+
+  DoFHandler<dim> dof_handler(triangulation);
+  FE_Q<dim> fe(fe_degree);
+  dof_handler.distribute_dofs (fe);
+
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                           locally_relevant_dofs);
+  ConstraintMatrix constraints;
+  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 ();
+
+  std::shared_ptr<MatrixFree<dim,double> > mf_data(new MatrixFree<dim,double> ());
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,double>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,double>::AdditionalData::partition_color;
+    data.mapping_update_flags = update_values | update_gradients | update_JxW_values;
+    mf_data->reinit (dof_handler, constraints, quad, data);
+  }
+
+  std::vector<LinearAlgebra::distributed::Vector<double> > eigenfunctions;
+  std::vector<double>                                      eigenvalues;
+  MatrixFreeOperators::MassOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > mass;
+  MatrixFreeOperators::LaplaceOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > laplace;
+  mass.initialize(mf_data);
+  laplace.initialize(mf_data);
+
+  eigenfunctions.resize (number_of_eigenvalues);
+  eigenvalues.resize (number_of_eigenvalues);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    mf_data->initialize_dof_vector (eigenfunctions[i]);
+
+  // test PArpack with matrix-free
+  {
+    std::vector<std::complex<double> > lambda(number_of_eigenvalues);
+
+    // set up iterative inverse
+    static ReductionControl inner_control_c(dof_handler.n_dofs(), 0.0, 1.e-14);
+
+    typedef LinearAlgebra::distributed::Vector<double> VectorType;
+    SolverCG<VectorType> solver_c(inner_control_c);
+    PreconditionIdentity preconditioner;
+    const auto invert =
+      inverse_operator(linear_operator<VectorType>(mass),
+                       solver_c,
+                       preconditioner);
+
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 10;
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::AdditionalData
+    additional_data(num_arnoldi_vectors,
+                    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::largest_magnitude,
+                    true,
+                    2);
+
+    SolverControl solver_control(
+      dof_handler.n_dofs(), 1e-9, /*log_history*/ false, /*log_results*/ false);
+
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> > eigensolver(
+      solver_control, mpi_communicator, additional_data);
+
+    eigensolver.reinit(eigenfunctions[0]);
+    // make sure initial vector is orthogonal to the space due to constraints
+    {
+      LinearAlgebra::distributed::Vector<double> init_vector;
+      mf_data->initialize_dof_vector(init_vector);
+      for (auto it = init_vector.begin(); it != init_vector.end(); ++it)
+        *it = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+      constraints.set_zero(init_vector);
+      eigensolver.set_initial_vector(init_vector);
+    }
+    // avoid output of iterative solver:
+    const unsigned int previous_depth = deallog.depth_file(0);
+    eigensolver.solve (laplace,
+                       mass,
+                       invert,
+                       lambda,
+                       eigenfunctions,
+                       eigenvalues.size());
+    deallog.depth_file(previous_depth);
+
+    for (unsigned int i = 0; i < lambda.size(); i++)
+      eigenvalues[i] = lambda[i].real();
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      deallog << eigenvalues[i] << std::endl;
+
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_j*B*x_i=\delta_{ij}
+    {
+      const double precision = 1e-7;
+      LinearAlgebra::distributed::Vector<double> Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          mass.vmult(Bx,eigenfunctions[i]);
+
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            {
+              const double err = std::abs( eigenfunctions[j] * Bx - (i==j));
+              Assert( err< precision,
+                      ExcMessage("Eigenvectors " +
+                                 Utilities::int_to_string(i) +
+                                 " and " +
+                                 Utilities::int_to_string(j) +
+                                 " are not orthonormal: " +
+                                 std::to_string(err)));
+            }
+
+          laplace.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*eigenvalues[i],Bx);
+          const double err = Ax.l2_norm();
+          Assert (err < precision,
+                  ExcMessage("Returned vector " +
+                             Utilities::int_to_string(i) +
+                             " is not an eigenvector: " +
+                             std::to_string(err)));
+        }
+    }
+  }
+
+
+  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/step-36_parpack_mf_02.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output b/tests/arpack/step-36_parpack_mf_02.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..cd93cf8
--- /dev/null
@@ -0,0 +1,6 @@
+DEAL::5931.69
+DEAL::5970.26
+DEAL::6035.05
+DEAL::6035.05
+DEAL::6099.84
+DEAL::Ok
diff --git a/tests/arpack/step-36_parpack_mf_03.cc b/tests/arpack/step-36_parpack_mf_03.cc
new file mode 100644 (file)
index 0000000..0bf3806
--- /dev/null
@@ -0,0 +1,247 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2017 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.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Same as step-36_parpack_mf_2 but solve for largest eigenvalues of SHEP Laplace.
+ * For the same problem slepc/step-36_parallel_03.cc produces:
+ * DEAL::3.98719
+ * DEAL::3.98719
+ * DEAL::3.97768
+ * DEAL::3.97768
+ * DEAL::3.96194
+ */
+
+#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/distributed/tria.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/matrix_free/operators.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/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;
+
+using namespace dealii;
+
+const double eps = 1e-10;
+
+const unsigned int fe_degree = 1;
+
+void test ()
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 5;
+
+  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);
+
+  parallel::distributed::Triangulation<dim> triangulation (mpi_communicator);
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (global_mesh_refinement_steps);
+
+
+  DoFHandler<dim> dof_handler(triangulation);
+  FE_Q<dim> fe(fe_degree);
+  dof_handler.distribute_dofs (fe);
+
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                           locally_relevant_dofs);
+  ConstraintMatrix constraints;
+  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 ();
+
+  std::shared_ptr<MatrixFree<dim,double> > mf_data(new MatrixFree<dim,double> ());
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,double>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,double>::AdditionalData::partition_color;
+    data.mapping_update_flags = update_values | update_gradients | update_JxW_values;
+    mf_data->reinit (dof_handler, constraints, quad, data);
+  }
+
+  std::vector<LinearAlgebra::distributed::Vector<double> > eigenfunctions;
+  std::vector<double>                                      eigenvalues;
+  MatrixFreeOperators::MassOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > mass;
+  MatrixFreeOperators::LaplaceOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > laplace;
+  mass.initialize(mf_data);
+  laplace.initialize(mf_data);
+
+  eigenfunctions.resize (number_of_eigenvalues);
+  eigenvalues.resize (number_of_eigenvalues);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    mf_data->initialize_dof_vector (eigenfunctions[i]);
+
+  // test PArpack with matrix-free
+  {
+    std::vector<std::complex<double> > lambda(number_of_eigenvalues);
+
+    // set up iterative inverse
+    static ReductionControl inner_control_c(dof_handler.n_dofs(), 0.0, 1.e-14);
+
+    typedef LinearAlgebra::distributed::Vector<double> VectorType;
+    SolverCG<VectorType> solver_c(inner_control_c);
+    PreconditionIdentity preconditioner;
+    const auto invert =
+      inverse_operator(linear_operator<VectorType>(mass),
+                       solver_c,
+                       preconditioner);
+
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 10;
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::AdditionalData
+    additional_data(num_arnoldi_vectors,
+                    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::largest_magnitude,
+                    true,
+                    1);
+
+    SolverControl solver_control(
+      dof_handler.n_dofs(), 1e-9, /*log_history*/ false, /*log_results*/ false);
+
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> > eigensolver(
+      solver_control, mpi_communicator, additional_data);
+
+    eigensolver.reinit(eigenfunctions[0]);
+    // make sure initial vector is orthogonal to the space due to constraints
+    {
+      LinearAlgebra::distributed::Vector<double> init_vector;
+      mf_data->initialize_dof_vector(init_vector);
+      for (auto it = init_vector.begin(); it != init_vector.end(); ++it)
+        *it = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+      constraints.set_zero(init_vector);
+      eigensolver.set_initial_vector(init_vector);
+    }
+    // avoid output of iterative solver:
+    const unsigned int previous_depth = deallog.depth_file(0);
+    eigensolver.solve (laplace,
+                       mass,
+                       invert,
+                       lambda,
+                       eigenfunctions,
+                       eigenvalues.size());
+    deallog.depth_file(previous_depth);
+
+    for (unsigned int i = 0; i < lambda.size(); i++)
+      eigenvalues[i] = lambda[i].real();
+
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      deallog << eigenvalues[i] << std::endl;
+
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*x_i).L2() == 0
+    // b) x_j*x_i=\delta_{ij}
+    {
+      const double precision = 1e-7;
+      LinearAlgebra::distributed::Vector<double> Ax(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            {
+              const double err = std::abs( eigenfunctions[j] * eigenfunctions[i] - (i==j));
+              Assert( err< precision,
+                      ExcMessage("Eigenvectors " +
+                                 Utilities::int_to_string(i) +
+                                 " and " +
+                                 Utilities::int_to_string(j) +
+                                 " are not orthonormal: " +
+                                 std::to_string(err)));
+            }
+
+          laplace.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*eigenvalues[i],eigenfunctions[i]);
+          const double err = Ax.l2_norm();
+          Assert (err < precision,
+                  ExcMessage("Returned vector " +
+                             Utilities::int_to_string(i) +
+                             " is not an eigenvector: " +
+                             std::to_string(err)));
+        }
+    }
+  }
+
+
+  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/step-36_parpack_mf_03.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output b/tests/arpack/step-36_parpack_mf_03.with_mpi=true.with_p4est=true.with_arpack_with_parpack=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..56aa2cf
--- /dev/null
@@ -0,0 +1,6 @@
+DEAL::3.96194
+DEAL::3.97768
+DEAL::3.97768
+DEAL::3.98719
+DEAL::3.98719
+DEAL::Ok
index 3c332bae669c73c258b1aee434e3461862c8a213..42a390997ac125e62e0200a7b32751b28be31641 100644 (file)
@@ -1,3 +1,20 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2017 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.
+//
+// ---------------------------------------------------------------------
+
+// Solve for smallest eigenvalues of Laplace operator in Step-36 tutorial
+
 #include "../tests.h"
 
 #include <deal.II/base/logstream.h>
diff --git a/tests/slepc/step-36_parallel_02.cc b/tests/slepc/step-36_parallel_02.cc
new file mode 100644 (file)
index 0000000..4ef4090
--- /dev/null
@@ -0,0 +1,383 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2017 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.
+//
+// ---------------------------------------------------------------------
+
+// same as step-36_parallel, but solve for largest eigenvalues of GHEP
+
+#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/dynamic_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/lac/slepc_solver.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 <fstream>
+#include <iostream>
+
+// test parallel (MPI) version of Step-36
+
+const unsigned int dim = 2;//run in 2d to save time
+
+const double eps = 1e-10;
+
+void test (std::string solver_name,
+           std::string preconditioner_name)
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 5;
+
+  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<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
+    = DoFTools::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::DynamicSparsityPattern 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);
+
+  // initialize 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);
+
+  eigenfunctions.resize (5);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    {
+      eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+      for (unsigned int j=0; j<locally_owned_dofs.n_elements(); ++j)
+        eigenfunctions[i][locally_owned_dofs.nth_index_in_set(j)] = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+      eigenfunctions[i].compress(dealii::VectorOperation::insert);
+    }
+
+  eigenvalues.resize (eigenfunctions.size ());
+
+
+  // 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)
+          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)
+                   ) * 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 SLEPc by
+  {
+    PETScWrappers::PreconditionerBase *preconditioner;
+
+    dealii::deallog<<preconditioner_name<<std::endl;
+    if (preconditioner_name == "Jacobi")
+      {
+        preconditioner = new PETScWrappers::PreconditionJacobi(mpi_communicator);
+      }
+    else if (preconditioner_name == "Boomer")
+      {
+        PETScWrappers::PreconditionBoomerAMG::AdditionalData data;
+        data.symmetric_operator = true;
+
+        preconditioner = new PETScWrappers::PreconditionBoomerAMG(mpi_communicator,
+                                                                  data);
+      }
+    else if (preconditioner_name == "BlockJacobi")
+      {
+        preconditioner = new PETScWrappers::PreconditionBlockJacobi(mpi_communicator);
+      }
+    else
+      {
+        AssertThrow (false, ExcMessage ("not supported preconditioner"));
+
+        // make compiler happy
+        preconditioner = new PETScWrappers::PreconditionJacobi(mpi_communicator);
+      }
+
+    dealii::SolverControl linear_solver_control (dof_handler.n_dofs(), 1e-15,/*log_history*/false,/*log_results*/false);
+    PETScWrappers::SolverCG  linear_solver(linear_solver_control,mpi_communicator);
+    linear_solver.initialize(*preconditioner);
+
+    dealii::SolverControl solver_control (100, 1e-11,/*log_history*/false,/*log_results*/false);
+
+    dealii::SLEPcWrappers::SolverBase *eigensolver;
+
+    dealii::deallog<<solver_name<<std::endl;
+    // Get a handle on the wanted eigenspectrum solver
+    if (solver_name == "KrylovSchur")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control,
+                                                                    mpi_communicator);
+      }
+
+    else if (solver_name == "GeneralizedDavidson")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverGeneralizedDavidson (solver_control,
+            mpi_communicator);
+      }
+    else if (solver_name == "JacobiDavidson")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverJacobiDavidson (solver_control,
+            mpi_communicator);
+      }
+    else if (solver_name == "Lanczos")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverLanczos (solver_control,
+                                                                mpi_communicator);
+      }
+    else
+      {
+        AssertThrow (false, ExcMessage ("not supported eigensolver"));
+
+        // Make compiler happy and not complaining about non
+        // uninitialized variables
+        eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control,
+                                                                    mpi_communicator);
+      }
+
+    // Set the initial vector. This is optional, if not done the initial vector is set to random values
+    eigensolver->set_initial_space(eigenfunctions);
+
+    eigensolver->set_which_eigenpairs (EPS_LARGEST_REAL);
+    eigensolver->set_problem_type (EPS_GHEP);
+
+    eigensolver->solve (stiffness_matrix,
+                        mass_matrix,
+                        eigenvalues,
+                        eigenfunctions,
+                        eigenfunctions.size());
+
+    // TODO make this robust on different platforms. Seems related to GHEP
+    // as solve_04 works ok.
+    //dealii::deallog << "outer iterations: "<< solver_control.last_step ()<<std::endl;
+    //dealii::deallog << "last inner iterations: "<<linear_solver_control.last_step()<<std::endl;
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      dealii::deallog << eigenvalues[i] << std::endl;
+
+    delete preconditioner;
+    delete eigensolver;
+
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*B*x_i).L2() == 0
+    // b) x_j*B*x_i=\delta_{ij}
+    {
+      const double precision = 1e-5;
+      PETScWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          mass_matrix.vmult(Bx,eigenfunctions[i]);
+
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< precision,
+                    ExcMessage("Eigenvectors " +
+                               Utilities::int_to_string(i) +
+                               " and " +
+                               Utilities::int_to_string(j) +
+                               " are not orthonormal!"));
+
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*eigenvalues[i],Bx);
+          Assert (Ax.l2_norm() < precision,
+                  ExcMessage(Utilities::to_string(Ax.l2_norm())));
+        }
+    }
+  }
+
+
+  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.depth_console(0);
+  dealii::deallog.threshold_double(eps);
+
+  try
+    {
+      dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        test ("KrylovSchur","Jacobi");
+        test ("KrylovSchur","BlockJacobi");
+        test ("KrylovSchur","Boomer");
+//        test ("GeneralizedDavidson");
+//        test ("JacobiDavidson");
+      }
+    }
+  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/slepc/step-36_parallel_02.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output b/tests/slepc/step-36_parallel_02.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..2a73664
--- /dev/null
@@ -0,0 +1,24 @@
+DEAL::Jacobi
+DEAL::KrylovSchur
+DEAL::6099.84
+DEAL::6035.05
+DEAL::6035.05
+DEAL::5970.26
+DEAL::5931.69
+DEAL::Ok
+DEAL::BlockJacobi
+DEAL::KrylovSchur
+DEAL::6099.84
+DEAL::6035.05
+DEAL::6035.05
+DEAL::5970.26
+DEAL::5931.69
+DEAL::Ok
+DEAL::Boomer
+DEAL::KrylovSchur
+DEAL::6099.84
+DEAL::6035.05
+DEAL::6035.05
+DEAL::5970.26
+DEAL::5931.69
+DEAL::Ok
diff --git a/tests/slepc/step-36_parallel_03.cc b/tests/slepc/step-36_parallel_03.cc
new file mode 100644 (file)
index 0000000..6adf4a9
--- /dev/null
@@ -0,0 +1,380 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2017 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.
+//
+// ---------------------------------------------------------------------
+
+// same as step-36_parallel_02, but solve SHEP
+
+#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/dynamic_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/lac/slepc_solver.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 <fstream>
+#include <iostream>
+
+// test parallel (MPI) version of Step-36
+
+const unsigned int dim = 2;//run in 2d to save time
+
+const double eps = 1e-10;
+
+void test (std::string solver_name,
+           std::string preconditioner_name)
+{
+  const unsigned int global_mesh_refinement_steps = 5;
+  const unsigned int number_of_eigenvalues        = 5;
+
+  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<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
+    = DoFTools::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::DynamicSparsityPattern 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);
+
+  // initialize 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);
+
+  eigenfunctions.resize (5);
+  for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+    {
+      eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+      for (unsigned int j=0; j<locally_owned_dofs.n_elements(); ++j)
+        eigenfunctions[i][locally_owned_dofs.nth_index_in_set(j)] = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+      eigenfunctions[i].compress(dealii::VectorOperation::insert);
+    }
+
+  eigenvalues.resize (eigenfunctions.size ());
+
+
+  // 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)
+          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)
+                   ) * 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 SLEPc by
+  {
+    PETScWrappers::PreconditionerBase *preconditioner;
+
+    dealii::deallog<<preconditioner_name<<std::endl;
+    if (preconditioner_name == "Jacobi")
+      {
+        preconditioner = new PETScWrappers::PreconditionJacobi(mpi_communicator);
+      }
+    else if (preconditioner_name == "Boomer")
+      {
+        PETScWrappers::PreconditionBoomerAMG::AdditionalData data;
+        data.symmetric_operator = true;
+
+        preconditioner = new PETScWrappers::PreconditionBoomerAMG(mpi_communicator,
+                                                                  data);
+      }
+    else if (preconditioner_name == "BlockJacobi")
+      {
+        preconditioner = new PETScWrappers::PreconditionBlockJacobi(mpi_communicator);
+      }
+    else
+      {
+        AssertThrow (false, ExcMessage ("not supported preconditioner"));
+
+        // make compiler happy
+        preconditioner = new PETScWrappers::PreconditionJacobi(mpi_communicator);
+      }
+
+    dealii::SolverControl linear_solver_control (dof_handler.n_dofs(), 1e-15,/*log_history*/false,/*log_results*/false);
+    PETScWrappers::SolverCG  linear_solver(linear_solver_control,mpi_communicator);
+    linear_solver.initialize(*preconditioner);
+
+    dealii::SolverControl solver_control (100, 1e-11,/*log_history*/false,/*log_results*/false);
+
+    dealii::SLEPcWrappers::SolverBase *eigensolver;
+
+    dealii::deallog<<solver_name<<std::endl;
+    // Get a handle on the wanted eigenspectrum solver
+    if (solver_name == "KrylovSchur")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control,
+                                                                    mpi_communicator);
+      }
+
+    else if (solver_name == "GeneralizedDavidson")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverGeneralizedDavidson (solver_control,
+            mpi_communicator);
+      }
+    else if (solver_name == "JacobiDavidson")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverJacobiDavidson (solver_control,
+            mpi_communicator);
+      }
+    else if (solver_name == "Lanczos")
+      {
+        eigensolver = new dealii::SLEPcWrappers::SolverLanczos (solver_control,
+                                                                mpi_communicator);
+      }
+    else
+      {
+        AssertThrow (false, ExcMessage ("not supported eigensolver"));
+
+        // Make compiler happy and not complaining about non
+        // uninitialized variables
+        eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control,
+                                                                    mpi_communicator);
+      }
+
+    // Set the initial vector. This is optional, if not done the initial vector is set to random values
+    eigensolver->set_initial_space(eigenfunctions);
+
+    eigensolver->set_which_eigenpairs (EPS_LARGEST_REAL);
+    eigensolver->set_problem_type (EPS_HEP);
+
+    eigensolver->solve (stiffness_matrix,
+                        eigenvalues,
+                        eigenfunctions,
+                        eigenfunctions.size());
+
+    // TODO make this robust on different platforms. Seems related to GHEP
+    // as solve_04 works ok.
+    //dealii::deallog << "outer iterations: "<< solver_control.last_step ()<<std::endl;
+    //dealii::deallog << "last inner iterations: "<<linear_solver_control.last_step()<<std::endl;
+    for (unsigned int i=0; i < eigenvalues.size(); i++)
+      dealii::deallog << eigenvalues[i] << std::endl;
+
+    delete preconditioner;
+    delete eigensolver;
+
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*x_i).L2() == 0
+    // b) x_j*x_i=\delta_{ij}
+    {
+      const double precision = 1e-5;
+      PETScWrappers::MPI::Vector Ax(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            Assert( std::abs( eigenfunctions[j] * eigenfunctions[i] - (i==j))< precision,
+                    ExcMessage("Eigenvectors " +
+                               Utilities::int_to_string(i) +
+                               " and " +
+                               Utilities::int_to_string(j) +
+                               " are not orthonormal!"));
+
+          stiffness_matrix.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*eigenvalues[i],eigenfunctions[i]);
+          Assert (Ax.l2_norm() < precision,
+                  ExcMessage(Utilities::to_string(Ax.l2_norm())));
+        }
+    }
+  }
+
+
+  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.depth_console(0);
+  dealii::deallog.threshold_double(eps);
+
+  try
+    {
+      dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        test ("KrylovSchur","Jacobi");
+        test ("KrylovSchur","BlockJacobi");
+        test ("KrylovSchur","Boomer");
+//        test ("GeneralizedDavidson");
+//        test ("JacobiDavidson");
+      }
+    }
+  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/slepc/step-36_parallel_03.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output b/tests/slepc/step-36_parallel_03.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..3fdcfca
--- /dev/null
@@ -0,0 +1,24 @@
+DEAL::Jacobi
+DEAL::KrylovSchur
+DEAL::3.98719
+DEAL::3.98719
+DEAL::3.97768
+DEAL::3.97768
+DEAL::3.96194
+DEAL::Ok
+DEAL::BlockJacobi
+DEAL::KrylovSchur
+DEAL::3.98719
+DEAL::3.98719
+DEAL::3.97768
+DEAL::3.97768
+DEAL::3.96194
+DEAL::Ok
+DEAL::Boomer
+DEAL::KrylovSchur
+DEAL::3.98719
+DEAL::3.98719
+DEAL::3.97768
+DEAL::3.97768
+DEAL::3.96194
+DEAL::Ok

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.