From e711f0cfa120058e84ee16fc4ad935dd11b1b777 Mon Sep 17 00:00:00 2001 From: vishalkenchan Date: Tue, 11 Jul 2017 16:51:00 +0200 Subject: [PATCH] new SolverFIRE inertial relaxation algorithm --- doc/news/changes/minor/20170711VishalBoddu | 4 + include/deal.II/lac/solver_fire.h | 392 +++++++++++++++++++++ tests/lac/solver.cc | 12 + tests/lac/solver.output | 26 ++ tests/lac/solver_fire_01.cc | 79 +++++ tests/lac/solver_fire_01.output | 10 + tests/lac/solver_fire_rosenbrock_01.cc | 117 ++++++ tests/lac/solver_fire_rosenbrock_01.output | 10 + 8 files changed, 650 insertions(+) create mode 100644 doc/news/changes/minor/20170711VishalBoddu create mode 100644 include/deal.II/lac/solver_fire.h create mode 100644 tests/lac/solver_fire_01.cc create mode 100644 tests/lac/solver_fire_01.output create mode 100644 tests/lac/solver_fire_rosenbrock_01.cc create mode 100644 tests/lac/solver_fire_rosenbrock_01.output diff --git a/doc/news/changes/minor/20170711VishalBoddu b/doc/news/changes/minor/20170711VishalBoddu new file mode 100644 index 0000000000..c31f544988 --- /dev/null +++ b/doc/news/changes/minor/20170711VishalBoddu @@ -0,0 +1,4 @@ +New: SolverFIRE implements FIRE (Fast Inertial Relaxation Engine) for solving +the problem of minimization of a given objective function. +
+(Vishal Boddu, 2017/07/11) diff --git a/include/deal.II/lac/solver_fire.h b/include/deal.II/lac/solver_fire.h new file mode 100644 index 0000000000..76e5b021de --- /dev/null +++ b/include/deal.II/lac/solver_fire.h @@ -0,0 +1,392 @@ +#ifndef dealii__solver_fire_h +#define dealii__solver_fire_h + + +#include +#include +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup Solvers */ +/*@{*/ + +/** + * FIRE (Fast Inertial Relaxation Engine) for minimization of (potentially + * non-linear) objective function $E(U)$, $U$ is a vector of $n$ variables + * ($n$ is the number of variables of the objective function). + * Like all other solver classes, it can work on any kind of vector and matrix + * as long as they satisfy certain requirements (for the requirements on + * matrices and vectors in order to work with this class, see the documentation + * of the Solver base class). The type of the solution vector must be passed as + * template argument, and defaults to dealii::Vector. + * + * FIRE is a damped dynamics method described in + * Structural + * Relaxation Made Simple by Bitzek et al. 2006, typically used to find + * stable equilibrium configurations of atomistic systems in computational + * material science. Starting from a given initial configuration of the + * atomistic system, the algorithm relies on inertia to obtain (nearest) + * configuration with least potential energy. + * + * The problem of solving the system of linear equations $AU = F$ for $U$ + * can be re-casted into the problem of finding $U$ for which the quadratic + * function $\frac{1}{2} U^{T}AU - U^{T}F$ is minimized. + * + * Notation: + * The global vector of unknowns: $U$.
+ * Objective function: $E(U)$.
+ * Rate of change of unknowns: $V$.
+ * Gradient of the objective + * function w.r.t unknowns: $G$.
+ * Mass matrix: $M$.
+ * Initial guess of unknowns: $U_0$.
+ * Time step: $\Delta t$.
+ * + * Given initial values for $\Delta t$, $\alpha = \alpha_0$, $U = U_0$ and $V=0$ + * along with a given mass matrix $M$, + * FIRE algorithm is as follows, + * 1. Calculate $G = \nabla E(U)$ and check for convergence. + * Calculate $U$ and $V$ using any common Molecular Dynamics time integrator. + * Here we use simple Euler integration step,
+ * $U = U + \Delta t V$,
+ * $V = V + \Delta t M^{-1} \cdot G$. + * 2. Calculate $P = G \cdot V$. + * 3. Set $V = (1-\alpha) V + \alpha \frac{|V|}{|G|} G$. + * 4. If $P \leq 0$ and number of steps since P was non-negative is larger than + * certain value, then increase time step $\Delta t$ and decrease $\alpha$. + * 5. If $P>0$, then decrease the time step, freeze the system i.e., V = 0 and + * reset $\alpha = \alpha_0$. + * 6. Return to 1. + * + * Alse see + * + * Energy-Minimization in Atomic-to-Continuum Scale-Bridging Methods by + * Eidel et al. 2011. + * + * @author Vishal Boddu 2017 + */ +template > +class SolverFIRE : public Solver +{ + +public: + + struct AdditionalData + { + explicit + AdditionalData (const double initial_timestep = 0.1, + const double maximum_timestep = 1, + const double maximum_linfty_norm = 1); + + /** + * Initial time step. + */ + const double initial_timestep; + + /** + * Maximum time step. + */ + const double maximum_timestep; + + /** + * Maximum change allowed in any degree of freedom. + */ + const double maximum_linfty_norm; + + }; + + /** + * Constructor. + */ + SolverFIRE (SolverControl &solver_control, + VectorMemory &vector_memory, + const AdditionalData &data ); + + /** + * Constructor. Use an object of type GrowingVectorMemory as a default to + * allocate memory. + */ + SolverFIRE (SolverControl &solver_control, + const AdditionalData &data ); + + /** + * Virtual destructor. + */ + virtual ~SolverFIRE(); + + /** + * Obtain a set of #p u (variables) that minimize an objective function + * described by the polymorphic function wrapper @p compute, with a given + * preconditioner @p inverse_masses and initial @p u values. + * The function @p compute takes in the state of the (u) variables as + * argument and returns a pair of objective function's value and + * objective function's gradient (with respect to the variables). + */ + template> + void solve + (std::function compute, + VectorType &u, + const PreconditionerType &inverse_masses); + + /** + * Solve for x that minimizes the quadratic function + * $E(x) = \frac{1}{2} x^{T}Ax - x^{T}b$. + */ + template + void solve (const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &precondition); + +protected: + + /** + * Interface for derived class. This function gets the current iteration + * u, u's time derivative and the gradient in each step. It can be used + * for a graphical output of the convergence history. + */ + void print_vectors (const unsigned int, + const VectorType &, + const VectorType &, + const VectorType &) const; + + /** + * Additional parameters. + */ + const AdditionalData additional_data; + +}; + +/*@}*/ + +/*------------------------- Implementation ----------------------------*/ + +#ifndef DOXYGEN + +template +SolverFIRE::AdditionalData:: +AdditionalData (const double initial_timestep, + const double maximum_timestep, + const double maximum_linfty_norm) + : + initial_timestep(initial_timestep), + maximum_timestep(maximum_timestep), + maximum_linfty_norm(maximum_linfty_norm) +{ + AssertThrow (initial_timestep > 0. && + maximum_timestep > 0. && + maximum_linfty_norm > 0., + ExcMessage("Expected positive values for initial_timestep, " + "maximum_timestep and maximum_linfty_norm but one " + "or more of the these values are not positive.")); +} + + + +template +SolverFIRE:: +SolverFIRE (SolverControl &solver_control, + VectorMemory &vector_memory, + const AdditionalData &data ) + : + Solver(solver_control, vector_memory), + additional_data(data) +{} + + + +template +SolverFIRE:: +SolverFIRE (SolverControl &solver_control, + const AdditionalData &data ) + : + Solver(solver_control), + additional_data(data) +{} + + + +template +SolverFIRE::~SolverFIRE() +{} + + + +template +template +void +SolverFIRE::solve +(std::function compute, + VectorType &u, + const PreconditionerType &inverse_masses) +{ + deallog.push("FIRE"); + + // FIRE algorithm constants + const double DELAYSTEP = 5; + const double TIMESTEP_GROW = 1.1; + const double TIMESTEP_SHRINK = 0.5; + const double ALPHA_0 = 0.1; + const double ALPHA_SHRINK = 0.99; + + using real_type = typename VectorType::real_type; + + typename VectorMemory::Pointer v(this->memory); + typename VectorMemory::Pointer g(this->memory); + + // Set velocities to zero but not gradients + // as we are going to compute them soon. + v->reinit(u,false); + g->reinit(u,true); + + // Refer to v and g with some readable names. + VectorType &velocities = *v; + VectorType &gradients = *g; + + // Update gradients for the new u. + compute(gradients, u); + + unsigned int iter = 0; + + SolverControl::State conv = SolverControl::iterate; + conv = this->iteration_status (iter, gradients * gradients, u); + if (conv != SolverControl::iterate) + return; + + // Refer to additional data members with some readable names. + const auto &maximum_timestep = additional_data.maximum_timestep; + double timestep = additional_data.initial_timestep; + + // First scaling factor. + double alpha = ALPHA_0; + + unsigned int previous_iter_with_positive_v_dot_g = 0; + + while (conv == SolverControl::iterate) + { + ++iter; + // Euler integration step. + u.add (timestep, velocities); // U += dt * V + inverse_masses.vmult(gradients, gradients); // G = M^{-1} * G + velocities.add (-timestep, gradients); // V -= dt * G + + // Compute gradients for the new u. + compute(gradients, u); + + const real_type gradient_norm_squared = gradients * gradients; + conv = this->iteration_status(iter, gradient_norm_squared, u); + if (conv != SolverControl::iterate) + break; + + // v_dot_g = V * G + const real_type v_dot_g = velocities * gradients; + + if (v_dot_g < 0.) + { + const real_type velocities_norm_squared = + velocities * velocities; + + // Check if we divide by zero in DEBUG mode. + Assert (gradient_norm_squared > 0., ExcInternalError()); + + // beta = - alpha |V|/|G| + const real_type beta = -alpha * + std::sqrt (velocities_norm_squared + / + gradient_norm_squared); + + // V = (1-alpha) V + beta G. + velocities.sadd (1. - alpha, beta, gradients); + + if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP) + { + // Increase timestep and decrease alpha. + timestep = std::min (timestep*TIMESTEP_GROW, maximum_timestep); + alpha *= ALPHA_SHRINK; + } + } + else + { + // Decrease timestep, reset alpha and set V = 0. + previous_iter_with_positive_v_dot_g = iter; + timestep *= TIMESTEP_SHRINK; + alpha = ALPHA_0; + velocities = 0.; + } + + real_type vmax = velocities.linfty_norm(); + + // Change timestep if any dof would move more than maximum_linfty_norm. + if (vmax > 0.) + { + const double minimal_timestep = additional_data.maximum_linfty_norm + / + vmax; + if (minimal_timestep < timestep) + timestep = minimal_timestep; + } + + print_vectors(iter, u, velocities, gradients); + + } // While we need to iterate. + + deallog.pop(); + + // In the case of failure: throw exception. + if (conv != SolverControl::success) + AssertThrow (false, + SolverControl::NoConvergence (iter, gradients * gradients)); + +} + + + +template +template +void SolverFIRE::solve (const MatrixType &A, + VectorType &x, + const VectorType &b, + const PreconditionerType &precondition) +{ + + std::function compute_func = + [&] (decltype(x) &G, decltype(b) &x) -> double + { + // Residual of the quadratic form $ \frac{1}{2} xAx - xb $. + // G = b - Ax + A.residual(G, x, b); + + // Gradient G = Ax -b. + G *= -1.; + + // The quadratic form $ \frac{1}{2} x^t A x - x^{t} b $. + return 0.5*A.matrix_norm_square(x) - x*b; + }; + + this->solve (compute_func, x, precondition); +} + + + +template +void +SolverFIRE:: +print_vectors (const unsigned int, + const VectorType &, + const VectorType &, + const VectorType &) const +{} + + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index 582bc7a4a4..ef88d74f65 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -93,6 +94,8 @@ int main() SolverBicgstab<> bicgstab(control, mem); SolverRichardson<> rich(control, mem); SolverQMRS<> qmrs(control, mem); + SolverFIRE<>::AdditionalData fire_data(0.1, 1, 1); + SolverFIRE<> fire(control, mem, fire_data); for (unsigned int size=4; size <= 30; size *= 3) { @@ -165,6 +168,10 @@ int main() check_solve(gmresright,A,u,f,prec_no); // check_solve(minres,A,u,f,prec_no); check_solve(qmrs,A,u,f,prec_no); + + control.set_max_steps(50); + check_solve(fire,A,u,f,prec_no); + control.set_max_steps(100); deallog.pop(); @@ -178,6 +185,7 @@ int main() check_solve(gmres,A,u,f,prec_no); check_solve(gmresright,A,u,f,prec_no); check_solve(qmrs,A,u,f,prec_no); + check_solve(fire,A,u,f,prec_no); rich.set_omega(1.); deallog.pop(); @@ -191,6 +199,7 @@ int main() check_solve(gmres,A,u,f,prec_richardson); check_solve(gmresright,A,u,f,prec_richardson); check_solve(qmrs,A,u,f,prec_richardson); + check_solve(fire,A,u,f,prec_richardson); rich.set_omega(1.); deallog.pop(); @@ -204,6 +213,7 @@ int main() check_solve(gmres,A,u,f,prec_ssor); check_solve(gmresright,A,u,f,prec_ssor); check_solve(qmrs,A,u,f,prec_ssor); + check_solve(fire,A,u,f,prec_ssor); deallog.pop(); @@ -215,6 +225,7 @@ int main() check_solve(bicgstab,A,u,f,prec_sor); check_solve(gmres,A,u,f,prec_sor); check_solve(gmresright,A,u,f,prec_sor); + check_solve(fire,A,u,f,prec_sor); deallog.pop(); @@ -226,6 +237,7 @@ int main() check_solve(bicgstab,A,u,f,prec_psor); check_solve(gmres,A,u,f,prec_psor); check_solve(gmresright,A,u,f,prec_psor); + check_solve(fire,A,u,f,prec_psor); deallog.pop(); } diff --git a/tests/lac/solver.output b/tests/lac/solver.output index d77ff53976..ef1e643dcb 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -11,6 +11,8 @@ DEAL:no-fail:GMRES::Starting value 3.000 DEAL:no-fail:GMRES::Convergence step 3 value 0 DEAL:no-fail:QMRS::Starting value 3.000 DEAL:no-fail:QMRS::Convergence step 3 value 0 +DEAL:no-fail:FIRE::Starting value 9.000 +DEAL:no-fail:FIRE::Convergence step 38 value 0.0002184 DEAL:no:Richardson::Starting value 3.000 DEAL:no:Richardson::Convergence step 24 value 0.0007118 DEAL:no:cg::Starting value 3.000 @@ -23,6 +25,8 @@ DEAL:no:GMRES::Starting value 3.000 DEAL:no:GMRES::Convergence step 3 value 0 DEAL:no:QMRS::Starting value 3.000 DEAL:no:QMRS::Convergence step 3 value 0 +DEAL:no:FIRE::Starting value 9.000 +DEAL:no:FIRE::Convergence step 38 value 0.0002184 DEAL:rich:Richardson::Starting value 3.000 DEAL:rich:Richardson::Convergence step 42 value 0.0008696 DEAL:rich:cg::Starting value 3.000 @@ -35,6 +39,8 @@ DEAL:rich:GMRES::Starting value 3.000 DEAL:rich:GMRES::Convergence step 3 value 0 DEAL:rich:QMRS::Starting value 3.000 DEAL:rich:QMRS::Convergence step 3 value 0 +DEAL:rich:FIRE::Starting value 9.000 +DEAL:rich:FIRE::Convergence step 32 value 0.0009708 DEAL:ssor:RichardsonT::Starting value 3.000 DEAL:ssor:RichardsonT::Convergence step 7 value 0.0006128 DEAL:ssor:Richardson::Starting value 3.000 @@ -49,6 +55,8 @@ DEAL:ssor:GMRES::Starting value 3.000 DEAL:ssor:GMRES::Convergence step 4 value 5.994e-05 DEAL:ssor:QMRS::Starting value 3.000 DEAL:ssor:QMRS::Convergence step 4 value 0.0001345 +DEAL:ssor:FIRE::Starting value 9.000 +DEAL:ssor:FIRE::Convergence step 38 value 0.0008314 DEAL:sor:RichardsonT::Starting value 3.000 DEAL:sor:RichardsonT::Convergence step 7 value 0.0004339 DEAL:sor:Richardson::Starting value 3.000 @@ -64,6 +72,8 @@ DEAL:sor:GMRES::Convergence step 5 value 0 DEAL:sor:GMRES::Starting value 3.000 DEAL:sor:GMRES::Re-orthogonalization enabled at step 5 DEAL:sor:GMRES::Convergence step 5 value 0 +DEAL:sor:FIRE::Starting value 9.000 +DEAL:sor:FIRE::Convergence step 51 value 0.0004993 DEAL:psor:RichardsonT::Starting value 3.000 DEAL:psor:RichardsonT::Convergence step 8 value 0.0004237 DEAL:psor:Richardson::Starting value 3.000 @@ -77,6 +87,8 @@ DEAL:psor:GMRES::Starting value 1.467 DEAL:psor:GMRES::Convergence step 5 value 0.0006649 DEAL:psor:GMRES::Starting value 3.000 DEAL:psor:GMRES::Convergence step 6 value 0.0004455 +DEAL:psor:FIRE::Starting value 9.000 +DEAL:psor:FIRE::Convergence step 45 value 0.0009217 DEAL::Size 12 Unknowns 121 DEAL::SOR-diff:0 DEAL:no-fail:cg::Starting value 11.00 @@ -95,6 +107,9 @@ DEAL:no-fail::Exception: SolverControl::NoConvergence (accumulated_iterations, l DEAL:no-fail:QMRS::Starting value 11.00 DEAL:no-fail:QMRS::Failure step 10 value 0.4215 DEAL:no-fail::Exception: SolverControl::NoConvergence (step, state.last_residual) +DEAL:no-fail:FIRE::Starting value 121.0 +DEAL:no-fail:FIRE::Failure step 50 value 2.971 +DEAL:no-fail::Exception: SolverControl::NoConvergence (iter, gradients * gradients) DEAL:no:Richardson::Starting value 11.00 DEAL:no:Richardson::Failure step 100 value 0.3002 DEAL:no::Exception: SolverControl::NoConvergence (iter, last_criterion) @@ -108,6 +123,9 @@ DEAL:no:GMRES::Starting value 11.00 DEAL:no:GMRES::Convergence step 43 value 0.0009788 DEAL:no:QMRS::Starting value 11.00 DEAL:no:QMRS::Convergence step 16 value 0.0002583 +DEAL:no:FIRE::Starting value 121.0 +DEAL:no:FIRE::Failure step 100 value 0.003317 +DEAL:no::Exception: SolverControl::NoConvergence (iter, gradients * gradients) DEAL:rich:Richardson::Starting value 11.00 DEAL:rich:Richardson::Failure step 100 value 1.219 DEAL:rich::Exception: SolverControl::NoConvergence (iter, last_criterion) @@ -121,6 +139,8 @@ DEAL:rich:GMRES::Starting value 11.00 DEAL:rich:GMRES::Convergence step 43 value 0.0009788 DEAL:rich:QMRS::Starting value 11.00 DEAL:rich:QMRS::Convergence step 16 value 0.0002583 +DEAL:rich:FIRE::Starting value 121.0 +DEAL:rich:FIRE::Convergence step 98 value 0.0009234 DEAL:ssor:RichardsonT::Starting value 11.00 DEAL:ssor:RichardsonT::Convergence step 48 value 0.0009502 DEAL:ssor:Richardson::Starting value 11.00 @@ -135,6 +155,8 @@ DEAL:ssor:GMRES::Starting value 11.00 DEAL:ssor:GMRES::Convergence step 8 value 0.0004929 DEAL:ssor:QMRS::Starting value 11.00 DEAL:ssor:QMRS::Convergence step 9 value 0.0004140 +DEAL:ssor:FIRE::Starting value 121.0 +DEAL:ssor:FIRE::Convergence step 65 value 0.0006547 DEAL:sor:RichardsonT::Starting value 11.00 DEAL:sor:RichardsonT::Convergence step 88 value 0.0009636 DEAL:sor:Richardson::Starting value 11.00 @@ -148,6 +170,8 @@ DEAL:sor:GMRES::Starting value 7.322 DEAL:sor:GMRES::Convergence step 23 value 0.0006981 DEAL:sor:GMRES::Starting value 11.00 DEAL:sor:GMRES::Convergence step 24 value 0.0007120 +DEAL:sor:FIRE::Starting value 121.0 +DEAL:sor:FIRE::Convergence step 92 value 0.0009080 DEAL:psor:RichardsonT::Starting value 11.00 DEAL:psor:RichardsonT::Convergence step 89 value 0.0009736 DEAL:psor:Richardson::Starting value 11.00 @@ -161,3 +185,5 @@ DEAL:psor:GMRES::Starting value 7.345 DEAL:psor:GMRES::Convergence step 20 value 0.0008491 DEAL:psor:GMRES::Starting value 11.00 DEAL:psor:GMRES::Convergence step 23 value 0.0007600 +DEAL:psor:FIRE::Starting value 121.0 +DEAL:psor:FIRE::Convergence step 93 value 0.0006976 diff --git a/tests/lac/solver_fire_01.cc b/tests/lac/solver_fire_01.cc new file mode 100644 index 0000000000..f2f1268826 --- /dev/null +++ b/tests/lac/solver_fire_01.cc @@ -0,0 +1,79 @@ + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +// Test to verify correctness of SolverFIRE::solve() +// The objective function is f(x,y) = x^2 + y^2. + + +using vector_t = typename dealii::Vector; + + +double compute (vector_t &G, const vector_t &X) +{ + AssertThrow (X.size() == 2 && G.size() == 2, + ExcInternalError()); + + G(0) = 2*X(0); + G(1) = 2*X(1); + + return X.norm_sqr(); +} + + + +void check_value (const double x, + const double y, + const double tol) +{ + vector_t X; + + X.reinit(2, true); + + // Use this to initialize DiagonalMatrix + X = 1.; + + // Create inverse diagonal matrix. + DiagonalMatrix inv_mass; + inv_mass.reinit(X); + + // Set initial iterate. + X(0) = x; + X(1) = y; + + auto additional_data = + SolverFIRE::AdditionalData(0.1, 1., 1); + + SolverControl solver_control (200, tol); + + SolverFIRE fire (solver_control, additional_data); + + fire.solve(compute, X, inv_mass); + + deallog << "FIRE::Solution vector: "; + + X.print(deallog); +} + +int main () +{ + std::ofstream logfile("output"); +// logfile.setf(std::ios::fixed); + deallog << std::setprecision(4); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + check_value( 10, -2, 1e-15); + check_value(-0.1, 0.1, 1e-15); + check_value( 9.1,-6.1, 1e-15); + +} diff --git a/tests/lac/solver_fire_01.output b/tests/lac/solver_fire_01.output new file mode 100644 index 0000000000..ea7264d829 --- /dev/null +++ b/tests/lac/solver_fire_01.output @@ -0,0 +1,10 @@ + +DEAL:FIRE::Starting value 416.0 +DEAL:FIRE::Convergence step 101 value 0 +DEAL::FIRE::Solution vector: 0 0 +DEAL:FIRE::Starting value 0.08000 +DEAL:FIRE::Convergence step 89 value 0 +DEAL::FIRE::Solution vector: 0 0 +DEAL:FIRE::Starting value 480.1 +DEAL:FIRE::Convergence step 86 value 0 +DEAL::FIRE::Solution vector: 0 0 diff --git a/tests/lac/solver_fire_rosenbrock_01.cc b/tests/lac/solver_fire_rosenbrock_01.cc new file mode 100644 index 0000000000..51c8b3fefc --- /dev/null +++ b/tests/lac/solver_fire_rosenbrock_01.cc @@ -0,0 +1,117 @@ + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +// Test to verify correctness of SolverFIRE::sovle() +// The objective function is the extended Rosenbrock function. +// The Rosenbrock function is a non-convex function used as a test problem +// for optimization algorithms introduced by Howard H. Rosenbrock. +// +// f(X) = f(x_0, x_1, ..., x_{N-1}) +// +// = \sum_{i=0}^{\frac{N}{2} -1} +// +// \left[ +// a ( x_{2i}^2 - x_{2i+1} )^2 +// + +// b ( x_{2i} - 1 )^2 +// \right], +// +// where N is even and a = 100 and b = 1. +// +// DOI: 10.1007/BF02196600 + + +using vector_t = typename dealii::Vector; + + +double compute (vector_t &G, const vector_t &X) +{ + AssertThrow (X.size() % 2 == 0, + ExcInternalError()); + + double value = 0.; + + // Value of the objective function. + for (unsigned int i = 0; i < X.size()/2; ++i) + value += 100 * + dealii::Utilities::fixed_power<2>( X(2*i) * X(2*i) - X(2*i+1) ) + + + dealii::Utilities::fixed_power<2>( X(2*i) - 1 ); + + // Gradient of the objective function. + for (unsigned int i = 0; i < X.size()/2; ++i) + { + G(2*i) = ( X(2*i) * X(2*i) - X(2*i+1) ) * X(2*i) * 400 + + + ( X(2*i) - 1 ) * 2; + + G(2*i+1) = ( X(2*i) * X(2*i) - X(2*i+1) ) * -200; + } + + return value; +} + + + +void check_value (const unsigned int N, + const double tol) +{ + AssertThrow (N % 2 == 0, + ExcInternalError()); + + vector_t X (N); + + // Use this to initialize DiagonalMatrix + X = 1.; + + // Create inverse diagonal matrix. + DiagonalMatrix inv_mass; + inv_mass.reinit(X); + + // Set initial guess. + for (unsigned int i=0; i < N/2; i++) + { + X(2*i) = -1.2; + X(2*i+1) = 1.0; + } + + auto additional_data = + SolverFIRE::AdditionalData(0.1, 1, 1); + + SolverControl solver_control (1e5, tol); + + SolverFIRE fire (solver_control, additional_data); + + fire.solve(compute, X, inv_mass); + + deallog << "FIRE::Solution vector: "; + + X.print(deallog); +} + + + +int main () +{ + std::ofstream logfile("output"); +// logfile.setf(std::ios::fixed); + deallog << std::setprecision(4); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + check_value( 2, 1e-14); + check_value(10, 1e-14); + check_value(20, 1e-14); + +} diff --git a/tests/lac/solver_fire_rosenbrock_01.output b/tests/lac/solver_fire_rosenbrock_01.output new file mode 100644 index 0000000000..ee4a049472 --- /dev/null +++ b/tests/lac/solver_fire_rosenbrock_01.output @@ -0,0 +1,10 @@ + +DEAL:FIRE::Starting value 5.423e+04 +DEAL:FIRE::Convergence step 17364 value 0 +DEAL::FIRE::Solution vector: 1.000 1.000 +DEAL:FIRE::Starting value 2.711e+05 +DEAL:FIRE::Convergence step 18427 value 0 +DEAL::FIRE::Solution vector: 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 +DEAL:FIRE::Starting value 5.423e+05 +DEAL:FIRE::Convergence step 18843 value 0 +DEAL::FIRE::Solution vector: 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -- 2.39.5