+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
#ifndef dealii__solver_fire_h
#define dealii__solver_fire_h
/**
* 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).
+ * non-linear) objective function $E(\mathbf x)$, $\mathbf x$ 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
* 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>
+ * - The global vector of unknown variables: $\mathbf x$.
+ * - Objective function: $E(\mathbf x)$.
+ * - Rate of change of unknowns: $\mathbf v$.
+ * - Gradient of the objective
+ * function w.r.t unknowns: $\mathbf g = \nabla E(\mathbf x)$.
+ * - Mass matrix: $\mathbf M$.
+ * - Initial guess of unknowns: $\mathbf x_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, <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.
+ * Given initial values for $\Delta t$, $\alpha = \alpha_0$, $\epsilon$,
+ * $\mathbf x = \mathbf x_0$ and $\mathbf v= \mathbf 0$ along with a given mass
+ * matrix $\mathbf M$, FIRE algorithm is as follows,
+ * 1. Calculate $\mathbf g = \nabla E(\mathbf x)$ and check for convergence
+ * ($\mathbf g \cdot \mathbf g < \epsilon^2 $).
+ * 2. Update $\mathbf x$ and $V$ using simple (forward) Euler integration step,
+ * <BR>
+ * $\mathbf x = \mathbf x + \Delta t \mathbf v$, <BR>
+ * $\mathbf v = \mathbf v + \Delta t \mathbf M^{-1} \cdot \mathbf g$.
+ * 3. Calculate $p = \mathbf g \cdot \mathbf v$.
+ * 4. Set $\mathbf v = (1-\alpha) \mathbf v
+ * + \alpha \frac{|\mathbf v|}{|\mathbf g|} \mathbf g$.
+ * 5. If $p<0$ and number of steps since $p$ was last negative is larger
+ * than certain value, then increase time step $\Delta t$ and decrease
+ * $\alpha$.
+ * 6. If $p>0$, then decrease the time step, freeze the system i.e.,
+ * $\mathbf v = \mathbf 0$ and reset $\alpha = \alpha_0$.
+ * 7. Return to 1.
*
- * Alse see
+ * Also 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.
public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
struct AdditionalData
{
+ /**
+ * Constructor. By default, set the initial time step for the (forward)
+ * Euler integration step to 0.1, the maximum time step to 1 and the
+ * maximum change allowed in any variable (per iteration) to 1.
+ */
explicit
AdditionalData (const double initial_timestep = 0.1,
const double maximum_timestep = 1,
const double maximum_linfty_norm = 1);
/**
- * Initial time step.
+ * Initial time step for the (forward) Euler integration step.
*/
const double initial_timestep;
/**
- * Maximum time step.
+ * Maximum time step for the (forward) Euler integration step.
*/
const double maximum_timestep;
/**
- * Maximum change allowed in any degree of freedom.
+ * Maximum change allowed in any variable of the objective function.
*/
const double maximum_linfty_norm;
*/
SolverFIRE (SolverControl &solver_control,
VectorMemory<VectorType> &vector_memory,
- const AdditionalData &data );
+ const AdditionalData &data = AdditionalData());
/**
* Constructor. Use an object of type GrowingVectorMemory as a default to
virtual ~SolverFIRE();
/**
- * Obtain a set of #p u (variables) that minimize an objective function
+ * Obtain a set of variables @p x 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).
+ * preconditioner @p inverse_mass_matrix and initial @p x values.
+ * The function @p compute returns the objective function's value and updates
+ * the objective function's gradient (with respect to the variables) when
+ * passed in as first argument based on the second argument-- the state of
+ * variables.
+ *
+ * @author Vishal Boddu, Denis Davydov, 2017
*/
template<typename PreconditionerType = DiagonalMatrix<VectorType>>
void solve
- (std::function<double(VectorType &, const VectorType &)> compute,
- VectorType &u,
- const PreconditionerType &inverse_masses);
+ (const std::function<double(VectorType &, const VectorType &)> &compute,
+ VectorType &x,
+ const PreconditionerType &inverse_mass_matrix);
/**
- * Solve for x that minimizes the quadratic function
- * $E(x) = \frac{1}{2} x^{T}Ax - x^{T}b$.
+ * Solve for x that minimizes $E(\mathbf x)$ for the <EM>special case</EM>
+ * when $E(\mathbf x)
+ * = \frac{1}{2} \mathbf x^{T} \mathbf A \mathbf x - \mathbf x^{T} \mathbf b$.
*/
template<typename MatrixType, typename PreconditionerType>
void solve (const MatrixType &A,
/**
* 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.
+ * @p x (variables), @p v (x's time derivative) and @p g (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;
+ virtual void print_vectors (const unsigned int,
+ const VectorType &x,
+ const VectorType &v,
+ const VectorType &g) const;
/**
- * Additional parameters.
+ * Additional data to the solver.
*/
const AdditionalData additional_data;
template<typename PreconditionerType>
void
SolverFIRE<VectorType>::solve
-(std::function<double(VectorType &, const VectorType &)> compute,
- VectorType &u,
- const PreconditionerType &inverse_masses)
+(const std::function<double(VectorType &, const VectorType &)> &compute,
+ VectorType &x,
+ const PreconditionerType &inverse_mass_matrix)
{
deallog.push("FIRE");
// Set velocities to zero but not gradients
// as we are going to compute them soon.
- v->reinit(u,false);
- g->reinit(u,true);
+ v->reinit(x,false);
+ g->reinit(x,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);
+ // Update gradients for the new x.
+ compute(gradients, x);
unsigned int iter = 0;
SolverControl::State conv = SolverControl::iterate;
- conv = this->iteration_status (iter, gradients * gradients, u);
+ conv = this->iteration_status (iter, gradients * gradients, x);
if (conv != SolverControl::iterate)
return;
{
++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
+ x.add (timestep, velocities); // x += dt * v
+ inverse_mass_matrix.vmult(gradients, gradients); // g = M^{-1} * g
+ velocities.add (-timestep, gradients); // v -= dt * h
- // Compute gradients for the new u.
- compute(gradients, u);
+ // Compute gradients for the new x.
+ compute(gradients, x);
const real_type gradient_norm_squared = gradients * gradients;
- conv = this->iteration_status(iter, gradient_norm_squared, u);
+ conv = this->iteration_status(iter, gradient_norm_squared, x);
if (conv != SolverControl::iterate)
break;
timestep = minimal_timestep;
}
- print_vectors(iter, u, velocities, gradients);
+ print_vectors(iter, x, velocities, gradients);
} // While we need to iterate.
{
std::function<double(VectorType &, const VectorType &)> compute_func =
- [&] (decltype(x) &G, decltype(b) &x) -> double
+ [&] (VectorType &g, const VectorType &x) -> double
{
// Residual of the quadratic form $ \frac{1}{2} xAx - xb $.
// G = b - Ax
- A.residual(G, x, b);
+ A.residual(g, x, b);
// Gradient G = Ax -b.
- G *= -1.;
+ g *= -1.;
- // The quadratic form $ \frac{1}{2} x^t A x - x^{t} b $.
+ // The quadratic form $\frac{1}{2} xAx - xb $.
return 0.5*A.matrix_norm_square(x) - x*b;
};
template <typename VectorType>
void
-SolverFIRE<VectorType>::
-print_vectors (const unsigned int,
- const VectorType &,
- const VectorType &,
- const VectorType &) const
+SolverFIRE<VectorType>::print_vectors (const unsigned int,
+ const VectorType &,
+ const VectorType &,
+ const VectorType &) const
{}