// ---------------------------------------------------------------------
//
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
* way to give these additional data to the @p SolverSelector object for each
* solver which it may use.
*
+ *
+ * <h3>Observing the progress of linear solver iterations</h3>
+ *
+ * The Solver class, being the base class for all of the iterative solvers such
+ * as SolverCG, SolverGMRES, etc, provides the facilities by which actual solver
+ * implementations determine whether the iteration is converged, not yet
+ * converged, or has failed. Typically, this is done using an object of type
+ * SolverControl that is passed to the solver classes's constructors and from
+ * them down to the constructor of this base class. Every one of the tutorial
+ * programs that solves a linear problem (starting with step-3) uses this
+ * method and it is described in detail there. However, the underlying mechanism
+ * is more general and allows for many other uses to observe how the linear
+ * solver iterations progress.
+ *
+ * The basic approach is that the iterative solvers invoke a <i>signal</i> at
+ * the end of each iteration to determine whether the solution is converged.
+ * A signal is a class that has, conceptually, a list of pointers to functions
+ * and every time the signal is invoked, each of these functions are called.
+ * In the language of signals, the functions called are called <i>slots</i>
+ * and one can attach any number of slots to a signal.
+ * (The implementation of signals and slots we use here is the one from the
+ * BOOST.signals2 library.) A number of details may clarify what is happening
+ * underneath:
+ * - In reality, the signal object does not store pointers to functions, but
+ * function objects as slots. Each slot must conform to a particular
+ * signature: here, it is an object that can be called with three arguments
+ * (the number of the current linear iteration, the current residual, and
+ * the current iterate; more specifics are discussed in the documentation of
+ * the connect() function). A pointer to a function with this argument list
+ * satisfies the requirements, but you can also pass a member function
+ * whose <code>this</code> argument has been bound using the
+ * <code>std_cxx11::bind</code> mechanism (see the example below).
+ * - Each of the slots will return a value that indicates whether
+ * the iteration should continue, should stop because it has succeeded,
+ * or stop because it has failed. The return type of slots is therefore
+ * of type SolverControl::State. The returned values from all of the
+ * slots will then have to be combined before they are returned to the
+ * iterative solver that invoked the signal. The way this works is that
+ * if at least one slot returned SolverControl::failure, then the
+ * combined value is SolverControl::failure; otherwise, if at least
+ * one slot returned SolverControl::iterate, then this is going to be
+ * the return value of the signal; finally, only if all slots return
+ * SolverControl::success will the signal's return value be
+ * SolverControl::success.
+ * - It may of course be that a particular
+ * slot has been connected to the signal only to observe how the
+ * solution or a specific part of it converges, but has no particular
+ * opinion on whether the iteration should continue or not. In such
+ * cases, the slot should just return SolverControl::success, which
+ * is the weakest of all return values according to the rules laid
+ * out above.
+ *
+ * Given all this, it should now be obvious how the SolverControl
+ * object fits into this scheme: when a SolverControl object is passed
+ * to the constructor of the current class, we simply connect the
+ * SolverControl::check() function of that object as a slot to the
+ * signal we maintain here. In other words, since a Solver object
+ * is always constructed using a SolverControl object, there is
+ * always at least one slot associated with the signal, namely the
+ * one that determines convergence.
+ *
+ * On the other hand, using the connect() member function, it is
+ * possible to connect any number of other slots to the signal to
+ * observe whatever it is you want to observe. The connect() function
+ * also returns an object that describes the connection from the
+ * signal to the slot, and the corresponding BOOST functions then allow
+ * you to disconnect the slot if you want.
+ *
+ * An example may illuminate these issues. In the step-3 tutorial
+ * program, let us add a member function as follows to the main class:
+ * @code
+ * SolverControl::State
+ * Step3::write_intermediate_solution (const unsigned int iteration,
+ * const double , //check_value
+ * const Vector<double> ¤t_iterate) const
+ * {
+ * DataOut<2> data_out;
+ * data_out.attach_dof_handler (dof_handler);
+ * data_out.add_data_vector (current_iterate, "solution");
+ * data_out.build_patches ();
+ *
+ * std::ofstream output ((std::string("solution-")
+ * + Utilities::int_to_string(iteration,4) + ".vtu").c_str());
+ * data_out.write_vtu (output);
+ *
+ * return SolverControl::success;
+ * }
+ * @endcode
+ * The function satisfies the signature necessary to be a slot for the
+ * signal discussed above, with the exception that it is a member
+ * function and consequently requires a <code>this</code> pointer.
+ * What the function does is to take the vector given as last
+ * argument and write it into a file in VTU format with a file
+ * name derived from the number of the iteration.
+ *
+ * This function can then be hooked into the CG solver by modifying
+ * the <code>Step3::solve()</code> function as follows:
+ * @code
+ * void Step3::solve ()
+ * {
+ * SolverControl solver_control (1000, 1e-12);
+ * SolverCG<> solver (solver_control);
+ *
+ * solver.connect (std_cxx11::bind (&Step3::write_intermediate_solution,
+ * this,
+ * std_cxx11::_1,
+ * std_cxx11::_2,
+ * std_cxx11::_3));
+ * solver.solve (system_matrix, solution, system_rhs,
+ * PreconditionIdentity());
+ * }
+ * @endcode
+ * The use of <code>std_cxx11::bind</code> here ensures that we convert
+ * the member function with its three arguments plus the <code>this</code>
+ * pointer, to a function that only takes three arguments, by fixing
+ * the implicit <code>this</code> argument of the function to the
+ * <code>this</code> pointer in the current function.
+ *
+ * It is well understood that the CG method is a smoothing iteration (in
+ * the same way as the more commonly used Jacobi or SSOR iterations are
+ * smoothers). The code above therefore allows to observe how the solution
+ * becomes smoother and smoother in every iteration. This is best
+ * observed by initializing the solution vector with randomly distributed
+ * numbers in $[-1,1]$, using code such as
+ * @code
+ * for (unsigned int i=0; i<solution.size(); ++i)
+ * solution(i) = 2.*rand()/RAND_MAX-1;
+ * @endcode
+ * Using this, the slot will then generate files that when visualized
+ * look like this over the course of iterations zero to five:
+ * <table>
+ * <tr>
+ * <td>
+ * @image html "cg-monitor-smoothing-0.png"
+ * </td>
+ * <td>
+ * @image html "cg-monitor-smoothing-1.png"
+ * </td>
+ * <td>
+ * @image html "cg-monitor-smoothing-2.png"
+ * </td>
+ * </tr>
+ * <tr>
+ * <td>
+ * @image html "cg-monitor-smoothing-3.png"
+ * </td>
+ * <td>
+ * @image html "cg-monitor-smoothing-4.png"
+ * </td>
+ * <td>
+ * @image html "cg-monitor-smoothing-5.png"
+ * </td>
+ * </tr>
+ * </table>
+ *
* @ingroup Solvers
- * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997-2001, 2005
+ * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997-2001, 2005, 2014
*/
template <class VECTOR = Vector<double> >
class Solver : public Subscriptor
* as long as that of the solver
* object.
*/
- Solver (SolverControl &solver_control);
+ Solver (SolverControl &solver_control);
+
+ /**
+ * Connect a function object that will be called periodically within
+ * iterative solvers. This function is used to attach monitors to
+ * iterative solvers, either to determine when convergence has happened,
+ * or simply to observe the progress of an iteration. See the documentation
+ * of this class for more information.
+ *
+ * @param slot A function object specified here will, with
+ * each call, receive the number of the current iteration,
+ * the value that is used to check for convergence (typically the residual
+ * of the current iterate with respect to the linear system to be solved)
+ * and the currently best available guess for the current iterate.
+ * Note that some solvers do not update the approximate solution in every
+ * iteration but only after convergence or failure has been determined
+ * (GMRES is an example); in such cases, the vector passed as the last
+ * argument to the signal is simply the best approximate at the time
+ * the signal is called, but not the vector that will be returned if
+ * the signal's return value indicates that the iteration should be
+ * terminated. The function object must return a SolverControl::State
+ * value that indicates whether the iteration should continue, has
+ * failed, or has succeeded. The results of all connected functions
+ * will then be combined to determine what should happen with the
+ * iteration.
+ *
+ * @return A connection object that represents the connection from the
+ * signal to the function object. It can be used to disconnect the
+ * function object again from the signal. See the documentation of
+ * the BOOST Signals2 library for more information on connection
+ * management.
+ */
+ boost::signals2::connection
+ connect (const std_cxx11::function<SolverControl::State (const unsigned int iteration,
+ const double check_value,
+ const VECTOR ¤t_iterate)> &slot);
+
+
protected:
/**
// only takes two arguments, the iteration and the check_value, and so
// we simply ignore the third argument that is passed in whenever the
// signal is executed
- iteration_status.connect (std_cxx11::bind(&SolverControl::check,
- std_cxx11::ref(solver_control),
- std_cxx11::_1,
- std_cxx11::_2));
+ connect (std_cxx11::bind(&SolverControl::check,
+ std_cxx11::ref(solver_control),
+ std_cxx11::_1,
+ std_cxx11::_2));
}
// only takes two arguments, the iteration and the check_value, and so
// we simply ignore the third argument that is passed in whenever the
// signal is executed
- iteration_status.connect (std_cxx11::bind(&SolverControl::check,
- std_cxx11::ref(solver_control),
- std_cxx11::_1,
- std_cxx11::_2));
+ connect (std_cxx11::bind(&SolverControl::check,
+ std_cxx11::ref(solver_control),
+ std_cxx11::_1,
+ std_cxx11::_2));
}
+template<class VECTOR>
+inline
+boost::signals2::connection
+Solver<VECTOR>::
+connect (const std_cxx11::function<SolverControl::State (const unsigned int iteration,
+ const double check_value,
+ const VECTOR ¤t_iterate)> &slot)
+{
+ return iteration_status.connect (slot);
+}
+
+
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// use signals to monitor solutions converging
+
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#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/solver_minres.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/precondition.h>
+
+
+SolverControl::State monitor_norm (const unsigned int iteration,
+ const double check_value,
+ const Vector<double> ¤t_iterate)
+{
+ deallog << " -- " << iteration << ' ' << check_value << std::endl;
+ deallog << " Norm=" << current_iterate.l2_norm() << std::endl;
+ return SolverControl::success;
+}
+
+
+SolverControl::State monitor_mean (const unsigned int iteration,
+ const double check_value,
+ const Vector<double> ¤t_iterate)
+{
+ deallog << " Mean=" << current_iterate.mean_value() << std::endl;
+ return SolverControl::success;
+}
+
+
+
+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");
+// logfile.setf(std::ios::fixed);
+ deallog << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ GrowingVectorMemory<> mem;
+
+ SolverControl control(100, 1.e-3);
+
+ // create CG and GMRES solvers and attach monitors to it
+ SolverCG<> cg(control, mem);
+ cg.connect (&monitor_norm);
+ cg.connect (&monitor_mean);
+
+ SolverGMRES<> gmres(control, mem, 8);
+ gmres.connect (&monitor_norm);
+ gmres.connect (&monitor_mean);
+
+ for (unsigned int size=4; size <= 30; size *= 3)
+ {
+ unsigned int dim = (size-1)*(size-1);
+
+ deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+ // 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);
+
+ PreconditionSSOR<> prec_ssor;
+ prec_ssor.initialize(A, 1.2);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+ Vector<double> res(dim);
+
+ f = 1.;
+ u = 1.;
+
+ A.residual(res,u,f);
+ A.SOR(res);
+ res.add(1.,u);
+ A.SOR_step(u,f);
+ res.add(-1.,u);
+
+ deallog << "SOR-diff:" << res *res << std::endl;
+
+ try
+ {
+ check_solve(cg,A,u,f,prec_ssor);
+ check_solve(gmres,A,u,f,prec_ssor);
+ }
+ catch (std::exception &e)
+ {
+ std::cerr << "Exception: " << e.what() << std::endl;
+ }
+ }
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// use signals to monitor solutions converging. test that we can
+// disconnect from a signal
+
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#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/solver_minres.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/precondition.h>
+
+
+SolverControl::State monitor_norm (const unsigned int iteration,
+ const double check_value,
+ const Vector<double> ¤t_iterate)
+{
+ deallog << " -- " << iteration << ' ' << check_value << std::endl;
+ deallog << " Norm=" << current_iterate.l2_norm() << std::endl;
+ return SolverControl::success;
+}
+
+
+SolverControl::State monitor_mean (const unsigned int iteration,
+ const double check_value,
+ const Vector<double> ¤t_iterate)
+{
+ deallog << " Mean=" << current_iterate.mean_value() << std::endl;
+ return SolverControl::success;
+}
+
+
+
+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");
+// logfile.setf(std::ios::fixed);
+ deallog << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ GrowingVectorMemory<> mem;
+
+ SolverControl control(100, 1.e-3);
+
+ // create CG and GMRES solvers and attach monitors to it
+ SolverCG<> cg(control, mem);
+ boost::signals2::connection cg_c1 = cg.connect (&monitor_norm);
+ boost::signals2::connection cg_c2 = cg.connect (&monitor_mean);
+
+ for (unsigned int size=4; size <= 12; size *= 3)
+ {
+ unsigned int dim = (size-1)*(size-1);
+
+ deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+ // 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);
+
+ PreconditionSSOR<> prec_ssor;
+ prec_ssor.initialize(A, 1.2);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+ Vector<double> res(dim);
+
+ f = 1.;
+ u = 1.;
+
+ A.residual(res,u,f);
+ A.SOR(res);
+ res.add(1.,u);
+ A.SOR_step(u,f);
+ res.add(-1.,u);
+
+ deallog << "SOR-diff:" << res *res << std::endl;
+
+ try
+ {
+ check_solve(cg,A,u,f,prec_ssor);
+ }
+ catch (std::exception &e)
+ {
+ std::cerr << "Exception: " << e.what() << std::endl;
+ }
+
+ // disconnect one of the signals
+ cg_c1.disconnect();
+ }
+}
+
projects / dealii.git / commitdiff
