PlayaNonlinearOperatorBase.hpp

/* @HEADER@ */
//
 /* @HEADER@ */

#ifndef PLAYA_NONLINEAROPERATORBASE_HPP
#define PLAYA_NONLINEAROPERATORBASE_HPP

#include "PlayaDefs.hpp"
#include "PlayaHandleable.hpp"
#include "PlayaOut.hpp"
#include "PlayaObjectWithVerbosity.hpp"
#include "PlayaVectorDecl.hpp"
#include "PlayaLinearOperatorDecl.hpp"
#include "PlayaLinearCombinationDecl.hpp"

namespace Playa
{
using namespace Teuchos;

template <class Scalar>
class NonlinearOperatorBase 
  : public Handleable<NonlinearOperatorBase<Scalar> >,
    public ObjectWithVerbosity
{
public:
  NonlinearOperatorBase() 
    : domain_(), range_(), 
      jacobianIsValid_(false),
      residualIsValid_(false),
      currentEvalPt_(),
      currentFunctionValue_(),
      currentJ_()
    {;}

  NonlinearOperatorBase(const VectorSpace<Scalar>& domain,
    const VectorSpace<Scalar>& range) 
    : domain_(domain.ptr()), range_(range.ptr()), 
      jacobianIsValid_(false),
      residualIsValid_(false),
      currentEvalPt_(),
      currentFunctionValue_(),
      currentJ_()
    {;}
                            
  const RCP<const VectorSpaceBase<Scalar> >& domain() const 
    {return domain_;}

  const RCP<const VectorSpaceBase<Scalar> >& range() const 
    {return range_;}

  void setEvalPt(const Vector<Scalar>& x) const 
    {
      if (this->verb() >= 1)
      {
        Out::os() << "NonlinearOperatorBase Setting new eval pt";
        if (this->verb() > 3)
        {
          Out::os() << " to " << std::endl ;
          x.print(Out::os());
        }
        Out::os() << std::endl;
      }
      jacobianIsValid_ = false;
      residualIsValid_ = false;

      currentEvalPt_.acceptCopyOf(x);
    }

  const Vector<double>& currentEvalPt() const {return currentEvalPt_;}

  LinearOperator<double> getJacobian() const 
    {
      if (this->verb() > 1)
      {
        Out::os() << "NonlinearOperatorBase getting Jacobian" << std::endl;
      }
      if (!jacobianIsValid_)
      {
        if (this->verb() > 3)
        {
          Out::os() << "...computing new J and F" << std::endl;
        }
        currentJ_ 
          = computeJacobianAndFunction(currentFunctionValue_);
        jacobianIsValid_ = true;
        residualIsValid_ = true;
      }
      else
      {
        if (this->verb() > 1)
        {
          Out::os() << "...reusing valid J" << std::endl;
        }
      }
      if (this->verb() > 3)
      {
        Out::os() << "J is " << std::endl;
        currentJ_.print(Out::os());
        Out::os() << std::endl;
      }
      return currentJ_;
    }

      

  Vector<double> getFunctionValue() const 
    {
      if (this->verb() > 1)
      {
        Out::os() << "NonlinearOperatorBase getting function value" << std::endl;
      }
      if (!residualIsValid_)
      {
        if (this->verb() > 1)
        {
          Out::os() << "...computing new F" << std::endl;
        }
        currentFunctionValue_ = computeFunctionValue();
        residualIsValid_ = true;
      }
      else
      {
        if (this->verb() > 1)
        {
          Out::os() << "...reusing valid F" << std::endl;
        }
      }

      if (this->verb() > 3)
      {
        Out::os() << "F is " << std::endl;
        currentFunctionValue_.print(Out::os());
        Out::os() << std::endl;
      }
      return currentFunctionValue_;
    }


  virtual Vector<double> getInitialGuess() const = 0 ;


protected:

  virtual LinearOperator<Scalar> computeJacobianAndFunction(Vector<double>& functionValue) const = 0 ;

  virtual Vector<Scalar> computeFunctionValue() const 
    {
      computeJacobianAndFunction(currentFunctionValue_);
      return currentFunctionValue_;
    }

      
  void setDomainAndRange(const VectorSpace<Scalar>& domain,
    const VectorSpace<Scalar>& range)
    {
      domain_ = domain.ptr();
      range_ = range.ptr();
    }


private:
  RCP<const VectorSpaceBase<Scalar> > domain_;

  RCP<const VectorSpaceBase<Scalar> > range_;

  mutable bool jacobianIsValid_;

  mutable bool residualIsValid_;

  mutable Vector<double> currentEvalPt_;

  mutable Vector<double> currentFunctionValue_;

  mutable LinearOperator<double> currentJ_;
};



 
}


#endif