hessian_precond.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/nlopt_precond.hpp>
#include <kernel/solver/test_aux/analytic_function_operator.hpp>
#include <kernel/util/dist.hpp>
 
namespace FEAT
{
  namespace Solver
  {
#ifndef DOXYGEN
    template<typename Operator_, typename Filter_>
    class HessianPrecond;
#else
    template<typename Operator_, typename Filter_>
    class HessianPrecond : public NLOptPrecond<typename Operator_::VectorType, Filter_>
    {};
#endif
 
 
    template<typename DT_, typename IT_, typename Function_, typename Filter_>
    class HessianPrecond<AnalyticFunctionOperator<DT_, IT_, Function_>, Filter_>
    : public NLOptPrecond<typename AnalyticFunctionOperator<DT_, IT_, Function_>::VectorTypeL, Filter_>
    {
      public:
        typedef AnalyticFunctionOperator<DT_, IT_, Function_> OperatorType;
        typedef Filter_ FilterType;
 
        typedef typename OperatorType::VectorTypeL VectorType;
        typedef NLOptPrecond<VectorType, Filter_> BaseClass;
 
        typedef typename OperatorType::DataType DataType;
        typedef typename OperatorType::HessianType HessianType;
 
      protected:
        OperatorType& _op;
        const FilterType& _filter;
        HessianType _hessian;
        HessianType _inv_hessian;
 
      public:
        explicit HessianPrecond(OperatorType& op_, const FilterType& filter_) :
          _op(op_),
          _filter(filter_),
          _hessian(DT_(0)),
          _inv_hessian(DT_(0))
        {
        }
 
        virtual ~HessianPrecond()
        {
        }
 
        virtual String name() const override
        {
          return "HessianPrecond";
        }
 
        virtual void prepare(const VectorType& DOXY(vec_state), FilterType& DOXY(filter)) override
        {
          _op.compute_hess(_hessian);
          _inv_hessian.set_inverse(_hessian);
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
          Statistics::add_solver_expression(std::make_shared<ExpressionCallPrecond>(this->name(), _op.name()));
 
          vec_cor(0, _inv_hessian*vec_def(0));
          this->_filter.filter_cor(vec_cor);
 
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::success, 1));
 
          return Status::success;
        }
    };
 
    template<typename Operator_, typename Filter_>
    inline std::shared_ptr<HessianPrecond<Operator_, Filter_>> new_hessian_precond(
      Operator_& op_, const Filter_& filter_)
      {
        return std::make_shared<HessianPrecond<Operator_, Filter_>>(op_, filter_);
      }
 
#ifndef DOXYGEN
    template<typename Operator_, typename Filter_>
    class ApproximateHessianPrecond;
#else
    template<typename Operator_, typename Filter_>
    class ApproximateHessianPrecond : public NLOptPrecond<typename OperatorType_::VectorType, Filter>
    {};
#endif
 
 
    template<typename DT_, typename IT_, typename Function_, typename Filter_>
    class ApproximateHessianPrecond<AnalyticFunctionOperator<DT_, IT_, Function_>, Filter_>
    : public NLOptPrecond<typename AnalyticFunctionOperator<DT_, IT_, Function_>::VectorTypeL, Filter_>
    {
      public:
        typedef AnalyticFunctionOperator<DT_, IT_, Function_> OperatorType;
        typedef Filter_ FilterType;
 
        typedef typename OperatorType::VectorTypeL VectorType;
        typedef NLOptPrecond<VectorType, Filter_> BaseClass;
 
        typedef typename OperatorType::DataType DataType;
        typedef typename OperatorType::HessianType HessianType;
 
      protected:
        OperatorType& _op;
        const FilterType& _filter;
        HessianType _hessian;
        HessianType _inv_hessian;
 
      public:
        explicit ApproximateHessianPrecond(OperatorType& op_, const FilterType& filter_) :
          _op(op_),
          _filter(filter_),
          _hessian(DT_(0)),
          _inv_hessian(DT_(0))
          {
          }
 
        virtual ~ApproximateHessianPrecond()
        {
        }
 
        virtual String name() const override
        {
          return "ApproximateHessianPrecond";
        }
 
        //virtual void init_symbolic() override
        //{
        //  _inv_diag = _op.create_vector_r();
        //}
 
        //virtual void done_symbolic() override
        //{
        //  _inv_diag.clear();
        //}
 
        //virtual void init_numeric() override
        //{
        //  // extract matrix diagonal
        //  _op.extract_diag(_inv_diag);
 
        //  // invert diagonal elements
        //  _inv_diag.component_invert(_inv_diag, _omega);
        //}
 
 
        virtual void prepare(const VectorType& DOXY(vec_state), FilterType& DOXY(filter)) override
        {
          _op.compute_approx_hess(_hessian);
          _inv_hessian.set_inverse(_hessian);
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
          Statistics::add_solver_expression(std::make_shared<ExpressionCallPrecond>(this->name(), _op.name()));
 
          vec_cor(0, _inv_hessian*vec_def(0));
          this->_filter.filter_cor(vec_cor);
 
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::success, 1));
 
          return Status::success;
        }
    };
 
    template<typename Operator_, typename Filter_>
    inline std::shared_ptr<ApproximateHessianPrecond<Operator_, Filter_>> new_approximate_hessian_precond(
      Operator_& op_, const Filter_& filter_)
      {
        return std::make_shared<ApproximateHessianPrecond<Operator_, Filter_>>(op_, filter_);
      }
 
  } // namespace Solver
} // namespace FEAT