42#include "Sacado_DynamicArrayTraits.hpp"
43#include "Teuchos_TimeMonitor.hpp"
49template <
typename T,
typename Storage>
56 expansion_ = const_expansion_;
59template <
typename T,
typename Storage>
66 expansion_ = const_expansion_;
69template <
typename T,
typename Storage>
71OrthogPoly(
const Teuchos::RCP<expansion_type>& expansion) :
72 expansion_(expansion),
78template <
typename T,
typename Storage>
80OrthogPoly(
const Teuchos::RCP<expansion_type>& expansion,
82 expansion_(expansion),
88template <
typename T,
typename Storage>
91 expansion_(x.expansion_),
97template <
typename T,
typename Storage>
103template <
typename T,
typename Storage>
106reset(
const Teuchos::RCP<expansion_type>& expansion)
108 expansion_ = expansion;
109 th->reset(expansion_->getBasis());
112template <
typename T,
typename Storage>
115reset(
const Teuchos::RCP<expansion_type>& expansion, ordinal_type sz)
117 expansion_ = expansion;
118 th->reset(expansion_->getBasis(), sz);
121template <
typename T,
typename Storage>
126 return th->evaluate(point);
129template <
typename T,
typename Storage>
136 return th->evaluate(point, bvals);
139template <
typename T,
typename Storage>
143 typedef IsEqual<value_type> IE;
144 if (x.size() != this->size())
return false;
147 if (expansion_ != x.expansion_) {
150 if ((expansion_ != const_expansion_) &&
151 (x.expansion_ != x.const_expansion_))
155 for (
int i=0; i<this->size(); i++)
156 eq = eq && IE::eval(x.coeff(i), this->coeff(i));
160template <
typename T,
typename Storage>
170template <
typename T,
typename Storage>
175 expansion_ = x.expansion_;
180template <
typename T,
typename Storage>
188template <
typename T,
typename Storage>
194 expansion_->unaryMinus(*(x.th), *th);
198template <
typename T,
typename Storage>
204 expansion_->plusEqual(*th, v);
208template <
typename T,
typename Storage>
214 expansion_->minusEqual(*th, v);
218template <
typename T,
typename Storage>
224 expansion_->timesEqual(*th, v);
228template <
typename T,
typename Storage>
234 expansion_->divideEqual(*th, v);
238template <
typename T,
typename Storage>
244 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
245 if (x.size() > size()) {
249 e->plusEqual(*th, *x.th);
253template <
typename T,
typename Storage>
259 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
260 if (x.size() > size()) {
264 e->minusEqual(*th, *x.th);
268template <
typename T,
typename Storage>
274 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
275 if (x.size() > size()) {
279 e->timesEqual(*th, *x.th);
283template <
typename T,
typename Storage>
289 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
290 if (x.size() > size()) {
294 e->divideEqual(*th, *x.th);
298template <
typename T,
typename Storage>
305 ordinal_type da = a.size();
306 ordinal_type db = b.size();
307 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
308 if (da == db || da > 1)
314 e->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
315 b.getOrthogPolyApprox());
320template <
typename T,
typename Storage>
326 b.expansion()->plus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
330template <
typename T,
typename Storage>
336 a.expansion()->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
340template <
typename T,
typename Storage>
347 ordinal_type da = a.size();
348 ordinal_type db = b.size();
349 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
350 if (da == db || da > 1)
356 e->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
357 b.getOrthogPolyApprox());
362template <
typename T,
typename Storage>
368 b.expansion()->minus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
372template <
typename T,
typename Storage>
378 a.expansion()->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
382template <
typename T,
typename Storage>
389 ordinal_type da = a.size();
390 ordinal_type db = b.size();
391 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
392 if (da == db || da > 1)
398 e->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
399 b.getOrthogPolyApprox());
404template <
typename T,
typename Storage>
410 b.expansion()->times(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
414template <
typename T,
typename Storage>
420 a.expansion()->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
424template <
typename T,
typename Storage>
431 ordinal_type da = a.size();
432 ordinal_type db = b.size();
433 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
434 if (da == db || da > 1)
440 e->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
441 b.getOrthogPolyApprox());
446template <
typename T,
typename Storage>
452 b.expansion()->divide(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
456template <
typename T,
typename Storage>
462 a.expansion()->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
466template <
typename T,
typename Storage>
471 a.expansion()->exp(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
475template <
typename T,
typename Storage>
479 TEUCHOS_FUNC_TIME_MONITOR(
"LOG");
482 TEUCHOS_FUNC_TIME_MONITOR(
"OPA LOG");
483 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
490template <
typename T,
typename Storage>
495 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
498template <
typename T,
typename Storage>
503 a.expansion()->log10(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
507template <
typename T,
typename Storage>
512 a.expansion()->sqrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
516template <
typename T,
typename Storage>
521 a.expansion()->cbrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
525template <
typename T,
typename Storage>
532 ordinal_type da = a.size();
533 ordinal_type db = b.size();
534 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
535 if (da == db || da > 1)
541 e->pow(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
542 b.getOrthogPolyApprox());
547template <
typename T,
typename Storage>
553 b.expansion()->pow(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
557template <
typename T,
typename Storage>
563 a.expansion()->pow(c.getOrthogPolyApprox(),a.getOrthogPolyApprox(), b);
567template <
typename T,
typename Storage>
572 a.expansion()->sin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
576template <
typename T,
typename Storage>
581 a.expansion()->cos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
585template <
typename T,
typename Storage>
590 a.expansion()->tan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
594template <
typename T,
typename Storage>
599 a.expansion()->sinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
603template <
typename T,
typename Storage>
608 a.expansion()->cosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
612template <
typename T,
typename Storage>
617 a.expansion()->tanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
621template <
typename T,
typename Storage>
626 a.expansion()->acos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
630template <
typename T,
typename Storage>
635 a.expansion()->asin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
639template <
typename T,
typename Storage>
644 a.expansion()->atan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
648template <
typename T,
typename Storage>
653 a.expansion()->acosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
657template <
typename T,
typename Storage>
662 a.expansion()->asinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
666template <
typename T,
typename Storage>
671 a.expansion()->atanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
675template <
typename T,
typename Storage>
680 a.expansion()->fabs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
684template <
typename T,
typename Storage>
689 a.expansion()->abs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
693template <
typename T,
typename Storage>
700 ordinal_type da = a.size();
701 ordinal_type db = b.size();
702 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
703 if (da == db || da > 1)
709 e->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
710 b.getOrthogPolyApprox());
715template <
typename T,
typename Storage>
721 b.expansion()->max(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
725template <
typename T,
typename Storage>
731 a.expansion()->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
735template <
typename T,
typename Storage>
742 ordinal_type da = a.size();
743 ordinal_type db = b.size();
744 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = a.expansion();
745 if (da == db || da > 1)
751 e->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
752 b.getOrthogPolyApprox());
757template <
typename T,
typename Storage>
763 b.expansion()->min(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
767template <
typename T,
typename Storage>
773 a.expansion()->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
777template <
typename T,
typename Storage>
782 int n = std::max(a.size(), b.size());
783 for (
int i=0; i<n; i++)
784 if (a.coeff(i) != b.coeff(i))
789template <
typename T,
typename Storage>
796 for (
int i=1; i<b.size(); i++)
797 if (b.coeff(i) != T(0.0))
802template <
typename T,
typename Storage>
809 for (
int i=1; i<a.size(); i++)
810 if (a.coeff(i) != T(0.0))
815template <
typename T,
typename Storage>
823template <
typename T,
typename Storage>
831template <
typename T,
typename Storage>
839template <
typename T,
typename Storage>
844 return a.two_norm() <= b.two_norm();
847template <
typename T,
typename Storage>
852 return a <= b.two_norm();
855template <
typename T,
typename Storage>
860 return a.two_norm() <= b;
863template <
typename T,
typename Storage>
868 return a.two_norm() >= b.two_norm();
871template <
typename T,
typename Storage>
876 return a >= b.two_norm();
879template <
typename T,
typename Storage>
884 return a.two_norm() >= b;
887template <
typename T,
typename Storage>
892 return a.two_norm() < b.two_norm();
895template <
typename T,
typename Storage>
900 return a < b.two_norm();
903template <
typename T,
typename Storage>
908 return a.two_norm() < b;
911template <
typename T,
typename Storage>
916 return a.two_norm() > b.two_norm();
919template <
typename T,
typename Storage>
924 return a > b.two_norm();
927template <
typename T,
typename Storage>
932 return a.two_norm() > b;
935template <
typename T,
typename Storage>
938 for (
int i=0; i<x.size(); i++)
939 is_zero = is_zero && (x.coeff(i) == 0.0);
943template <
typename T,
typename Storage>
950template <
typename T,
typename Storage>
958template <
typename T,
typename Storage>
966template <
typename T,
typename Storage>
973template <
typename T,
typename Storage>
981template <
typename T,
typename Storage>
989template <
typename T,
typename Storage>
997 for (ordinal_type i=0; i<a.size(); i++) {
998 os << a.coeff(i) <<
" ";
1005template <
typename T,
typename Storage>
1015 for (ordinal_type i=0; i<a.size(); i++) {
1016 is >> a.fastAccessCoeff(i);
Orthogonal polynomial expansion class for constant (size 1) expansions.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
std::istream & operator>>(std::istream &is, OrthogPoly< T, Storage > &a)
bool operator>=(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > operator+(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > operator/(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator!=(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > operator-(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator==(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator&&(const OrthogPoly< T, Storage > &x1, const OrthogPoly< T, Storage > &x2)
bool operator<=(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
std::ostream & operator<<(std::ostream &os, const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > min(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > max(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator<(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > pow(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool toBool(const OrthogPoly< T, Storage > &x)
OrthogPoly< T, Storage > operator*(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator>(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
bool operator||(const OrthogPoly< T, Storage > &x1, const OrthogPoly< T, Storage > &x2)