1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_solver_control_h
17 #define dealii_solver_control_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/subscriptor.h>
23
24 #include <vector>
25
26 DEAL_II_NAMESPACE_OPEN
27
28 #ifndef DOXYGEN
29 class ParameterHandler;
30 #endif
31
32 /*!@addtogroup Solvers */
33 /*@{*/
34
35 /**
36 * Control class to determine convergence of iterative solvers.
37 *
38 * Used by iterative methods to determine whether the iteration should be
39 * continued. To this end, the virtual function <tt>check()</tt> is called in
40 * each iteration with the current iteration step and the value indicating
41 * convergence (usually the residual).
42 *
43 * After the iteration has terminated, the functions last_value() and
44 * last_step() can be used to obtain information about the final state of the
45 * iteration.
46 *
47 * check() can be replaced in derived classes to allow for more sophisticated
48 * tests.
49 *
50 *
51 * <h3>State</h3> The return states of the check function are of type #State,
52 * which is an enum local to this class. It indicates the state the solver is
53 * in.
54 *
55 * The possible values of State are
56 * <ul>
57 * <li> <tt>iterate = 0</tt>: continue the iteration.
58 * <li> @p success: the goal is reached, the iterative method can terminate
59 * successfully.
60 * <li> @p failure: the iterative method should stop because convergence could
61 * not be achieved or at least was not achieved within the given maximal
62 * number of iterations.
63 * </ul>
64 */
65 class SolverControl : public Subscriptor
66 {
67 public:
68 /**
69 * Enum denoting the different states a solver can be in. See the general
70 * documentation of this class for more information.
71 */
72 enum State
73 {
74 /// Continue iteration
75 iterate = 0,
76 /// Stop iteration, goal reached
77 success,
78 /// Stop iteration, goal not reached
79 failure
80 };
81
82
83
84 /**
85 * Class to be thrown upon failing convergence of an iterative solver, when
86 * either the number of iterations exceeds the limit or the residual fails
87 * to reach the desired limit, e.g. in the case of a break-down.
88 *
89 * The residual in the last iteration, as well as the iteration number of
90 * the last step are stored in this object and can be recovered upon
91 * catching an exception of this class.
92 */
93
94 class NoConvergence : public dealii::ExceptionBase
95 {
96 public:
NoConvergence(const unsigned int last_step,const double last_residual)97 NoConvergence(const unsigned int last_step, const double last_residual)
98 : last_step(last_step)
99 , last_residual(last_residual)
100 {}
101
102 virtual ~NoConvergence() noexcept override = default;
103
104 virtual void
print_info(std::ostream & out)105 print_info(std::ostream &out) const override
106 {
107 out
108 << "Iterative method reported convergence failure in step " << last_step
109 << ". The residual in the last step was " << last_residual << ".\n\n"
110 << "This error message can indicate that you have simply not allowed "
111 << "a sufficiently large number of iterations for your iterative solver "
112 << "to converge. This often happens when you increase the size of your "
113 << "problem. In such cases, the last residual will likely still be very "
114 << "small, and you can make the error go away by increasing the allowed "
115 << "number of iterations when setting up the SolverControl object that "
116 << "determines the maximal number of iterations you allow."
117 << "\n\n"
118 << "The other situation where this error may occur is when your matrix "
119 << "is not invertible (e.g., your matrix has a null-space), or if you "
120 << "try to apply the wrong solver to a matrix (e.g., using CG for a "
121 << "matrix that is not symmetric or not positive definite). In these "
122 << "cases, the residual in the last iteration is likely going to be large."
123 << std::endl;
124 }
125
126 /**
127 * Iteration number of the last step.
128 */
129 const unsigned int last_step;
130
131 /**
132 * Residual in the last step.
133 */
134 const double last_residual;
135 };
136
137
138
139 /**
140 * Constructor. The parameters @p n and @p tol are the maximum number of
141 * iteration steps before failure and the tolerance to determine success of
142 * the iteration.
143 *
144 * @p log_history specifies whether the history (i.e. the value to be
145 * checked and the number of the iteration step) shall be printed to @p
146 * deallog stream. Default is: do not print. Similarly, @p log_result
147 * specifies the whether the final result is logged to @p deallog. Default
148 * is yes.
149 */
150 SolverControl(const unsigned int n = 100,
151 const double tol = 1.e-10,
152 const bool log_history = false,
153 const bool log_result = true);
154
155 /**
156 * Virtual destructor is needed as there are virtual functions in this
157 * class.
158 */
159 virtual ~SolverControl() override = default;
160
161 /**
162 * Interface to parameter file.
163 */
164 static void
165 declare_parameters(ParameterHandler ¶m);
166
167 /**
168 * Read parameters from file.
169 */
170 void
171 parse_parameters(ParameterHandler ¶m);
172
173 /**
174 * Decide about success or failure of an iteration. This function gets the
175 * current iteration step to determine, whether the allowed number of steps
176 * has been exceeded and returns @p failure in this case. If @p check_value
177 * is below the prescribed tolerance, it returns @p success. In all other
178 * cases @p iterate is returned to suggest continuation of the iterative
179 * procedure.
180 *
181 * The iteration is also aborted if the residual becomes a denormalized
182 * value (@p NaN).
183 *
184 * <tt>check()</tt> additionally preserves @p step and @p check_value. These
185 * values are accessible by <tt>last_value()</tt> and <tt>last_step()</tt>.
186 *
187 * Derived classes may overload this function, e.g. to log the convergence
188 * indicators (@p check_value) or to do other computations.
189 */
190 virtual State
191 check(const unsigned int step, const double check_value);
192
193 /**
194 * Return the result of the last check operation.
195 */
196 State
197 last_check() const;
198
199 /**
200 * Return the initial convergence criterion.
201 */
202 double
203 initial_value() const;
204
205 /**
206 * Return the convergence value of last iteration step for which @p check
207 * was called by the solver.
208 */
209 double
210 last_value() const;
211
212 /**
213 * Number of last iteration step.
214 */
215 unsigned int
216 last_step() const;
217
218 /**
219 * Maximum number of steps.
220 */
221 unsigned int
222 max_steps() const;
223
224 /**
225 * Change maximum number of steps.
226 */
227 unsigned int
228 set_max_steps(const unsigned int);
229
230 /**
231 * Enables the failure check. Solving is stopped with @p ReturnState @p
232 * failure if <tt>residual>failure_residual</tt> with
233 * <tt>failure_residual := rel_failure_residual*first_residual</tt>.
234 */
235 void
236 set_failure_criterion(const double rel_failure_residual);
237
238 /**
239 * Disables failure check and resets @p relative_failure_residual and @p
240 * failure_residual to zero.
241 */
242 void
243 clear_failure_criterion();
244
245 /**
246 * Tolerance.
247 */
248 double
249 tolerance() const;
250
251 /**
252 * Change tolerance.
253 */
254 double
255 set_tolerance(const double);
256
257 /**
258 * Enables writing residuals of each step into a vector for later analysis.
259 */
260 void
261 enable_history_data();
262
263 /**
264 * Provide read access to the collected residual data.
265 */
266 const std::vector<double> &
267 get_history_data() const;
268
269 /**
270 * Average error reduction over all steps.
271 *
272 * Requires enable_history_data()
273 */
274 double
275 average_reduction() const;
276 /**
277 * Error reduction of the last step; for stationary iterations, this
278 * approximates the norm of the iteration matrix.
279 *
280 * Requires enable_history_data()
281 */
282 double
283 final_reduction() const;
284
285 /**
286 * Error reduction of any iteration step.
287 *
288 * Requires enable_history_data()
289 */
290 double
291 step_reduction(unsigned int step) const;
292
293 /**
294 * Log each iteration step. Use @p log_frequency for skipping steps.
295 */
296 void
297 log_history(const bool);
298
299 /**
300 * Return the @p log_history flag.
301 */
302 bool
303 log_history() const;
304
305 /**
306 * Set logging frequency.
307 */
308 unsigned int
309 log_frequency(unsigned int);
310
311 /**
312 * Log start and end step.
313 */
314 void
315 log_result(const bool);
316
317 /**
318 * Return the @p log_result flag.
319 */
320 bool
321 log_result() const;
322
323 /**
324 * This exception is thrown if a function operating on the vector of history
325 * data of a SolverControl object id called, but storage of history data was
326 * not enabled by enable_history_data().
327 */
328 DeclException0(ExcHistoryDataRequired);
329
330 protected:
331 /**
332 * Maximum number of steps.
333 */
334 unsigned int maxsteps;
335
336 /**
337 * Prescribed tolerance to be achieved.
338 */
339 double tol;
340
341 /**
342 * Result of last check operation.
343 */
344 State lcheck;
345
346 /**
347 * Initial value.
348 */
349 double initial_val;
350
351 /**
352 * Last value of the convergence criterion.
353 */
354 double lvalue;
355
356 /**
357 * Last step.
358 */
359 unsigned int lstep;
360
361 /**
362 * Is set to @p true by @p set_failure_criterion and enables failure
363 * checking.
364 */
365 bool check_failure;
366
367 /**
368 * Stores the @p rel_failure_residual set by @p set_failure_criterion
369 */
370 double relative_failure_residual;
371
372 /**
373 * @p failure_residual equals the first residual multiplied by @p
374 * relative_crit set by @p set_failure_criterion (see there).
375 *
376 * Until the first residual is known it is 0.
377 */
378 double failure_residual;
379
380 /**
381 * Log convergence history to @p deallog.
382 */
383 bool m_log_history;
384
385 /**
386 * Log only every nth step.
387 */
388 unsigned int m_log_frequency;
389
390 /**
391 * Log iteration result to @p deallog. If true, after finishing the
392 * iteration, a statement about failure or success together with @p lstep
393 * and @p lvalue are logged.
394 */
395 bool m_log_result;
396
397 /**
398 * Control over the storage of history data. Set by enable_history_data().
399 */
400 bool history_data_enabled;
401
402 /**
403 * Vector storing the result after each iteration step for later statistical
404 * analysis.
405 *
406 * Use of this vector is enabled by enable_history_data().
407 */
408 std::vector<double> history_data;
409 };
410
411
412 /**
413 * Specialization of @p SolverControl which returns @p success if either the
414 * specified tolerance is achieved or if the initial residual (or whatever
415 * criterion was chosen by the solver class) is reduced by a given factor.
416 * This is useful in cases where you don't want to solve exactly, but rather
417 * want to gain two digits or if the maximal number of iterations is achieved.
418 * For example: The maximal number of iterations is 20, the reduction factor
419 * is 1% and the tolerance is 0.1%. The initial residual is 2.5. The process
420 * will break if 20 iteration are completed or the new residual is less then
421 * 2.5*1% or if it is less then 0.1%.
422 */
423 class ReductionControl : public SolverControl
424 {
425 public:
426 /**
427 * Constructor. Provide the reduction factor in addition to arguments that
428 * have the same meaning as those of the constructor of the SolverControl
429 * constructor.
430 */
431 ReductionControl(const unsigned int maxiter = 100,
432 const double tolerance = 1.e-10,
433 const double reduce = 1.e-2,
434 const bool log_history = false,
435 const bool log_result = true);
436
437 /**
438 * Initialize with a SolverControl object. The result will emulate
439 * SolverControl by setting @p reduce to zero.
440 */
441 ReductionControl(const SolverControl &c);
442
443 /**
444 * Assign a SolverControl object to ReductionControl. The result of the
445 * assignment will emulate SolverControl by setting @p reduce to zero.
446 */
447 ReductionControl &
448 operator=(const SolverControl &c);
449
450 /**
451 * Virtual destructor is needed as there are virtual functions in this
452 * class.
453 */
454 virtual ~ReductionControl() override = default;
455
456 /**
457 * Interface to parameter file.
458 */
459 static void
460 declare_parameters(ParameterHandler ¶m);
461
462 /**
463 * Read parameters from file.
464 */
465 void
466 parse_parameters(ParameterHandler ¶m);
467
468 /**
469 * Decide about success or failure of an iteration. This function calls the
470 * one in the base class, but sets the tolerance to <tt>reduction * initial
471 * value</tt> upon the first iteration.
472 */
473 virtual State
474 check(const unsigned int step, const double check_value) override;
475
476 /**
477 * Reduction factor.
478 */
479 double
480 reduction() const;
481
482 /**
483 * Change reduction factor.
484 */
485 double
486 set_reduction(const double);
487
488 protected:
489 /**
490 * Desired reduction factor.
491 */
492 double reduce;
493
494 /**
495 * Reduced tolerance. Stop iterations if either this value is achieved or if
496 * the base class indicates success.
497 */
498 double reduced_tol;
499 };
500
501 /**
502 * Specialization of @p SolverControl which returns @p success if a given
503 * number of iteration was performed, irrespective of the actual residual.
504 * This is useful in cases where you don't want to solve exactly, but rather
505 * want to perform a fixed number of iterations, e.g. in an inner solver. The
506 * arguments given to this class are exactly the same as for the SolverControl
507 * class and the solver terminates similarly when one of the given tolerance
508 * or the maximum iteration count were reached. The only difference to
509 * SolverControl is that the solver returns success in the latter case.
510 */
511 class IterationNumberControl : public SolverControl
512 {
513 public:
514 /**
515 * Constructor. Provide exactly the same arguments as the constructor of
516 * the SolverControl class.
517 */
518 IterationNumberControl(const unsigned int maxiter = 100,
519 const double tolerance = 1e-12,
520 const bool log_history = false,
521 const bool log_result = true);
522
523 /**
524 * Initialize with a SolverControl object. The result will emulate
525 * SolverControl by setting the reduction target to zero.
526 */
527 IterationNumberControl(const SolverControl &c);
528
529 /**
530 * Assign a SolverControl object to ReductionControl. The result of the
531 * assignment will emulate SolverControl by setting the reduction target to
532 * zero.
533 */
534 IterationNumberControl &
535 operator=(const SolverControl &c);
536
537 /**
538 * Virtual destructor is needed as there are virtual functions in this
539 * class.
540 */
541 virtual ~IterationNumberControl() override = default;
542
543 /**
544 * Decide about success or failure of an iteration. This function bases
545 * success solely on the fact if a given number of iterations was reached or
546 * the check value reached exactly zero.
547 */
548 virtual State
549 check(const unsigned int step, const double check_value) override;
550 };
551
552
553 /**
554 * Specialization of @p SolverControl which returns SolverControl::State::success if
555 * and only if a certain positive number of consecutive iterations satisfy the
556 * specified tolerance. This is useful in cases when solving nonlinear problems
557 * using inexact Hessian.
558 *
559 * For example: The requested number of consecutively converged iterations is 2,
560 * the tolerance is 0.2. The ConsecutiveControl will return
561 * SolverControl::State::success only at the last step in the sequence 0.5,
562 * 0.0005, 1.0, 0.05, 0.01.
563 */
564 class ConsecutiveControl : public SolverControl
565 {
566 public:
567 /**
568 * Constructor. @p n_consecutive_iterations is the number of
569 * consecutive iterations which should satisfy the prescribed tolerance for
570 * convergence. Other arguments have the same meaning as those of the
571 * constructor of the SolverControl.
572 */
573 ConsecutiveControl(const unsigned int maxiter = 100,
574 const double tolerance = 1.e-10,
575 const unsigned int n_consecutive_iterations = 2,
576 const bool log_history = false,
577 const bool log_result = false);
578
579 /**
580 * Initialize with a SolverControl object. The result will emulate
581 * SolverControl by setting @p n_consecutive_iterations to one.
582 */
583 ConsecutiveControl(const SolverControl &c);
584
585 /**
586 * Assign a SolverControl object to ConsecutiveControl. The result of the
587 * assignment will emulate SolverControl by setting @p n_consecutive_iterations
588 * to one.
589 */
590 ConsecutiveControl &
591 operator=(const SolverControl &c);
592
593 /**
594 * Virtual destructor is needed as there are virtual functions in this
595 * class.
596 */
597 virtual ~ConsecutiveControl() override = default;
598
599 /**
600 * Decide about success or failure of an iteration, see the class description
601 * above.
602 */
603 virtual State
604 check(const unsigned int step, const double check_value) override;
605
606 protected:
607 /**
608 * The number of consecutive iterations which should satisfy the prescribed
609 * tolerance for convergence.
610 */
611 unsigned int n_consecutive_iterations;
612
613 /**
614 * Counter for the number of consecutively converged iterations.
615 */
616 unsigned int n_converged_iterations;
617 };
618
619 /*@}*/
620 //---------------------------------------------------------------------------
621
622 #ifndef DOXYGEN
623
624 inline unsigned int
max_steps()625 SolverControl::max_steps() const
626 {
627 return maxsteps;
628 }
629
630
631
632 inline unsigned int
set_max_steps(const unsigned int newval)633 SolverControl::set_max_steps(const unsigned int newval)
634 {
635 unsigned int old = maxsteps;
636 maxsteps = newval;
637 return old;
638 }
639
640
641
642 inline void
set_failure_criterion(const double rel_failure_residual)643 SolverControl::set_failure_criterion(const double rel_failure_residual)
644 {
645 relative_failure_residual = rel_failure_residual;
646 check_failure = true;
647 }
648
649
650
651 inline void
clear_failure_criterion()652 SolverControl::clear_failure_criterion()
653 {
654 relative_failure_residual = 0;
655 failure_residual = 0;
656 check_failure = false;
657 }
658
659
660
661 inline double
tolerance()662 SolverControl::tolerance() const
663 {
664 return tol;
665 }
666
667
668
669 inline double
set_tolerance(const double t)670 SolverControl::set_tolerance(const double t)
671 {
672 double old = tol;
673 tol = t;
674 return old;
675 }
676
677
678 inline void
log_history(const bool newval)679 SolverControl::log_history(const bool newval)
680 {
681 m_log_history = newval;
682 }
683
684
685
686 inline bool
log_history()687 SolverControl::log_history() const
688 {
689 return m_log_history;
690 }
691
692
693 inline void
log_result(const bool newval)694 SolverControl::log_result(const bool newval)
695 {
696 m_log_result = newval;
697 }
698
699
700 inline bool
log_result()701 SolverControl::log_result() const
702 {
703 return m_log_result;
704 }
705
706
707 inline double
reduction()708 ReductionControl::reduction() const
709 {
710 return reduce;
711 }
712
713
714 inline double
set_reduction(const double t)715 ReductionControl::set_reduction(const double t)
716 {
717 double old = reduce;
718 reduce = t;
719 return old;
720 }
721
722 #endif // DOXYGEN
723
724 DEAL_II_NAMESPACE_CLOSE
725
726 #endif
727