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>
161OrthogPoly<T,Storage>&
170template <
typename T,
typename Storage>
171OrthogPoly<T,Storage>&
175 expansion_ = x.expansion_;
180template <
typename T,
typename Storage>
188template <
typename T,
typename Storage>
193 OrthogPoly<T,Storage> x(expansion_);
194 expansion_->unaryMinus(*(x.th), *th);
198template <
typename T,
typename Storage>
199OrthogPoly<T,Storage>&
204 expansion_->plusEqual(*th, v);
208template <
typename T,
typename Storage>
209OrthogPoly<T,Storage>&
214 expansion_->minusEqual(*th, v);
218template <
typename T,
typename Storage>
219OrthogPoly<T,Storage>&
224 expansion_->timesEqual(*th, v);
228template <
typename T,
typename Storage>
229OrthogPoly<T,Storage>&
234 expansion_->divideEqual(*th, v);
238template <
typename T,
typename Storage>
239OrthogPoly<T,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>
254OrthogPoly<T,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>
269OrthogPoly<T,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>
284OrthogPoly<T,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>
301 const OrthogPoly<T,Storage>& b)
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)
313 OrthogPoly<T,Storage> c(e, 0);
314 e->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
315 b.getOrthogPolyApprox());
320template <
typename T,
typename Storage>
323 const OrthogPoly<T,Storage>& b)
325 OrthogPoly<T,Storage> c(b.expansion(), 0);
326 b.expansion()->plus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
330template <
typename T,
typename Storage>
335 OrthogPoly<T,Storage> c(a.expansion(), 0);
336 a.expansion()->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
340template <
typename T,
typename Storage>
343 const OrthogPoly<T,Storage>& b)
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)
355 OrthogPoly<T,Storage> c(e, 0);
356 e->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
357 b.getOrthogPolyApprox());
362template <
typename T,
typename Storage>
365 const OrthogPoly<T,Storage>& b)
367 OrthogPoly<T,Storage> c(b.expansion(), 0);
368 b.expansion()->minus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
372template <
typename T,
typename Storage>
377 OrthogPoly<T,Storage> c(a.expansion(), 0);
378 a.expansion()->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
382template <
typename T,
typename Storage>
385 const OrthogPoly<T,Storage>& b)
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)
397 OrthogPoly<T,Storage> c(e, 0);
398 e->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
399 b.getOrthogPolyApprox());
404template <
typename T,
typename Storage>
407 const OrthogPoly<T,Storage>& b)
409 OrthogPoly<T,Storage> c(b.expansion(), 0);
410 b.expansion()->times(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
414template <
typename T,
typename Storage>
419 OrthogPoly<T,Storage> c(a.expansion(), 0);
420 a.expansion()->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
424template <
typename T,
typename Storage>
427 const OrthogPoly<T,Storage>& b)
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)
439 OrthogPoly<T,Storage> c(e, 0);
440 e->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
441 b.getOrthogPolyApprox());
446template <
typename T,
typename Storage>
449 const OrthogPoly<T,Storage>& b)
451 OrthogPoly<T,Storage> c(b.expansion(), 0);
452 b.expansion()->divide(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
456template <
typename T,
typename Storage>
461 OrthogPoly<T,Storage> c(a.expansion(), 0);
462 a.expansion()->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
466template <
typename T,
typename Storage>
468exp(
const OrthogPoly<T,Storage>& a)
470 OrthogPoly<T,Storage> c(a.expansion(), 0);
471 a.expansion()->exp(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
475template <
typename T,
typename Storage>
477log(
const OrthogPoly<T,Storage>& a)
479 TEUCHOS_FUNC_TIME_MONITOR(
"LOG");
480 OrthogPoly<T,Storage> c(a.expansion(), 0);
482 TEUCHOS_FUNC_TIME_MONITOR(
"OPA LOG");
483 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
486 OrthogPoly<T,Storage> d(c);
490template <
typename T,
typename Storage>
492log(OrthogPoly<T,Storage>& c,
const OrthogPoly<T,Storage>& a)
494 OrthogPoly<T,Storage> d(a.expansion(), 0);
495 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
498template <
typename T,
typename Storage>
500log10(
const OrthogPoly<T,Storage>& a)
502 OrthogPoly<T,Storage> c(a.expansion(), 0);
503 a.expansion()->log10(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
507template <
typename T,
typename Storage>
509sqrt(
const OrthogPoly<T,Storage>& a)
511 OrthogPoly<T,Storage> c(a.expansion(), 0);
512 a.expansion()->sqrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
516template <
typename T,
typename Storage>
518cbrt(
const OrthogPoly<T,Storage>& a)
520 OrthogPoly<T,Storage> c(a.expansion(), 0);
521 a.expansion()->cbrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
525template <
typename T,
typename Storage>
527pow(
const OrthogPoly<T,Storage>& a,
528 const OrthogPoly<T,Storage>& b)
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)
540 OrthogPoly<T,Storage> c(e, 0);
541 e->pow(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
542 b.getOrthogPolyApprox());
547template <
typename T,
typename Storage>
550 const OrthogPoly<T,Storage>& b)
552 OrthogPoly<T,Storage> c(b.expansion(), 0);
553 b.expansion()->pow(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
557template <
typename T,
typename Storage>
559pow(
const OrthogPoly<T,Storage>& a,
562 OrthogPoly<T,Storage> c(a.expansion(), 0);
563 a.expansion()->pow(c.getOrthogPolyApprox(),a.getOrthogPolyApprox(), b);
567template <
typename T,
typename Storage>
569sin(
const OrthogPoly<T,Storage>& a)
571 OrthogPoly<T,Storage> c(a.expansion(), 0);
572 a.expansion()->sin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
576template <
typename T,
typename Storage>
578cos(
const OrthogPoly<T,Storage>& a)
580 OrthogPoly<T,Storage> c(a.expansion(), 0);
581 a.expansion()->cos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
585template <
typename T,
typename Storage>
587tan(
const OrthogPoly<T,Storage>& a)
589 OrthogPoly<T,Storage> c(a.expansion(), 0);
590 a.expansion()->tan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
594template <
typename T,
typename Storage>
596sinh(
const OrthogPoly<T,Storage>& a)
598 OrthogPoly<T,Storage> c(a.expansion(), 0);
599 a.expansion()->sinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
603template <
typename T,
typename Storage>
605cosh(
const OrthogPoly<T,Storage>& a)
607 OrthogPoly<T,Storage> c(a.expansion(), 0);
608 a.expansion()->cosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
612template <
typename T,
typename Storage>
614tanh(
const OrthogPoly<T,Storage>& a)
616 OrthogPoly<T,Storage> c(a.expansion(), 0);
617 a.expansion()->tanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
621template <
typename T,
typename Storage>
623acos(
const OrthogPoly<T,Storage>& a)
625 OrthogPoly<T,Storage> c(a.expansion(), 0);
626 a.expansion()->acos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
630template <
typename T,
typename Storage>
632asin(
const OrthogPoly<T,Storage>& a)
634 OrthogPoly<T,Storage> c(a.expansion(), 0);
635 a.expansion()->asin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
639template <
typename T,
typename Storage>
641atan(
const OrthogPoly<T,Storage>& a)
643 OrthogPoly<T,Storage> c(a.expansion(), 0);
644 a.expansion()->atan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
648template <
typename T,
typename Storage>
650acosh(
const OrthogPoly<T,Storage>& a)
652 OrthogPoly<T,Storage> c(a.expansion(), 0);
653 a.expansion()->acosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
657template <
typename T,
typename Storage>
659asinh(
const OrthogPoly<T,Storage>& a)
661 OrthogPoly<T,Storage> c(a.expansion(), 0);
662 a.expansion()->asinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
666template <
typename T,
typename Storage>
668atanh(
const OrthogPoly<T,Storage>& a)
670 OrthogPoly<T,Storage> c(a.expansion(), 0);
671 a.expansion()->atanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
675template <
typename T,
typename Storage>
677fabs(
const OrthogPoly<T,Storage>& a)
679 OrthogPoly<T,Storage> c(a.expansion(), 0);
680 a.expansion()->fabs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
684template <
typename T,
typename Storage>
686abs(
const OrthogPoly<T,Storage>& a)
688 OrthogPoly<T,Storage> c(a.expansion(), 0);
689 a.expansion()->abs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
693template <
typename T,
typename Storage>
695max(
const OrthogPoly<T,Storage>& a,
696 const OrthogPoly<T,Storage>& b)
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)
708 OrthogPoly<T,Storage> c(e, 0);
709 e->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
710 b.getOrthogPolyApprox());
715template <
typename T,
typename Storage>
718 const OrthogPoly<T,Storage>& b)
720 OrthogPoly<T,Storage> c(b.expansion(), 0);
721 b.expansion()->max(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
725template <
typename T,
typename Storage>
727max(
const OrthogPoly<T,Storage>& a,
730 OrthogPoly<T,Storage> c(a.expansion(), 0);
731 a.expansion()->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
735template <
typename T,
typename Storage>
737min(
const OrthogPoly<T,Storage>& a,
738 const OrthogPoly<T,Storage>& b)
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)
750 OrthogPoly<T,Storage> c(e, 0);
751 e->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
752 b.getOrthogPolyApprox());
757template <
typename T,
typename Storage>
760 const OrthogPoly<T,Storage>& b)
762 OrthogPoly<T,Storage> c(b.expansion(), 0);
763 b.expansion()->min(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
767template <
typename T,
typename Storage>
769min(
const OrthogPoly<T,Storage>& a,
772 OrthogPoly<T,Storage> c(a.expansion(), 0);
773 a.expansion()->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
777template <
typename T,
typename Storage>
780 const OrthogPoly<T,Storage>& b)
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>
792 const OrthogPoly<T,Storage>& b)
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>
818 const OrthogPoly<T,Storage>& b)
823template <
typename T,
typename Storage>
826 const OrthogPoly<T,Storage>& b)
831template <
typename T,
typename Storage>
839template <
typename T,
typename Storage>
842 const OrthogPoly<T,Storage>& b)
844 return a.two_norm() <= b.two_norm();
847template <
typename T,
typename Storage>
850 const OrthogPoly<T,Storage>& b)
852 return a <= b.two_norm();
855template <
typename T,
typename Storage>
860 return a.two_norm() <= b;
863template <
typename T,
typename Storage>
866 const OrthogPoly<T,Storage>& b)
868 return a.two_norm() >= b.two_norm();
871template <
typename T,
typename Storage>
874 const OrthogPoly<T,Storage>& b)
876 return a >= b.two_norm();
879template <
typename T,
typename Storage>
884 return a.two_norm() >= b;
887template <
typename T,
typename Storage>
890 const OrthogPoly<T,Storage>& b)
892 return a.two_norm() < b.two_norm();
895template <
typename T,
typename Storage>
898 const OrthogPoly<T,Storage>& b)
900 return a < b.two_norm();
903template <
typename T,
typename Storage>
908 return a.two_norm() < b;
911template <
typename T,
typename Storage>
914 const OrthogPoly<T,Storage>& b)
916 return a.two_norm() > b.two_norm();
919template <
typename T,
typename Storage>
922 const OrthogPoly<T,Storage>& b)
924 return a > b.two_norm();
927template <
typename T,
typename Storage>
932 return a.two_norm() > b;
935template <
typename T,
typename Storage>
936bool toBool(
const OrthogPoly<T,Storage>& x) {
938 for (
int i=0; i<x.size(); i++)
939 is_zero = is_zero && (x.coeff(i) == 0.0);
943template <
typename T,
typename Storage>
945operator && (
const OrthogPoly<T,Storage>& x1,
const OrthogPoly<T,Storage>& x2)
950template <
typename T,
typename Storage>
953 const OrthogPoly<T,Storage>& x2)
958template <
typename T,
typename Storage>
966template <
typename T,
typename Storage>
968operator || (
const OrthogPoly<T,Storage>& x1,
const OrthogPoly<T,Storage>& x2)
973template <
typename T,
typename Storage>
976 const OrthogPoly<T,Storage>& x2)
981template <
typename T,
typename Storage>
989template <
typename T,
typename Storage>
991operator << (std::ostream& os,
const OrthogPoly<T,Storage>& a)
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)
OrthogPoly< T, Storage > cosh(const 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)
OrthogPoly< T, Storage > tan(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > asin(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > sqrt(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > abs(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > exp(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > cos(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > fabs(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > tanh(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > cbrt(const OrthogPoly< T, Storage > &a)
bool operator!=(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > acos(const OrthogPoly< T, Storage > &a)
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 > log(const OrthogPoly< T, Storage > &a)
bool operator&&(const OrthogPoly< T, Storage > &x1, const OrthogPoly< T, Storage > &x2)
OrthogPoly< T, Storage > log10(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > atan(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > acosh(const OrthogPoly< T, Storage > &a)
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 > sin(const OrthogPoly< T, Storage > &a)
OrthogPoly< T, Storage > pow(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
OrthogPoly< T, Storage > sinh(const OrthogPoly< T, Storage > &a)
bool toBool(const OrthogPoly< T, Storage > &x)
OrthogPoly< T, Storage > atanh(const OrthogPoly< T, Storage > &a)
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 > asinh(const OrthogPoly< T, Storage > &a)
bool operator||(const OrthogPoly< T, Storage > &x1, const OrthogPoly< T, Storage > &x2)