Thyra_ModelEvaluatorBase_decl.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
// 
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//                 Copyright (2004) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef THYRA_MODEL_EVALUATOR_BASE_DECL_HPP
#define THYRA_MODEL_EVALUATOR_BASE_DECL_HPP


#include "Thyra_LinearOpWithSolveBase.hpp"
#include "Teuchos_Describable.hpp"
#include "Teuchos_Assert.hpp"

#ifdef HAVE_THYRA_ME_POLYNOMIAL
#  include "Teuchos_Polynomial.hpp"
#endif


namespace Thyra {


class ModelEvaluatorBase
  : virtual public Teuchos::Describable,
    virtual public Teuchos::VerboseObject<ModelEvaluatorBase>
{
public:


  enum EInArgsMembers {
    IN_ARG_x_dot 
    ,IN_ARG_x 
    ,IN_ARG_x_dot_poly 
    ,IN_ARG_x_poly 
    ,IN_ARG_t 
    ,IN_ARG_alpha 
    ,IN_ARG_beta 
  };
  static const int NUM_E_IN_ARGS_MEMBERS=7;

  template<class Scalar>
  class InArgs : public Teuchos::Describable {
  public:
    typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
    InArgs();
    int Np() const;
    bool supports(EInArgsMembers arg) const;
    void set_x_dot( const RCP<const VectorBase<Scalar> > &x_dot );
    RCP<const VectorBase<Scalar> > get_x_dot() const;
    void set_x( const RCP<const VectorBase<Scalar> > &x );
    RCP<const VectorBase<Scalar> > get_x() const;
#ifdef HAVE_THYRA_ME_POLYNOMIAL

    void set_x_poly( 
      const RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > &x_poly );
    RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > get_x_poly() const;
    void set_x_dot_poly(
      const RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > &x_dot_poly );
    RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > get_x_dot_poly() const;
#endif // HAVE_THYRA_ME_POLYNOMIAL

    void set_p( int l, const RCP<const VectorBase<Scalar> > &p_l );
    RCP<const VectorBase<Scalar> > get_p(int l) const;
    void set_t( ScalarMag t );
    ScalarMag get_t() const;
    void set_alpha( Scalar alpha );
    Scalar get_alpha() const;
    void set_beta( Scalar beta );
    Scalar get_beta() const;
    void setArgs(
      const InArgs<Scalar>& inArgs, bool ignoreUnsupported = false,
      bool cloneObjects = false
      );
    void assertSameSupport( const InArgs<Scalar> &inArgs ) const;
    std::string modelEvalDescription() const;
    std::string description() const;
    void describe(
      Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel
      ) const;
  protected:
    void _setModelEvalDescription( const std::string &modelEvalDescription );
    void _set_Np(int Np);
    void _setSupports( EInArgsMembers arg, bool supports );
    void _setSupports( const InArgs<Scalar>& inputInArgs, const int Np );
    void _setUnsupportsAndRelated( EInArgsMembers arg );
  private:
    // types
    typedef Teuchos::Array<RCP<const VectorBase<Scalar> > > p_t;
    // data
    std::string modelEvalDescription_;
    RCP<const VectorBase<Scalar> > x_dot_;
    RCP<const VectorBase<Scalar> > x_;
#ifdef HAVE_THYRA_ME_POLYNOMIAL
    RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > x_dot_poly_;
    RCP<const Teuchos::Polynomial< VectorBase<Scalar> > > x_poly_;
#endif // HAVE_THYRA_ME_POLYNOMIAL
    p_t p_;
    ScalarMag t_;
    Scalar alpha_;
    Scalar beta_;
    bool supports_[NUM_E_IN_ARGS_MEMBERS];
    // functions
    void assert_supports(EInArgsMembers arg) const;
    void assert_l(int l) const;
  };

  enum EDerivativeMultiVectorOrientation {
    DERIV_MV_JACOBIAN_FORM, 
    DERIV_MV_GRADIENT_FORM, 
    DERIV_MV_BY_COL = DERIV_MV_JACOBIAN_FORM, 
    DERIV_TRANS_MV_BY_ROW = DERIV_MV_GRADIENT_FORM 
  };

  enum EDerivativeLinearOp {
    DERIV_LINEAR_OP 
  };

  class DerivativeSupport {
  public:
    DerivativeSupport()
      :supportsLinearOp_(false), supportsMVByCol_(false), supportsTransMVByRow_(false)
      {}
    DerivativeSupport( EDerivativeLinearOp )
      :supportsLinearOp_(true), supportsMVByCol_(false), supportsTransMVByRow_(false)
      {}
    DerivativeSupport( EDerivativeMultiVectorOrientation mvOrientation )
      :supportsLinearOp_(false), supportsMVByCol_(mvOrientation==DERIV_MV_BY_COL)
      ,supportsTransMVByRow_(mvOrientation==DERIV_TRANS_MV_BY_ROW)
      {}
    DerivativeSupport& plus(EDerivativeLinearOp)
      { supportsLinearOp_ = true; return *this; }
    DerivativeSupport& plus(EDerivativeMultiVectorOrientation mvOrientation)
      {
        switch(mvOrientation) {
          case DERIV_MV_BY_COL: supportsMVByCol_ = true; break;
          case DERIV_TRANS_MV_BY_ROW: supportsTransMVByRow_ = true; break;
          default: TEST_FOR_EXCEPT(true);
        }
        return *this;
      }
    bool none() const
      { return ( !supportsLinearOp_ && !supportsMVByCol_ && !supportsTransMVByRow_ ); }
    bool supports(EDerivativeLinearOp) const
      { return supportsLinearOp_; }
    bool supports(EDerivativeMultiVectorOrientation mvOrientation) const
      {
        switch(mvOrientation) {
          case DERIV_MV_BY_COL: return supportsMVByCol_;
          case DERIV_TRANS_MV_BY_ROW: return supportsTransMVByRow_;
          default: TEST_FOR_EXCEPT(true);
        }
        return false; // Will never be called!
      }
    bool isSameSupport(const DerivativeSupport &derivSupport) const
      {
        return (
          supportsLinearOp_ == derivSupport.supportsLinearOp_
          && supportsMVByCol_ == derivSupport.supportsMVByCol_
          && supportsTransMVByRow_ == derivSupport.supportsTransMVByRow_
          );
      } 
    std::string description() const;
  private:
    bool supportsLinearOp_;
    bool supportsMVByCol_;
    bool supportsTransMVByRow_;
  public:
  };
  
  enum EDerivativeLinearity {
    DERIV_LINEARITY_UNKNOWN      
    ,DERIV_LINEARITY_CONST       
    ,DERIV_LINEARITY_NONCONST    
  };

  enum ERankStatus {
    DERIV_RANK_UNKNOWN       
    ,DERIV_RANK_FULL         
    ,DERIV_RANK_DEFICIENT    
  };

  struct DerivativeProperties {
    EDerivativeLinearity     linearity;
    ERankStatus              rank;
    bool                     supportsAdjoint;
    DerivativeProperties()
      :linearity(DERIV_LINEARITY_UNKNOWN),
       rank(DERIV_RANK_UNKNOWN),supportsAdjoint(false)
      {}
    DerivativeProperties(
      EDerivativeLinearity in_linearity, ERankStatus in_rank,
      bool in_supportsAdjoint
      )
      :linearity(in_linearity),rank(in_rank),
       supportsAdjoint(in_supportsAdjoint)
      {}
  };

  template<class Scalar>
  class DerivativeMultiVector {
  public:
    DerivativeMultiVector()
      :orientation_(DERIV_MV_BY_COL)
      {}
    DerivativeMultiVector(
      const RCP<MultiVectorBase<Scalar> > &mv
      ,const EDerivativeMultiVectorOrientation orientation = DERIV_MV_BY_COL
      ) : mv_(mv.assert_not_null()), orientation_(orientation) {}
    void changeOrientation( const EDerivativeMultiVectorOrientation orientation )
      { orientation_ = orientation; };
    const DerivativeMultiVector<Scalar>& assert_not_null() const
      { mv_.assert_not_null(); return *this; }
    RCP<MultiVectorBase<Scalar> > getMultiVector() const
      { return mv_; }
    EDerivativeMultiVectorOrientation getOrientation() const
      { return orientation_; }
    std::string description() const;
    void describe( 
      Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel
      ) const;
  private:
    RCP<MultiVectorBase<Scalar> > mv_;
    EDerivativeMultiVectorOrientation orientation_;
  };

  template<class Scalar>
  class Derivative {
  public:
    Derivative() {}
    Derivative( const RCP<LinearOpBase<Scalar> > &lo )
      : lo_(lo.assert_not_null()) {}
    Derivative(
      const RCP<MultiVectorBase<Scalar> > &mv,
      const EDerivativeMultiVectorOrientation orientation = DERIV_MV_BY_COL
      ) : dmv_(mv,orientation) {}
    Derivative( const DerivativeMultiVector<Scalar> &dmv )
      : dmv_(dmv) {}
    bool isEmpty() const
      { return ( lo_.get()==NULL && dmv_.getMultiVector().get()==NULL ); }
    const Derivative<Scalar>& assert_not_null() const
      { dmv_.assert_not_null(); lo_.assert_not_null(); return *this; }
    RCP<LinearOpBase<Scalar> > getLinearOp() const
      { return lo_; }
    RCP<MultiVectorBase<Scalar> > getMultiVector() const
      { return dmv_.getMultiVector(); }
    EDerivativeMultiVectorOrientation getMultiVectorOrientation() const
      { return dmv_.getOrientation(); }
    DerivativeMultiVector<Scalar> getDerivativeMultiVector() const
      { return dmv_; }
    bool isSupportedBy( const DerivativeSupport &derivSupport ) const
      {
        // If there is not derivative support then we will return false!
        if (derivSupport.none())
          return false;
        if (!is_null(getMultiVector())) {
          return derivSupport.supports(getMultiVectorOrientation());
        }
        else if(!is_null(getLinearOp())) {
          return derivSupport.supports(DERIV_LINEAR_OP);
        }
        // If nothing is set then of course we support that!
        return true;
      }
    std::string description() const;
    void describe( 
      Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel
      ) const;
  private:
    RCP<LinearOpBase<Scalar> > lo_;
    DerivativeMultiVector<Scalar> dmv_;
  };

  enum EOutArgsMembers {
    OUT_ARG_f       
    ,OUT_ARG_W      
    ,OUT_ARG_W_op   
    ,OUT_ARG_f_poly 
  };
  static const int NUM_E_OUT_ARGS_MEMBERS=4;

  enum EOutArgsDfDp {
    OUT_ARG_DfDp   
  };

  enum EOutArgsDgDx_dot {
    OUT_ARG_DgDx_dot   
  };

  enum EOutArgsDgDx {
    OUT_ARG_DgDx   
  };

  enum EOutArgsDgDp {
    OUT_ARG_DgDp   
  };
  
  template<class Scalar>
  class OutArgs : public Teuchos::Describable {
  public:
    OutArgs();
    int Np() const;
    int Ng() const;
    bool supports(EOutArgsMembers arg) const;
    const DerivativeSupport& supports(EOutArgsDfDp arg, int l) const;
    const DerivativeSupport& supports(EOutArgsDgDx_dot arg, int j) const;
    const DerivativeSupport& supports(EOutArgsDgDx arg, int j) const;
    const DerivativeSupport& supports(EOutArgsDgDp arg, int j, int l) const;
    void set_f( const RCP<VectorBase<Scalar> > &f );
    RCP<VectorBase<Scalar> > get_f() const;
    void set_g( int j, const RCP<VectorBase<Scalar> > &g_j );
    RCP<VectorBase<Scalar> > get_g(int j) const;
    void set_W( const RCP<LinearOpWithSolveBase<Scalar> > &W );
    RCP<LinearOpWithSolveBase<Scalar> > get_W() const;
    void set_W_op( const RCP<LinearOpBase<Scalar> > &W_op );
    RCP<LinearOpBase<Scalar> > get_W_op() const;
    DerivativeProperties get_W_properties() const;
    void set_DfDp(int l,  const Derivative<Scalar> &DfDp_l);
    Derivative<Scalar> get_DfDp(int l) const;
    DerivativeProperties get_DfDp_properties(int l) const;
    void set_DgDx_dot(int j, const Derivative<Scalar> &DgDx_dot_j);
    Derivative<Scalar> get_DgDx_dot(int j) const;
    DerivativeProperties get_DgDx_dot_properties(int j) const;
    void set_DgDx(int j, const Derivative<Scalar> &DgDx_j);
    Derivative<Scalar> get_DgDx(int j) const;
    DerivativeProperties get_DgDx_properties(int j) const;
    void set_DgDp( int j, int l, const Derivative<Scalar> &DgDp_j_l );
    Derivative<Scalar> get_DgDp(int j, int l) const;
    DerivativeProperties get_DgDp_properties(int j, int l) const;
#ifdef HAVE_THYRA_ME_POLYNOMIAL

    void set_f_poly( const RCP<Teuchos::Polynomial< VectorBase<Scalar> > > &f_poly );
    RCP<Teuchos::Polynomial< VectorBase<Scalar> > > get_f_poly() const;
#endif // HAVE_THYRA_ME_POLYNOMIAL

    void setArgs( const OutArgs<Scalar>& outArgs, bool ignoreUnsupported = false );
    void setFailed() const;
    bool isFailed() const;
    bool isEmpty() const;
    void assertSameSupport( const OutArgs<Scalar> &outArgs ) const;
    std::string modelEvalDescription() const;
    std::string description() const;
    void describe(
      Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel
      ) const;
  protected:
    void _setModelEvalDescription( const std::string &modelEvalDescription );
    void _set_Np_Ng(int Np, int Ng);
    void _setSupports( EOutArgsMembers arg, bool supports );
    void _setSupports( EOutArgsDfDp arg, int l, const DerivativeSupport& );
    void _setSupports( EOutArgsDgDx_dot arg, int j, const DerivativeSupport& );
    void _setSupports( EOutArgsDgDx arg, int j, const DerivativeSupport& );
    void _setSupports( EOutArgsDgDp arg, int j, int l, const DerivativeSupport& );
    void _set_W_properties( const DerivativeProperties &properties );
    void _set_DfDp_properties( int l, const DerivativeProperties &properties );
    void _set_DgDx_dot_properties( int j, const DerivativeProperties &properties );
    void _set_DgDx_properties( int j, const DerivativeProperties &properties );
    void _set_DgDp_properties( int j, int l, const DerivativeProperties &properties );
    void _setSupports( const OutArgs<Scalar>& inputOutArgs );
    void _setUnsupportsAndRelated( EInArgsMembers arg );
    void _setUnsupportsAndRelated( EOutArgsMembers arg );
  private:
    // types
    typedef Teuchos::Array<RCP<VectorBase<Scalar> > > g_t;
    typedef Teuchos::Array<Derivative<Scalar> > deriv_t;
    typedef Teuchos::Array<DerivativeProperties> deriv_properties_t;
    typedef Teuchos::Array<DerivativeSupport> supports_t;
    // data
    std::string modelEvalDescription_;
    bool supports_[NUM_E_OUT_ARGS_MEMBERS];
    supports_t supports_DfDp_; // Np
    supports_t supports_DgDx_dot_; // Ng
    supports_t supports_DgDx_; // Ng
    supports_t supports_DgDp_; // Ng x Np
    RCP<VectorBase<Scalar> > f_;
    g_t g_; // Ng
    RCP<LinearOpWithSolveBase<Scalar> > W_;
    RCP<LinearOpBase<Scalar> > W_op_;
    DerivativeProperties W_properties_;
    deriv_t DfDp_; // Np
    deriv_properties_t DfDp_properties_; // Np
    deriv_t DgDx_dot_; // Ng
    deriv_t DgDx_; // Ng
    deriv_properties_t DgDx_dot_properties_; // Ng
    deriv_properties_t DgDx_properties_; // Ng
    deriv_t DgDp_; // Ng x Np
    deriv_properties_t DgDp_properties_; // Ng x Np
#ifdef HAVE_THYRA_ME_POLYNOMIAL
   RCP<Teuchos::Polynomial< VectorBase<Scalar> > > f_poly_;
#endif // HAVE_THYRA_ME_POLYNOMIAL
    mutable bool isFailed_;
    // functions
    void assert_supports(EOutArgsMembers arg) const;
    void assert_supports(
      EOutArgsDfDp arg, int l,
      const Derivative<Scalar> &deriv = Derivative<Scalar>()
      ) const;
    void assert_supports(
      EOutArgsDgDx_dot arg, int j,
      const Derivative<Scalar> &deriv = Derivative<Scalar>()
      ) const;
    void assert_supports(
      EOutArgsDgDx arg, int j,
      const Derivative<Scalar> &deriv = Derivative<Scalar>()
      ) const;
    void assert_supports(
      EOutArgsDgDp arg, int j, int l,
      const Derivative<Scalar> &deriv = Derivative<Scalar>()
      ) const;
    void assert_l(int l) const;
    void assert_j(int j) const;
  };


// Added since at least gcc 3.3.4 does not do the right thing here!
#ifdef HAVE_PROTECTED_NESTED_TEMPLATE_CLASS_ACCESS
protected:
#endif


  template<class Scalar>
  class InArgsSetup : public InArgs<Scalar> {
  public:
    InArgsSetup();
    InArgsSetup( const InArgs<Scalar>& );
    void setModelEvalDescription( const std::string &modelEvalDescription );
    void set_Np(int Np);
    void setSupports( EInArgsMembers arg, bool supports = true );
    void setSupports( const InArgs<Scalar>& inputInArgs, const int Np = -1 );
    void setUnsupportsAndRelated( EInArgsMembers arg );
  };

  template<class Scalar>
  class OutArgsSetup : public OutArgs<Scalar> {
  public:
    OutArgsSetup();
    OutArgsSetup( const OutArgs<Scalar>& );
    void setModelEvalDescription( const std::string &modelEvalDescription );
    void set_Np_Ng(int Np, int Ng);
    void setSupports( EOutArgsMembers arg, bool supports = true );
    void setSupports(EOutArgsDfDp arg, int l, const DerivativeSupport& );
    void setSupports(EOutArgsDgDx_dot arg, int j, const DerivativeSupport& );
    void setSupports(EOutArgsDgDx arg, int j, const DerivativeSupport& );
    void setSupports(EOutArgsDgDp arg, int j, int l, const DerivativeSupport& );
    void set_W_properties( const DerivativeProperties &properties );
    void set_DfDp_properties( int l, const DerivativeProperties &properties );
    void set_DgDx_dot_properties( int j, const DerivativeProperties &properties );
    void set_DgDx_properties( int j, const DerivativeProperties &properties );
    void set_DgDp_properties( int j, int l, const DerivativeProperties &properties );
    void setSupports( const OutArgs<Scalar>& inputOutArgs );
    void setUnsupportsAndRelated( EInArgsMembers arg );
    void setUnsupportsAndRelated( EOutArgsMembers arg );
   };


};


std::string toString(ModelEvaluatorBase::EInArgsMembers);


std::string toString(ModelEvaluatorBase::EOutArgsMembers);


std::string toString(
  ModelEvaluatorBase::EDerivativeMultiVectorOrientation orientation
  );


ModelEvaluatorBase::EDerivativeMultiVectorOrientation
getOtherDerivativeMultiVectorOrientation(
  ModelEvaluatorBase::EDerivativeMultiVectorOrientation orientation
  );


} // namespace Thyra


// //////////////////////////////////
// Inline Defintions


//
// Thyra_MEB_helper_functions_grp
//


inline
std::string Thyra::toString(ModelEvaluatorBase::EInArgsMembers arg)
{
  switch(arg) {
    case ModelEvaluatorBase::IN_ARG_x_dot:
      return "IN_ARG_x_dot";
    case ModelEvaluatorBase::IN_ARG_x:
      return "IN_ARG_x";
    case ModelEvaluatorBase::IN_ARG_x_dot_poly:
      return "IN_ARG_x_dot_poly";
    case ModelEvaluatorBase::IN_ARG_x_poly:
      return "IN_ARG_x_poly";
    case ModelEvaluatorBase::IN_ARG_t:
      return "IN_ARG_t";
    case ModelEvaluatorBase::IN_ARG_alpha:
      return "IN_ARG_alpha";
    case ModelEvaluatorBase::IN_ARG_beta:
      return "IN_ARG_beta";
#ifdef TEUCHOS_DEBUG
    default:
      TEST_FOR_EXCEPT(true);
#endif
  }
  return ""; // Will never be executed!
}


inline
std::string Thyra::toString(ModelEvaluatorBase::EOutArgsMembers arg)
{
  switch(arg) {
    case ModelEvaluatorBase::OUT_ARG_f:
      return "OUT_ARG_f";
    case ModelEvaluatorBase::OUT_ARG_W:
      return "OUT_ARG_W";
    case ModelEvaluatorBase::OUT_ARG_W_op:
      return "OUT_ARG_W_op";
    case ModelEvaluatorBase::OUT_ARG_f_poly:
      return "OUT_ARG_f_poly";
#ifdef TEUCHOS_DEBUG
    default:
      TEST_FOR_EXCEPT(true);
#endif
  }
  return ""; // Will never be executed!
}


inline
std::string Thyra::toString(
  ModelEvaluatorBase::EDerivativeMultiVectorOrientation orientation
  )
{
  switch(orientation) {
    case ModelEvaluatorBase::DERIV_MV_BY_COL:
      return "DERIV_MV_BY_COL";
    case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW:
      return "DERIV_TRANS_MV_BY_ROW";
#ifdef TEUCHOS_DEBUG
    default:
      TEST_FOR_EXCEPT(true);
#endif
  }
  return ""; // Should never execute this!
}


inline
Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation
Thyra::getOtherDerivativeMultiVectorOrientation(
  ModelEvaluatorBase::EDerivativeMultiVectorOrientation orientation
  )
{
  switch(orientation) {
    case ModelEvaluatorBase::DERIV_MV_BY_COL:
      return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
    case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW:
      return ModelEvaluatorBase::DERIV_MV_BY_COL;
#ifdef TEUCHOS_DEBUG
    default:
      TEST_FOR_EXCEPT(true);
#endif
  }
  return ModelEvaluatorBase::DERIV_MV_BY_COL; // Should never execute this!
}


#endif // THYRA_MODEL_EVALUATOR_BASE_DECL_HPP