From: Denis Davydov <davydden@gmail.com>
Date: Mon, 9 Nov 2015 14:41:02 +0000 (+0100)
Subject: GHEP and SHEP unit tests for SLEPc
X-Git-Tag: v8.4.0-rc2~233^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bb703601544b65918abb31afb54e0f9452a86ee6;p=dealii.git

GHEP and SHEP unit tests for SLEPc
---

diff --git a/tests/slepc/solve_01.cc b/tests/slepc/solve_01.cc
new file mode 100644
index 0000000000..41e874f3da
--- /dev/null
+++ b/tests/slepc/solve_01.cc
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test the SLEPc solvers with GHEP
+// the unit tests is a mirror of
+// SLEPc-3.6.1/src/eps/examples/tests/test1.c
+
+
+#include "../tests.h"
+#include "../lac/testmatrix.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/slepc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/vector_memory.h>
+#include <typeinfo>
+
+template<class SOLVER, class MATRIX, class VECTOR>
+void
+check_solve( SOLVER &solver,
+             const SolverControl &solver_control,
+             const MATRIX &A, const MATRIX &B,
+             std::vector<VECTOR> &u, std::vector<PetscScalar > &v)
+{
+  deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+  try
+    {
+      solver.set_problem_type(EPS_GHEP);
+      // reset vectors and set them as initial space
+      // to avoid dependency on random numbers:
+      for (unsigned int i = 0; i < u.size(); i++)
+        {
+          u[i] = 0.0;
+          u[i][i] = 1.0;
+        }
+      solver.set_initial_space(u);
+
+      solver.solve(A,B,v,u,v.size());
+    }
+  catch (dealii::SolverControl::NoConvergence &e)
+    {
+      deallog << "Exception: " << e.get_exc_name() << std::endl;
+    }
+
+  deallog << "Solver stopped after " << solver_control.last_step()
+          << " iterations" << std::endl;
+
+  deallog << "Eigenvalues:";
+  for (unsigned int i = 0; i < v.size(); i++)
+    {
+      deallog << " "<<v[i];
+      if ( i != (v.size()-1) )
+        deallog <<",";
+    }
+  deallog << std::endl << std::endl;
+}
+
+
+int main(int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog << std::setprecision(7);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  {
+    SolverControl control(1000, 1000*PETSC_MACHINE_EPSILON);
+
+    const unsigned int size = 46;
+    unsigned int dim = (size-1)*(size-1);
+
+    const unsigned n_eigenvalues = 4;
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl << std::endl;
+
+    // Make matrix
+    FDMatrix testproblem(size, size);
+    FDDiagMatrix diagonal(size,size);
+    PETScWrappers::SparseMatrix  A(dim, dim, 5);
+    PETScWrappers::SparseMatrix  B(dim, dim, 5);
+    testproblem.five_point(A);
+    A.compress (VectorOperation::insert);
+    diagonal.diag(B);
+    B.compress (VectorOperation::insert);
+
+    std::vector<PETScWrappers::Vector>  u(n_eigenvalues,
+                                          PETScWrappers::Vector(dim));
+    std::vector<PetscScalar> v(n_eigenvalues);
+
+    {
+      SLEPcWrappers::SolverKrylovSchur solver(control);
+      check_solve (solver, control, A,B,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverArnoldi solver(control);
+      check_solve (solver, control, A,B,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverLanczos solver(control);
+      check_solve (solver, control, A,B,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverGeneralizedDavidson solver(control);
+      check_solve (solver, control, A,B,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverJacobiDavidson solver(control);
+      check_solve (solver, control, A,B,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverGeneralizedDavidson::AdditionalData data(true);
+      SLEPcWrappers::SolverGeneralizedDavidson solver(control,PETSC_COMM_SELF,data);
+      check_solve (solver, control, A,B,u,v);
+    }
+  }
+
+}
+
diff --git a/tests/slepc/solve_01.output b/tests/slepc/solve_01.output
new file mode 100644
index 0000000000..635bea8a23
--- /dev/null
+++ b/tests/slepc/solve_01.output
@@ -0,0 +1,30 @@
+
+DEAL::Size 46 Unknowns 2025
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers17SolverKrylovSchurE
+DEAL::Convergence step 23 value 0
+DEAL::Solver stopped after 23 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers13SolverArnoldiE
+DEAL::Convergence step 80 value 0
+DEAL::Solver stopped after 80 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers13SolverLanczosE
+DEAL::Convergence step 79 value 0
+DEAL::Solver stopped after 79 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers25SolverGeneralizedDavidsonE
+DEAL::Solver stopped after 278 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers20SolverJacobiDavidsonE
+DEAL::Solver stopped after 40 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers25SolverGeneralizedDavidsonE
+DEAL::Solver stopped after 581 iterations
+DEAL::Eigenvalues: 29.73524, 29.67536, 29.58873, 29.46785
+DEAL::
diff --git a/tests/slepc/solve_04.cc b/tests/slepc/solve_04.cc
new file mode 100644
index 0000000000..5514930dbf
--- /dev/null
+++ b/tests/slepc/solve_04.cc
@@ -0,0 +1,142 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test the SLEPc solvers with HEP
+// the unit tests is a mirror of
+// SLEPc-3.6.1/src/eps/examples/tests/test4.c
+
+
+#include "../tests.h"
+#include "../lac/testmatrix.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/slepc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/vector_memory.h>
+#include <typeinfo>
+
+template<class SOLVER, class MATRIX, class VECTOR>
+void
+check_solve( SOLVER &solver,
+             const SolverControl &solver_control,
+             const MATRIX &A,
+             std::vector<VECTOR> &u, std::vector<PetscScalar > &v)
+{
+  deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+  try
+    {
+      solver.set_problem_type(EPS_HEP);
+      // reset vectors and set them as initial space
+      // to avoid dependency on random numbers:
+      for (unsigned int i = 0; i < u.size(); i++)
+        {
+          u[i] = 0.0;
+          u[i][i] = 1.0;
+        }
+      solver.set_initial_space(u);
+      // solve
+      solver.solve(A,v,u,v.size());
+    }
+  catch (dealii::SolverControl::NoConvergence &e)
+    {
+      deallog << "Exception: " << e.get_exc_name() << std::endl;
+    }
+
+  deallog << "Solver stopped after " << solver_control.last_step()
+          << " iterations" << std::endl;
+
+  deallog << "Eigenvalues:";
+  for (unsigned int i = 0; i < v.size(); i++)
+    {
+      deallog << " "<<v[i];
+      if ( i != (v.size()-1) )
+        deallog <<",";
+    }
+  deallog << std::endl << std::endl;
+}
+
+
+int main(int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog << std::setprecision(6);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  {
+    SolverControl control(500, 1000*PETSC_MACHINE_EPSILON);
+
+    const unsigned int size = 31;
+    unsigned int dim = (size-1);
+
+    const unsigned n_eigenvalues = 4;
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl << std::endl;
+
+    // Make matrix
+    FD1DLaplaceMatrix testproblem(size);
+    PETScWrappers::SparseMatrix  A(dim, dim, 3);
+    testproblem.three_point(A);
+    A.compress (VectorOperation::insert);
+
+    std::vector<PETScWrappers::Vector>  u(n_eigenvalues,
+                                          PETScWrappers::Vector(dim));
+    std::vector<PetscScalar> v(n_eigenvalues);
+
+    {
+      SLEPcWrappers::SolverKrylovSchur solver(control);
+      check_solve (solver, control, A,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverArnoldi solver(control);
+      check_solve (solver, control, A,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverLanczos solver(control);
+      check_solve (solver, control, A,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverGeneralizedDavidson solver(control);
+      check_solve (solver, control, A,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverJacobiDavidson solver(control);
+      check_solve (solver, control, A,u,v);
+    }
+
+    {
+      SLEPcWrappers::SolverGeneralizedDavidson::AdditionalData data(true);
+      SLEPcWrappers::SolverGeneralizedDavidson solver(control,PETSC_COMM_SELF,data);
+      check_solve (solver, control, A,u,v);
+    }
+  }
+
+}
+
diff --git a/tests/slepc/solve_04.output b/tests/slepc/solve_04.output
new file mode 100644
index 0000000000..2f8efd4c1f
--- /dev/null
+++ b/tests/slepc/solve_04.output
@@ -0,0 +1,33 @@
+
+DEAL::Size 31 Unknowns 30
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers17SolverKrylovSchurE
+DEAL::Convergence step 5 value 0
+DEAL::Solver stopped after 5 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers13SolverArnoldiE
+DEAL::Convergence step 21 value 0
+DEAL::Solver stopped after 21 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers13SolverLanczosE
+DEAL::Convergence step 23 value 0
+DEAL::Solver stopped after 23 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers25SolverGeneralizedDavidsonE
+DEAL::Convergence step 120 value 0
+DEAL::Solver stopped after 120 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers20SolverJacobiDavidsonE
+DEAL::Convergence step 399 value 0
+DEAL::Solver stopped after 399 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
+DEAL::Solver type: N6dealii13SLEPcWrappers25SolverGeneralizedDavidsonE
+DEAL::Convergence step 145 value 0
+DEAL::Solver stopped after 145 iterations
+DEAL::Eigenvalues: 3.98974, 3.95906, 3.90828, 3.83792
+DEAL::
diff --git a/tests/slepc/testmatrix.h b/tests/slepc/testmatrix.h
new file mode 100644
index 0000000000..969b4c29ec
--- /dev/null
+++ b/tests/slepc/testmatrix.h
@@ -0,0 +1,185 @@
+
+#include "../tests.h"
+#include <deal.II/base/exceptions.h>
+#include <iostream>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/sparse_matrix_ez.h>
+#include <deal.II/lac/vector.h>
+
+
+/**
+ * Finite difference diagonal matrix on uniform grid.
+ */
+
+class FDDiagMatrix
+{
+public:
+  /**
+   * Constructor specifying grid resolution.
+   */
+  FDDiagMatrix(unsigned int nx, unsigned int ny);
+
+  /**
+   * Generate the matrix structure.
+   */
+  template <typename SP>
+  void diag_structure(SP &structure) const;
+
+  /**
+   * Fill the matrix with values.
+   */
+  template <typename MATRIX>
+  void diag(MATRIX &) const;
+
+  template <typename number>
+  void gnuplot_print(std::ostream &, const Vector<number> &) const;
+
+private:
+  /**
+   * Number of gridpoints in x-direction.
+   */
+  unsigned int nx;
+
+  /**
+   * Number of gridpoints in y-direction.
+   */
+  unsigned int ny;
+};
+
+
+// --------------- inline and template functions -----------------
+
+inline
+FDDiagMatrix::FDDiagMatrix(unsigned int nx, unsigned int ny)
+  :
+  nx(nx), ny(ny)
+{}
+
+
+
+template <typename SP>
+inline
+void
+FDDiagMatrix::diag_structure(SP &structure) const
+{
+  for (unsigned int i=0; i<=ny-2; i++)
+    {
+      for (unsigned int j=0; j<=nx-2; j++)
+        {
+          // Number of the row to be entered
+          unsigned int row = j+(nx-1)*i;
+          structure.add(row, row);
+        }
+    }
+}
+
+
+template<typename MATRIX>
+inline
+void
+FDDiagMatrix::diag(MATRIX &A) const
+{
+  for (unsigned int i=0; i<=ny-2; i++)
+    {
+      for (unsigned int j=0; j<=nx-2; j++)
+        {
+          // Number of the row to be entered
+          unsigned int row = j+(nx-1)*i;
+          A.set(row, row, 2./std::log(2.0+row));
+        }
+    }
+}
+
+template<typename number>
+inline
+void
+FDDiagMatrix::gnuplot_print(std::ostream &s, const Vector<number> &V) const
+{
+  for (unsigned int i=0; i<=ny-2; i++)
+    {
+      for (unsigned int j=0; j<=nx-2; j++)
+        {
+          // Number of the row to be entered
+          unsigned int row = j+(nx-1)*i;
+          s << (j+1)/(float)nx << '\t' << (i+1)/(float)ny << '\t' << V(row) << std::endl;
+        }
+      s << std::endl;
+    }
+  s << std::endl;
+}
+
+
+/**
+ * Finite difference matrix for laplace operator in 1D on uniform grid.
+ */
+
+class FD1DLaplaceMatrix
+{
+public:
+  /**
+   * Constructor specifying grid resolution.
+   */
+  FD1DLaplaceMatrix(unsigned int n);
+
+  /**
+   * Generate the matrix structure.
+   */
+  template <typename SP>
+  void three_point_structure(SP &structure) const;
+
+  /**
+   * Fill the matrix with values.
+   */
+  template <typename MATRIX>
+  void three_point(MATRIX &) const;
+
+private:
+  /**
+   * Number of gridpoints in x-direction.
+   */
+  unsigned int n;
+};
+
+
+// --------------- inline and template functions -----------------
+
+inline
+FD1DLaplaceMatrix::FD1DLaplaceMatrix(unsigned int n)
+  :
+  n(n)
+{}
+
+
+
+template <typename SP>
+inline
+void
+FD1DLaplaceMatrix::three_point_structure(SP &structure) const
+{
+  for (unsigned int i=0; i<=n-2; i++)
+    {
+      structure.add(i, i);
+      if (i >0)
+        structure.add(i,i-1);
+      if (i < n-1)
+        structure.add(i,i+1);
+    }
+}
+
+
+template<typename MATRIX>
+inline
+void
+FD1DLaplaceMatrix::three_point(MATRIX &A) const
+{
+  for (unsigned int i=0; i<=n-2; i++)
+    {
+      A.set(i,i,2.0);
+      if (i > 0)
+        A.set(i,i-1,-1.0);
+
+      if (i < n-2)
+        A.set(i,i+1,-1.0);
+    }
+}