From 2dc1c04092e5d8a85fcffef9b9f82a4f3fdeb442 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
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<IndexSet> &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<double> 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<VectorType>::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<VectorType>::set_shift(const double s )
   shift_value = s;
 }
 
+template <typename VectorType>
+void PArpackSolver<VectorType>::
+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 <typename VectorType>
 void PArpackSolver<VectorType>::
 internal_reinit(const IndexSet &locally_owned_dofs)
@@ -684,7 +709,7 @@ void PArpackSolver<VectorType>::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