From 2dc1c04092e5d8a85fcffef9b9f82a4f3fdeb442 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 14 Jun 2016 11:00:47 +0200 Subject: [PATCH] parpack: add set_initial_vector() for Krylov subspace --- include/deal.II/lac/parpack_solver.h | 29 ++++++++++++++++++++++-- tests/arpack/step-36_parpack.cc | 2 ++ tests/arpack/step-36_parpack_trilinos.cc | 3 ++- 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/include/deal.II/lac/parpack_solver.h b/include/deal.II/lac/parpack_solver.h index 972de1286a..5cfbf42591 100644 --- a/include/deal.II/lac/parpack_solver.h +++ b/include/deal.II/lac/parpack_solver.h @@ -259,7 +259,13 @@ public: const std::vector &partitioning); /** - * Set desired shift value. + * Set initial vector for building Krylov space. + */ + void set_initial_vector(const VectorType &vec); + + /** + * Set desired shift value. If this function is not called, the + * shift is assumed to be zero. */ void set_shift(const double s ); @@ -344,6 +350,11 @@ protected: */ std::vector v; + /** + * An auxiliary flag which is set to true when initial vector is provided. + */ + bool initial_vector_provided; + /** * The initial residual vector, possibly from a previous run. On output, it * contains the final residual vector. @@ -491,6 +502,7 @@ PArpackSolver::PArpackSolver (SolverControl &control, additional_data (data), mpi_communicator( mpi_communicator ), mpi_communicator_fortran ( MPI_Comm_c2f( mpi_communicator ) ), + initial_vector_provided(false), shift_value(0.0) {} @@ -501,6 +513,19 @@ void PArpackSolver::set_shift(const double s ) shift_value = s; } +template +void PArpackSolver:: +set_initial_vector(const VectorType &vec) +{ + initial_vector_provided = true; + Assert (resid.size() == local_indices.size(), + ExcDimensionMismatch(resid.size(),local_indices.size())); + vec.extract_subvector_to (local_indices.begin(), + local_indices.end(), + &resid[0]); +} + + template void PArpackSolver:: internal_reinit(const IndexSet &locally_owned_dofs) @@ -684,7 +709,7 @@ void PArpackSolver::solve // possibly from a previous run. // Typical choices in this situation might be to use the final value // of the starting vector from the previous eigenvalue calculation - int info = 1; + int info = initial_vector_provided? 1 : 0; // Number of eigenvalues of OP to be computed. 0 < NEV < N. int nev = n_eigenvalues; diff --git a/tests/arpack/step-36_parpack.cc b/tests/arpack/step-36_parpack.cc index 7e4df13cb8..b9cef33e15 100644 --- a/tests/arpack/step-36_parpack.cc +++ b/tests/arpack/step-36_parpack.cc @@ -286,6 +286,8 @@ void test () mpi_communicator, additional_data); eigensolver.reinit(locally_owned_dofs); + eigenfunctions[0] = 1.; + eigensolver.set_initial_vector(eigenfunctions[0]); eigensolver.solve (stiffness_matrix, mass_matrix, inverse, diff --git a/tests/arpack/step-36_parpack_trilinos.cc b/tests/arpack/step-36_parpack_trilinos.cc index 90a2a65ccf..ee2ff903b8 100644 --- a/tests/arpack/step-36_parpack_trilinos.cc +++ b/tests/arpack/step-36_parpack_trilinos.cc @@ -276,7 +276,8 @@ void test () additional_data); eigensolver.reinit(locally_owned_dofs); eigensolver.set_shift(shift); - + eigenfunctions[0] = 1.; + eigensolver.set_initial_vector(eigenfunctions[0]); // avoid output of iterative solver: const unsigned int previous_depth = deallog.depth_file(0); eigensolver.solve (stiffness_matrix, -- 2.39.5