Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_PCE_OrthogPolyImp.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#include "Sacado_DynamicArrayTraits.hpp"
43#include "Teuchos_TimeMonitor.hpp"
45
46namespace Sacado {
47namespace PCE {
48
49template <typename T, typename Storage>
51OrthogPoly() :
52 expansion_(),
54{
55 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
56 expansion_ = const_expansion_;
57}
58
59template <typename T, typename Storage>
62 expansion_(),
63 th(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(Teuchos::null, 1, &x))
64{
65 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
66 expansion_ = const_expansion_;
67}
68
69template <typename T, typename Storage>
71OrthogPoly(const Teuchos::RCP<expansion_type>& expansion) :
72 expansion_(expansion),
73 th(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(expansion_->getBasis()))
74{
75 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
76}
77
78template <typename T, typename Storage>
80OrthogPoly(const Teuchos::RCP<expansion_type>& expansion,
81 ordinal_type sz) :
82 expansion_(expansion),
83 th(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(expansion_->getBasis(), sz))
84{
85 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
86}
87
88template <typename T, typename Storage>
90OrthogPoly(const OrthogPoly<T,Storage>& x) :
91 expansion_(x.expansion_),
92 th(x.th)
93{
94 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
95}
96
97template <typename T, typename Storage>
100{
101}
102
103template <typename T, typename Storage>
104void
106reset(const Teuchos::RCP<expansion_type>& expansion)
107{
108 expansion_ = expansion;
109 th->reset(expansion_->getBasis());
110}
111
112template <typename T, typename Storage>
113void
115reset(const Teuchos::RCP<expansion_type>& expansion, ordinal_type sz)
116{
117 expansion_ = expansion;
118 th->reset(expansion_->getBasis(), sz);
119}
120
121template <typename T, typename Storage>
124evaluate(const Teuchos::Array<typename OrthogPoly<T,Storage>::value_type>& point) const
125{
126 return th->evaluate(point);
127}
128
129template <typename T, typename Storage>
133 const Teuchos::Array<typename OrthogPoly<T,Storage>::value_type>& point,
134 const Teuchos::Array<typename OrthogPoly<T,Storage>::value_type>& bvals) const
135{
136 return th->evaluate(point, bvals);
137}
138
139template <typename T, typename Storage>
140bool
142isEqualTo(const OrthogPoly& x) const {
143 typedef IsEqual<value_type> IE;
144 if (x.size() != this->size()) return false;
145 // Allow expansions to be different if their size is 1 and one
146 // of them is a constant expansion
147 if (expansion_ != x.expansion_) {
148 if (x.size() != 1)
149 return false;
150 if ((expansion_ != const_expansion_) &&
151 (x.expansion_ != x.const_expansion_))
152 return false;
153 }
154 bool eq = true;
155 for (int i=0; i<this->size(); i++)
156 eq = eq && IE::eval(x.coeff(i), this->coeff(i));
157 return eq;
158}
159
160template <typename T, typename Storage>
161OrthogPoly<T,Storage>&
164{
165 th.makeOwnCopy();
166 *th = v;
167 return *this;
168}
169
170template <typename T, typename Storage>
171OrthogPoly<T,Storage>&
173operator=(const OrthogPoly<T,Storage>& x)
174{
175 expansion_ = x.expansion_;
176 th = x.th;
177 return *this;
178}
179
180template <typename T, typename Storage>
181OrthogPoly<T,Storage>
183operator+() const
184{
185 return *this;
186}
187
188template <typename T, typename Storage>
189OrthogPoly<T,Storage>
191operator-() const
192{
193 OrthogPoly<T,Storage> x(expansion_);
194 expansion_->unaryMinus(*(x.th), *th);
195 return x;
196}
197
198template <typename T, typename Storage>
199OrthogPoly<T,Storage>&
202{
203 th.makeOwnCopy();
204 expansion_->plusEqual(*th, v);
205 return *this;
206}
207
208template <typename T, typename Storage>
209OrthogPoly<T,Storage>&
212{
213 th.makeOwnCopy();
214 expansion_->minusEqual(*th, v);
215 return *this;
216}
217
218template <typename T, typename Storage>
219OrthogPoly<T,Storage>&
222{
223 th.makeOwnCopy();
224 expansion_->timesEqual(*th, v);
225 return *this;
226}
227
228template <typename T, typename Storage>
229OrthogPoly<T,Storage>&
232{
233 th.makeOwnCopy();
234 expansion_->divideEqual(*th, v);
235 return *this;
236}
237
238template <typename T, typename Storage>
239OrthogPoly<T,Storage>&
241operator+=(const OrthogPoly<T,Storage>& x)
242{
243 th.makeOwnCopy();
244 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
245 if (x.size() > size()) {
246 e = x.expansion();
247 reset(e, size());
248 }
249 e->plusEqual(*th, *x.th);
250 return *this;
251}
252
253template <typename T, typename Storage>
254OrthogPoly<T,Storage>&
256operator-=(const OrthogPoly<T,Storage>& x)
257{
258 th.makeOwnCopy();
259 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
260 if (x.size() > size()) {
261 e = x.expansion();
262 reset(e, size());
263 }
264 e->minusEqual(*th, *x.th);
265 return *this;
266}
267
268template <typename T, typename Storage>
269OrthogPoly<T,Storage>&
271operator*=(const OrthogPoly<T,Storage>& x)
272{
273 th.makeOwnCopy();
274 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
275 if (x.size() > size()) {
276 e = x.expansion();
277 reset(e, size());
278 }
279 e->timesEqual(*th, *x.th);
280 return *this;
281}
282
283template <typename T, typename Storage>
284OrthogPoly<T,Storage>&
286operator/=(const OrthogPoly<T,Storage>& x)
287{
288 th.makeOwnCopy();
289 Teuchos::RCP<typename OrthogPoly<T,Storage>::expansion_type> e = expansion_;
290 if (x.size() > size()) {
291 e = x.expansion();
292 reset(e, size());
293 }
294 e->divideEqual(*th, *x.th);
295 return *this;
296}
297
298template <typename T, typename Storage>
299OrthogPoly<T,Storage>
300operator+(const OrthogPoly<T,Storage>& a,
301 const OrthogPoly<T,Storage>& b)
302{
303 // Get expansion
304 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
309 e = a.expansion();
310 else
311 e = b.expansion();
312
313 OrthogPoly<T,Storage> c(e, 0);
314 e->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
315 b.getOrthogPolyApprox());
316
317 return c;
318}
319
320template <typename T, typename Storage>
321OrthogPoly<T,Storage>
323 const OrthogPoly<T,Storage>& b)
324{
325 OrthogPoly<T,Storage> c(b.expansion(), 0);
326 b.expansion()->plus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
327 return c;
328}
329
330template <typename T, typename Storage>
331OrthogPoly<T,Storage>
332operator+(const OrthogPoly<T,Storage>& a,
333 const typename OrthogPoly<T,Storage>::value_type& b)
334{
335 OrthogPoly<T,Storage> c(a.expansion(), 0);
336 a.expansion()->plus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
337 return c;
338}
339
340template <typename T, typename Storage>
341OrthogPoly<T,Storage>
342operator-(const OrthogPoly<T,Storage>& a,
343 const OrthogPoly<T,Storage>& b)
344{
345 // Get expansion
346 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
351 e = a.expansion();
352 else
353 e = b.expansion();
354
355 OrthogPoly<T,Storage> c(e, 0);
356 e->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
357 b.getOrthogPolyApprox());
358
359 return c;
360}
361
362template <typename T, typename Storage>
363OrthogPoly<T,Storage>
365 const OrthogPoly<T,Storage>& b)
366{
367 OrthogPoly<T,Storage> c(b.expansion(), 0);
368 b.expansion()->minus(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
369 return c;
370}
371
372template <typename T, typename Storage>
373OrthogPoly<T,Storage>
374operator-(const OrthogPoly<T,Storage>& a,
375 const typename OrthogPoly<T,Storage>::value_type& b)
376{
377 OrthogPoly<T,Storage> c(a.expansion(), 0);
378 a.expansion()->minus(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
379 return c;
380}
381
382template <typename T, typename Storage>
383OrthogPoly<T,Storage>
384operator*(const OrthogPoly<T,Storage>& a,
385 const OrthogPoly<T,Storage>& b)
386{
387 // Get expansion
388 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
393 e = a.expansion();
394 else
395 e = b.expansion();
396
397 OrthogPoly<T,Storage> c(e, 0);
398 e->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
399 b.getOrthogPolyApprox());
400
401 return c;
402}
403
404template <typename T, typename Storage>
405OrthogPoly<T,Storage>
407 const OrthogPoly<T,Storage>& b)
408{
409 OrthogPoly<T,Storage> c(b.expansion(), 0);
410 b.expansion()->times(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
411 return c;
412}
413
414template <typename T, typename Storage>
415OrthogPoly<T,Storage>
416operator*(const OrthogPoly<T,Storage>& a,
417 const typename OrthogPoly<T,Storage>::value_type& b)
418{
419 OrthogPoly<T,Storage> c(a.expansion(), 0);
420 a.expansion()->times(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
421 return c;
422}
423
424template <typename T, typename Storage>
425OrthogPoly<T,Storage>
426operator/(const OrthogPoly<T,Storage>& a,
427 const OrthogPoly<T,Storage>& b)
428{
429 // Get expansion
430 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
435 e = a.expansion();
436 else
437 e = b.expansion();
438
439 OrthogPoly<T,Storage> c(e, 0);
440 e->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
441 b.getOrthogPolyApprox());
442
443 return c;
444}
445
446template <typename T, typename Storage>
447OrthogPoly<T,Storage>
449 const OrthogPoly<T,Storage>& b)
450{
451 OrthogPoly<T,Storage> c(b.expansion(), 0);
452 b.expansion()->divide(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
453 return c;
454}
455
456template <typename T, typename Storage>
457OrthogPoly<T,Storage>
458operator/(const OrthogPoly<T,Storage>& a,
459 const typename OrthogPoly<T,Storage>::value_type& b)
460{
461 OrthogPoly<T,Storage> c(a.expansion(), 0);
462 a.expansion()->divide(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
463 return c;
464}
465
466template <typename T, typename Storage>
467OrthogPoly<T,Storage>
468exp(const OrthogPoly<T,Storage>& a)
469{
470 OrthogPoly<T,Storage> c(a.expansion(), 0);
471 a.expansion()->exp(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
472 return c;
473}
474
475template <typename T, typename Storage>
476OrthogPoly<T,Storage>
477log(const OrthogPoly<T,Storage>& a)
478{
479 TEUCHOS_FUNC_TIME_MONITOR("LOG");
480 OrthogPoly<T,Storage> c(a.expansion(), 0);
481 {
482 TEUCHOS_FUNC_TIME_MONITOR("OPA LOG");
483 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
484 }
485
486 OrthogPoly<T,Storage> d(c);
487 return d;
488}
489
490template <typename T, typename Storage>
491void
492log(OrthogPoly<T,Storage>& c, const OrthogPoly<T,Storage>& a)
493{
494 OrthogPoly<T,Storage> d(a.expansion(), 0);
495 a.expansion()->log(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
496}
497
498template <typename T, typename Storage>
499OrthogPoly<T,Storage>
500log10(const OrthogPoly<T,Storage>& a)
501{
502 OrthogPoly<T,Storage> c(a.expansion(), 0);
503 a.expansion()->log10(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
504 return c;
505}
506
507template <typename T, typename Storage>
508OrthogPoly<T,Storage>
509sqrt(const OrthogPoly<T,Storage>& a)
510{
511 OrthogPoly<T,Storage> c(a.expansion(), 0);
512 a.expansion()->sqrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
513 return c;
514}
515
516template <typename T, typename Storage>
517OrthogPoly<T,Storage>
518cbrt(const OrthogPoly<T,Storage>& a)
519{
520 OrthogPoly<T,Storage> c(a.expansion(), 0);
521 a.expansion()->cbrt(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
522 return c;
523}
524
525template <typename T, typename Storage>
526OrthogPoly<T,Storage>
527pow(const OrthogPoly<T,Storage>& a,
528 const OrthogPoly<T,Storage>& b)
529{
530 // Get expansion
531 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
536 e = a.expansion();
537 else
538 e = b.expansion();
539
540 OrthogPoly<T,Storage> c(e, 0);
541 e->pow(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
542 b.getOrthogPolyApprox());
543
544 return c;
545}
546
547template <typename T, typename Storage>
548OrthogPoly<T,Storage>
549pow(const T& a,
550 const OrthogPoly<T,Storage>& b)
551{
552 OrthogPoly<T,Storage> c(b.expansion(), 0);
553 b.expansion()->pow(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
554 return c;
555}
556
557template <typename T, typename Storage>
558OrthogPoly<T,Storage>
559pow(const OrthogPoly<T,Storage>& a,
560 const T& b)
561{
562 OrthogPoly<T,Storage> c(a.expansion(), 0);
563 a.expansion()->pow(c.getOrthogPolyApprox(),a.getOrthogPolyApprox(), b);
564 return c;
565}
566
567template <typename T, typename Storage>
568OrthogPoly<T,Storage>
569sin(const OrthogPoly<T,Storage>& a)
570{
571 OrthogPoly<T,Storage> c(a.expansion(), 0);
572 a.expansion()->sin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
573 return c;
574}
575
576template <typename T, typename Storage>
577OrthogPoly<T,Storage>
578cos(const OrthogPoly<T,Storage>& a)
579{
580 OrthogPoly<T,Storage> c(a.expansion(), 0);
581 a.expansion()->cos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
582 return c;
583}
584
585template <typename T, typename Storage>
586OrthogPoly<T,Storage>
587tan(const OrthogPoly<T,Storage>& a)
588{
589 OrthogPoly<T,Storage> c(a.expansion(), 0);
590 a.expansion()->tan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
591 return c;
592}
593
594template <typename T, typename Storage>
595OrthogPoly<T,Storage>
596sinh(const OrthogPoly<T,Storage>& a)
597{
598 OrthogPoly<T,Storage> c(a.expansion(), 0);
599 a.expansion()->sinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
600 return c;
601}
602
603template <typename T, typename Storage>
604OrthogPoly<T,Storage>
605cosh(const OrthogPoly<T,Storage>& a)
606{
607 OrthogPoly<T,Storage> c(a.expansion(), 0);
608 a.expansion()->cosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
609 return c;
610}
611
612template <typename T, typename Storage>
613OrthogPoly<T,Storage>
614tanh(const OrthogPoly<T,Storage>& a)
615{
616 OrthogPoly<T,Storage> c(a.expansion(), 0);
617 a.expansion()->tanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
618 return c;
619}
620
621template <typename T, typename Storage>
622OrthogPoly<T,Storage>
623acos(const OrthogPoly<T,Storage>& a)
624{
625 OrthogPoly<T,Storage> c(a.expansion(), 0);
626 a.expansion()->acos(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
627 return c;
628}
629
630template <typename T, typename Storage>
631OrthogPoly<T,Storage>
632asin(const OrthogPoly<T,Storage>& a)
633{
634 OrthogPoly<T,Storage> c(a.expansion(), 0);
635 a.expansion()->asin(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
636 return c;
637}
638
639template <typename T, typename Storage>
640OrthogPoly<T,Storage>
641atan(const OrthogPoly<T,Storage>& a)
642{
643 OrthogPoly<T,Storage> c(a.expansion(), 0);
644 a.expansion()->atan(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
645 return c;
646}
647
648template <typename T, typename Storage>
649OrthogPoly<T,Storage>
650acosh(const OrthogPoly<T,Storage>& a)
651{
652 OrthogPoly<T,Storage> c(a.expansion(), 0);
653 a.expansion()->acosh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
654 return c;
655}
656
657template <typename T, typename Storage>
658OrthogPoly<T,Storage>
659asinh(const OrthogPoly<T,Storage>& a)
660{
661 OrthogPoly<T,Storage> c(a.expansion(), 0);
662 a.expansion()->asinh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
663 return c;
664}
665
666template <typename T, typename Storage>
667OrthogPoly<T,Storage>
668atanh(const OrthogPoly<T,Storage>& a)
669{
670 OrthogPoly<T,Storage> c(a.expansion(), 0);
671 a.expansion()->atanh(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
672 return c;
673}
674
675template <typename T, typename Storage>
676OrthogPoly<T,Storage>
677fabs(const OrthogPoly<T,Storage>& a)
678{
679 OrthogPoly<T,Storage> c(a.expansion(), 0);
680 a.expansion()->fabs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
681 return c;
682}
683
684template <typename T, typename Storage>
685OrthogPoly<T,Storage>
686abs(const OrthogPoly<T,Storage>& a)
687{
688 OrthogPoly<T,Storage> c(a.expansion(), 0);
689 a.expansion()->abs(c.getOrthogPolyApprox(), a.getOrthogPolyApprox());
690 return c;
691}
692
693template <typename T, typename Storage>
694OrthogPoly<T,Storage>
695max(const OrthogPoly<T,Storage>& a,
696 const OrthogPoly<T,Storage>& b)
697{
698 // Get expansion
699 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
704 e = a.expansion();
705 else
706 e = b.expansion();
707
708 OrthogPoly<T,Storage> c(e, 0);
709 e->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
710 b.getOrthogPolyApprox());
711
712 return c;
713}
714
715template <typename T, typename Storage>
716OrthogPoly<T,Storage>
718 const OrthogPoly<T,Storage>& b)
719{
720 OrthogPoly<T,Storage> c(b.expansion(), 0);
721 b.expansion()->max(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
722 return c;
723}
724
725template <typename T, typename Storage>
726OrthogPoly<T,Storage>
727max(const OrthogPoly<T,Storage>& a,
728 const typename OrthogPoly<T,Storage>::value_type& b)
729{
730 OrthogPoly<T,Storage> c(a.expansion(), 0);
731 a.expansion()->max(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
732 return c;
733}
734
735template <typename T, typename Storage>
736OrthogPoly<T,Storage>
737min(const OrthogPoly<T,Storage>& a,
738 const OrthogPoly<T,Storage>& b)
739{
740 // Get expansion
741 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
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)
746 e = a.expansion();
747 else
748 e = b.expansion();
749
750 OrthogPoly<T,Storage> c(e, 0);
751 e->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(),
752 b.getOrthogPolyApprox());
753
754 return c;
755}
756
757template <typename T, typename Storage>
758OrthogPoly<T,Storage>
760 const OrthogPoly<T,Storage>& b)
761{
762 OrthogPoly<T,Storage> c(b.expansion(), 0);
763 b.expansion()->min(c.getOrthogPolyApprox(), a, b.getOrthogPolyApprox());
764 return c;
765}
766
767template <typename T, typename Storage>
768OrthogPoly<T,Storage>
769min(const OrthogPoly<T,Storage>& a,
770 const typename OrthogPoly<T,Storage>::value_type& b)
771{
772 OrthogPoly<T,Storage> c(a.expansion(), 0);
773 a.expansion()->min(c.getOrthogPolyApprox(), a.getOrthogPolyApprox(), b);
774 return c;
775}
776
777template <typename T, typename Storage>
778bool
779operator==(const OrthogPoly<T,Storage>& a,
780 const OrthogPoly<T,Storage>& b)
781{
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))
785 return false;
786 return true;
787}
788
789template <typename T, typename Storage>
790bool
792 const OrthogPoly<T,Storage>& b)
793{
794 if (a != b.coeff(0))
795 return false;
796 for (int i=1; i<b.size(); i++)
797 if (b.coeff(i) != T(0.0))
798 return false;
799 return true;
800}
801
802template <typename T, typename Storage>
803bool
804operator==(const OrthogPoly<T,Storage>& a,
805 const typename OrthogPoly<T,Storage>::value_type& b)
806{
807 if (a.coeff(0) != b)
808 return false;
809 for (int i=1; i<a.size(); i++)
810 if (a.coeff(i) != T(0.0))
811 return false;
812 return true;
813}
814
815template <typename T, typename Storage>
816bool
817operator!=(const OrthogPoly<T,Storage>& a,
818 const OrthogPoly<T,Storage>& b)
819{
820 return !(a == b);
821}
822
823template <typename T, typename Storage>
824bool
826 const OrthogPoly<T,Storage>& b)
827{
828 return !(a == b);
829}
830
831template <typename T, typename Storage>
832bool
833operator!=(const OrthogPoly<T,Storage>& a,
834 const typename OrthogPoly<T,Storage>::value_type& b)
835{
836 return !(a == b);
837}
838
839template <typename T, typename Storage>
840bool
841operator<=(const OrthogPoly<T,Storage>& a,
842 const OrthogPoly<T,Storage>& b)
843{
844 return a.two_norm() <= b.two_norm();
845}
846
847template <typename T, typename Storage>
848bool
850 const OrthogPoly<T,Storage>& b)
851{
852 return a <= b.two_norm();
853}
854
855template <typename T, typename Storage>
856bool
857operator<=(const OrthogPoly<T,Storage>& a,
858 const typename OrthogPoly<T,Storage>::value_type& b)
859{
860 return a.two_norm() <= b;
861}
862
863template <typename T, typename Storage>
864bool
865operator>=(const OrthogPoly<T,Storage>& a,
866 const OrthogPoly<T,Storage>& b)
867{
868 return a.two_norm() >= b.two_norm();
869}
870
871template <typename T, typename Storage>
872bool
874 const OrthogPoly<T,Storage>& b)
875{
876 return a >= b.two_norm();
877}
878
879template <typename T, typename Storage>
880bool
881operator>=(const OrthogPoly<T,Storage>& a,
882 const typename OrthogPoly<T,Storage>::value_type& b)
883{
884 return a.two_norm() >= b;
885}
886
887template <typename T, typename Storage>
888bool
889operator<(const OrthogPoly<T,Storage>& a,
890 const OrthogPoly<T,Storage>& b)
891{
892 return a.two_norm() < b.two_norm();
893}
894
895template <typename T, typename Storage>
896bool
898 const OrthogPoly<T,Storage>& b)
899{
900 return a < b.two_norm();
901}
902
903template <typename T, typename Storage>
904bool
905operator<(const OrthogPoly<T,Storage>& a,
906 const typename OrthogPoly<T,Storage>::value_type& b)
907{
908 return a.two_norm() < b;
909}
910
911template <typename T, typename Storage>
912bool
913operator>(const OrthogPoly<T,Storage>& a,
914 const OrthogPoly<T,Storage>& b)
915{
916 return a.two_norm() > b.two_norm();
917}
918
919template <typename T, typename Storage>
920bool
922 const OrthogPoly<T,Storage>& b)
923{
924 return a > b.two_norm();
925}
926
927template <typename T, typename Storage>
928bool
929operator>(const OrthogPoly<T,Storage>& a,
930 const typename OrthogPoly<T,Storage>::value_type& b)
931{
932 return a.two_norm() > b;
933}
934
935template <typename T, typename Storage>
936bool toBool(const OrthogPoly<T,Storage>& x) {
937 bool is_zero = true;
938 for (int i=0; i<x.size(); i++)
939 is_zero = is_zero && (x.coeff(i) == 0.0);
940 return !is_zero;
941}
942
943template <typename T, typename Storage>
944inline bool
945operator && (const OrthogPoly<T,Storage>& x1, const OrthogPoly<T,Storage>& x2)
946{
947 return toBool(x1) && toBool(x2);
948}
949
950template <typename T, typename Storage>
951inline bool
953 const OrthogPoly<T,Storage>& x2)
954{
955 return a && toBool(x2);
956}
957
958template <typename T, typename Storage>
959inline bool
960operator && (const OrthogPoly<T,Storage>& x1,
961 const typename OrthogPoly<T,Storage>::value_type& b)
962{
963 return toBool(x1) && b;
964}
965
966template <typename T, typename Storage>
967inline bool
968operator || (const OrthogPoly<T,Storage>& x1, const OrthogPoly<T,Storage>& x2)
969{
970 return toBool(x1) || toBool(x2);
971}
972
973template <typename T, typename Storage>
974inline bool
976 const OrthogPoly<T,Storage>& x2)
977{
978 return a || toBool(x2);
979}
980
981template <typename T, typename Storage>
982inline bool
983operator || (const OrthogPoly<T,Storage>& x1,
984 const typename OrthogPoly<T,Storage>::value_type& b)
985{
986 return toBool(x1) || b;
987}
988
989template <typename T, typename Storage>
990std::ostream&
991operator << (std::ostream& os, const OrthogPoly<T,Storage>& a)
992{
993 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
994
995 os << "[ ";
996
997 for (ordinal_type i=0; i<a.size(); i++) {
998 os << a.coeff(i) << " ";
999 }
1000
1001 os << "]\n";
1002 return os;
1003}
1004
1005template <typename T, typename Storage>
1006std::istream&
1007operator >> (std::istream& is, OrthogPoly<T,Storage>& a)
1008{
1009 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
1010
1011 // Read in the opening "["
1012 char bracket;
1013 is >> bracket;
1014
1015 for (ordinal_type i=0; i<a.size(); i++) {
1016 is >> a.fastAccessCoeff(i);
1017 }
1018
1019 // Read in the closing "]"
1020
1021 is >> bracket;
1022 return is;
1023}
1024
1025} // namespace PCE
1026} // namespace Sacado
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)