Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Sacado_ETPCE_OrthogPoly.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef SACADO_ETPCE_ORTHOGPOLY_HPP
43 #define SACADO_ETPCE_ORTHOGPOLY_HPP
44 
45 #include "Stokhos_ConfigDefs.h"
46 
47 #ifdef HAVE_STOKHOS_SACADO
48 
49 #include "Teuchos_RCP.hpp"
50 
51 #include "Sacado_Traits.hpp"
52 #include "Sacado_Handle.hpp"
53 
57 
58 #include "Sacado_mpl_apply.hpp"
59 
60 #include <ostream> // for std::ostream
61 
62 #ifdef HAVE_STOKHOS_THRUST
63 #include "thrust/tuple.h"
64 #endif
65 
66 #ifndef KERNEL_PREFIX
67 #define KERNEL_PREFIX
68 #endif
69 
70 namespace Sacado {
71 
73  namespace ETPCE {
74 
75  template <int k, typename T> KERNEL_PREFIX T&
76  get(T* a) { return a[k]; }
77  template <int k, typename T> KERNEL_PREFIX const T&
78  get(const T* a) { return a[k]; }
79 
80  template <int k, int N, typename T> KERNEL_PREFIX T&
81  get(T a[N]) { return a[k]; }
82  template <int k, int N, typename T> KERNEL_PREFIX const T&
83  get(const T a[N]) { return a[k]; }
84 
86 
90  template <typename ExprT> class Expr {};
91 
92  // Forward declaration
93  template <typename T, typename S> class OrthogPoly;
94 
96  template <typename T, typename Storage >
97  class OrthogPolyImpl {
98  public:
99 
101  typedef T value_type;
102 
104  typedef typename ScalarType<T>::type scalar_type;
105 
107  typedef int ordinal_type;
108 
110  typedef Storage storage_type;
111 
114 
117 
120 
123 
124  typedef typename approx_type::pointer pointer;
125  typedef typename approx_type::const_pointer const_pointer;
126  typedef typename approx_type::reference reference;
127  typedef typename approx_type::const_reference const_reference;
128 
129 
131 
134  OrthogPolyImpl();
135 
137 
140  OrthogPolyImpl(const value_type& x);
141 
143 
146  OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion);
147 
149 
152  OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion,
153  ordinal_type sz);
154 
156  OrthogPolyImpl(const OrthogPolyImpl& x);
157 
159  template <typename S> OrthogPolyImpl(const Expr<S>& x);
160 
162  ~OrthogPolyImpl() {}
163 
165  void init(const T& v) { th_->init(v); }
166 
168  void init(const T* v) { th_->init(v); }
169 
171  template <typename S>
172  void init(const OrthogPolyImpl<T,S>& v) { th_->init(v.getOrthogPolyApprox()); }
173 
175  void load(T* v) { th_->load(v); }
176 
178  template <typename S>
179  void load(OrthogPolyImpl<T,S>& v) { th_->load(v.getOrthogPolyApprox()); }
180 
182 
185  void reset(const Teuchos::RCP<expansion_type>& expansion);
186 
188 
191  void reset(const Teuchos::RCP<expansion_type>& expansion,
192  ordinal_type sz);
193 
195 
204  void copyForWrite() { th_.makeOwnCopy(); }
205 
207  value_type evaluate(const Teuchos::Array<value_type>& point) const;
208 
210  value_type evaluate(const Teuchos::Array<value_type>& point,
211  const Teuchos::Array<value_type>& bvals) const;
212 
214  value_type mean() const {return th_->mean(); }
215 
217  value_type standard_deviation() const { return th_->standard_deviation(); }
218 
220  value_type two_norm() const { return th_->two_norm(); }
221 
223  value_type two_norm_squared() const { return th_->two_norm_squared(); }
224 
226  value_type inner_product(const OrthogPolyImpl& b) const {
227  return th_->inner_product(b.getOrthogPolyApprox()); }
228 
230  std::ostream& print(std::ostream& os) const { return th_->print(os); }
231 
233  template <typename S> bool isEqualTo(const Expr<S>& x) const;
234 
239 
241  OrthogPolyImpl& operator=(const value_type& val);
242 
244  OrthogPolyImpl& operator=(const OrthogPolyImpl& x);
245 
247  template <typename S>
248  OrthogPolyImpl& operator=(const Expr<S>& x);
249 
251 
256 
258  Teuchos::RCP<basis_type> basis() const { return th_->basis(); }
259 
261  Teuchos::RCP<expansion_type> expansion() const { return expansion_; }
262 
264  Teuchos::RCP<quad_expansion_type> quad_expansion() const { return quad_expansion_; }
265 
267 
272 
274  const_reference val() const { return (*th_)[0]; }
275 
277  reference val() { return (*th_)[0]; }
278 
280 
285 
287  ordinal_type size() const { return th_->size();}
288 
290  bool hasFastAccess(ordinal_type sz) const { return th_->size()>=sz;}
291 
293  const_pointer coeff() const { return th_->coeff();}
294 
296  pointer coeff() { return th_->coeff();}
297 
299  value_type coeff(ordinal_type i) const {
300  return i<th_->size() ? (*th_)[i]:value_type(0.); }
301 
303  reference fastAccessCoeff(ordinal_type i) { return (*th_)[i];}
304 
306  value_type fastAccessCoeff(ordinal_type i) const { return (*th_)[i];}
307 
309  reference term(ordinal_type dimension, ordinal_type order) {
310  return th_->term(dimension, order); }
311 
313  const_reference term(ordinal_type dimension, ordinal_type order) const {
314  return th_->term(dimension, order); }
315 
317  Teuchos::Array<ordinal_type> order(ordinal_type term) const {
318  return th_->order(term); }
319 
321 
326 
328  OrthogPolyImpl& operator += (const value_type& x);
329 
331  OrthogPolyImpl& operator -= (const value_type& x);
332 
334  OrthogPolyImpl& operator *= (const value_type& x);
335 
337  OrthogPolyImpl& operator /= (const value_type& x);
338 
340 
342  const approx_type& getOrthogPolyApprox() const { return *th_; }
343 
345  approx_type& getOrthogPolyApprox() { return *th_; }
346 
347  protected:
348 
350  template <typename S> void expressionCopy(const Expr<S>& x);
351 
352  protected:
353 
355  Teuchos::RCP<expansion_type> expansion_;
356 
358  Teuchos::RCP<quad_expansion_type> quad_expansion_;
359 
361  static Teuchos::RCP<expansion_type> const_expansion_;
362 
364  Sacado::Handle< Stokhos::OrthogPolyApprox<int,value_type,Storage> > th_;
365 
366  }; // class OrthogPolyImpl
367 
369 
374  template <typename T, typename Storage>
375  class Expr< OrthogPolyImpl<T,Storage> > :
376  public OrthogPolyImpl<T,Storage> {
377 
378  public:
379 
383  typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;
385  typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;
386 
387  typedef OrthogPoly<T,Storage> base_expr_type;
388 
390  static const int num_args = 1;
391 
393  Expr() :
394  OrthogPolyImpl<T,Storage>() {}
395 
397  Expr(const T & x) :
398  OrthogPolyImpl<T,Storage>(x) {}
399 
401  Expr(const Teuchos::RCP<typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion) :
402  OrthogPolyImpl<T,Storage>(expansion) {}
403 
405  Expr(const Teuchos::RCP<typename OrthogPolyImpl<T,Storage>::expansion_type>& expansion,
407  OrthogPolyImpl<T,Storage>(expansion, sz) {}
408 
410  Expr(const Expr& x) :
411  OrthogPolyImpl<T,Storage>(static_cast<const OrthogPolyImpl<T,Storage>&>(x)) {}
412 
414  Expr(const OrthogPolyImpl<T,Storage>& x) :
415  OrthogPolyImpl<T,Storage>(x) {}
416 
418  template <typename S> Expr(const Expr<S>& x) :
419  OrthogPolyImpl<T,Storage>(x) {}
420 
422  ~Expr() {}
423 
424  const approx_type& getArg(int i) const {
425  return this->getOrthogPolyApprox(); }
426 
427  bool has_fast_access(int sz) const { return this->size() >= sz; }
428 
429  bool has_nonconst_expansion() const {
430  return this->expansion_ != this->const_expansion_;
431  }
432 
433  int order() const { return this->size() == 1 ? 0 : 1; }
434 
435  value_type fast_higher_order_coeff(int i) const {
436  return this->fastAccessCoeff(i);
437  }
438 
439  value_type higher_order_coeff(int i) const {
440  return this->coeff(i);
441  }
442 
443  template <int offset, typename tuple_type>
444  KERNEL_PREFIX
445  value_type eval_sample(tuple_type x) const {
446  return get<offset>(x);
447  }
448 
449  std::string name() const { return "x"; }
450 
451  }; // class Expr<OrthogPolyImpl>
452 
453  } // namespace ETPCE
454 
455 } // namespace Sacado
456 
461 
462 namespace Sacado {
463 
464  namespace ETPCE {
465 
467  template <typename T, typename Storage>
468  class OrthogPoly : public Expr< OrthogPolyImpl<T,Storage> > {
469  public:
470 
473 
476 
479 
482 
485 
487  typedef typename OrthogPolyImpl<T,Storage>::expansion_type expansion_type;
488 
490  typedef typename OrthogPolyImpl<T,Storage>::approx_type approx_type;
491 
492  typedef typename OrthogPolyImpl<T,Storage>::pointer pointer;
493  typedef typename OrthogPolyImpl<T,Storage>::const_pointer const_pointer;
494  typedef typename OrthogPolyImpl<T,Storage>::reference reference;
495  typedef typename OrthogPolyImpl<T,Storage>::const_reference const_reference;
496 
498  template <typename S>
499  struct apply {
500  typedef typename Sacado::mpl::apply<Storage,ordinal_type,S>::type storage_type;
501  typedef OrthogPoly<S,storage_type> type;
502  };
503 
505 
508  OrthogPoly() :
509  Expr< OrthogPolyImpl<T,Storage> >() {}
510 
512 
515  OrthogPoly(const value_type& x) :
516  Expr< OrthogPolyImpl<T,Storage> >(x) {}
517 
519 
522  OrthogPoly(const Teuchos::RCP<expansion_type>& expansion) :
523  Expr< OrthogPolyImpl<T,Storage> >(expansion) {}
524 
526 
529  OrthogPoly(const Teuchos::RCP<expansion_type>& expansion,
530  ordinal_type sz) :
531  Expr< OrthogPolyImpl<T,Storage> >(expansion, sz) {}
532 
534  OrthogPoly(const OrthogPoly& x) :
535  Expr< OrthogPolyImpl<T,Storage> >(x) {}
536 
538  template <typename S> OrthogPoly(const Expr<S>& x) :
539  Expr< OrthogPolyImpl<T,Storage> >(x) {}
540 
542  ~OrthogPoly() {}
543 
545  OrthogPoly& operator=(const value_type& val) {
546  OrthogPolyImpl<T,Storage>::operator=(val);
547  return *this;
548  }
549 
551  OrthogPoly& operator=(const OrthogPoly& x) {
552  OrthogPolyImpl<T,Storage>::operator=(static_cast<const OrthogPolyImpl<T,Storage>&>(x));
553  return *this;
554  }
555 
557  OrthogPoly& operator=(const Expr< OrthogPolyImpl<T,Storage> >& x) {
558  OrthogPolyImpl<T,Storage>::operator=(static_cast<const OrthogPolyImpl<T,Storage>&>(x));
559  return *this;
560  }
561 
563  template <typename S>
564  OrthogPoly& operator=(const Expr<S>& x) {
565  OrthogPolyImpl<T,Storage>::operator=(x);
566  return *this;
567  }
568 
570 
572  OrthogPoly& operator += (const value_type& x) {
573  OrthogPolyImpl<T,Storage>::operator+=(x);
574  return *this;
575  }
576 
578  OrthogPoly& operator -= (const value_type& x) {
579  OrthogPolyImpl<T,Storage>::operator-=(x);
580  return *this;
581  }
582 
584  OrthogPoly& operator *= (const value_type& x) {
585  OrthogPolyImpl<T,Storage>::operator*=(x);
586  return *this;
587  }
588 
590  OrthogPoly& operator /= (const value_type& x) {
591  OrthogPolyImpl<T,Storage>::operator/=(x);
592  return *this;
593  }
594 
596  template <typename S>
597  OrthogPoly& operator += (const Expr<S>& x) {
598  *this = *this + x;
599  return *this;
600  }
601 
603  template <typename S>
604  OrthogPoly& operator -= (const Expr<S>& x) {
605  *this = *this - x;
606  return *this;
607  }
608 
610  template <typename S>
611  OrthogPoly& operator *= (const Expr<S>& x) {
612  *this = *this * x;
613  return *this;
614  }
615 
617  template <typename S>
618  OrthogPoly& operator /= (const Expr<S>& x) {
619  *this = *this / x;
620  return *this;
621  }
622 
624 
625  }; // class OrthogPoly
626 
627  } // namespace ETPCE
628 
629  template <typename T>
630  struct IsExpr< ETPCE::Expr<T> > {
631  static const bool value = true;
632  };
633 
634  template <typename T>
635  struct BaseExprType< ETPCE::Expr<T> > {
636  typedef typename ETPCE::Expr<T>::base_expr_type type;
637  };
638 
639  template <typename T, typename S>
640  struct IsExpr< ETPCE::OrthogPoly<T,S> > {
641  static const bool value = true;
642  };
643 
644  template <typename T, typename S>
645  struct BaseExprType< ETPCE::OrthogPoly<T,S> > {
646  typedef ETPCE::OrthogPoly<T,S> type;
647  };
648 
649 } // namespace Sacado
650 
651 #endif // HAVE_STOKHOS_SACADO
652 
653 #endif // SACADO_ETPCE_ORTHOGPOLY_HPP
Stokhos::StandardStorage< int, double > storage_type
expr val()
Stokhos::LegendreBasis< int, double > basis_type
Abstract base class for orthogonal polynomial-based expansions.
Abstract base class for multivariate orthogonal polynomials.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Orthogonal polynomial expansions based on numerical quadrature.