nloptls.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
#include <kernel/base_header.hpp>
#include <kernel/solver/base.hpp>
#include <kernel/solver/iterative.hpp>
#include <kernel/solver/linesearch.hpp>
#include <kernel/solver/nlopt_precond.hpp>
 
namespace FEAT
{
  namespace Solver
  {
 
    template<typename Functional_, typename Filter_>
    class NLOptLS : public PreconditionedIterativeSolver<typename Functional_::VectorTypeR>
    {
      public:
        typedef Functional_ FunctionalType;
        typedef Filter_ FilterType;
        typedef Solver::Linesearch<Functional_, Filter_> LinesearchType;
 
        typedef typename Functional_::GradientType GradientType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename Functional_::DataType DataType;
 
        typedef PreconditionedIterativeSolver<VectorType> BaseClass;
        typedef NLOptPrecond<typename Functional_::VectorTypeL, Filter_> PrecondType;
 
      protected:
        Functional_& _functional;
        Filter_& _filter;
 
        DataType _tol_fval;
        DataType _tol_step;
 
        DataType _fval_init;
        DataType _fval;
        DataType _fval_prev;
        DataType _steplength;
        Index _ls_iter_digits;
        Index _ls_its;
 
      public:
        explicit NLOptLS(const String& plot_name_, Functional_& functional, Filter_& filter,
          std::shared_ptr<PrecondType> precond = nullptr) :
          BaseClass(plot_name_, precond),
          _functional(functional),
          _filter(filter),
          _tol_fval(DataType(0)),
          _tol_step(Math::eps<DataType>()),
          _fval_init(-Math::huge<DataType>()),
          _fval(Math::huge<DataType>()),
          _fval_prev(-Math::huge<DataType>()),
          _steplength(1),
          _ls_iter_digits(2),
          _ls_its(~Index(0))
        {
          // set communicator by functional (same interface as matrix)
          this->_set_comm_by_matrix(functional);
        }
 
        explicit NLOptLS(const String& plot_name, const String& section_name, const PropertyMap* section,
          Functional_& functional, Filter_& filter, std::shared_ptr<PrecondType> precond = nullptr) :
          BaseClass(plot_name, section_name, section, precond),
          _functional(functional),
          _filter(filter),
          _tol_fval(DataType(0)),
          _tol_step(Math::eps<DataType>()),
          _fval_init(-Math::huge<DataType>()),
          _fval(Math::huge<DataType>()),
          _fval_prev(-Math::huge<DataType>()),
          _steplength(1),
          _ls_iter_digits(2),
          _ls_its(~Index(0))
        {
          // set communicator by functional (same interface as matrix)
          this->_set_comm_by_matrix(functional);
 
          auto tol_fval_p = section->get_entry("tol_fval");
          if (tol_fval_p.second)
          {
            set_tol_fval((DataType)std::stod(tol_fval_p.first));
          }
 
          auto tol_step_p = section->get_entry("tol_step");
          if (tol_step_p.second)
          {
            set_tol_step((DataType)std::stod(tol_step_p.first));
          }
        }
 
        virtual ~NLOptLS()
        {
        }
 
        virtual String get_summary() const override
        {
          String msg(this->get_plot_name()+ ": its: "+stringify(this->get_num_iter())
              +" ("+ stringify(this->get_status())+")"
              +", evals: "+stringify(_functional.get_num_func_evals())+" (func) "
              + stringify(_functional.get_num_grad_evals()) + " (grad) "
              + stringify(_functional.get_num_hess_evals()) + " (hess)"
              +" last step: "+stringify_fp_sci(_steplength)+"\n");
          msg +=this->get_plot_name()+": fval: "+stringify_fp_sci(_fval_init)
            + " -> "+stringify_fp_sci(_fval)
            + ", factor "+stringify_fp_sci(_fval/_fval_init)
            + ", last reduction "+stringify_fp_sci(_fval_prev - _fval)+"\n";
          msg += this->get_plot_name()  +": grad: "+stringify_fp_sci(this->_def_init)
            + " -> "+stringify_fp_sci(this->_def_cur)
            + ", factor " +stringify_fp_sci(this->_def_cur/this->_def_init);
 
          return msg;
        }
 
        virtual DataType get_tol_fval()
        {
          return _tol_fval;
        }
 
        virtual DataType get_tol_step()
        {
          return _tol_step;
        }
 
        virtual void set_tol_fval(DataType tol_fval)
        {
          _tol_fval = tol_fval;
        }
 
        virtual void set_tol_step(DataType tol_step)
        {
          _tol_step = tol_step;
        }
 
        virtual void set_ls_iter_digits(Index digits)
        {
          _ls_iter_digits = digits;
        }
 
      protected:
        virtual Status _set_initial_defect(const VectorType& vec_def, const VectorType& vec_sol) override
        {
          // store new defect
          this->_def_init = this->_def_cur = this->_calc_def_norm(vec_def, vec_sol);
          this->_num_iter = Index(0);
          this->_num_stag_iter = Index(0);
 
          _fval_init = _fval;
          _fval_prev = Math::huge<DataType>();
          _steplength = DataType(1);
          _ls_its = Index(0);
 
          Statistics::add_solver_expression(
            std::make_shared<ExpressionDefect>(this->name(), this->_def_init, this->get_num_iter()));
 
          if(this->_plot_iter())
          {
            String msg = this->_plot_name
              +  ": " + stringify(this->_num_iter).pad_front(this->_iter_digits)
              +  " (" + stringify(this->_ls_its).pad_front(_ls_iter_digits) + ")"
              + " : " + stringify_fp_sci(this->_def_cur)
              + " / " + stringify_fp_sci(this->_def_cur / this->_def_init)
              + " : " + stringify_fp_sci(this->_fval);
 
            this->_print_line(msg);
          }
 
          // Ensure that the initial fval and defect are neither NaN nor infinity
          if( !(Math::isfinite(this->_def_init) && Math::isfinite(this->_fval_init)) )
          {
            return Status::aborted;
          }
 
          // Check if the initial defect is already low enough
          if(this->_def_init <= Math::sqr(Math::eps<DataType>()))
          {
            return Status::success;
          }
 
          // continue iterating
          return Status::progress;
        }
 
        virtual Status _set_new_defect(const VectorType& vec_r, const VectorType& vec_sol) override
        {
          // increase iteration count
          ++this->_num_iter;
 
          // first, let's see if we have to compute the defect at all
          bool calc_def = false;
          calc_def = calc_def || (this->_min_iter < this->_max_iter);
          calc_def = calc_def || this->_plot_iter();
          calc_def = calc_def || (this->_min_stag_iter > Index(0));
 
          // compute new defect
          if(calc_def)
          {
            this->_def_cur = this->_calc_def_norm(vec_r, vec_sol);
            Statistics::add_solver_expression(
              std::make_shared<ExpressionDefect>(this->name(), this->_def_cur, this->get_num_iter()));
          }
 
          // plot?
          if(this->_plot_iter())
          {
            String msg = this->_plot_name
              +  ": " + stringify(this->_num_iter).pad_front(this->_iter_digits)
              +  " (" + stringify(this->_ls_its).pad_front(_ls_iter_digits) + ")"
              + " : " + stringify_fp_sci(this->_def_cur)
              + " / " + stringify_fp_sci(this->_def_cur / this->_def_init)
              + " : " + stringify_fp_sci(this->_fval)
              + " : " + stringify_fp_sci(this->_steplength);
 
            // print message line via comm (if available)
            this->_print_line(msg);
          }
 
          // ensure that the defect is neither NaN nor infinity
          if(!Math::isfinite(this->_def_cur))
          {
            return Status::aborted;
          }
 
          // is diverged?
          if(this->is_diverged())
          {
            return Status::diverged;
          }
 
          // minimum number of iterations performed?
          if(this->_num_iter < this->_min_iter)
          {
            return Status::progress;
          }
 
          // maximum number of iterations performed?
          if(this->_num_iter >= this->_max_iter)
          {
            return Status::max_iter;
          }
 
          // Check for convergence of the gradient norm (relative AND absolute criterion)
          if(this->is_converged())
          {
            return Status::success;
          }
 
          // Check for convergence wrt. the function value improvement if _tol_fval says so
          if(_tol_fval > DataType(0))
          {
            // This is the factor for the relative function value
            DataType scale(Math::max(_fval, _fval_prev));
            // Make sure it is at least 1
            scale = Math::max(scale, DataType(1));
            // Check for success
            if(Math::abs(_fval_prev - _fval) <= _tol_fval*scale)
            {
              return Status::success;
            }
          }
 
          if(_steplength <= _tol_step)
          {
            return Status::success;
          }
 
          // If there were too many stagnated iterations, the solver is stagnated
          if(this->_min_stag_iter > 0 && this->_num_stag_iter > this->_min_stag_iter)
          {
            return Status::stagnated;
          }
 
          // continue iterating
          return Status::progress;
        }
 
    };
  } // namespace Solver
} //namespace FEAT