Sacado_PCE_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_PCE_ORTHOGPOLY_HPP
#define SACADO_PCE_ORTHOGPOLY_HPP

#include "Stokhos_ConfigDefs.h"

#ifdef HAVE_STOKHOS_SACADO

#include "Teuchos_RCP.hpp"

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

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

#include <cmath>
#include <algorithm>  // for std::min and std::max
#include <ostream>  // for std::ostream

namespace Sacado {

  namespace PCE {


    template <typename T, typename Storage > 
    class OrthogPoly {
    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::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;

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


      OrthogPoly();


      OrthogPoly(const value_type& x);


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


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

      OrthogPoly(const OrthogPoly& x);

      ~OrthogPoly();

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

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

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

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

      template <typename S>
      void load(OrthogPoly<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 OrthogPoly& b) const { 
  return th->inner_product(b.getOrthogPolyApprox()); }

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

      bool isEqualTo(const OrthogPoly& x) const;


      OrthogPoly<T,Storage>& operator=(const value_type& val);

      OrthogPoly<T,Storage>& operator=(const OrthogPoly<T,Storage>& x);



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

      Teuchos::RCP<expansion_type> expansion() const { return 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 { 
  value_type tmp= i<th->size() ? (*th)[i]:value_type(0.); return tmp;}
    
      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); }
    


      OrthogPoly<T,Storage> operator + () const;

      OrthogPoly<T,Storage> operator - () const;

      OrthogPoly<T,Storage>& operator += (const value_type& x);

      OrthogPoly<T,Storage>& operator -= (const value_type& x);

      OrthogPoly<T,Storage>& operator *= (const value_type& x);

      OrthogPoly<T,Storage>& operator /= (const value_type& x);

      OrthogPoly<T,Storage>& operator += (const OrthogPoly<T,Storage>& x);

      OrthogPoly<T,Storage>& operator -= (const OrthogPoly<T,Storage>& x);
  
      OrthogPoly<T,Storage>& operator *= (const OrthogPoly<T,Storage>& x);

      OrthogPoly<T,Storage>& operator /= (const OrthogPoly<T,Storage>& x);


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

      approx_type& getOrthogPolyApprox() { return *th; }

    protected:

      Teuchos::RCP<expansion_type> expansion_;

      static Teuchos::RCP<expansion_type> const_expansion_;

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

    }; // class Hermite

    template <typename T, typename Storage>
    Teuchos::RCP<typename OrthogPoly<T, Storage>::expansion_type>
    OrthogPoly<T,Storage>::const_expansion_ =
      Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);

    // Operations
    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator+(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator+(const typename OrthogPoly<T,Storage>::value_type& a, 
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator+(const OrthogPoly<T,Storage>& a, 
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator-(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator-(const typename OrthogPoly<T,Storage>::value_type& a, 
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator-(const OrthogPoly<T,Storage>& a, 
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator*(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator*(const typename OrthogPoly<T,Storage>::value_type& a, 
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator*(const OrthogPoly<T,Storage>& a, 
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator/(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator/(const typename OrthogPoly<T,Storage>::value_type& a, 
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    operator/(const OrthogPoly<T,Storage>& a, 
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    exp(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    log(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> void
    log(OrthogPoly<T,Storage>& c, const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    log10(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    sqrt(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    cbrt(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    pow(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    pow(const T& a, 
  const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    pow(const OrthogPoly<T,Storage>& a, 
  const T& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    cos(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    sin(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    tan(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    cosh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage>
    sinh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    tanh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    acos(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    asin(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    atan(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    atan2(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    atan2(const typename OrthogPoly<T,Storage>::value_type& a, 
    const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    atan2(const OrthogPoly<T,Storage>& a, 
    const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    acosh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    asinh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    atanh(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    abs(const OrthogPoly<T,Storage>& a);
    
    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    fabs(const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    max(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    max(const typename OrthogPoly<T,Storage>::value_type& a, 
  const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    max(const OrthogPoly<T,Storage>& a, 
  const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    min(const OrthogPoly<T,Storage>& a, const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    min(const typename OrthogPoly<T,Storage>::value_type& a, 
  const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> OrthogPoly<T,Storage> 
    min(const OrthogPoly<T,Storage>& a, 
  const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator==(const OrthogPoly<T,Storage>& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator==(const typename OrthogPoly<T,Storage>::value_type& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator==(const OrthogPoly<T,Storage>& a,
         const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator!=(const OrthogPoly<T,Storage>& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator!=(const typename OrthogPoly<T,Storage>::value_type& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator!=(const OrthogPoly<T,Storage>& a,
         const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator<=(const OrthogPoly<T,Storage>& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator<=(const typename OrthogPoly<T,Storage>::value_type& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator<=(const OrthogPoly<T,Storage>& a,
         const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator>=(const OrthogPoly<T,Storage>& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator>=(const typename OrthogPoly<T,Storage>::value_type& a,
         const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator>=(const OrthogPoly<T,Storage>& a,
         const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator<(const OrthogPoly<T,Storage>& a,
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator<(const typename OrthogPoly<T,Storage>::value_type& a,
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator<(const OrthogPoly<T,Storage>& a,
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> bool 
    operator>(const OrthogPoly<T,Storage>& a,
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator>(const typename OrthogPoly<T,Storage>::value_type& a,
        const OrthogPoly<T,Storage>& b);

    template <typename T, typename Storage> bool 
    operator>(const OrthogPoly<T,Storage>& a,
        const typename OrthogPoly<T,Storage>::value_type& b);

    template <typename T, typename Storage> std::ostream& 
    operator << (std::ostream& os, const OrthogPoly<T,Storage>& a);

    template <typename T, typename Storage> std::istream& 
    operator >> (std::istream& os, OrthogPoly<T,Storage>& a);

  } // namespace PCE

} // namespace Sacado

#include "Sacado_PCE_OrthogPolyTraits.hpp"
#include "Sacado_PCE_OrthogPolyImp.hpp"

#endif // HAVE_STOKHOS_SACADO

#endif // SACADO_PCE_ORTHOGPOLY_HPP