42#ifndef SACADO_ETPCE_ORTHOGPOLY_HPP
43#define SACADO_ETPCE_ORTHOGPOLY_HPP
47#ifdef HAVE_STOKHOS_SACADO
49#include "Teuchos_RCP.hpp"
51#include "Sacado_Traits.hpp"
52#include "Sacado_Handle.hpp"
59#include "Sacado_mpl_apply.hpp"
63#ifdef HAVE_STOKHOS_THRUST
64#include "thrust/tuple.h"
76 template <
int k,
typename T> KERNEL_PREFIX T&
77 get(T* a) {
return a[k]; }
78 template <
int k,
typename T> KERNEL_PREFIX
const T&
79 get(
const T* a) {
return a[k]; }
81 template <
int k,
int N,
typename T> KERNEL_PREFIX T&
82 get(T a[N]) {
return a[k]; }
83 template <
int k,
int N,
typename T> KERNEL_PREFIX
const T&
84 get(
const T a[N]) {
return a[k]; }
91 template <
typename ExprT>
class Expr {};
94 template <
typename T,
typename S>
class OrthogPoly;
97 template <
typename T,
typename Storage >
98 class OrthogPolyImpl {
125 typedef typename approx_type::pointer pointer;
126 typedef typename approx_type::const_pointer const_pointer;
127 typedef typename approx_type::reference reference;
128 typedef typename approx_type::const_reference const_reference;
141 OrthogPolyImpl(
const value_type& x);
147 OrthogPolyImpl(
const Teuchos::RCP<expansion_type>& expansion);
153 OrthogPolyImpl(
const Teuchos::RCP<expansion_type>& expansion,
157 OrthogPolyImpl(
const OrthogPolyImpl& x);
160 template <
typename S> OrthogPolyImpl(
const Expr<S>& x);
166 void init(
const T& v) { th_->init(v); }
169 void init(
const T* v) { th_->init(v); }
172 template <
typename S>
173 void init(
const OrthogPolyImpl<T,S>& v) { th_->init(v.getOrthogPolyApprox()); }
176 void load(T* v) { th_->load(v); }
179 template <
typename S>
180 void load(OrthogPolyImpl<T,S>& v) { th_->load(v.getOrthogPolyApprox()); }
186 void reset(
const Teuchos::RCP<expansion_type>& expansion);
192 void reset(
const Teuchos::RCP<expansion_type>& expansion,
205 void copyForWrite() { th_.makeOwnCopy(); }
208 value_type evaluate(
const Teuchos::Array<value_type>& point)
const;
211 value_type evaluate(
const Teuchos::Array<value_type>& point,
212 const Teuchos::Array<value_type>& bvals)
const;
215 value_type mean()
const {
return th_->mean(); }
218 value_type standard_deviation()
const {
return th_->standard_deviation(); }
221 value_type two_norm()
const {
return th_->two_norm(); }
224 value_type two_norm_squared()
const {
return th_->two_norm_squared(); }
227 value_type inner_product(
const OrthogPolyImpl& b)
const {
228 return th_->inner_product(b.getOrthogPolyApprox()); }
231 std::ostream& print(std::ostream& os)
const {
return th_->print(os); }
234 template <
typename S>
bool isEqualTo(
const Expr<S>& x)
const;
242 OrthogPolyImpl& operator=(
const value_type&
val);
245 OrthogPolyImpl& operator=(
const OrthogPolyImpl& x);
248 template <
typename S>
249 OrthogPolyImpl& operator=(
const Expr<S>& x);
259 Teuchos::RCP<basis_type> basis()
const {
return th_->basis(); }
262 Teuchos::RCP<expansion_type> expansion()
const {
return expansion_; }
265 Teuchos::RCP<quad_expansion_type> quad_expansion()
const {
return quad_expansion_; }
275 const_reference
val()
const {
return (*th_)[0]; }
278 reference
val() {
return (*th_)[0]; }
291 bool hasFastAccess(ordinal_type sz)
const {
return th_->size()>=sz;}
294 const_pointer coeff()
const {
return th_->coeff();}
297 pointer coeff() {
return th_->coeff();}
301 return i<th_->size() ? (*th_)[i]:
value_type(0.); }
310 reference term(ordinal_type dimension, ordinal_type order) {
311 return th_->term(dimension, order); }
314 const_reference term(ordinal_type dimension, ordinal_type order)
const {
315 return th_->term(dimension, order); }
318 Teuchos::Array<ordinal_type> order(ordinal_type term)
const {
319 return th_->order(term); }
329 OrthogPolyImpl& operator += (
const value_type& x);
332 OrthogPolyImpl& operator -= (
const value_type& x);
335 OrthogPolyImpl& operator *= (
const value_type& x);
338 OrthogPolyImpl& operator /= (
const value_type& x);
343 const approx_type& getOrthogPolyApprox()
const {
return *th_; }
346 approx_type& getOrthogPolyApprox() {
return *th_; }
351 template <
typename S>
void expressionCopy(
const Expr<S>& x);
356 Teuchos::RCP<expansion_type> expansion_;
359 Teuchos::RCP<quad_expansion_type> quad_expansion_;
362 Teuchos::RCP<expansion_type> const_expansion_;
365 Sacado::Handle< Stokhos::OrthogPolyApprox<int,value_type,Storage> > th_;
375 template <
typename T,
typename Storage>
376 class Expr< OrthogPolyImpl<T,
Storage> > :
377 public OrthogPolyImpl<T,Storage> {
382 typedef typename OrthogPolyImpl<T,Storage>::value_type
value_type;
383 typedef typename OrthogPolyImpl<T,Storage>::scalar_type
scalar_type;
384 typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;
385 typedef typename OrthogPolyImpl<T,Storage>::storage_type
storage_type;
386 typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;
388 typedef OrthogPoly<T,Storage> base_expr_type;
391 static const int num_args = 1;
402 Expr(
const Teuchos::RCP<
typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion) :
403 OrthogPolyImpl<T,
Storage>(expansion) {}
406 Expr(
const Teuchos::RCP<
typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion,
407 typename OrthogPolyImpl<T,Storage>::ordinal_type sz) :
408 OrthogPolyImpl<T,
Storage>(expansion, sz) {}
411 Expr(
const Expr& x) :
412 OrthogPolyImpl<T,
Storage>(static_cast<const OrthogPolyImpl<T,
Storage>&>(
x)) {}
415 Expr(
const OrthogPolyImpl<T,Storage>& x) :
419 template <
typename S> Expr(
const Expr<S>& x) :
425 const approx_type& getArg(
int i)
const {
426 return this->getOrthogPolyApprox(); }
428 bool has_fast_access(
int sz)
const {
return this->size() >= sz; }
430 bool has_nonconst_expansion()
const {
431 return this->expansion_ != this->const_expansion_;
434 int order()
const {
return this->size() == 1 ? 0 : 1; }
436 value_type fast_higher_order_coeff(
int i)
const {
441 return this->coeff(i);
444 template <
int offset,
typename tuple_type>
447 return get<offset>(x);
450 std::string name()
const {
return "x"; }
468 template <
typename T,
typename Storage>
469 class OrthogPoly :
public Expr< OrthogPolyImpl<T,Storage> > {
473 typedef typename OrthogPolyImpl<T,Storage>::value_type
value_type;
476 typedef typename OrthogPolyImpl<T,Storage>::scalar_type
scalar_type;
479 typedef typename OrthogPolyImpl<T,Storage>::ordinal_type
ordinal_type;
482 typedef typename OrthogPolyImpl<T,Storage>::storage_type
storage_type;
485 typedef typename OrthogPolyImpl<T,Storage>::basis_type
basis_type;
488 typedef typename OrthogPolyImpl<T,Storage>::expansion_type expansion_type;
491 typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;
493 typedef typename OrthogPolyImpl<T,Storage>::pointer pointer;
494 typedef typename OrthogPolyImpl<T,Storage>::const_pointer const_pointer;
495 typedef typename OrthogPolyImpl<T,Storage>::reference reference;
496 typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;
499 template <
typename S>
501 typedef typename Sacado::mpl::apply<Storage,ordinal_type,S>::type
storage_type;
502 typedef OrthogPoly<S,storage_type> type;
510 Expr< OrthogPolyImpl<T,
Storage> >() {}
516 OrthogPoly(
const value_type& x) :
517 Expr< OrthogPolyImpl<T,
Storage> >(
x) {}
523 OrthogPoly(
const Teuchos::RCP<expansion_type>& expansion) :
524 Expr< OrthogPolyImpl<T,
Storage> >(expansion) {}
530 OrthogPoly(
const Teuchos::RCP<expansion_type>& expansion,
532 Expr< OrthogPolyImpl<T,
Storage> >(expansion, sz) {}
535 OrthogPoly(
const OrthogPoly& x) :
536 Expr< OrthogPolyImpl<T,
Storage> >(
x) {}
539 template <
typename S> OrthogPoly(
const Expr<S>& x) :
540 Expr< OrthogPolyImpl<T,
Storage> >(
x) {}
546 OrthogPoly& operator=(
const value_type&
val) {
547 OrthogPolyImpl<T,Storage>::operator=(
val);
552 OrthogPoly& operator=(
const OrthogPoly& x) {
553 OrthogPolyImpl<T,Storage>::operator=(
static_cast<const OrthogPolyImpl<T,Storage>&
>(x));
558 OrthogPoly& operator=(
const Expr< OrthogPolyImpl<T,Storage> >& x) {
559 OrthogPolyImpl<T,Storage>::operator=(
static_cast<const OrthogPolyImpl<T,Storage>&
>(x));
564 template <
typename S>
565 OrthogPoly& operator=(
const Expr<S>& x) {
566 OrthogPolyImpl<T,Storage>::operator=(x);
573 OrthogPoly& operator += (
const value_type& x) {
574 OrthogPolyImpl<T,Storage>::operator+=(x);
579 OrthogPoly& operator -= (
const value_type& x) {
580 OrthogPolyImpl<T,Storage>::operator-=(x);
585 OrthogPoly& operator *= (
const value_type& x) {
586 OrthogPolyImpl<T,Storage>::operator*=(x);
591 OrthogPoly& operator /= (
const value_type& x) {
592 OrthogPolyImpl<T,Storage>::operator/=(x);
597 template <
typename S>
598 OrthogPoly& operator += (
const Expr<S>& x) {
604 template <
typename S>
605 OrthogPoly& operator -= (
const Expr<S>& x) {
611 template <
typename S>
612 OrthogPoly& operator *= (
const Expr<S>& x) {
618 template <
typename S>
619 OrthogPoly& operator /= (
const Expr<S>& x) {
630 template <
typename T>
631 struct IsExpr< ETPCE::Expr<T> > {
632 static const bool value =
true;
635 template <
typename T>
636 struct BaseExprType< ETPCE::Expr<T> > {
637 typedef typename ETPCE::Expr<T>::base_expr_type type;
640 template <
typename T,
typename S>
641 struct IsExpr< ETPCE::OrthogPoly<T,S> > {
642 static const bool value =
true;
645 template <
typename T,
typename S>
646 struct BaseExprType< ETPCE::OrthogPoly<T,S> > {
647 typedef ETPCE::OrthogPoly<T,S> type;
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
Stokhos::StandardStorage< int, double > storage_type
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Abstract base class for multivariate orthogonal polynomials.
Abstract base class for orthogonal polynomial-based expansions.
Orthogonal polynomial expansions based on numerical quadrature.
Stokhos::LegendreBasis< int, double > basis_type
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x