31 typename Solver::float_type delta_parameter_prev,
32 const typename Solver::value_type& uh_prev,
33 const typename Solver::value_type& duh_dparameter,
34 const typename Solver::float_type& duh_dparameter_sign,
37 typename Solver::value_type& uh_guess)
39 typedef typename Solver::value_type value_type;
40 typedef typename Solver::float_type float_type;
41 std::string name = F.parameter_name();
42 float_type delta_parameter = delta_parameter_prev;
43 if (p_err) *p_err <<
"#[continuation] delta_"<<name<<
"(0) = " << delta_parameter << std::endl;
44 for (
size_t k = 1, k_max = 10;
true; ++k) {
45 F.set_parameter (F.parameter() + delta_parameter);
46 value_type uh0 = uh_prev;
48 uh0 = uh0 + (duh_dparameter_sign*delta_parameter)*duh_dparameter;
50 F.update_derivative (uh0);
51 value_type m_r0h = F.residue (uh0);
52 value_type delta_uh0 = - F.derivative_solve (m_r0h);
53 value_type uh1 = uh0 + delta_uh0;
54 F.update_derivative (uh1);
55 value_type m_r1h = F.residue (uh1);
56 value_type delta_uh1 = - F.derivative_solve (m_r1h);
57 value_type uh2 = uh1 + delta_uh1;
58 value_type m_r2h = F.residue (uh2);
59 float_type r0 = F.dual_space_norm (m_r0h);
60 float_type r1 = F.dual_space_norm (m_r1h);
61 float_type r2 = F.dual_space_norm (m_r2h);
62 value_type uh_possible_guess = (r2 < r1) ? uh2 : uh1;
63 r0 = max(r0, std::numeric_limits<float_type>::epsilon());
64 if (sqr(r0) < opts.
tol && sqr(r1) < opts.
tol) {
67 float_type new_delta_parameter = delta_parameter*theta;
70 delta_parameter = new_delta_parameter;
71 uh_guess = uh_possible_guess;
72 if (p_err) *p_err <<
"#[continuation] prediction: too small residues (r0="<<r0<<
", r1="<<r1<<
"): increase delta_"<<name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
73 return delta_parameter;
75 float_type kappa0 = r1/r0;
76 float_type theta = sqrt(opts.
kappa/kappa0);
81 if (p_err) *p_err <<
"#[continuation] prediction: not-a-number in residues (r0="<<r0<<
", r1="<<r1<<
"): set theta="<<theta<<std::endl;
85 if (p_err) *p_err <<
"#[continuation] prediction: optimal rate (kappa="<<kappa0<<
", theta="<<theta<<
") with delta_"<<name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
86 uh_guess = uh_possible_guess;
87 return delta_parameter;
92 float_type new_delta_parameter = delta_parameter*theta;
97 new_delta_parameter = min (new_delta_parameter, opts.
theta_incr*delta_parameter_prev);
99 if (k >= k_max || fabs(new_delta_parameter - delta_parameter) < std::numeric_limits<float_type>::epsilon()) {
101 if (p_err) *p_err <<
"#[continuation] prediction: avoid infinite loop with delta_"<<name<<
"("<<k<<
") = " << delta_parameter<<std::endl;
102 uh_guess = uh_possible_guess;
103 return delta_parameter;
105 F.set_parameter (F.parameter() - delta_parameter);
106 delta_parameter = new_delta_parameter;
107 if (p_err) *p_err <<
"#[continuation] prediction: loop with delta_"<<name<<
"("<<k<<
") = " << new_delta_parameter<<std::endl;
Solver::float_type step_adjust(Solver &F, size_t n, typename Solver::float_type delta_parameter_prev, const typename Solver::value_type &uh_prev, const typename Solver::value_type &duh_dparameter, const typename Solver::float_type &duh_dparameter_sign, const continuation_option &opts, odiststream *p_err, typename Solver::value_type &uh_guess)