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>
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,
253 if (line_search_signal.empty())
259 const Number
g0 =
g * p;
262 "Function does not decrease along the current direction"));
272 ExcMessage(
"Function value is not decreasing"));
273 df =
std::max(
df, 100. * std::numeric_limits<Number>::epsilon());
282 [&](
const Number &
x_line) -> std::pair<Number, Number> {
287 return std::make_pair(f,
g_line);
291 const auto res = LineMinimization::line_search<Number>(
295 LineMinimization::poly_fit<Number>,
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();
384 if (additional_data.debug_output)
387 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)
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