reducednewtonalgorithmwithinverse.h

Go to the documentation of this file.

#ifndef REDUCEDNEWTON__ALGORITHM_INVERSE_H_
#define REDUCEDNEWTON__ALGORITHM_INVERSE_H_

#include "reducednewtonalgorithm.h"
#include "parameterreader.h"

#include <iostream>
#include <assert.h>
#include <iomanip>
namespace DOpE
{
  template <typename PROBLEM, typename VECTOR>
    class ReducedNewtonAlgorithmWithInverse : public ReducedNewtonAlgorithm<PROBLEM, VECTOR>
  {
  public:
    ReducedNewtonAlgorithmWithInverse(PROBLEM* OP, 
                                      ReducedProblemInterface<PROBLEM, VECTOR>* S,
                                      ParameterReader &param_reader,
                                      DOpEExceptionHandler<VECTOR>* Except=NULL,
                                      DOpEOutputHandler<VECTOR>* Output=NULL,
                                      int base_priority=0);
    ~ReducedNewtonAlgorithmWithInverse();
    
     static void declare_params(ParameterReader &param_reader);


    virtual int Solve(ControlVector<VECTOR>& q,double global_tol=-1.);

  protected:
    int SolveReducedLinearSystem(const ControlVector<VECTOR>& q, 
                                 const ControlVector<VECTOR>& gradient,
                                 const ControlVector<VECTOR>& gradient_transposed, 
                                 ControlVector<VECTOR>& dq);
    
    double Residual(const ControlVector<VECTOR>& gradient,
                    const ControlVector<VECTOR>& /*gradient_transposed*/)
    {return  gradient*gradient;}
  private:
    unsigned int line_maxiter_;
    double       linesearch_rho_, linesearch_c_;
    
  };

  /***************************************************************************************/
  /****************************************IMPLEMENTATION*********************************/
  /***************************************************************************************/
  using namespace dealii;

  /******************************************************/

template <typename PROBLEM, typename VECTOR>
void ReducedNewtonAlgorithmWithInverse<PROBLEM, VECTOR>::declare_params(ParameterReader &param_reader)
  {
    ReducedNewtonAlgorithm<PROBLEM, VECTOR>::declare_params(param_reader);
  }
/******************************************************/

template <typename PROBLEM, typename VECTOR>
ReducedNewtonAlgorithmWithInverse<PROBLEM, VECTOR>
  ::ReducedNewtonAlgorithmWithInverse(PROBLEM* OP, 
                                      ReducedProblemInterface<PROBLEM, VECTOR>* S,
                                      ParameterReader &param_reader,
                                      DOpEExceptionHandler<VECTOR>* Except,
                                      DOpEOutputHandler<VECTOR>* Output,
                                      int base_priority) 
  : ReducedNewtonAlgorithm<PROBLEM, VECTOR>(OP,S,param_reader,Except,Output,base_priority)
  {
    param_reader.SetSubsection("reducednewtonalgorithm parameters");
    line_maxiter_         = param_reader.get_integer ("line_maxiter");
    linesearch_rho_       = param_reader.get_double ("linesearch_rho");
    linesearch_c_         = param_reader.get_double ("linesearch_c");
  }

/******************************************************/

template <typename PROBLEM, typename VECTOR>
ReducedNewtonAlgorithmWithInverse<PROBLEM, VECTOR>::~ReducedNewtonAlgorithmWithInverse()
  {
    
  }
/******************************************************/

template <typename PROBLEM, typename VECTOR>
int ReducedNewtonAlgorithmWithInverse<PROBLEM, VECTOR>::Solve(ControlVector<VECTOR>& q,double global_tol)
{
  return ReducedNewtonAlgorithm<PROBLEM, VECTOR>::Solve(q,global_tol);
}

/******************************************************/

template <typename PROBLEM, typename VECTOR>
int ReducedNewtonAlgorithmWithInverse<PROBLEM, VECTOR>
  ::SolveReducedLinearSystem(const ControlVector<VECTOR>& q, 
                             const ControlVector<VECTOR>& gradient,
                             const ControlVector<VECTOR>& /*gradient_transposed*/, 
                             ControlVector<VECTOR>& dq)
{
  int iter = 0;
  this->GetReducedProblem()->ComputeReducedHessianInverseVector(q,gradient,dq);
//  {
//    //Check...
//    ControlVector<VECTOR> Hd(q), Hd_transposed(q);
//    this->GetReducedProblem()->ComputeReducedHessianVector(q,dq,Hd,Hd_transposed);
//    Hd.add(-1.,gradient);
//    std::cout<<" Linear residual "<<Hd*Hd<<std::endl;
//  }
  return iter;
}


/******************************************************/
}
#endif