]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add signals for retrieving additional data from solver_cg and solver_gmres. 703/head
authorSimon Sticko <simon.sticko@it.uu.se>
Sun, 29 Mar 2015 11:25:18 +0000 (13:25 +0200)
committerSimon Sticko <simon.sticko@it.uu.se>
Fri, 10 Jul 2015 19:44:30 +0000 (21:44 +0200)
SolverCG and SolverGMRES compute quantities such as eigenvalues, condition numbers
and coefficients during the solution process. Add boost::signals2 to be able to
communicate these estimates back to the caller.

include/deal.II/lac/solver_cg.h
include/deal.II/lac/solver_gmres.h
tests/lac/solver_signals.cc [new file with mode: 0644]
tests/lac/solver_signals.output [new file with mode: 0644]

index 72455efcc5635f33a22f4fcf736a3359939563a5..9e2f65abf68d728e23ea01f622a1556612803723 100644 (file)
@@ -66,9 +66,8 @@ class PreconditionIdentity;
  * are equal. But, even for small numbers of iteration steps, the condition
  * number of @p T is a good estimate for the one of @p PA.
  *
- * With the coefficients @p alpha and @p beta written to the log file if
- * <tt>AdditionalData::log_coefficients = true</tt>, the matrix @p T_m after
- * @p m steps is the tri-diagonal matrix with diagonal elements
+ * After @p m steps the matrix T_m can be written in terms of the coefficients
+ * @p alpha and @p beta as the tri-diagonal matrix with diagonal elements
  * <tt>1/alpha_0</tt>, <tt>1/alpha_1 + beta_0/alpha_0</tt>, ...,
  * <tt>1/alpha_{m-1</tt>+beta_{m-2}/alpha_{m-2}} and off-diagonal elements
  * <tt>sqrt(beta_0)/alpha_0</tt>, ..., <tt>sqrt(beta_{m-2</tt>)/alpha_{m-2}}.
@@ -77,6 +76,15 @@ class PreconditionIdentity;
  * @see Y. Saad: "Iterative methods for Sparse Linear Systems", section 6.7.3
  * for details.
  *
+ * The coefficients, eigenvalues and condition number (computed as the ratio of
+ * the largest over smallest eigenvalue) can be obtained by connecting a
+ * function as a slot to the solver using one of the functions
+ * @p connect_coefficients_slot, @p connect_eigenvalues_slot and
+ * @p connect_condition_number_slot. These slots will then be called from the
+ * solver with the estimates as argument.
+ *
+ * @deprecated Alternatively these estimates can be written to deallog by
+ * setting flags in @p AdditionalData.
  *
  * <h3>Observing the progress of linear solver iterations</h3>
  *
@@ -130,11 +138,18 @@ public:
 
     /**
      * Constructor. Initialize data fields.  Confer the description of those.
+     * @deprecated Instead use: connect_coefficients_slot,
+     * connect_condition_number_slot, and connect_eigenvalues_slot.
      */
-    AdditionalData (const bool log_coefficients = false,
+    AdditionalData (const bool log_coefficients,
                     const bool compute_condition_number = false,
                     const bool compute_all_condition_numbers = false,
-                    const bool compute_eigenvalues = false);
+                    const bool compute_eigenvalues = false) DEAL_II_DEPRECATED;
+
+    /**
+     * Constructor. Initializes all data fields to false.
+     */
+    AdditionalData();
   };
 
   /**
@@ -166,6 +181,38 @@ public:
          const VECTOR         &b,
          const PRECONDITIONER &precondition);
 
+  /**
+   * Connect a slot to retrieve the CG coefficients. The slot will be called
+   * with alpha as the first argument and with beta as the second argument,
+   * where alpha and beta follow the notation in
+   * Y. Saad: "Iterative methods for Sparse Linear Systems", section 6.7.
+   * Called once per iteration
+   */
+  boost::signals2::connection
+  connect_coefficients_slot(
+    const std_cxx11::function<void (double,double)> &slot);
+
+  /**
+   * Connect a slot to retrieve the estimated condition number.
+   * Called on each 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_condition_number_slot(const std_cxx11::function<void (double)> &slot,
+                                const bool every_iteration=false);
+
+  /**
+   * Connect a slot to retrieve the estimated eigenvalues.
+   * Called on each 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_eigenvalues_slot(
+    const std_cxx11::function<void (const std::vector<double> &)> &slot,
+    const bool every_iteration=false);
+
 protected:
   /**
    * Implementation of the computation of the norm of the residual. This can
@@ -183,6 +230,22 @@ protected:
                              const VECTOR &r,
                              const VECTOR &d) const;
 
+  /**
+   * Estimates the eigenvalues from diagonal and offdiagonal. Uses
+   * these estimate to compute the condition number. Calls the signals
+   * eigenvalues_signal and cond_signal with these estimates as arguments.
+   * Outputs the eigenvalues/condition-number to deallog if
+   * log_eigenvalues/log_cond is true.
+   */
+  static void
+  compute_eigs_and_cond(
+    const std::vector<double> &diagonal,
+    const std::vector<double> &offdiagonal,
+    const boost::signals2::signal<void (const std::vector<double> &)> &eigenvalues_signal,
+    const boost::signals2::signal<void (double)> &cond_signal,
+    const bool log_eigenvalues,
+    const bool log_cond);
+
   /**
    * Temporary vectors, allocated through the @p VectorMemory object at the
    * start of the actual solution process and deallocated at the end.
@@ -204,6 +267,36 @@ protected:
    */
   AdditionalData additional_data;
 
+  /**
+   * Signal used to retrieve the CG coefficients.
+   * Called on each iteration.
+   */
+  boost::signals2::signal<void (double,double)> coefficients_signal;
+
+  /**
+   * Signal used to retrieve the estimated condition number.
+   * Called once when all iterations are ended.
+   */
+  boost::signals2::signal<void (double)> condition_number_signal;
+
+  /**
+   * Signal used to retrieve the estimated condition numbers.
+   * Called on each iteration.
+   */
+  boost::signals2::signal<void (double)> all_condition_numbers_signal;
+
+  /**
+   * Signal used to retrieve the estimated eigenvalues.
+   * Called once when all iterations are ended.
+   */
+  boost::signals2::signal<void (const std::vector<double> &)> eigenvalues_signal;
+
+  /**
+   * Signal used to retrieve the estimated eigenvalues.
+   * Called on each iteration.
+   */
+  boost::signals2::signal<void (const std::vector<double> &)> all_eigenvalues_signal;
+
 private:
   void cleanup();
 };
@@ -230,6 +323,19 @@ AdditionalData (const bool log_coefficients,
 
 
 
+template <class VECTOR>
+inline
+SolverCG<VECTOR>::AdditionalData::
+AdditionalData ()
+  :
+  log_coefficients (false),
+  compute_condition_number(false),
+  compute_all_condition_numbers(false),
+  compute_eigenvalues(false)
+{}
+
+
+
 template <class VECTOR>
 SolverCG<VECTOR>::SolverCG (SolverControl        &cn,
                             VectorMemory<VECTOR> &mem,
@@ -288,6 +394,63 @@ SolverCG<VECTOR>::print_vectors(const unsigned int,
 
 
 
+template <class VECTOR>
+inline void
+SolverCG<VECTOR>::compute_eigs_and_cond(
+  const std::vector<double> &diagonal,
+  const std::vector<double> &offdiagonal,
+  const boost::signals2::signal<void (const std::vector<double> &)> &eigenvalues_signal,
+  const boost::signals2::signal<void (double)> &cond_signal,
+  const bool log_eigenvalues,
+  const bool log_cond)
+{
+  //Avoid computing eigenvalues unless they are needed.
+  if (!cond_signal.empty()|| !eigenvalues_signal.empty()  || log_cond ||
+      log_eigenvalues)
+    {
+      TridiagonalMatrix<double> T(diagonal.size(), true);
+      for (size_type i=0; i<diagonal.size(); ++i)
+        {
+          T(i,i) = diagonal[i];
+          if (i< diagonal.size()-1)
+            T(i,i+1) = offdiagonal[i];
+        }
+      T.compute_eigenvalues();
+      //Need two eigenvalues to estimate the condition number.
+      if (diagonal.size()>1)
+        {
+          double condition_number=T.eigenvalue(T.n()-1)/T.eigenvalue(0);
+          cond_signal(condition_number);
+          //Send to deallog
+          if (log_cond)
+            {
+              deallog << "Condition number estimate: " <<
+                      condition_number << std::endl;
+            }
+        }
+      //Avoid copying the eigenvalues of T to a vector unless a signal is
+      //connected.
+      if (!eigenvalues_signal.empty())
+        {
+          std::vector<double> eigenvalues(T.n());
+          for (unsigned int j = 0; j < T.n(); ++j)
+            {
+              eigenvalues.at(j)=T.eigenvalue(j);
+            }
+          eigenvalues_signal(eigenvalues);
+        }
+      if (log_eigenvalues)
+        {
+          for (size_type i=0; i<T.n(); ++i)
+            deallog << ' ' << T.eigenvalue(i);
+          deallog << std::endl;
+        }
+    }
+
+}
+
+
+
 template <class VECTOR>
 template <class MATRIX, class PRECONDITIONER>
 void
@@ -306,9 +469,13 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
   Vp = this->memory.alloc();
   // Should we build the matrix for
   // eigenvalue computations?
-  bool do_eigenvalues = additional_data.compute_condition_number
-                        | additional_data.compute_all_condition_numbers
-                        | additional_data.compute_eigenvalues;
+  const bool do_eigenvalues = !condition_number_signal.empty()
+                              |!all_condition_numbers_signal.empty()
+                              |!eigenvalues_signal.empty()
+                              |!all_eigenvalues_signal.empty()
+                              | additional_data.compute_condition_number
+                              | additional_data.compute_all_condition_numbers
+                              | additional_data.compute_eigenvalues;
   double eigen_beta_alpha = 0;
 
   // vectors used for eigenvalue
@@ -404,6 +571,7 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
               d.sadd(beta,-1.,g);
             }
 
+          this->coefficients_signal(alpha,beta);
           if (additional_data.log_coefficients)
             deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
           // set up the vectors
@@ -416,20 +584,9 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
               eigen_beta_alpha = beta/alpha;
               offdiagonal.push_back(std::sqrt(beta)/alpha);
             }
-
-          if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
-            {
-              TridiagonalMatrix<double> T(diagonal.size(), true);
-              for (size_type i=0; i<diagonal.size(); ++i)
-                {
-                  T(i,i) = diagonal[i];
-                  if (i< diagonal.size()-1)
-                    T(i,i+1) = offdiagonal[i];
-                }
-              T.compute_eigenvalues();
-              deallog << "Condition number estimate: " <<
-                      T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
-            }
+          compute_eigs_and_cond(diagonal,offdiagonal,all_eigenvalues_signal,
+                                all_condition_numbers_signal,false,
+                                additional_data.compute_all_condition_numbers);
         }
     }
   catch (...)
@@ -437,30 +594,11 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
       cleanup();
       throw;
     }
-
-  // Write eigenvalues or condition number
-  if (do_eigenvalues)
-    {
-      TridiagonalMatrix<double> T(diagonal.size(), true);
-      for (size_type i=0; i<diagonal.size(); ++i)
-        {
-          T(i,i) = diagonal[i];
-          if (i< diagonal.size()-1)
-            T(i,i+1) = offdiagonal[i];
-        }
-      T.compute_eigenvalues();
-      if (additional_data.compute_condition_number
-          && ! additional_data.compute_all_condition_numbers
-          && (diagonal.size() > 1))
-        deallog << "Condition number estimate: " <<
-                T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
-      if (additional_data.compute_eigenvalues)
-        {
-          for (size_type i=0; i<T.n(); ++i)
-            deallog << ' ' << T.eigenvalue(i);
-          deallog << std::endl;
-        }
-    }
+  compute_eigs_and_cond(diagonal,offdiagonal,eigenvalues_signal,
+                        condition_number_signal,
+                        additional_data.compute_eigenvalues,
+                        (additional_data.compute_condition_number &&
+                         !additional_data.compute_all_condition_numbers));
 
   // Deallocate Memory
   cleanup();
@@ -470,6 +608,54 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
   // otherwise exit as normal
 }
 
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_coefficients_slot(
+  const std_cxx11::function<void(double,double)> &slot)
+{
+  return coefficients_signal.connect(slot);
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_condition_number_slot(
+  const std_cxx11::function<void(double)> &slot,
+  const bool every_iteration)
+{
+  if (every_iteration)
+    {
+      return all_condition_numbers_signal.connect(slot);
+    }
+  else
+    {
+      return condition_number_signal.connect(slot);
+    }
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_eigenvalues_slot(
+  const std_cxx11::function<void (const std::vector<double> &)> &slot,
+  const bool every_iteration)
+{
+  if (every_iteration)
+    {
+      return all_eigenvalues_signal.connect(slot);
+    }
+  else
+    {
+      return eigenvalues_signal.connect(slot);
+    }
+}
+
+
+
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
index dbe037f2d2a53798e34d14df1c88a07b20d3fd9a..e7fe89d9df5f66ca1bea629966d704efaccb3b95 100644 (file)
@@ -30,6 +30,7 @@
 
 #include <vector>
 #include <cmath>
+#include <algorithm>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -154,6 +155,18 @@ namespace internal
  * to observe the progress of the iteration.
  *
  *
+ * <h3>Eigenvalue and condition number estimates</h3>
+ *
+ * This class can estimate eigenvalues and condition number during the solution
+ * process. This is done by creating the Hessenberg matrix during the inner
+ * iterations. The eigenvalues are estimated as the eigenvalues of the
+ * Hessenberg matrix and the condition number is estimated as the ratio of the
+ * largest and smallest singular value of the Hessenberg matrix. The estimates
+ * can be obtained by connecting a function as a slot using
+ * @p connect_condition_number_slot and @p connect_eigenvalues_slot. These slots
+ * will then be called from the solver with the estimates as argument.
+ *
+ *
  * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann.
  */
 template <class VECTOR = Vector<double> >
@@ -174,8 +187,18 @@ public:
     AdditionalData (const unsigned int max_n_tmp_vectors = 30,
                     const bool right_preconditioning = false,
                     const bool use_default_residual = true,
-                    const bool force_re_orthogonalization = false,
-                    const bool compute_eigenvalues = false);
+                    const bool force_re_orthogonalization = false);
+
+    /**
+     * Constructor.
+     * @deprecated To obtain the estimated eigenvalues instead use:
+     * connect_eigenvalues_slot
+     */
+    AdditionalData (const unsigned int max_n_tmp_vectors,
+                    const bool right_preconditioning,
+                    const bool use_default_residual,
+                    const bool force_re_orthogonalization,
+                    const bool compute_eigenvalues) DEAL_II_DEPRECATED;
 
     /**
      * Maximum number of temporary vectors. This parameter controls the size
@@ -242,6 +265,28 @@ public:
          const VECTOR         &b,
          const PRECONDITIONER &precondition);
 
+  /**
+   * Connect a slot to retrieve the estimated condition number.
+   * 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_condition_number_slot(const std_cxx11::function<void (double)> &slot,
+                                const bool every_iteration=false);
+
+  /**
+   * Connect a slot to retrieve the estimated eigenvalues.
+   * 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_eigenvalues_slot(
+    const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
+    const bool every_iteration=false);
+
+
   DeclException1 (ExcTooFewTmpVectors,
                   int,
                   << "The number of temporary vectors you gave ("
@@ -249,11 +294,37 @@ public:
                   << "any results, and much more for reasonable ones.");
 
 protected:
+
   /**
    * Includes the maximum number of tmp vectors.
    */
   AdditionalData additional_data;
 
+  /**
+   * Signal used to retrieve the estimated condition number.
+   * Called once when all iterations are ended.
+   */
+  boost::signals2::signal<void (double)> condition_number_signal;
+
+  /**
+   * Signal used to retrieve the estimated condition numbers.
+   * Called on each outer iteration.
+   */
+  boost::signals2::signal<void (double)> all_condition_numbers_signal;
+
+  /**
+   * Signal used to retrieve the estimated eigenvalues.
+   * Called once when all iterations are ended.
+   */
+  boost::signals2::signal<void (const std::vector<std::complex<double> > &)> eigenvalues_signal;
+
+  /**
+   * Signal used to retrieve the estimated eigenvalues.
+   * Called on each outer iteration.
+   */
+  boost::signals2::signal<void (const std::vector<std::complex<double> > &)> all_eigenvalues_signal;
+
+
   /**
    * Implementation of the computation of the norm of the residual.
    */
@@ -284,6 +355,21 @@ protected:
                          Vector<double>     &h,
                          bool               &re_orthogonalize);
 
+  /**
+    * Estimates the eigenvalues from the Hessenberg matrix, H_orig, generated
+    * during the inner iterations. Uses these estimate to compute the condition
+    * number. Calls the signals eigenvalues_signal and cond_signal with these
+    * estimates as arguments. Outputs the eigenvalues to deallog if
+    * log_eigenvalues is true.
+    */
+  static void
+  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(double)> &cond_signal,
+    const bool log_eigenvalues);
+
   /**
    * Projected system matrix
    */
@@ -457,6 +543,23 @@ namespace internal
 
 
 
+template <class VECTOR>
+inline
+SolverGMRES<VECTOR>::AdditionalData::
+AdditionalData (const unsigned int max_n_tmp_vectors,
+                const bool         right_preconditioning,
+                const bool         use_default_residual,
+                const bool         force_re_orthogonalization)
+  :
+  max_n_tmp_vectors(max_n_tmp_vectors),
+  right_preconditioning(right_preconditioning),
+  use_default_residual(use_default_residual),
+  force_re_orthogonalization(force_re_orthogonalization),
+  compute_eigenvalues(false)
+{}
+
+
+
 template <class VECTOR>
 inline
 SolverGMRES<VECTOR>::AdditionalData::
@@ -470,10 +573,11 @@ AdditionalData (const unsigned int max_n_tmp_vectors,
   right_preconditioning(right_preconditioning),
   use_default_residual(use_default_residual),
   force_re_orthogonalization(force_re_orthogonalization),
-  compute_eigenvalues (compute_eigenvalues)
+  compute_eigenvalues(compute_eigenvalues)
 {}
 
 
+
 template <class VECTOR>
 SolverGMRES<VECTOR>::SolverGMRES (SolverControl        &cn,
                                   VectorMemory<VECTOR> &mem,
@@ -584,6 +688,57 @@ SolverGMRES<VECTOR>::modified_gram_schmidt (const internal::SolverGMRES::TmpVect
 
 
 
+template<class VECTOR>
+inline void
+SolverGMRES<VECTOR>::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 (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 )
+    {
+      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);
+      //Avoid computing eigenvalues if they are not needed.
+      if (!eigenvalues_signal.empty() || log_eigenvalues )
+        {
+          //Copy mat so that we can compute svd below. Necessary since
+          //compute_eigenvalues will leave mat in state LAPACKSupport::unusable.
+          LAPACKFullMatrix<double> mat_eig(mat);
+          mat_eig.compute_eigenvalues();
+          std::vector<std::complex<double> > eigenvalues(dim);
+          for (unsigned int i=0; i<mat_eig.n(); ++i)
+            eigenvalues[i] = mat_eig.eigenvalue(i);
+          eigenvalues_signal(eigenvalues);
+          if (log_eigenvalues)
+            {
+              //Sort eigenvalues for nicer output.
+              std::sort(eigenvalues.begin(), eigenvalues.end(),
+                        internal::SolverGMRES::complex_less_pred);
+              deallog << "Eigenvalue estimate: ";
+              for (unsigned int i=0; i<mat_eig.n(); ++i)
+                deallog << ' ' << eigenvalues[i];
+              deallog << std::endl;
+            }
+        }
+      //Calculate condition number, avoid calculating the svd if a slot
+      //isn't connected. Need at least a 2-by-2 matrix to do the estimate.
+      if (!cond_signal.empty() && (mat.n()>1))
+        {
+          mat.compute_svd();
+          double condition_number=mat.singular_value(0)/mat.singular_value(mat.n()-1);
+          cond_signal(condition_number);
+        }
+    }
+}
+
+
+
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
@@ -610,10 +765,16 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
   // restart
   unsigned int accumulated_iterations = 0;
 
+  const bool do_eigenvalues=
+    !condition_number_signal.empty()
+    |!all_condition_numbers_signal.empty()
+    |!eigenvalues_signal.empty()
+    |!all_eigenvalues_signal.empty()
+    |additional_data.compute_eigenvalues;
   // for eigenvalue computation, need to collect the Hessenberg matrix (before
   // applying Givens rotations)
   FullMatrix<double> H_orig;
-  if (additional_data.compute_eigenvalues)
+  if (do_eigenvalues)
     H_orig.reinit(n_tmp_vectors, n_tmp_vectors-1);
 
   // matrix used for the orthogonalization process later
@@ -765,7 +926,7 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
 
           // for eigenvalues, get the resulting coefficients from the
           // orthogonalization process
-          if (additional_data.compute_eigenvalues)
+          if (do_eigenvalues)
             for (unsigned int i=0; i<dim+1; ++i)
               H_orig(i,inner_iteration) = h(i);
 
@@ -840,27 +1001,9 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
         for (unsigned int j=0; j<dim; ++j)
           H1(i,j) = H(i,j);
 
-      // compute eigenvalues from the Hessenberg matrix generated during the
-      // inner iterations
-      if (additional_data.compute_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);
-          mat.compute_eigenvalues();
-          std::vector<std::complex<double> > eigenvalues(dim);
-          for (unsigned int i=0; i<mat.n(); ++i)
-            eigenvalues[i] = mat.eigenvalue(i);
-
-          std::sort(eigenvalues.begin(), eigenvalues.end(),
-                    internal::SolverGMRES::complex_less_pred);
-
-          deallog << "Eigenvalue estimate: ";
-          for (unsigned int i=0; i<mat.n(); ++i)
-            deallog << ' ' << eigenvalues[i];
-          deallog << std::endl;
-        }
+      compute_eigs_and_cond(H_orig,dim,all_eigenvalues_signal,
+                            all_condition_numbers_signal,
+                            additional_data.compute_eigenvalues);
 
       H1.backward(h,gamma);
 
@@ -880,6 +1023,8 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
     }
   while (iteration_state == SolverControl::iterate);
 
+  compute_eigs_and_cond(H_orig,dim,eigenvalues_signal,condition_number_signal,
+                        false);
   if (!use_default_residual)
     {
       this->memory.free(r);
@@ -898,6 +1043,42 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
 
 
 
+template<class VECTOR>
+boost::signals2::connection
+SolverGMRES<VECTOR>::connect_condition_number_slot(
+  const std_cxx11::function<void(double)> &slot,
+  const bool every_iteration)
+{
+  if (every_iteration)
+    {
+      return all_condition_numbers_signal.connect(slot);
+    }
+  else
+    {
+      return condition_number_signal.connect(slot);
+    }
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverGMRES<VECTOR>::connect_eigenvalues_slot(
+  const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
+  const bool every_iteration)
+{
+  if (every_iteration)
+    {
+      return all_eigenvalues_signal.connect(slot);
+    }
+  else
+    {
+      return eigenvalues_signal.connect(slot);
+    }
+}
+
+
+
 template<class VECTOR>
 double
 SolverGMRES<VECTOR>::criterion ()
diff --git a/tests/lac/solver_signals.cc b/tests/lac/solver_signals.cc
new file mode 100644 (file)
index 0000000..a60dd2e
--- /dev/null
@@ -0,0 +1,137 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Connects slots to all signals in solver_cg and solver_gmres and writes all
+//output to deallog.
+
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <complex>
+#include <fstream>
+#include <iomanip>
+#include <string>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+
+
+void output_double_number(double input,const std::string &text)
+{
+  deallog<<text<< input<<std::endl;
+}
+
+void output_coefficients(double alpha,double beta)
+{
+  deallog<<"alpha: "<< alpha << " beta: "<< beta << std::endl;
+}
+
+template<class NUMBER>
+void output_eigenvalues(const std::vector<NUMBER> &eigenvalues,const std::string &text)
+{
+  deallog<< text;
+  for (unsigned int j = 0; j < eigenvalues.size(); ++j)
+    {
+      deallog<< ' ' << eigenvalues.at(j);
+    }
+  deallog << std::endl;
+}
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_solve( SOLVER &solver, const MATRIX &A,
+             VECTOR &u, VECTOR &f, const PRECONDITION &P)
+{
+  u = 0.;
+  f = 1.;
+  try
+    {
+      solver.solve(A,u,f,P);
+    }
+  catch (dealii::SolverControl::NoConvergence &e)
+    {
+      deallog << "Exception: " << e.get_exc_name() << std::endl;
+    }
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(4);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  SolverControl solver_control(100, 1.e-3);
+
+  unsigned int size=30;
+  unsigned int dim = (size-1)*(size-1);
+
+  // Make matrix
+  FDMatrix testproblem(size, size);
+  SparsityPattern structure(dim, dim, 5);
+  testproblem.five_point_structure(structure);
+  structure.compress();
+  SparseMatrix<double>  A(structure);
+  testproblem.five_point(A);
+
+  Vector<double>  f(dim);
+  Vector<double>  u(dim);
+  f = 1.;
+
+  try
+    {
+      SolverCG<> solver_cg(solver_control);
+      //Attach all possible slots.
+      solver_cg.connect_coefficients_slot(&output_coefficients);
+      solver_cg.connect_condition_number_slot(
+        std_cxx11::bind(output_double_number,std_cxx11::_1,"Condition number: "),true);
+      solver_cg.connect_condition_number_slot(
+        std_cxx11::bind(output_double_number,std_cxx11::_1,"Final condition number: "));
+      solver_cg.connect_eigenvalues_slot(
+        std_cxx11::bind(output_eigenvalues<double>,std_cxx11::_1,"Eigenvalues: "),true);
+      solver_cg.connect_eigenvalues_slot(
+        std_cxx11::bind(output_eigenvalues<double>,std_cxx11::_1,"Final Eigenvalues: "));
+      solver_cg.solve(A,u,f,PreconditionIdentity());
+
+      u=0;
+      SolverGMRES<> solver_gmres(solver_control);
+      //Attach all possible slots.
+      solver_gmres.connect_condition_number_slot(
+        std_cxx11::bind(output_double_number,std_cxx11::_1,"Condition number: "),true);
+      solver_gmres.connect_condition_number_slot(
+        std_cxx11::bind(output_double_number,std_cxx11::_1,"Final condition number: "));
+      solver_gmres.connect_eigenvalues_slot(
+        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.solve(A,u,f,PreconditionIdentity());
+    }
+  catch (std::exception &e)
+    {
+      std::cerr << "Exception: " << e.what() << std::endl;
+    }
+
+}
+
diff --git a/tests/lac/solver_signals.output b/tests/lac/solver_signals.output
new file mode 100644 (file)
index 0000000..4a591e2
--- /dev/null
@@ -0,0 +1,131 @@
+
+DEAL:cg::Starting value 29.00
+DEAL:cg::alpha: 7.250 beta: 6.750
+DEAL:cg::Eigenvalues:  0.1379
+DEAL:cg::alpha: 0.8749 beta: 0.8595
+DEAL:cg::Condition number: 29.00
+DEAL:cg::Eigenvalues:  0.07373 2.138
+DEAL:cg::alpha: 0.7611 beta: 1.026
+DEAL:cg::Condition number: 61.17
+DEAL:cg::Eigenvalues:  0.05344 1.186 3.269
+DEAL:cg::alpha: 0.4874 beta: 0.9388
+DEAL:cg::Condition number: 99.58
+DEAL:cg::Eigenvalues:  0.04473 0.8257 2.584 4.454
+DEAL:cg::alpha: 0.4926 beta: 0.7872
+DEAL:cg::Condition number: 154.6
+DEAL:cg::Eigenvalues:  0.03835 0.5853 1.901 3.408 5.931
+DEAL:cg::alpha: 0.4826 beta: 0.9174
+DEAL:cg::Condition number: 190.9
+DEAL:cg::Eigenvalues:  0.03430 0.4500 1.488 2.865 4.152 6.546
+DEAL:cg::alpha: 0.4614 beta: 0.7765
+DEAL:cg::Condition number: 222.3
+DEAL:cg::Eigenvalues:  0.03119 0.3550 1.172 2.327 3.506 5.277 6.935
+DEAL:cg::alpha: 0.4739 beta: 0.8809
+DEAL:cg::Condition number: 248.2
+DEAL:cg::Eigenvalues:  0.02893 0.2914 0.9558 1.935 3.052 4.057 5.895 7.180
+DEAL:cg::alpha: 0.4500 beta: 0.7691
+DEAL:cg::Condition number: 270.6
+DEAL:cg::Eigenvalues:  0.02716 0.2444 0.7908 1.617 2.606 3.576 4.940 6.324 7.352
+DEAL:cg::alpha: 0.4645 beta: 0.8436
+DEAL:cg::Condition number: 289.5
+DEAL:cg::Eigenvalues:  0.02581 0.2101 0.6685 1.374 2.250 3.186 4.018 5.501 6.632 7.474
+DEAL:cg::alpha: 0.4416 beta: 0.7523
+DEAL:cg::Condition number: 305.7
+DEAL:cg::Eigenvalues:  0.02475 0.1837 0.5730 1.179 1.949 2.805 3.628 4.739 5.912 6.863 7.565
+DEAL:cg::alpha: 0.4535 beta: 0.8013
+DEAL:cg::Condition number: 319.1
+DEAL:cg::Eigenvalues:  0.02393 0.1636 0.4988 1.025 1.704 2.483 3.287 4.000 5.239 6.229 7.040 7.634
+DEAL:cg::alpha: 0.4317 beta: 0.7212
+DEAL:cg::Condition number: 330.1
+DEAL:cg::Eigenvalues:  0.02329 0.1480 0.4398 0.8995 1.500 2.203 2.954 3.669 4.608 5.620 6.480 7.180 7.688
+DEAL:cg::alpha: 0.4389 beta: 0.7451
+DEAL:cg::Condition number: 338.9
+DEAL:cg::Eigenvalues:  0.02281 0.1359 0.3930 0.7982 1.333 1.967 2.663 3.368 3.993 5.055 5.929 6.682 7.291 7.731
+DEAL:cg::alpha: 0.4165 beta: 0.6656
+DEAL:cg::Condition number: 345.7
+DEAL:cg::Eigenvalues:  0.02246 0.1267 0.3557 0.7156 1.193 1.766 2.406 3.073 3.704 4.519 5.405 6.181 6.846 7.381 7.766
+DEAL:cg::alpha: 0.4162 beta: 0.6562
+DEAL:cg::Condition number: 350.8
+DEAL:cg::Eigenvalues:  0.02222 0.1198 0.3265 0.6488 1.078 1.597 2.184 2.810 3.436 3.994 4.921 5.699 6.391 6.983 7.456 7.795
+DEAL:cg::alpha: 0.3885 beta: 0.5568
+DEAL:cg::Condition number: 354.3
+DEAL:cg::Eigenvalues:  0.02207 0.1151 0.3044 0.5959 0.9836 1.455 1.993 2.576 3.174 3.738 4.462 5.246 5.946 6.568 7.097 7.518 7.818
+DEAL:cg::alpha: 0.3719 beta: 0.4876
+DEAL:cg::Condition number: 356.6
+DEAL:cg::Eigenvalues:  0.02198 0.1121 0.2894 0.5567 0.9102 1.341 1.835 2.376 2.942 3.504 4.009 4.829 5.526 6.159 6.720 7.194 7.571 7.838
+DEAL:cg::alpha: 0.3309 beta: 0.3590
+DEAL:cg::Condition number: 358.0
+DEAL:cg::Eigenvalues:  0.02194 0.1107 0.2807 0.5314 0.8581 1.254 1.708 2.209 2.739 3.277 3.780 4.443 5.138 5.770 6.345 6.851 7.279 7.616 7.856
+DEAL:cg::alpha: 0.3081 beta: 0.3265
+DEAL:cg::Condition number: 358.9
+DEAL:cg::Eigenvalues:  0.02193 0.1101 0.2769 0.5181 0.8261 1.193 1.612 2.073 2.567 3.078 3.582 4.046 4.774 5.401 5.979 6.503 6.963 7.350 7.655 7.870
+DEAL:cg::alpha: 0.3281 beta: 0.4916
+DEAL:cg::Condition number: 359.5
+DEAL:cg::Eigenvalues:  0.02193 0.1099 0.2751 0.5103 0.8018 1.136 1.507 1.920 2.374 2.854 3.342 3.801 4.403 5.032 5.612 6.147 6.630 7.053 7.407 7.686 7.882
+DEAL:cg::alpha: 0.4942 beta: 1.139
+DEAL:cg::Condition number: 360.0
+DEAL:cg::Eigenvalues:  0.02192 0.1097 0.2731 0.4961 0.7283 0.9649 1.297 1.703 2.151 2.624 3.106 3.579 4.005 4.676 5.257 5.798 6.296 6.744 7.134 7.459 7.714 7.893
+DEAL:cg::alpha: 0.5077 beta: 0.7756
+DEAL:cg::Condition number: 360.7
+DEAL:cg::Eigenvalues:  0.02192 0.1094 0.2667 0.3974 0.5549 0.8536 1.214 1.617 2.047 2.496 2.953 3.406 3.828 4.389 4.962 5.495 5.993 6.452 6.863 7.219 7.515 7.745 7.905
+DEAL:cg::alpha: 0.3174 beta: 0.3195
+DEAL:cg::Condition number: 361.3
+DEAL:cg::Eigenvalues:  0.02191 0.1092 0.2555 0.3298 0.5368 0.8410 1.198 1.589 1.999 2.417 2.838 3.263 3.681 4.082 4.676 5.200 5.695 6.157 6.581 6.962 7.292 7.565 7.775 7.918
+DEAL:cg::alpha: 0.3428 beta: 0.5580
+DEAL:cg::Condition number: 361.8
+DEAL:cg::Eigenvalues:  0.02191 0.1091 0.2468 0.3123 0.5328 0.8364 1.188 1.561 1.909 2.227 2.599 3.015 3.437 3.832 4.346 4.877 5.373 5.840 6.274 6.673 7.030 7.341 7.599 7.797 7.929
+DEAL:cg::alpha: 0.4717 beta: 0.7491
+DEAL:cg::Condition number: 362.3
+DEAL:cg::Eigenvalues:  0.02191 0.1090 0.2350 0.3000 0.5286 0.8253 1.084 1.255 1.619 2.030 2.452 2.872 3.286 3.684 4.057 4.611 5.100 5.561 5.996 6.399 6.766 7.095 7.384 7.627 7.815 7.940
+DEAL:cg::alpha: 0.3388 beta: 0.3950
+DEAL:cg::Condition number: 362.7
+DEAL:cg::Eigenvalues:  0.02191 0.1090 0.2247 0.2935 0.5238 0.7650 0.8817 1.210 1.599 2.006 2.410 2.796 3.168 3.535 3.880 4.367 4.845 5.300 5.733 6.141 6.519 6.859 7.159 7.422 7.648 7.827 7.948
+DEAL:cg::alpha: 0.3941 beta: 0.7196
+DEAL:cg::Condition number: 363.0
+DEAL:cg::Eigenvalues:  0.02191 0.1089 0.2179 0.2902 0.5179 0.6680 0.8533 1.202 1.587 1.962 2.224 2.520 2.908 3.304 3.685 4.034 4.554 5.012 5.447 5.860 6.248 6.609 6.934 7.216 7.456 7.665 7.836 7.954
+DEAL:cg::alpha: 0.4369 beta: 0.5564
+DEAL:cg::Condition number: 363.2
+DEAL:cg::Eigenvalues:  0.02191 0.1089 0.2098 0.2865 0.4941 0.5708 0.8420 1.185 1.408 1.625 2.023 2.432 2.829 3.206 3.564 3.891 4.348 4.792 5.216 5.621 6.005 6.364 6.700 7.008 7.277 7.497 7.683 7.844 7.959
+DEAL:cg::alpha: 0.2956 beta: 0.4380
+DEAL:cg::Condition number: 363.4
+DEAL:cg::Eigenvalues:  0.02191 0.1089 0.2053 0.2845 0.4628 0.5484 0.8375 1.147 1.267 1.608 2.013 2.416 2.789 3.117 3.433 3.760 4.096 4.563 4.984 5.387 5.774 6.138 6.475 6.783 7.069 7.327 7.536 7.702 7.850 7.962
+DEAL:cg::alpha: 0.4195 beta: 0.5047
+DEAL:cg::Condition number: 363.6
+DEAL:cg::Eigenvalues:  0.02191 0.1089 0.2013 0.2827 0.4285 0.5392 0.8312 1.030 1.217 1.596 1.966 2.127 2.453 2.844 3.219 3.571 3.891 4.326 4.752 5.158 5.548 5.919 6.273 6.603 6.899 7.157 7.393 7.599 7.745 7.863 7.968
+DEAL:cg::alpha: 0.3754 beta: 0.4222
+DEAL:cg::Condition number: 363.9
+DEAL:cg::Eigenvalues:  0.02191 0.1089 0.1988 0.2813 0.4042 0.5349 0.8174 0.9179 1.204 1.555 1.680 2.025 2.429 2.816 3.166 3.479 3.772 4.060 4.476 4.863 5.244 5.616 5.975 6.317 6.642 6.937 7.192 7.419 7.624 7.774 7.875 7.973
+DEAL:cg::alpha: 0.3396 beta: 0.4622
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1974 0.2804 0.3902 0.5326 0.7846 0.8676 1.195 1.432 1.619 2.015 2.410 2.694 2.903 3.243 3.585 3.895 4.315 4.721 5.100 5.446 5.759 6.058 6.362 6.668 6.957 7.207 7.429 7.631 7.784 7.881 7.975
+DEAL:cg::alpha: 0.3734 beta: 0.3948
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1965 0.2797 0.3787 0.5303 0.7245 0.8487 1.175 1.288 1.603 1.978 2.106 2.440 2.828 3.186 3.506 3.791 4.058 4.447 4.825 5.205 5.575 5.928 6.256 6.532 6.750 6.990 7.225 7.437 7.637 7.790 7.885 7.977
+DEAL:cg::alpha: 0.3032 beta: 0.3507
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1961 0.2793 0.3728 0.5288 0.6842 0.8432 1.132 1.234 1.594 1.825 2.032 2.430 2.810 3.115 3.339 3.621 3.910 4.317 4.694 5.015 5.308 5.633 5.971 6.303 6.614 6.873 7.063 7.252 7.446 7.640 7.794 7.887 7.977
+DEAL:cg::alpha: 0.3564 beta: 0.3733
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1959 0.2791 0.3690 0.5273 0.6520 0.8393 1.065 1.214 1.560 1.663 2.018 2.393 2.525 2.840 3.197 3.519 3.798 4.052 4.427 4.803 5.177 5.530 5.828 6.068 6.345 6.647 6.927 7.153 7.317 7.461 7.643 7.796 7.889 7.978
+DEAL:cg::alpha: 0.2812 beta: 0.2570
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1958 0.2790 0.3672 0.5263 0.6336 0.8366 1.012 1.205 1.475 1.618 2.005 2.190 2.442 2.825 3.171 3.450 3.687 3.940 4.335 4.684 4.956 5.250 5.593 5.932 6.235 6.447 6.680 6.950 7.188 7.382 7.493 7.647 7.798 7.889 7.978
+DEAL:cg::alpha: 0.3218 beta: 0.3291
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1958 0.2789 0.3664 0.5257 0.6233 0.8344 0.9748 1.199 1.393 1.607 1.947 2.052 2.430 2.791 2.941 3.216 3.532 3.807 4.050 4.415 4.787 5.148 5.445 5.682 5.979 6.297 6.580 6.762 6.975 7.206 7.413 7.545 7.653 7.799 7.890 7.978
+DEAL:cg::alpha: 0.2821 beta: 0.2217
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues:  0.02191 0.1088 0.1958 0.2789 0.3660 0.5253 0.6175 0.8323 0.9462 1.192 1.322 1.598 1.804 2.024 2.402 2.518 2.834 3.186 3.485 3.728 3.962 4.346 4.672 4.910 5.218 5.562 5.863 6.072 6.333 6.629 6.870 7.024 7.222 7.427 7.590 7.663 7.801 7.891 7.978
+DEAL:cg::Convergence step 40 value 0.0009064
+DEAL:cg::Final condition number: 364.1
+DEAL:cg::Final Eigenvalues:  0.02191 0.1088 0.1958 0.2789 0.3660 0.5253 0.6175 0.8323 0.9462 1.192 1.322 1.598 1.804 2.024 2.402 2.518 2.834 3.186 3.485 3.728 3.962 4.346 4.672 4.910 5.218 5.562 5.863 6.072 6.333 6.629 6.870 7.024 7.222 7.427 7.590 7.663 7.801 7.891 7.978
+DEAL:GMRES::Starting value 29.00
+DEAL:GMRES::Eigenvalues:  (7.954,0.000) (7.836,0.000) (7.665,0.000) (7.456,0.000) (7.216,0.000) (6.934,0.000) (6.609,0.000) (6.248,0.000) (5.860,0.000) (5.447,0.000) (5.012,0.000) (4.554,0.000) (0.02191,0.000) (0.1089,0.000) (0.2179,0.000) (0.2902,0.000) (0.5179,0.000) (0.6680,0.000) (0.8533,0.000) (1.202,0.000) (4.034,0.000) (3.685,0.000) (3.304,0.000) (1.587,0.000) (2.908,0.000) (2.520,0.000) (1.962,0.000) (2.224,0.000)
+DEAL:GMRES::Condition number: 363.0
+DEAL:GMRES::Eigenvalues:  (7.978,0.000) (7.891,0.000) (7.797,0.000) (7.647,0.000) (7.393,0.000) (7.129,0.000) (6.901,0.000) (6.553,0.000) (6.146,0.000) (5.741,0.000) (5.333,0.000) (4.910,0.000) (4.475,0.000) (4.075,0.000) (3.694,0.000) (3.205,0.000) (2.760,0.000) (0.02765,0.000) (0.1418,0.000) (0.2014,0.000) (0.3484,0.000) (0.5334,0.000) (0.6447,0.000) (1.032,0.000) (1.319,0.000) (2.310,0.000) (1.634,0.000) (1.937,0.000)
+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) (7.726,0.000) (7.170,0.000) (6.307,0.000)
+DEAL:GMRES::Condition number: 288.2
+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) (7.726,0.000) (7.170,0.000) (6.307,0.000)
+DEAL:GMRES::Final condition number: 288.2

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.