--- /dev/null
+New: SolverFIRE implements FIRE (Fast Inertial Relaxation Engine) for solving
+the problem of minimization of a given objective function.
+<br>
+(Vishal Boddu, 2017/07/11)
--- /dev/null
+#ifndef dealii__solver_fire_h
+#define dealii__solver_fire_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/solver.h>
+
+#include <functional>
+
+
+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<double>.
+ *
+ * FIRE is a damped dynamics method described in
+ * <a href="https://doi.org/10.1103/PhysRevLett.97.170201">Structural
+ * Relaxation Made Simple</a> 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$. <BR>
+ * Objective function: $E(U)$. <BR>
+ * Rate of change of unknowns: $V$. <BR>
+ * Gradient of the objective
+ * function w.r.t unknowns: $G$. <BR>
+ * Mass matrix: $M$. <BR>
+ * Initial guess of unknowns: $U_0$. <BR>
+ * Time step: $\Delta t$. <BR>
+ *
+ * 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, <BR>
+ * $U = U + \Delta t V$, <BR>
+ * $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
+ * <a href="http://onlinelibrary.wiley.com/doi/10.1002/pamm.201110246/full">
+ * Energy-Minimization in Atomic-to-Continuum Scale-Bridging Methods </a> by
+ * Eidel et al. 2011.
+ *
+ * @author Vishal Boddu 2017
+ */
+template<typename VectorType = Vector<double> >
+class SolverFIRE : public Solver<VectorType>
+{
+
+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<VectorType> &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<typename PreconditionerType = DiagonalMatrix<VectorType>>
+ void solve
+ (std::function<double(VectorType &, const VectorType &)> 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<typename MatrixType, typename PreconditionerType>
+ 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<typename VectorType>
+SolverFIRE<VectorType>::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<typename VectorType>
+SolverFIRE<VectorType>::
+SolverFIRE (SolverControl &solver_control,
+ VectorMemory<VectorType> &vector_memory,
+ const AdditionalData &data )
+ :
+ Solver<VectorType>(solver_control, vector_memory),
+ additional_data(data)
+{}
+
+
+
+template<typename VectorType>
+SolverFIRE<VectorType>::
+SolverFIRE (SolverControl &solver_control,
+ const AdditionalData &data )
+ :
+ Solver<VectorType>(solver_control),
+ additional_data(data)
+{}
+
+
+
+template<typename VectorType>
+SolverFIRE<VectorType>::~SolverFIRE()
+{}
+
+
+
+template<typename VectorType>
+template<typename PreconditionerType>
+void
+SolverFIRE<VectorType>::solve
+(std::function<double(VectorType &, const VectorType &)> 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<VectorType>::Pointer v(this->memory);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+template<typename MatrixType, typename PreconditionerType>
+void SolverFIRE<VectorType>::solve (const MatrixType &A,
+ VectorType &x,
+ const VectorType &b,
+ const PreconditionerType &precondition)
+{
+
+ std::function<double(VectorType &, const VectorType &)> 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 <typename VectorType>
+void
+SolverFIRE<VectorType>::
+print_vectors (const unsigned int,
+ const VectorType &,
+ const VectorType &,
+ const VectorType &) const
+{}
+
+
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_fire.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/solver_bicgstab.h>
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)
{
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();
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();
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();
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();
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();
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();
}
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
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
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
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
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
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
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)
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)
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
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
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
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
--- /dev/null
+
+#include "../tests.h"
+#include <fstream>
+#include <iomanip>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_fire.h>
+
+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>;
+
+
+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<vector_t> inv_mass;
+ inv_mass.reinit(X);
+
+ // Set initial iterate.
+ X(0) = x;
+ X(1) = y;
+
+ auto additional_data =
+ SolverFIRE<vector_t>::AdditionalData(0.1, 1., 1);
+
+ SolverControl solver_control (200, tol);
+
+ SolverFIRE<vector_t> 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);
+
+}
--- /dev/null
+
+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
--- /dev/null
+
+#include "../tests.h"
+#include <fstream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_fire.h>
+
+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>;
+
+
+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<vector_t> 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<vector_t>::AdditionalData(0.1, 1, 1);
+
+ SolverControl solver_control (1e5, tol);
+
+ SolverFIRE<vector_t> 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);
+
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff