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