15#ifndef dealii_solver_bfgs_h
16#define dealii_solver_bfgs_h
56template <
typename VectorType>
63 using Number =
typename VectorType::value_type;
110 const std::function<
Number(
const VectorType &
x, VectorType &
g)> &compute,
121 boost::signals2::connection
124 Number(
Number &f, VectorType &
x, VectorType &
g,
const VectorType &p)>
153 boost::signals2::connection
155 const std::function<
void(VectorType &
g,
169 boost::signals2::signal<
170 Number(
Number &f, VectorType &
x, VectorType &
g,
const VectorType &p)>
176 boost::signals2::signal<void(VectorType &
g,
186template <
typename VectorType>
196template <
typename VectorType>
198 const AdditionalData &data)
200 , additional_data(data)
205template <
typename VectorType>
206boost::signals2::connection
209 Number(Number &f, VectorType &
x, VectorType &
g,
const VectorType &p)> &
slot)
211 Assert(line_search_signal.empty(),
212 ExcMessage(
"One should not attach more than one line search signal."));
213 return line_search_signal.connect(
slot);
218template <
typename VectorType>
219boost::signals2::connection
221 const std::function<
void(VectorType &
g,
225 Assert(preconditioner_signal.empty(),
227 "One should not attach more than one preconditioner signal."));
228 return preconditioner_signal.connect(
slot);
233template <
typename VectorType>
236 const std::function<
typename VectorType::value_type(
const VectorType &
x,
237 VectorType &f)> &compute,
251 if (line_search_signal.empty())
255 [&](Number &f, VectorType &
x, VectorType &
g,
const VectorType &
p) {
257 const Number
g0 =
g *
p;
260 "Function does not decrease along the current direction"));
270 ExcMessage(
"Function value is not decreasing"));
271 df =
std::max(
df, 100. * std::numeric_limits<Number>::epsilon());
280 [&](
const Number &
x_line) -> std::pair<Number, Number> {
285 return std::make_pair(f,
g_line);
289 const auto res = LineMinimization::line_search<Number>(
293 LineMinimization::poly_fit<Number>,
314 std::vector<Number>
c1;
315 c1.reserve(additional_data.max_history_size);
330 conv = this->iteration_status(
k,
g.l2_norm(),
x);
336 if (additional_data.debug_output)
337 deallog <<
"Iteration " <<
k <<
" history " << m << std::endl
338 <<
"f=" << f << std::endl;
344 for (
unsigned int i = 0; i < m; ++i)
346 c1[i] =
rho[i] * (s[i] *
p);
350 if (!preconditioner_signal.empty())
351 preconditioner_signal(p, s,
y);
354 for (
int i = m - 1; i >= 0; --i)
357 const Number
c2 =
rho[i] * (
y[i] *
p);
358 p.add(
c1[i] -
c2, s[i]);
365 const Number
alpha = line_search_signal(f,
x,
g, p)
370 if (additional_data.debug_output)
371 deallog <<
"Line search a=" <<
alpha <<
" f=" << f << std::endl;
375 const Number
g_l2 =
g.l2_norm();
382 if (additional_data.debug_output)
385 if (
curvature > 0. && additional_data.max_history_size > 0)
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)
::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