]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new SolverFIRE inertial relaxation algorithm
authorvishalkenchan <vishal.boddu@fau.de>
Tue, 11 Jul 2017 14:51:00 +0000 (16:51 +0200)
committerVishal Boddu <vishal.boddu@fau.de>
Fri, 14 Jul 2017 22:19:11 +0000 (00:19 +0200)
doc/news/changes/minor/20170711VishalBoddu [new file with mode: 0644]
include/deal.II/lac/solver_fire.h [new file with mode: 0644]
tests/lac/solver.cc
tests/lac/solver.output
tests/lac/solver_fire_01.cc [new file with mode: 0644]
tests/lac/solver_fire_01.output [new file with mode: 0644]
tests/lac/solver_fire_rosenbrock_01.cc [new file with mode: 0644]
tests/lac/solver_fire_rosenbrock_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170711VishalBoddu b/doc/news/changes/minor/20170711VishalBoddu
new file mode 100644 (file)
index 0000000..c31f544
--- /dev/null
@@ -0,0 +1,4 @@
+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)
diff --git a/include/deal.II/lac/solver_fire.h b/include/deal.II/lac/solver_fire.h
new file mode 100644 (file)
index 0000000..76e5b02
--- /dev/null
@@ -0,0 +1,392 @@
+#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
index 582bc7a4a41f07cde34c927753d6595362e11736..ef88d74f6537c5e89621bfd83a34bfd3b81a8758 100644 (file)
@@ -26,6 +26,7 @@
 #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>
@@ -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();
         }
index d77ff53976cdbacd6bc53f9e2dcbba87731a05b8..ef1e643dcb263fa7af3ca9f450366eea66ca3381 100644 (file)
@@ -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 (file)
index 0000000..f2f1268
--- /dev/null
@@ -0,0 +1,79 @@
+
+#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);
+
+}
diff --git a/tests/lac/solver_fire_01.output b/tests/lac/solver_fire_01.output
new file mode 100644 (file)
index 0000000..ea7264d
--- /dev/null
@@ -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 (file)
index 0000000..51c8b3f
--- /dev/null
@@ -0,0 +1,117 @@
+
+#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);
+
+}
diff --git a/tests/lac/solver_fire_rosenbrock_01.output b/tests/lac/solver_fire_rosenbrock_01.output
new file mode 100644 (file)
index 0000000..ee4a049
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.