Sacado_ETPCE_OrthogPoly.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
// 
//                           Stokhos Package
//                 Copyright (2009) 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.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
// 
// ***********************************************************************
// @HEADER

#ifndef SACADO_ETPCE_ORTHOGPOLY_HPP
#define SACADO_ETPCE_ORTHOGPOLY_HPP

#include "Stokhos_ConfigDefs.h"

#ifdef HAVE_STOKHOS_SACADO

#include "Teuchos_RCP.hpp"

#include "Sacado_Traits.hpp"
#include "Sacado_Handle.hpp"

#include "Stokhos_OrthogPolyExpansion.hpp"
#include "Stokhos_QuadOrthogPolyExpansion.hpp"
#include "Stokhos_ConstantOrthogPolyExpansion.hpp"
#include "Stokhos_OrthogPolyApprox.hpp"

#include "Sacado_mpl_apply.hpp"

#include <ostream>      // for std::ostream

#ifdef HAVE_STOKHOS_THRUST
#include "thrust/tuple.h"
#endif

#ifndef KERNEL_PREFIX
#define KERNEL_PREFIX
#endif

namespace Sacado {

  namespace ETPCE {

    template <int k, typename T> KERNEL_PREFIX T&
    get(T* a) { return a[k]; }
    template <int k, typename T> KERNEL_PREFIX const T&
    get(const T* a) { return a[k]; }

    template <int k, int N, typename T> KERNEL_PREFIX T&
    get(T a[N]) { return a[k]; }
    template <int k, int N, typename T> KERNEL_PREFIX const T&
    get(const T a[N]) { return a[k]; }


    template <typename ExprT> class Expr {};

    // Forward declaration
    template <typename T, typename S> class OrthogPoly;

    template <typename T, typename Storage >
    class OrthogPolyImpl {
    public:

      typedef T value_type;

      typedef typename ScalarType<T>::type scalar_type;

      typedef int ordinal_type;

      typedef Storage storage_type;

      typedef Stokhos::OrthogPolyBasis<ordinal_type,T> basis_type;

      typedef Stokhos::OrthogPolyExpansion<ordinal_type,T,Storage> expansion_type;

      typedef Stokhos::QuadOrthogPolyExpansion<ordinal_type,T,Storage> quad_expansion_type;

      typedef Stokhos::OrthogPolyApprox<ordinal_type,T,Storage> approx_type;

      typedef typename approx_type::pointer pointer;
      typedef typename approx_type::const_pointer const_pointer;
      typedef typename approx_type::reference reference;
      typedef typename approx_type::const_reference const_reference;



      OrthogPolyImpl();


      OrthogPolyImpl(const value_type& x);


      OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion);


      OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion,
                     ordinal_type sz);

      OrthogPolyImpl(const OrthogPolyImpl& x);

      template <typename S> OrthogPolyImpl(const Expr<S>& x);

      ~OrthogPolyImpl() {}

      void init(const T& v) { th_->init(v); }

      void init(const T* v) { th_->init(v); }

      template <typename S>
      void init(const OrthogPolyImpl<T,S>& v) { th_->init(v.getOrthogPolyApprox()); }

      void load(T* v) { th_->load(v); }

      template <typename S>
      void load(OrthogPolyImpl<T,S>& v) { th_->load(v.getOrthogPolyApprox()); }


      void reset(const Teuchos::RCP<expansion_type>& expansion);


      void reset(const Teuchos::RCP<expansion_type>& expansion,
                 ordinal_type sz);


      void copyForWrite() { th_.makeOwnCopy(); }

      value_type evaluate(const Teuchos::Array<value_type>& point) const;

      value_type evaluate(const Teuchos::Array<value_type>& point,
                          const Teuchos::Array<value_type>& bvals) const;

      value_type mean() const {return th_->mean(); }

      value_type standard_deviation() const { return th_->standard_deviation(); }

      value_type two_norm() const { return th_->two_norm(); }

      value_type two_norm_squared() const { return th_->two_norm_squared(); }

      value_type inner_product(const OrthogPolyImpl& b) const {
        return th_->inner_product(b.getOrthogPolyApprox()); }

      std::ostream& print(std::ostream& os) const { return th_->print(os); }

      template <typename S> bool isEqualTo(const Expr<S>& x) const;


      OrthogPolyImpl& operator=(const value_type& val);

      OrthogPolyImpl& operator=(const OrthogPolyImpl& x);

      template <typename S>
      OrthogPolyImpl& operator=(const Expr<S>& x);



      Teuchos::RCP<basis_type> basis() const { return th_->basis(); }

      Teuchos::RCP<expansion_type> expansion() const { return expansion_; }

      Teuchos::RCP<quad_expansion_type> quad_expansion() const { return quad_expansion_; }



      const_reference val() const { return (*th_)[0]; }

      reference val() { return (*th_)[0]; }



      ordinal_type size() const { return th_->size();}

      bool hasFastAccess(ordinal_type sz) const { return th_->size()>=sz;}

      const_pointer coeff() const { return th_->coeff();}

      pointer coeff() { return th_->coeff();}

      value_type coeff(ordinal_type i) const {
        return i<th_->size() ? (*th_)[i]:value_type(0.); }

      reference fastAccessCoeff(ordinal_type i) { return (*th_)[i];}

      value_type fastAccessCoeff(ordinal_type i) const { return (*th_)[i];}

      reference term(ordinal_type dimension, ordinal_type order) {
        return th_->term(dimension, order); }

      const_reference term(ordinal_type dimension, ordinal_type order) const {
        return th_->term(dimension, order); }

      Teuchos::Array<ordinal_type> order(ordinal_type term) const {
        return th_->order(term); }



      OrthogPolyImpl& operator += (const value_type& x);

      OrthogPolyImpl& operator -= (const value_type& x);

      OrthogPolyImpl& operator *= (const value_type& x);

      OrthogPolyImpl& operator /= (const value_type& x);


      const approx_type& getOrthogPolyApprox() const { return *th_; }

      approx_type& getOrthogPolyApprox() { return *th_; }

    protected:

      template <typename S> void expressionCopy(const Expr<S>& x);

    protected:

      Teuchos::RCP<expansion_type> expansion_;

      Teuchos::RCP<quad_expansion_type> quad_expansion_;

      Teuchos::RCP<expansion_type> const_expansion_;

      Sacado::Handle< Stokhos::OrthogPolyApprox<int,value_type,Storage> > th_;

    }; // class OrthogPolyImpl


    template <typename T, typename Storage>
    class Expr< OrthogPolyImpl<T,Storage> > :
        public OrthogPolyImpl<T,Storage> {

    public:

      typedef typename OrthogPolyImpl<T,Storage>::value_type value_type;
      typedef typename OrthogPolyImpl<T,Storage>::scalar_type scalar_type;
      typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;
      typedef typename OrthogPolyImpl<T,Storage>::storage_type storage_type;
      typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;

      typedef OrthogPoly<T,Storage> base_expr_type;

      static const int num_args = 1;

      Expr() :
        OrthogPolyImpl<T,Storage>() {}

      Expr(const T & x) :
        OrthogPolyImpl<T,Storage>(x) {}

      Expr(const Teuchos::RCP<typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion) :
        OrthogPolyImpl<T,Storage>(expansion) {}

      Expr(const Teuchos::RCP<typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion,
           typename OrthogPolyImpl<T,Storage>::ordinal_type sz) :
        OrthogPolyImpl<T,Storage>(expansion, sz) {}

      Expr(const Expr& x) :
        OrthogPolyImpl<T,Storage>(static_cast<const OrthogPolyImpl<T,Storage>&>(x)) {}

      Expr(const OrthogPolyImpl<T,Storage>& x) :
        OrthogPolyImpl<T,Storage>(x) {}

      template <typename S> Expr(const Expr<S>& x) :
        OrthogPolyImpl<T,Storage>(x) {}

      ~Expr() {}

      const approx_type& getArg(int i) const {
        return this->getOrthogPolyApprox(); }

      bool has_fast_access(int sz) const { return this->size() >= sz; }

      bool has_nonconst_expansion() const {
        return this->expansion_ != this->const_expansion_;
      }

      int order() const { return this->size() == 1 ? 0 : 1; }

      value_type fast_higher_order_coeff(int i) const {
        return this->fastAccessCoeff(i);
      }

      value_type higher_order_coeff(int i) const {
        return this->coeff(i);
      }

      template <int offset, typename tuple_type>
      KERNEL_PREFIX
      value_type eval_sample(tuple_type x) const {
        return get<offset>(x);
      }

      std::string name() const { return "x"; }

    }; // class Expr<OrthogPolyImpl>

  } // namespace ETPCE

} // namespace Sacado

#include "Sacado_ETPCE_ExpressionTraits.hpp"
#include "Sacado_ETPCE_OrthogPolyTraits.hpp"
#include "Sacado_ETPCE_OrthogPolyImp.hpp"
#include "Sacado_ETPCE_OrthogPolyOps.hpp"

namespace Sacado {

  namespace ETPCE {

    template <typename T, typename Storage>
    class OrthogPoly : public Expr< OrthogPolyImpl<T,Storage> > {
    public:

      typedef typename OrthogPolyImpl<T,Storage>::value_type value_type;

      typedef typename OrthogPolyImpl<T,Storage>::scalar_type scalar_type;

      typedef typename OrthogPolyImpl<T,Storage>::ordinal_type ordinal_type;

      typedef typename OrthogPolyImpl<T,Storage>::storage_type storage_type;

      typedef typename OrthogPolyImpl<T,Storage>::basis_type basis_type;

      typedef typename OrthogPolyImpl<T,Storage>::expansion_type expansion_type;

      typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;

      typedef typename OrthogPolyImpl<T,Storage>::pointer pointer;
      typedef typename OrthogPolyImpl<T,Storage>::const_pointer const_pointer;
      typedef typename OrthogPolyImpl<T,Storage>::reference reference;
      typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;

      template <typename S>
      struct apply {
        typedef typename Sacado::mpl::apply<Storage,ordinal_type,S>::type storage_type;
        typedef OrthogPoly<S,storage_type> type;
      };


      OrthogPoly() :
        Expr< OrthogPolyImpl<T,Storage> >() {}


      OrthogPoly(const value_type& x) :
        Expr< OrthogPolyImpl<T,Storage> >(x) {}


      OrthogPoly(const Teuchos::RCP<expansion_type>& expansion) :
        Expr< OrthogPolyImpl<T,Storage> >(expansion) {}


      OrthogPoly(const Teuchos::RCP<expansion_type>& expansion,
                     ordinal_type sz) :
        Expr< OrthogPolyImpl<T,Storage> >(expansion, sz) {}

      OrthogPoly(const OrthogPoly& x) :
        Expr< OrthogPolyImpl<T,Storage> >(x) {}

      template <typename S> OrthogPoly(const Expr<S>& x) :
        Expr< OrthogPolyImpl<T,Storage> >(x) {}

      ~OrthogPoly() {}

      OrthogPoly& operator=(const value_type& val) {
        OrthogPolyImpl<T,Storage>::operator=(val);
        return *this;
      }

      OrthogPoly& operator=(const OrthogPoly& x) {
        OrthogPolyImpl<T,Storage>::operator=(static_cast<const OrthogPolyImpl<T,Storage>&>(x));
        return *this;
      }

      OrthogPoly& operator=(const Expr< OrthogPolyImpl<T,Storage> >& x) {
        OrthogPolyImpl<T,Storage>::operator=(static_cast<const OrthogPolyImpl<T,Storage>&>(x));
        return *this;
      }

      template <typename S>
      OrthogPoly& operator=(const Expr<S>& x) {
        OrthogPolyImpl<T,Storage>::operator=(x);
        return *this;
      }


      OrthogPoly& operator += (const value_type& x) {
        OrthogPolyImpl<T,Storage>::operator+=(x);
        return *this;
      }

      OrthogPoly& operator -= (const value_type& x) {
        OrthogPolyImpl<T,Storage>::operator-=(x);
        return *this;
      }

      OrthogPoly& operator *= (const value_type& x) {
        OrthogPolyImpl<T,Storage>::operator*=(x);
        return *this;
      }

      OrthogPoly& operator /= (const value_type& x) {
        OrthogPolyImpl<T,Storage>::operator/=(x);
        return *this;
      }

      template <typename S>
      OrthogPoly& operator += (const Expr<S>& x) {
        *this = *this + x;
        return *this;
      }

      template <typename S>
      OrthogPoly& operator -= (const Expr<S>& x) {
        *this = *this - x;
        return *this;
      }

      template <typename S>
      OrthogPoly& operator *= (const Expr<S>& x) {
        *this = *this * x;
        return *this;
      }

      template <typename S>
      OrthogPoly& operator /= (const Expr<S>& x) {
        *this = *this / x;
        return *this;
      }


    }; // class OrthogPoly

  } // namespace ETPCE

  template <typename T>
  struct IsExpr< ETPCE::Expr<T> > {
    static const bool value = true;
  };

  template <typename T>
  struct BaseExprType< ETPCE::Expr<T> > {
    typedef typename ETPCE::Expr<T>::base_expr_type type;
  };

  template <typename T, typename S>
  struct IsExpr< ETPCE::OrthogPoly<T,S> > {
    static const bool value = true;
  };

  template <typename T, typename S>
  struct BaseExprType< ETPCE::OrthogPoly<T,S> > {
    typedef ETPCE::OrthogPoly<T,S> type;
  };

} // namespace Sacado

#endif // HAVE_STOKHOS_SACADO

#endif // SACADO_ETPCE_ORTHOGPOLY_HPP