16 #ifndef dealii_solver_bfgs_h
17 #define dealii_solver_bfgs_h
57 template <
typename VectorType>
64 using Number =
typename VectorType::value_type;
122 boost::signals2::connection
154 boost::signals2::connection
170 boost::signals2::signal<
187 template <
typename VectorType>
189 const unsigned int max_history_size_,
190 const bool debug_output_)
191 : max_history_size(max_history_size_)
192 , debug_output(debug_output_)
197 template <
typename VectorType>
199 const AdditionalData &data)
201 , additional_data(data)
206 template <
class VectorType>
207 boost::signals2::connection
212 Assert(line_search_signal.empty(),
213 ExcMessage(
"One should not attach more than one line search signal."));
214 return line_search_signal.connect(slot);
219 template <
class VectorType>
220 boost::signals2::connection
226 Assert(preconditioner_signal.empty(),
228 "One should not attach more than one preconditioner signal."));
229 return preconditioner_signal.connect(slot);
234 template <
typename VectorType>
237 const std::function<
typename VectorType::value_type(
const VectorType &x,
248 bool first_step =
true;
252 if (line_search_signal.empty())
255 const auto default_line_min =
258 const Number g0 = g * p;
261 "Function does not decrease along the current direction"));
269 Number df = f_prev - f;
270 Assert(first_step || df >= 0.,
271 ExcMessage(
"Function value is not decreasing"));
275 (first_step ? 1. :
std::min(1., -1.01 * 2. * df / g0));
280 const auto line_func =
281 [&](
const Number &x_line) -> std::pair<Number, Number> {
285 const Number g_line = g * p;
286 return std::make_pair(f, g_line);
290 const auto res = LineMinimization::line_search<Number>(
294 LineMinimization::poly_fit<Number>,
304 this->connect_line_search_slot(default_line_min);
315 std::vector<Number> c1;
316 c1.reserve(additional_data.max_history_size);
331 conv = this->iteration_status(k, g.l2_norm(), x);
337 if (additional_data.debug_output)
338 deallog <<
"Iteration " << k <<
" history " << m << std::endl
339 <<
"f=" << f << std::endl;
345 for (
unsigned int i = 0; i < m; ++i)
347 c1[i] = rho[i] * (s[i] * p);
351 if (!preconditioner_signal.empty())
352 preconditioner_signal(p, s, y);
355 for (
int i = m - 1; i >= 0; --i)
358 const Number c2 = rho[i] * (y[i] * p);
359 p.add(c1[i] - c2, s[i]);
366 const Number alpha = line_search_signal(f, x, g, p)
371 if (additional_data.debug_output)
372 deallog <<
"Line search a=" << alpha <<
" f=" << f << std::endl;
376 const Number g_l2 = g.l2_norm();
377 conv = this->iteration_status(k, g_l2, x);
382 const Number curvature = s_k * y_k;
383 if (additional_data.debug_output)
384 deallog <<
"Curvature " << curvature << std::endl;
386 if (curvature > 0. && additional_data.max_history_size > 0)
390 rho.add(1. / curvature);