From 39a064f2f1c5ee3893057fdec918d8c02739350d Mon Sep 17 00:00:00 2001
From: Giuseppe Pitton <gpitton@sissa.it>
Date: Thu, 7 Apr 2016 16:01:56 -0500
Subject: [PATCH] Update a comment for a variable.

---
 doc/news/changes.h                 |  7 +++
 include/deal.II/lac/solver_gmres.h | 99 +++++++++++++++++++++++++++++-
 tests/lac/solver_signals.cc        | 24 ++++++++
 tests/lac/solver_signals.output    | 61 ++++++++++++++++++
 4 files changed, 189 insertions(+), 2 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index eb6858df67..6ae21c3be5 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -145,6 +145,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> New: added hessenberg_signal and krylov_space_signal to SolverGMRES.
+ These signals allow to retrieve the Hessenberg matrix and the basis vectors
+ generated by the Arnoldi algorithm.
+ <br>
+ (Giuseppe Pitton, Luca Heltai, 2016/04/11)
+ </li>
+
  <li> Fixed: Meshworker::Assembler::ResidualSimple now also works for
  multiple blocks if no constraints are given.
  <br>
diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h
index c3e031f821..f83ccd4538 100644
--- a/include/deal.II/lac/solver_gmres.h
+++ b/include/deal.II/lac/solver_gmres.h
@@ -83,6 +83,13 @@ namespace internal
       VectorType &operator() (const unsigned int i,
                               const VectorType   &temp);
 
+      /**
+       * Returns size of data vector. It is used in the solver to store
+       * the Arnoldi vectors.
+       */
+      unsigned int size() const;
+
+
     private:
       /**
        * Pool were vectors are obtained from.
@@ -288,6 +295,28 @@ public:
     const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
     const bool every_iteration=false);
 
+  /**
+   * Connect a slot to retrieve the Hessenberg matrix obtained by the
+   * projection of the initial matrix on the Krylov basis. Called on each
+   * outer iteration if every_iteration=true, otherwise called once when
+   * iterations are ended (i.e., either because convergence has been achieved,
+   * or because divergence has been detected).
+   */
+  boost::signals2::connection
+  connect_hessenberg_slot(
+    const std_cxx11::function<void (const FullMatrix<double> &)> &slot,
+    const bool every_iteration=true);
+
+  /**
+   * Connect a slot to retrieve the basis vectors of the Krylov space
+   * generated by the Arnoldi algorithm. Called at once when iterations
+   * are completed (i.e., either because convergence has been achieved,
+   * or because divergence has been detected).
+   */
+  boost::signals2::connection
+  connect_krylov_space_slot(
+    const std_cxx11::function<void (const internal::SolverGMRES::TmpVectors<VectorType> &)> &slot);
+
 
   DeclException1 (ExcTooFewTmpVectors,
                   int,
@@ -326,6 +355,23 @@ protected:
    */
   boost::signals2::signal<void (const std::vector<std::complex<double> > &)> all_eigenvalues_signal;
 
+  /**
+   * Signal used to retrieve the Hessenberg matrix. Called once when
+   * all iterations are ended.
+   */
+  boost::signals2::signal<void (const FullMatrix<double> &)> hessenberg_signal;
+
+  /**
+   * Signal used to retrieve the Hessenberg matrix. Called on each outer
+   * iteration.
+   */
+  boost::signals2::signal<void (const FullMatrix<double> &)> all_hessenberg_signal;
+
+  /**
+   * Signal used to retrieve the Krylov space basis vectors. Called once
+   * when all iterations are ended.
+   */
+  boost::signals2::signal<void (const internal::SolverGMRES::TmpVectors<VectorType> &)> krylov_space_signal;
 
   /**
    * Implementation of the computation of the norm of the residual.
@@ -369,6 +415,7 @@ protected:
     const FullMatrix<double> &H_orig ,
     const unsigned int dim,
     const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
+    const boost::signals2::signal<void (const FullMatrix<double> &)> &hessenberg_signal,
     const boost::signals2::signal<void(double)> &cond_signal,
     const bool log_eigenvalues);
 
@@ -537,6 +584,15 @@ namespace internal
       return *data[i-offset];
     }
 
+
+    template <class VectorType>
+    unsigned int
+    TmpVectors<VectorType>::size() const
+    {
+      return (data.size() > 0 ? data.size()-1 : 0);
+    }
+
+
     // A comparator for better printing eigenvalues
     inline
     bool complex_less_pred(const std::complex<double> &x,
@@ -701,16 +757,19 @@ SolverGMRES<VectorType>::compute_eigs_and_cond
 (const FullMatrix<double>                     &H_orig,
  const unsigned int                           dim,
  const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
+ const boost::signals2::signal<void (const FullMatrix<double> &)> &hessenberg_signal,
  const boost::signals2::signal<void (double)> &cond_signal,
  const bool                                   log_eigenvalues)
 {
   //Avoid copying the Hessenberg matrix if it isn't needed.
-  if (!eigenvalues_signal.empty() || !cond_signal.empty() || log_eigenvalues )
+  if (!eigenvalues_signal.empty() || !hessenberg_signal.empty()
+      || !cond_signal.empty() || log_eigenvalues )
     {
       LAPACKFullMatrix<double> mat(dim,dim);
       for (unsigned int i=0; i<dim; ++i)
         for (unsigned int j=0; j<dim; ++j)
           mat(i,j) = H_orig(i,j);
+      hessenberg_signal(H_orig);
       //Avoid computing eigenvalues if they are not needed.
       if (!eigenvalues_signal.empty() || log_eigenvalues )
         {
@@ -777,6 +836,8 @@ SolverGMRES<VectorType>::solve (const MatrixType         &A,
     |!all_condition_numbers_signal.empty()
     |!eigenvalues_signal.empty()
     |!all_eigenvalues_signal.empty()
+    |!hessenberg_signal.empty()
+    |!all_hessenberg_signal.empty()
     |additional_data.compute_eigenvalues;
   // for eigenvalue computation, need to collect the Hessenberg matrix (before
   // applying Givens rotations)
@@ -1003,6 +1064,7 @@ SolverGMRES<VectorType>::solve (const MatrixType         &A,
           H1(i,j) = H(i,j);
 
       compute_eigs_and_cond(H_orig,dim,all_eigenvalues_signal,
+                            all_hessenberg_signal,
                             all_condition_numbers_signal,
                             additional_data.compute_eigenvalues);
 
@@ -1024,8 +1086,13 @@ SolverGMRES<VectorType>::solve (const MatrixType         &A,
     }
   while (iteration_state == SolverControl::iterate);
 
-  compute_eigs_and_cond(H_orig,dim,eigenvalues_signal,condition_number_signal,
+  compute_eigs_and_cond(H_orig,dim,eigenvalues_signal,hessenberg_signal,
+                        condition_number_signal,
                         false);
+
+  if (!krylov_space_signal.empty())
+    krylov_space_signal(tmp_vectors);
+
   if (!use_default_residual)
     {
       this->memory.free(r);
@@ -1080,6 +1147,34 @@ SolverGMRES<VectorType>::connect_eigenvalues_slot
 
 
 
+template<class VectorType>
+boost::signals2::connection
+SolverGMRES<VectorType>::connect_hessenberg_slot
+(const std_cxx11::function<void (const FullMatrix<double> &)> &slot,
+ const bool every_iteration)
+{
+  if (every_iteration)
+    {
+      return all_hessenberg_signal.connect(slot);
+    }
+  else
+    {
+      return hessenberg_signal.connect(slot);
+    }
+}
+
+
+
+template<class VectorType>
+boost::signals2::connection
+SolverGMRES<VectorType>::connect_krylov_space_slot
+(const std_cxx11::function<void (const internal::SolverGMRES::TmpVectors<VectorType> &)> &slot)
+{
+  return krylov_space_signal.connect(slot);
+}
+
+
+
 template<class VectorType>
 double
 SolverGMRES<VectorType>::criterion ()
diff --git a/tests/lac/solver_signals.cc b/tests/lac/solver_signals.cc
index 18eedb9930..e2857462b9 100644
--- a/tests/lac/solver_signals.cc
+++ b/tests/lac/solver_signals.cc
@@ -56,6 +56,26 @@ void output_eigenvalues(const std::vector<NUMBER> &eigenvalues,const std::string
   deallog << std::endl;
 }
 
+template<class NUMBER>
+void output_hessenberg_matrix(const FullMatrix<NUMBER> &H,const std::string &text)
+{
+  deallog << text << std::endl;
+  for (unsigned int i=0; i<H.m(); ++i)
+    {
+      for (unsigned int j=0; j<H.n(); ++j)
+        deallog << H(i,j) << " ";
+      deallog << std::endl;
+    }
+}
+
+template<class NUMBER>
+void output_arnoldi_vectors_norms(const internal::SolverGMRES::TmpVectors<Vector<NUMBER> > &tmp_vector,const std::string &text)
+{
+  deallog << text << std::endl;
+  for (unsigned int i=0; i<tmp_vector.size(); ++i)
+    deallog << tmp_vector[i].l2_norm() << std::endl;
+}
+
 template<typename SolverType, typename MatrixType, typename VectorType, class PRECONDITION>
 void
 check_solve(SolverType         &solver,
@@ -128,6 +148,10 @@ int main()
         std_cxx11::bind(output_eigenvalues<std::complex<double> >,std_cxx11::_1,"Eigenvalues: "),true);
       solver_gmres.connect_eigenvalues_slot(
         std_cxx11::bind(output_eigenvalues<std::complex<double> >,std_cxx11::_1,"Final Eigenvalues: "));
+      solver_gmres.connect_hessenberg_slot(
+        std_cxx11::bind(output_hessenberg_matrix<double>,std_cxx11::_1,"Hessenberg matrix: "),false);
+      solver_gmres.connect_krylov_space_slot(
+        std_cxx11::bind(output_arnoldi_vectors_norms<double>,std_cxx11::_1,"Arnoldi vectors norms: "));
       solver_gmres.solve(A,u,f,PreconditionIdentity());
     }
   catch (std::exception &e)
diff --git a/tests/lac/solver_signals.output b/tests/lac/solver_signals.output
index caa51fc486..896dae4732 100644
--- a/tests/lac/solver_signals.output
+++ b/tests/lac/solver_signals.output
@@ -127,5 +127,66 @@ DEAL:GMRES::Condition number: 288.5
 DEAL:GMRES::Convergence step 65 value 0.0009799
 DEAL:GMRES::Eigenvalues:  (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000)
 DEAL:GMRES::Condition number: 288.2
+DEAL:GMRES::Hessenberg matrix: 
+DEAL:GMRES::0.1982 0.8628 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0.8628 6.098 1.971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 1.971 4.008 1.721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 1.721 3.583 2.004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 2.004 4.008 1.910 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 1.910 3.325 1.950 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 1.950 4.036 1.744 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 1.744 3.951 1.992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 1.992 4.080 1.994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 1.920 4.007 2.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 2.002 3.981 2.022 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 2.022 4.028 1.908 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 1.908 4.238 1.925 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 1.925 4.011 2.042 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 2.042 3.827 2.111 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.111 3.742 2.046 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 4.070 1.879 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.879 4.168 1.933 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.933 4.013 1.886 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.886 3.963 1.957 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.957 3.818 2.033 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.033 4.006 2.201 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.201 3.944 1.862 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.862 4.039 1.827 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.827 3.704 1.844 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.844 4.073 2.046 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 3.592 1.647 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.647 3.431 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.155 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:GMRES::Final Eigenvalues:  (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000)
 DEAL:GMRES::Final condition number: 288.2
+DEAL:GMRES::Arnoldi vectors norms: 
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
+DEAL:GMRES::1.000
-- 
2.39.5