newton_raphson_linesearch.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/linesearch.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Functional_, typename Filter_>
    class NewtonRaphsonLinesearch : public Linesearch<Functional_, Filter_>
    {
      public:
        typedef Filter_ FilterType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename Functional_::DataType DataType;
        typedef Linesearch<Functional_, Filter_> BaseClass;
 
      public:
        explicit NewtonRaphsonLinesearch(Functional_& functional, Filter_& filter, bool keep_iterates = false) :
          BaseClass("NR-LS", functional, filter, keep_iterates)
          {
          }
 
        explicit NewtonRaphsonLinesearch(const String& section_name, const PropertyMap* section,
        Functional_& functional, Filter_& filter) :
          BaseClass("NR-LS", section_name, section, functional, filter)
          {
          }
 
        virtual ~NewtonRaphsonLinesearch()
        {
        }
 
        virtual String name() const override
        {
          return "Newton-Raphson-Linesearch";
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_dir) override
        {
          // clear solution vector
          vec_cor.format();
          this->_functional.prepare(vec_cor, this->_filter);
 
          // apply
          this->_status = _apply_intern(vec_cor, vec_dir);
          this->plot_summary();
          return this->_status;
        }
 
        virtual Status correct(VectorType& vec_sol, const VectorType& vec_dir) override
        {
          this->_functional.prepare(vec_sol, this->_filter);
          // apply
          this->_status = _apply_intern(vec_sol, vec_dir);
          this->plot_summary();
          return this->_status;
        }
 
      protected:
        virtual Status _apply_intern(VectorType& vec_sol, const VectorType& vec_dir)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
          // The step length wrt. to the NORMALISED search direction
          DataType alpha(0);
          // The functional value
          DataType fval(0);
          DataType df(0);
 
          // Perform initializations and checks
          Status st = this->_startup(alpha, fval, df, vec_sol, vec_dir);
          //  Because of the Newton-Raphson update x[k+1] = x[k] + alpha[k+1] d and alpha[k+1] = alpha[k] - < g[k], d > / < d, H[k] d >,
          // || d ||_2 does not matter and the relative update is invariant
 
          // Compute new _alpha <- _alpha - grad.dot(vec_dir) / vec_dir.dot(Hess*vec_dir)
          this->_functional.apply_hess(this->_vec_tmp, this->_vec_pn);
          this->_filter.filter_def(this->_vec_tmp);
 
          alpha = -this->_vec_grad.dot(this->_vec_pn)/this->_vec_pn.dot(this->_vec_tmp);
 
          // start iterating
          while(st == Status::progress)
          {
            IterationStats stat(*this);
 
            // Increase iteration count
            ++this->_num_iter;
 
            this->_functional.prepare(vec_sol, this->_filter);
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            // DO NOT call trim() here, as this simple linesearch accepts bad steps too, which would then be trimmed
            // and used.
            //this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
 
            df = this->_vec_pn.dot(this->_vec_grad);
 
            if(fval <= this->_fval_min)
            {
              this->_fval_min = fval;
              this->_alpha_min = alpha;
            }
 
            st = this->_check_convergence(fval, df, alpha);
 
            if(st != Status::progress)
            {
              break;
            }
 
            // Compute new _alpha <- _alpha - grad.dot(this->_vec_pn) / this->_vec_pn.dot(Hess*this->_vec_pn)
            this->_functional.apply_hess(this->_vec_tmp, this->_vec_pn);
            this->_filter.filter_def(this->_vec_tmp);
 
            alpha -= this->_vec_grad.dot(this->_vec_pn)/this->_vec_pn.dot(this->_vec_tmp);
 
            // Update solution
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, alpha);
 
          }
 
          // If we are not successful, we update the best step length and need to re-evaluate everything for that
          // step
          if(st != Status::success)
          {
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, this->_alpha_min);
 
            // Prepare and evaluate
            this->_functional.prepare(vec_sol, this->_filter);
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            // DO NOT call trim() here, as this simple linesearch accepts bad steps too, which would then be trimmed
            // and used.
            this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
          }
 
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::undefined, this->get_num_iter()));
          return st;
        }
 
    }; // class NewtonRaphsonLinesearch
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<NewtonRaphsonLinesearch<Functional_, Filter_>> new_newton_raphson_linesearch(
      Functional_& functional, Filter_& filter, bool keep_iterates = false)
      {
        return std::make_shared<NewtonRaphsonLinesearch<Functional_, Filter_>>(functional, filter, keep_iterates);
      }
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<NewtonRaphsonLinesearch<Functional_, Filter_>> new_newton_raphson_linesearch(
      const String& section_name, const PropertyMap* section,
      Functional_& functional, Filter_& filter)
      {
        return std::make_shared<NewtonRaphsonLinesearch<Functional_, Filter_>>(section_name, section, functional, filter);
      }
  } // namespace Solver
} // namespace FEAT