From bb703601544b65918abb31afb54e0f9452a86ee6 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 9 Nov 2015 15:41:02 +0100 Subject: [PATCH] GHEP and SHEP unit tests for SLEPc --- tests/slepc/solve_01.cc | 146 ++++++++++++++++++++++++++++ tests/slepc/solve_01.output | 30 ++++++ tests/slepc/solve_04.cc | 142 +++++++++++++++++++++++++++ tests/slepc/solve_04.output | 33 +++++++ tests/slepc/testmatrix.h | 185 ++++++++++++++++++++++++++++++++++++ 5 files changed, 536 insertions(+) create mode 100644 tests/slepc/solve_01.cc create mode 100644 tests/slepc/solve_01.output create mode 100644 tests/slepc/solve_04.cc create mode 100644 tests/slepc/solve_04.output create mode 100644 tests/slepc/testmatrix.h 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void +check_solve( SOLVER &solver, + const SolverControl &solver_control, + const MATRIX &A, const MATRIX &B, + std::vector &u, std::vector &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 << " "< u(n_eigenvalues, + PETScWrappers::Vector(dim)); + std::vector 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void +check_solve( SOLVER &solver, + const SolverControl &solver_control, + const MATRIX &A, + std::vector &u, std::vector &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 << " "< u(n_eigenvalues, + PETScWrappers::Vector(dim)); + std::vector 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 +#include +#include +#include +#include +#include + + +/** + * 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 + void diag_structure(SP &structure) const; + + /** + * Fill the matrix with values. + */ + template + void diag(MATRIX &) const; + + template + void gnuplot_print(std::ostream &, const Vector &) 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 +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 +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 +inline +void +FDDiagMatrix::gnuplot_print(std::ostream &s, const Vector &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 + void three_point_structure(SP &structure) const; + + /** + * Fill the matrix with values. + */ + template + 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 +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 +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); + } +} -- 2.39.5