15#ifndef dealii_solver_bfgs_h
16#define dealii_solver_bfgs_h
26#include <boost/signals2.hpp>
58template <
typename VectorType>
65 using Number =
typename VectorType::value_type;
112 const std::function<
Number(
const VectorType &x, VectorType &g)> &compute,
123 boost::signals2::connection
126 Number(
Number &f, VectorType &x, VectorType &g,
const VectorType &p)>
155 boost::signals2::connection
157 const std::function<
void(VectorType &g,
171 boost::signals2::signal<
172 Number(
Number &f, VectorType &x, VectorType &g,
const VectorType &p)>
178 boost::signals2::signal<void(VectorType &g,
188template <
typename VectorType>
190 const unsigned int max_history_size_,
191 const bool debug_output_)
192 : max_history_size(max_history_size_)
193 , debug_output(debug_output_)
198template <
typename VectorType>
200 const AdditionalData &
data)
202 , additional_data(
data)
207template <
typename VectorType>
208boost::signals2::connection
211 Number(Number &f, VectorType &x, VectorType &g,
const VectorType &p)> &slot)
213 Assert(line_search_signal.empty(),
214 ExcMessage(
"One should not attach more than one line search signal."));
215 return line_search_signal.connect(slot);
220template <
typename VectorType>
221boost::signals2::connection
223 const std::function<
void(VectorType &g,
227 Assert(preconditioner_signal.empty(),
229 "One should not attach more than one preconditioner signal."));
230 return preconditioner_signal.connect(slot);
235template <
typename VectorType>
238 const std::function<
typename VectorType::value_type(
const VectorType &x,
239 VectorType &f)> &compute,
249 bool first_step =
true;
253 if (line_search_signal.empty())
256 const auto default_line_min =
259 const Number g0 = g * p;
262 "Function does not decrease along the current direction"));
270 Number df = f_prev - f;
271 Assert(first_step || df >= 0.,
272 ExcMessage(
"Function value is not decreasing"));
273 df =
std::max(df, 100. * std::numeric_limits<Number>::epsilon());
276 (first_step ? 1. :
std::min(1., -1.01 * 2. * df / g0));
281 const auto line_func =
282 [&](
const Number &x_line) -> std::pair<Number, Number> {
286 const Number g_line = g * p;
287 return std::make_pair(f, g_line);
291 const auto res = LineMinimization::line_search<Number>(
295 LineMinimization::poly_fit<Number>,
305 this->connect_line_search_slot(default_line_min);
316 std::vector<Number> c1;
317 c1.reserve(additional_data.max_history_size);
332 conv = this->iteration_status(k, g.l2_norm(), x);
338 if (additional_data.debug_output)
339 deallog <<
"Iteration " << k <<
" history " << m << std::endl
340 <<
"f=" << f << std::endl;
346 for (
unsigned int i = 0; i < m; ++i)
348 c1[i] = rho[i] * (s[i] * p);
352 if (!preconditioner_signal.empty())
353 preconditioner_signal(p, s, y);
356 for (
int i = m - 1; i >= 0; --i)
359 const Number c2 = rho[i] * (y[i] * p);
360 p.add(c1[i] - c2, s[i]);
367 const Number alpha = line_search_signal(f, x, g, p)
372 if (additional_data.debug_output)
373 deallog <<
"Line search a=" << alpha <<
" f=" << f << std::endl;
377 const Number g_l2 = g.l2_norm();
378 conv = this->iteration_status(k, g_l2, x);
383 const Number curvature = s_k * y_k;
384 if (additional_data.debug_output)
385 deallog <<
"Curvature " << curvature << std::endl;
387 if (curvature > 0. && additional_data.max_history_size > 0)
391 rho.add(1. / curvature);
void add(const T &element)
typename VectorType::value_type Number
boost::signals2::connection connect_line_search_slot(const std::function< Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> &slot)
boost::signals2::connection connect_preconditioner_slot(const std::function< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> &slot)
void solve(const std::function< Number(const VectorType &x, VectorType &g)> &compute, VectorType &x)
boost::signals2::signal< void(VectorType &g, const FiniteSizeHistory< VectorType > &s, const FiniteSizeHistory< VectorType > &y)> preconditioner_signal
const AdditionalData additional_data
boost::signals2::signal< Number(Number &f, VectorType &x, VectorType &g, const VectorType &p)> line_search_signal
SolverBFGS(SolverControl &residual_control, const AdditionalData &data=AdditionalData())
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_history_size=5, const bool debug_output=false)
unsigned int max_history_size