Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_UQ_PCE_Imp.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_SerialDenseMatrix.hpp"
11 
12 
13 namespace Sacado {
14 namespace UQ {
15 
16 /*
17 template <typename Storage>
18 KOKKOS_INLINE_FUNCTION
19 typename PCE<Storage>::value_type
20 PCE<Storage>::
21 evaluate(const Teuchos::Array<typename PCE<Storage>::value_type>& point) const
22 {
23  return s_.evaluate(point);
24 }
25 
26 template <typename Storage>
27 KOKKOS_INLINE_FUNCTION
28 typename PCE<Storage>::value_type
29 PCE<Storage>::
30 evaluate(
31  const Teuchos::Array<typename PCE<Storage>::value_type>& point,
32  const Teuchos::Array<typename PCE<Storage>::value_type>& bvals) const
33 {
34  return s_.evaluate(point, bvals);
35 }
36 */
37 
38 template <typename Storage>
39 KOKKOS_INLINE_FUNCTION
40 typename PCE<Storage>::value_type
41 PCE<Storage>::
42 standard_deviation() const {
43  value_type s = 0.0;
44  const ordinal_type sz = this->size();
45  const_pointer c = this->coeff();
46  for (ordinal_type i=1; i<sz; ++i)
47  s += c[i]*c[i];
48  return std::sqrt(s);
49 }
50 
51 template <typename Storage>
52 KOKKOS_INLINE_FUNCTION
53 typename PCE<Storage>::value_type
54 PCE<Storage>::
55 two_norm_squared() const {
56  value_type s = 0.0;
57  const ordinal_type sz = this->size();
58  const_pointer c = this->coeff();
59  for (ordinal_type i=0; i<sz; ++i)
60  s += c[i]*c[i];
61  return s;
62 }
63 
64 template <typename Storage>
65 KOKKOS_INLINE_FUNCTION
66 typename PCE<Storage>::value_type
67 PCE<Storage>::
68 inner_product(const PCE& x) const {
69  value_type s = 0.0;
70  const ordinal_type sz = this->size();
71  const ordinal_type xsz = x.size();
72  const ordinal_type n = sz < xsz ? sz : xsz;
73  const_pointer c = this->coeff();
74  const_pointer xc = x.coeff();
75  for (ordinal_type i=0; i<n; ++i)
76  s += c[i]*xc[i];
77  return s;
78 }
79 
80 template <typename Storage>
81 KOKKOS_INLINE_FUNCTION
82 bool
83 PCE<Storage>::
84 isEqualTo(const PCE& x) const {
85  typedef IsEqual<value_type> IE;
86  const ordinal_type sz = this->size();
87  if (x.size() != sz) return false;
88  bool eq = true;
89  for (ordinal_type i=0; i<sz; i++)
90  eq = eq && IE::eval(x.coeff(i), this->coeff(i));
91  return eq;
92 }
93 
94 template <typename Storage>
95 KOKKOS_INLINE_FUNCTION
96 bool
97 PCE<Storage>::
98 isEqualTo(const PCE& x) const volatile {
99  typedef IsEqual<value_type> IE;
100  const ordinal_type sz = this->size();
101  if (x.size() != sz) return false;
102  bool eq = true;
103  for (ordinal_type i=0; i<sz; i++)
104  eq = eq && IE::eval(x.coeff(i), this->coeff(i));
105  return eq;
106 }
107 
108 template <typename Storage>
109 KOKKOS_INLINE_FUNCTION
110 PCE<Storage>&
111 PCE<Storage>::
112 operator=(const typename PCE<Storage>::value_type v)
113 {
114  s_.init(value_type(0));
115  s_[0] = v;
116  return *this;
117 }
118 
119 template <typename Storage>
120 KOKKOS_INLINE_FUNCTION
121 /*volatile*/ PCE<Storage>&
122 PCE<Storage>::
123 operator=(const typename PCE<Storage>::value_type v) volatile
124 {
125  s_.init(value_type(0));
126  s_[0] = v;
127  return const_cast<PCE&>(*this);
128 }
129 
130 template <typename Storage>
131 KOKKOS_INLINE_FUNCTION
132 PCE<Storage>&
133 PCE<Storage>::
134 operator=(const PCE<Storage>& x)
135 {
136  if (this != &x) {
137  if (!s_.is_view())
138  cijk_ = x.cijk_;
139  if (!s_.is_view() && is_constant(x)) {
140  s_.resize(1);
141  s_[0] = x.s_[0];
142  }
143  else
144  s_ = x.s_;
145 
146  // For DyamicStorage as a view (is_owned=false), we need to set
147  // the trailing entries when assigning a constant vector (because
148  // the copy constructor in this case doesn't reset the size of this)
149  if (s_.size() > x.s_.size())
150  for (ordinal_type i=x.s_.size(); i<s_.size(); i++)
151  s_[i] = value_type(0);
152  }
153  return *this;
154 }
155 
156 template <typename Storage>
157 KOKKOS_INLINE_FUNCTION
158 PCE<Storage>&
159 PCE<Storage>::
160 operator=(const volatile PCE<Storage>& x)
161 {
162  if (this != &x) {
163  if (!s_.is_view())
164  cijk_ = const_cast<const my_cijk_type&>(x.cijk_);
165  if (!s_.is_view() && is_constant(x)) {
166  s_.resize(1);
167  s_[0] = x.s_[0];
168  }
169  else
170  s_ = x.s_;
171 
172  // For DyamicStorage as a view (is_owned=false), we need to set
173  // the trailing entries when assigning a constant vector (because
174  // the copy constructor in this case doesn't reset the size of this)
175  if (s_.size() > x.s_.size())
176  for (ordinal_type i=x.s_.size(); i<s_.size(); i++)
177  s_[i] = value_type(0);
178  }
179  return *this;
180 }
181 
182 template <typename Storage>
183 KOKKOS_INLINE_FUNCTION
184 /*volatile*/ PCE<Storage>&
185 PCE<Storage>::
186 operator=(const PCE<Storage>& x) volatile
187 {
188  if (this != &x) {
189  if (!s_.is_view())
190  const_cast<my_cijk_type&>(cijk_) = x.cijk_;
191  if (!s_.is_view() && is_constant(x)) {
192  s_.resize(1);
193  s_[0] = x.s_[0];
194  }
195  else
196  s_ = x.s_;
197 
198  // For DyamicStorage as a view (is_owned=false), we need to set
199  // the trailing entries when assigning a constant vector (because
200  // the copy constructor in this case doesn't reset the size of this)
201  if (s_.size() > x.s_.size())
202  for (ordinal_type i=x.s_.size(); i<s_.size(); i++)
203  s_[i] = value_type(0);
204  }
205  return const_cast<PCE&>(*this);
206 }
207 
208 template <typename Storage>
209 KOKKOS_INLINE_FUNCTION
210 /*volatile*/ PCE<Storage>&
211 PCE<Storage>::
212 operator=(const volatile PCE<Storage>& x) volatile
213 {
214  if (this != &x) {
215  if (!s_.is_view())
216  const_cast<my_cijk_type&>(cijk_) =
217  const_cast<const my_cijk_type&>(x.cijk_);
218  if (!s_.is_view() && is_constant(x)) {
219  s_.resize(1);
220  s_[0] = x.s_[0];
221  }
222  else
223  s_ = x.s_;
224 
225  // For DyamicStorage as a view (is_owned=false), we need to set
226  // the trailing entries when assigning a constant vector (because
227  // the copy constructor in this case doesn't reset the size of this)
228  if (s_.size() > x.s_.size())
229  for (ordinal_type i=x.s_.size(); i<s_.size(); i++)
230  s_[i] = value_type(0);
231  }
232  return const_cast<PCE&>(*this);
233 }
234 
235 template <typename Storage>
236 KOKKOS_INLINE_FUNCTION
237 PCE<Storage>
238 PCE<Storage>::
239 operator-() const
240 {
241  const ordinal_type sz = this->size();
242  PCE<Storage> x(cijk_, sz);
243  pointer xc = x.coeff();
244  const_pointer cc = this->coeff();
245  for (ordinal_type i=0; i<sz; ++i)
246  xc[i] = -cc[i];
247  return x;
248 }
249 
250 template <typename Storage>
251 KOKKOS_INLINE_FUNCTION
252 PCE<Storage>
253 PCE<Storage>::
254 operator-() const volatile
255 {
256  const ordinal_type sz = this->size();
257  PCE<Storage> x(const_cast<const my_cijk_type&>(cijk_), sz);
258  pointer xc = x.coeff();
259  const_volatile_pointer cc = this->coeff();
260  for (ordinal_type i=0; i<sz; ++i)
261  xc[i] = -cc[i];
262  return x;
263 }
264 
265 template <typename Storage>
266 KOKKOS_INLINE_FUNCTION
267 PCE<Storage>&
268 PCE<Storage>::
269 operator*=(const typename PCE<Storage>::value_type v)
270 {
271  pointer cc = this->coeff();
272  const ordinal_type sz = this->size();
273  for (ordinal_type i=0; i<sz; ++i)
274  cc[i] *= v;
275  return *this;
276 }
277 
278 template <typename Storage>
279 KOKKOS_INLINE_FUNCTION
280 /*volatile*/ PCE<Storage>&
281 PCE<Storage>::
282 operator*=(const typename PCE<Storage>::value_type v) volatile
283 {
284  volatile_pointer cc = this->coeff();
285  const ordinal_type sz = this->size();
286  for (ordinal_type i=0; i<sz; ++i)
287  cc[i] *= v;
288  return const_cast<PCE&>(*this);
289 }
290 
291 template <typename Storage>
292 KOKKOS_INLINE_FUNCTION
293 PCE<Storage>&
294 PCE<Storage>::
295 operator/=(const typename PCE<Storage>::value_type v)
296 {
297  pointer cc = this->coeff();
298  const ordinal_type sz = this->size();
299  for (ordinal_type i=0; i<sz; ++i)
300  cc[i] /= v;
301  return *this;
302 }
303 
304 template <typename Storage>
305 KOKKOS_INLINE_FUNCTION
306 /*volatile*/ PCE<Storage>&
307 PCE<Storage>::
308 operator/=(const typename PCE<Storage>::value_type v) volatile
309 {
310  volatile pointer cc = this->coeff();
311  const ordinal_type sz = this->size();
312  for (ordinal_type i=0; i<sz; ++i)
313  cc[i] /= v;
314  return const_cast<PCE&>(*this);
315 }
316 
317 template <typename Storage>
318 KOKKOS_INLINE_FUNCTION
319 PCE<Storage>&
320 PCE<Storage>::
321 operator+=(const PCE<Storage>& x)
322 {
323  const ordinal_type xsz = x.size();
324  const ordinal_type sz = this->size();
325  if (xsz > sz) {
326  this->reset(x.cijk_, xsz);
327  }
328  const_pointer xc = x.coeff();
329  pointer cc = this->coeff();
330  for (ordinal_type i=0; i<xsz; ++i)
331  cc[i] += xc[i];
332  return *this;
333 }
334 
335 template <typename Storage>
336 KOKKOS_INLINE_FUNCTION
337 PCE<Storage>&
338 PCE<Storage>::
339 operator-=(const PCE<Storage>& x)
340 {
341  const ordinal_type xsz = x.size();
342  const ordinal_type sz = this->size();
343  if (xsz > sz) {
344  this->reset(x.cijk_, xsz);
345  }
346  const_pointer xc = x.coeff();
347  pointer cc = this->coeff();
348  for (ordinal_type i=0; i<xsz; ++i)
349  cc[i] -= xc[i];
350  return *this;
351 }
352 
353 template <typename Storage>
354 KOKKOS_INLINE_FUNCTION
355 PCE<Storage>&
356 PCE<Storage>::
357 operator*=(const PCE<Storage>& x)
358 {
359  const ordinal_type xsz = x.size();
360  const ordinal_type sz = this->size();
361  const ordinal_type csz = sz > xsz ? sz : xsz;
362 
363 
364  KOKKOS_IF_ON_HOST(
365  TEUCHOS_TEST_FOR_EXCEPTION(
366  sz != xsz && sz != 1 && xsz != 1, std::logic_error,
367  "Sacado::UQ::PCE::operator*=(): input sizes do not match");
368  )
369 
370  if (cijk_.is_empty() && !x.cijk_.is_empty())
371  cijk_ = x.cijk_;
372 
373  KOKKOS_IF_ON_HOST(
374  TEUCHOS_TEST_FOR_EXCEPTION(
375  cijk_.is_empty() && csz != 1, std::logic_error,
376  "Sacado::UQ::PCE::operator*(): empty cijk but expansion size > 1");
377  )
378 
379  if (csz > sz)
380  s_.resize(csz);
381 
382  const_pointer xc = x.coeff();
383  pointer cc = this->coeff();
384  if (xsz == 1) {
385  const value_type xcz = xc[0];
386  for (ordinal_type i=0; i<sz; ++i)
387  cc[i] *= xcz;
388  }
389  else if (sz == 1) {
390  const value_type ccz = cc[0];
391  for (ordinal_type i=0; i<xsz; ++i)
392  cc[i] = ccz*xc[i];
393  }
394  else {
395  PCE<Storage> y(cijk_, csz);
396  pointer yc = y.coeff();
397  for (ordinal_type i=0; i<csz; ++i) {
398  const cijk_size_type num_entry = cijk_.num_entry(i);
399  const cijk_size_type entry_beg = cijk_.entry_begin(i);
400  const cijk_size_type entry_end = entry_beg + num_entry;
401  value_type ytmp = 0;
402  for (cijk_size_type entry = entry_beg; entry < entry_end; ++entry) {
403  const cijk_size_type j = cijk_.coord(entry,0);
404  const cijk_size_type k = cijk_.coord(entry,1);
405  ytmp += cijk_.value(entry) * ( cc[j] * xc[k] + cc[k] * xc[j] );
406  }
407  yc[i] = ytmp;
408  }
409  s_ = y.s_;
410  }
411  return *this;
412 }
413 
414 template <typename Storage>
415 KOKKOS_INLINE_FUNCTION
416 PCE<Storage>&
417 PCE<Storage>::
418 operator/=(const PCE<Storage>& x)
419 {
420  const ordinal_type xsz = x.size();
421  const ordinal_type sz = this->size();
422  const ordinal_type csz = sz > xsz ? sz : xsz;
423 
424  KOKKOS_IF_ON_HOST(
425  TEUCHOS_TEST_FOR_EXCEPTION(
426  sz != xsz && sz != 1 && xsz != 1, std::logic_error,
427  "Sacado::UQ::PCE::operator/=(): input sizes do not match");
428  )
429 
430  if (cijk_.is_empty() && !x.cijk_.is_empty())
431  cijk_ = x.cijk_;
432 
433  if (csz > sz)
434  s_.resize(csz);
435 
436  const_pointer xc = x.coeff();
437  pointer cc = this->coeff();
438 
439  KOKKOS_IF_ON_DEVICE(
440  const value_type xcz = xc[0];
441  for (ordinal_type i=0; i<sz; ++i)
442  cc[i] /= xcz;
443  )
444 
445  KOKKOS_IF_ON_HOST(
446  if (xsz == 1) {//constant denom
447  const value_type xcz = xc[0];
448  for (ordinal_type i=0; i<sz; ++i)
449  cc[i] /= xcz;
450  }
451  else {
452 
453  PCE<Storage> y(cijk_, csz);
454  CG_divide(*this, x, y);
455  s_ = y.s_;
456  }
457  )
458 
459  return *this;
460 
461 }
462 
463 template <typename Storage>
464 KOKKOS_INLINE_FUNCTION
465 PCE<Storage>
466 operator+(const PCE<Storage>& a, const PCE<Storage>& b)
467 {
468  typedef typename PCE<Storage>::my_cijk_type my_cijk_type;
469  typedef typename PCE<Storage>::pointer pointer;
470  typedef typename PCE<Storage>::const_pointer const_pointer;
471  typedef typename PCE<Storage>::ordinal_type ordinal_type;
472 
473  const ordinal_type asz = a.size();
474  const ordinal_type bsz = b.size();
475  const ordinal_type csz = asz > bsz ? asz : bsz;
476  my_cijk_type c_cijk = asz > bsz ? a.cijk() : b.cijk();
477 
478  PCE<Storage> c(c_cijk, csz);
479  const_pointer ac = a.coeff();
480  const_pointer bc = b.coeff();
481  pointer cc = c.coeff();
482  if (asz > bsz) {
483  for (ordinal_type i=0; i<bsz; ++i)
484  cc[i] = ac[i] + bc[i];
485  for (ordinal_type i=bsz; i<asz; ++i)
486  cc[i] = ac[i];
487  }
488  else {
489  for (ordinal_type i=0; i<asz; ++i)
490  cc[i] = ac[i] + bc[i];
491  for (ordinal_type i=asz; i<bsz; ++i)
492  cc[i] = bc[i];
493  }
494 
495  return c;
496 }
497 
498 template <typename Storage>
499 KOKKOS_INLINE_FUNCTION
500 PCE<Storage>
501 operator+(const typename PCE<Storage>::value_type& a,
502  const PCE<Storage>& b)
503 {
504  typedef typename PCE<Storage>::pointer pointer;
505  typedef typename PCE<Storage>::const_pointer const_pointer;
506  typedef typename PCE<Storage>::ordinal_type ordinal_type;
507 
508  const ordinal_type bsz = b.size();
509  PCE<Storage> c(b.cijk(), bsz);
510  const_pointer bc = b.coeff();
511  pointer cc = c.coeff();
512  cc[0] = a + bc[0];
513  for (ordinal_type i=1; i<bsz; ++i)
514  cc[i] = bc[i];
515  return c;
516 }
517 
518 template <typename Storage>
519 KOKKOS_INLINE_FUNCTION
520 PCE<Storage>
521 operator+(const PCE<Storage>& a,
522  const typename PCE<Storage>::value_type& b)
523 {
524  typedef typename PCE<Storage>::pointer pointer;
525  typedef typename PCE<Storage>::const_pointer const_pointer;
526  typedef typename PCE<Storage>::ordinal_type ordinal_type;
527 
528  const ordinal_type asz = a.size();
529  PCE<Storage> c(a.cijk(), asz);
530  const_pointer ac = a.coeff();
531  pointer cc = c.coeff();
532  cc[0] = ac[0] + b;
533  for (ordinal_type i=1; i<asz; ++i)
534  cc[i] = ac[i];
535  return c;
536 }
537 
538 template <typename Storage>
539 KOKKOS_INLINE_FUNCTION
540 PCE<Storage>
541 operator-(const PCE<Storage>& a, const PCE<Storage>& b)
542 {
543  typedef typename PCE<Storage>::my_cijk_type my_cijk_type;
544  typedef typename PCE<Storage>::pointer pointer;
545  typedef typename PCE<Storage>::const_pointer const_pointer;
546  typedef typename PCE<Storage>::ordinal_type ordinal_type;
547 
548  const ordinal_type asz = a.size();
549  const ordinal_type bsz = b.size();
550  const ordinal_type csz = asz > bsz ? asz : bsz;
551  my_cijk_type c_cijk = asz > bsz ? a.cijk() : b.cijk();
552 
553  PCE<Storage> c(c_cijk, csz);
554  const_pointer ac = a.coeff();
555  const_pointer bc = b.coeff();
556  pointer cc = c.coeff();
557  if (asz > bsz) {
558  for (ordinal_type i=0; i<bsz; ++i)
559  cc[i] = ac[i] - bc[i];
560  for (ordinal_type i=bsz; i<asz; ++i)
561  cc[i] = ac[i];
562  }
563  else {
564  for (ordinal_type i=0; i<asz; ++i)
565  cc[i] = ac[i] - bc[i];
566  for (ordinal_type i=asz; i<bsz; ++i)
567  cc[i] = -bc[i];
568  }
569 
570  return c;
571 }
572 
573 template <typename Storage>
574 KOKKOS_INLINE_FUNCTION
575 PCE<Storage>
576 operator-(const typename PCE<Storage>::value_type& a,
577  const PCE<Storage>& b)
578 {
579  typedef typename PCE<Storage>::pointer pointer;
580  typedef typename PCE<Storage>::const_pointer const_pointer;
581  typedef typename PCE<Storage>::ordinal_type ordinal_type;
582 
583  const ordinal_type bsz = b.size();
584  PCE<Storage> c(b.cijk(), bsz);
585  const_pointer bc = b.coeff();
586  pointer cc = c.coeff();
587  cc[0] = a - bc[0];
588  for (ordinal_type i=1; i<bsz; ++i)
589  cc[i] = -bc[i];
590  return c;
591 }
592 
593 template <typename Storage>
594 KOKKOS_INLINE_FUNCTION
595 PCE<Storage>
596 operator-(const PCE<Storage>& a,
597  const typename PCE<Storage>::value_type& b)
598 {
599  typedef typename PCE<Storage>::pointer pointer;
600  typedef typename PCE<Storage>::const_pointer const_pointer;
601  typedef typename PCE<Storage>::ordinal_type ordinal_type;
602 
603  const ordinal_type asz = a.size();
604  PCE<Storage> c(a.cijk(), asz);
605  const_pointer ac = a.coeff();
606  pointer cc = c.coeff();
607  cc[0] = ac[0] - b;
608  for (ordinal_type i=1; i<asz; ++i)
609  cc[i] = ac[i];
610  return c;
611 }
612 
613 template <typename Storage>
614 KOKKOS_INLINE_FUNCTION
615 PCE<Storage>
616 operator*(const PCE<Storage>& a, const PCE<Storage>& b)
617 {
618  typedef typename PCE<Storage>::my_cijk_type my_cijk_type;
619  typedef typename PCE<Storage>::pointer pointer;
620  typedef typename PCE<Storage>::const_pointer const_pointer;
621  typedef typename PCE<Storage>::ordinal_type ordinal_type;
622  typedef typename PCE<Storage>::value_type value_type;
623  typedef typename my_cijk_type::size_type cijk_size_type;
624 
625  const ordinal_type asz = a.size();
626  const ordinal_type bsz = b.size();
627  const ordinal_type csz = asz > bsz ? asz : bsz;
628 
629  KOKKOS_IF_ON_HOST((
630  TEUCHOS_TEST_FOR_EXCEPTION(
631  asz != bsz && asz != 1 && bsz != 1, std::logic_error,
632  "Sacado::UQ::PCE::operator*(): input sizes do not match");
633  ))
634 
635  my_cijk_type c_cijk = a.cijk().is_empty() ? b.cijk() : a.cijk();
636 
637  KOKKOS_IF_ON_HOST((
638  TEUCHOS_TEST_FOR_EXCEPTION(
639  c_cijk.is_empty() && csz != 1, std::logic_error,
640  "Sacado::UQ::PCE::operator*(): empty cijk but expansion size > 1");
641  ))
642 
643  PCE<Storage> c(c_cijk, csz);
644  const_pointer ac = a.coeff();
645  const_pointer bc = b.coeff();
646  pointer cc = c.coeff();
647 
648  if (asz == 1) {
649  const value_type acz = ac[0];
650  for (ordinal_type i=0; i<csz; ++i)
651  cc[i] = acz * bc[i];
652  }
653  else if (bsz == 1) {
654  const value_type bcz = bc[0];
655  for (ordinal_type i=0; i<csz; ++i)
656  cc[i] = ac[i] * bcz;
657  }
658  else {
659  for (ordinal_type i=0; i<csz; ++i) {
660  const cijk_size_type num_entry = c_cijk.num_entry(i);
661  const cijk_size_type entry_beg = c_cijk.entry_begin(i);
662  const cijk_size_type entry_end = entry_beg + num_entry;
663  value_type ytmp = 0;
664  for (cijk_size_type entry = entry_beg; entry < entry_end; ++entry) {
665  const cijk_size_type j = c_cijk.coord(entry,0);
666  const cijk_size_type k = c_cijk.coord(entry,1);
667  ytmp += c_cijk.value(entry) * ( ac[j] * bc[k] + ac[k] * bc[j] );
668  }
669  cc[i] = ytmp;
670  }
671  }
672 
673  return c;
674 }
675 
676 template <typename Storage>
677 KOKKOS_INLINE_FUNCTION
678 PCE<Storage>
679 operator*(const typename PCE<Storage>::value_type& a,
680  const PCE<Storage>& b)
681 {
682  typedef typename PCE<Storage>::pointer pointer;
683  typedef typename PCE<Storage>::const_pointer const_pointer;
684  typedef typename PCE<Storage>::ordinal_type ordinal_type;
685 
686  const ordinal_type bsz = b.size();
687  PCE<Storage> c(b.cijk(), bsz);
688  const_pointer bc = b.coeff();
689  pointer cc = c.coeff();
690  for (ordinal_type i=0; i<bsz; ++i)
691  cc[i] = a*bc[i];
692  return c;
693 }
694 
695 template <typename Storage>
696 KOKKOS_INLINE_FUNCTION
697 PCE<Storage>
698 operator*(const PCE<Storage>& a,
699  const typename PCE<Storage>::value_type& b)
700 {
701  typedef typename PCE<Storage>::pointer pointer;
702  typedef typename PCE<Storage>::const_pointer const_pointer;
703  typedef typename PCE<Storage>::ordinal_type ordinal_type;
704 
705  const ordinal_type asz = a.size();
706  PCE<Storage> c(a.cijk(), asz);
707  const_pointer ac = a.coeff();
708  pointer cc = c.coeff();
709  for (ordinal_type i=0; i<asz; ++i)
710  cc[i] = ac[i]*b;
711  return c;
712 }
713 
714 template <typename Storage>
715 KOKKOS_INLINE_FUNCTION
716 PCE<Storage>
717 operator/(const PCE<Storage>& a, const PCE<Storage>& b)
718 {
719  typedef typename PCE<Storage>::pointer pointer;
720  typedef typename PCE<Storage>::const_pointer const_pointer;
721  typedef typename PCE<Storage>::ordinal_type ordinal_type;
722  typedef typename PCE<Storage>::value_type value_type;
723  typedef typename PCE<Storage>::my_cijk_type my_cijk_type;
724 
725  const ordinal_type asz = a.size();
726  const ordinal_type bsz = b.size();
727  const ordinal_type csz = asz > bsz ? asz : bsz;
728 
729  KOKKOS_IF_ON_HOST((
730  TEUCHOS_TEST_FOR_EXCEPTION(
731  asz != bsz && asz != 1 && bsz != 1, std::logic_error,
732  "Sacado::UQ::PCE::operator/(): input sizes do not match");
733  ))
734  my_cijk_type c_cijk = asz == bsz || asz >1 ? a.cijk() : b.cijk();
735 
736  PCE<Storage> c(c_cijk, csz);
737 
738  KOKKOS_IF_ON_HOST((
739  const_pointer ac = a.coeff();
740  pointer cc = c.coeff();
741  value_type bcz = b.fastAccessCoeff(0);
742  for (ordinal_type i=0; i<asz; ++i)
743  cc[i] = ac[i]/bcz;
744  ))
745 
746  KOKKOS_IF_ON_HOST((
747  if (bsz == 1) {//constant denom
748  const_pointer ac = a.coeff();
749  const_pointer bc = b.coeff();
750  pointer cc = c.coeff();
751  const value_type bcz = bc[0];
752  for (ordinal_type i=0; i<csz; ++i)
753  cc[i] = ac[i] / bcz;
754  }
755  else {
756  CG_divide(a,b,c);
757  }
758  ))
759 
760  return c;
761 }
762 
763 template <typename Storage>
764 KOKKOS_INLINE_FUNCTION
765 PCE<Storage>
766 operator/(const typename PCE<Storage>::value_type& a,
767  const PCE<Storage>& b)
768 {
769  //Creat a 0-th order PCE for a
770  PCE<Storage> a_pce(a);
771  return operator/(a_pce,b);
772 }
773 
774 template <typename Storage>
775 KOKKOS_INLINE_FUNCTION
776 PCE<Storage>
777 operator/(const PCE<Storage>& a,
778  const typename PCE<Storage>::value_type& b)
779 {
780  typedef typename PCE<Storage>::pointer pointer;
781  typedef typename PCE<Storage>::const_pointer const_pointer;
782  typedef typename PCE<Storage>::ordinal_type ordinal_type;
783 
784  const ordinal_type asz = a.size();
785  PCE<Storage> c(a.cijk(), asz);
786  const_pointer ac = a.coeff();
787  pointer cc = c.coeff();
788  for (ordinal_type i=0; i<asz; ++i)
789  cc[i] = ac[i]/b;
790  return c;
791 }
792 
793 template <typename Storage>
794 KOKKOS_INLINE_FUNCTION
795 PCE<Storage>
796 exp(const PCE<Storage>& a)
797 {
798  KOKKOS_IF_ON_HOST(
799  TEUCHOS_TEST_FOR_EXCEPTION(
800  a.size() != 1, std::logic_error,
801  "Sacado::UQ::PCE::exp(): argument has size != 1");
802  )
803 
804  PCE<Storage> c(a.cijk(), 1);
805  c.fastAccessCoeff(0) = std::exp( a.fastAccessCoeff(0) );
806 
807  return c;
808 }
809 
810 template <typename Storage>
811 KOKKOS_INLINE_FUNCTION
812 PCE<Storage>
813 log(const PCE<Storage>& a)
814 {
815  KOKKOS_IF_ON_HOST(
816  TEUCHOS_TEST_FOR_EXCEPTION(
817  a.size() != 1, std::logic_error,
818  "Sacado::UQ::PCE::log(): argument has size != 1");
819  )
820 
821  PCE<Storage> c(a.cijk(), 1);
822  c.fastAccessCoeff(0) = std::log( a.fastAccessCoeff(0) );
823 
824  return c;
825 }
826 
827 template <typename Storage>
828 KOKKOS_INLINE_FUNCTION
829 PCE<Storage>
830 log10(const PCE<Storage>& a)
831 {
832  KOKKOS_IF_ON_HOST(
833  TEUCHOS_TEST_FOR_EXCEPTION(
834  a.size() != 1, std::logic_error,
835  "Sacado::UQ::PCE::log10(): argument has size != 1");
836  )
837 
838  PCE<Storage> c(a.cijk(), 1);
839  c.fastAccessCoeff(0) = std::log10( a.fastAccessCoeff(0) );
840 
841  return c;
842 }
843 
844 template <typename Storage>
845 KOKKOS_INLINE_FUNCTION
846 PCE<Storage>
847 sqrt(const PCE<Storage>& a)
848 {
849  KOKKOS_IF_ON_HOST(
850  TEUCHOS_TEST_FOR_EXCEPTION(
851  a.size() != 1, std::logic_error,
852  "Sacado::UQ::PCE::sqrt(): argument has size != 1");
853  )
854 
855  PCE<Storage> c(a.cijk(), 1);
856  c.fastAccessCoeff(0) = std::sqrt( a.fastAccessCoeff(0) );
857 
858  return c;
859 }
860 
861 template <typename Storage>
862 KOKKOS_INLINE_FUNCTION
863 PCE<Storage>
864 cbrt(const PCE<Storage>& a)
865 {
866  KOKKOS_IF_ON_HOST(
867  TEUCHOS_TEST_FOR_EXCEPTION(
868  a.size() != 1, std::logic_error,
869  "Sacado::UQ::PCE::cbrt(): argument has size != 1");
870  )
871 
872  PCE<Storage> c(a.cijk(), 1);
873  c.fastAccessCoeff(0) = std::cbrt( a.fastAccessCoeff(0) );
874 
875  return c;
876 }
877 
878 template <typename Storage>
879 KOKKOS_INLINE_FUNCTION
880 PCE<Storage>
881 pow(const PCE<Storage>& a, const PCE<Storage>& b)
882 {
883  KOKKOS_IF_ON_HOST(
884  TEUCHOS_TEST_FOR_EXCEPTION(
885  a.size() != 1 || b.size() != 1, std::logic_error,
886  "Sacado::UQ::PCE::pow(): arguments have size != 1");
887  )
888 
889  PCE<Storage> c(a.cijk(), 1);
890  c.fastAccessCoeff(0) = std::pow(a.fastAccessCoeff(0), b.fastAccessCoeff(0));
891 
892  return c;
893 }
894 
895 template <typename Storage>
896 KOKKOS_INLINE_FUNCTION
897 PCE<Storage>
898 pow(const typename PCE<Storage>::value_type& a,
899  const PCE<Storage>& b)
900 {
901  KOKKOS_IF_ON_HOST(
902  TEUCHOS_TEST_FOR_EXCEPTION(
903  b.size() != 1, std::logic_error,
904  "Sacado::UQ::PCE::pow(): arguments have size != 1");
905  )
906 
907  PCE<Storage> c(b.cijk(), 1);
908  c.fastAccessCoeff(0) = std::pow(a, b.fastAccessCoeff(0));
909 
910  return c;
911 }
912 
913 template <typename Storage>
914 KOKKOS_INLINE_FUNCTION
915 PCE<Storage>
916 pow(const PCE<Storage>& a,
917  const typename PCE<Storage>::value_type& b)
918 {
919  KOKKOS_IF_ON_HOST(
920  TEUCHOS_TEST_FOR_EXCEPTION(
921  a.size() != 1, std::logic_error,
922  "Sacado::UQ::PCE::pow(): arguments have size != 1");
923  )
924 
925  PCE<Storage> c(a.cijk(), 1);
926  c.fastAccessCoeff(0) = std::pow(a.fastAccessCoeff(0), b);
927 
928  return c;
929 }
930 
931 template <typename Storage>
932 KOKKOS_INLINE_FUNCTION
933 PCE<Storage>
934 sin(const PCE<Storage>& a)
935 {
936  KOKKOS_IF_ON_HOST(
937  TEUCHOS_TEST_FOR_EXCEPTION(
938  a.size() != 1, std::logic_error,
939  "Sacado::UQ::PCE::sin(): argument has size != 1");
940  )
941 
942  PCE<Storage> c(a.cijk(), 1);
943  c.fastAccessCoeff(0) = std::sin( a.fastAccessCoeff(0) );
944 
945  return c;
946 }
947 
948 template <typename Storage>
949 KOKKOS_INLINE_FUNCTION
950 PCE<Storage>
951 cos(const PCE<Storage>& a)
952 {
953  KOKKOS_IF_ON_HOST(
954  TEUCHOS_TEST_FOR_EXCEPTION(
955  a.size() != 1, std::logic_error,
956  "Sacado::UQ::PCE::cos(): argument has size != 1");
957  )
958 
959  PCE<Storage> c(a.cijk(), 1);
960  c.fastAccessCoeff(0) = std::cos( a.fastAccessCoeff(0) );
961 
962  return c;
963 }
964 
965 template <typename Storage>
966 KOKKOS_INLINE_FUNCTION
967 PCE<Storage>
968 tan(const PCE<Storage>& a)
969 {
970  KOKKOS_IF_ON_HOST(
971  TEUCHOS_TEST_FOR_EXCEPTION(
972  a.size() != 1, std::logic_error,
973  "Sacado::UQ::PCE::tan(): argument has size != 1");
974  )
975 
976  PCE<Storage> c(a.cijk(), 1);
977  c.fastAccessCoeff(0) = std::tan( a.fastAccessCoeff(0) );
978 
979  return c;
980 }
981 
982 template <typename Storage>
983 KOKKOS_INLINE_FUNCTION
984 PCE<Storage>
985 sinh(const PCE<Storage>& a)
986 {
987  KOKKOS_IF_ON_HOST(
988  TEUCHOS_TEST_FOR_EXCEPTION(
989  a.size() != 1, std::logic_error,
990  "Sacado::UQ::PCE::sinh(): argument has size != 1");
991  )
992 
993  PCE<Storage> c(a.cijk(), 1);
994  c.fastAccessCoeff(0) = std::sinh( a.fastAccessCoeff(0) );
995 
996  return c;
997 }
998 
999 template <typename Storage>
1000 KOKKOS_INLINE_FUNCTION
1001 PCE<Storage>
1002 cosh(const PCE<Storage>& a)
1003 {
1004  KOKKOS_IF_ON_HOST(
1005  TEUCHOS_TEST_FOR_EXCEPTION(
1006  a.size() != 1, std::logic_error,
1007  "Sacado::UQ::PCE::cosh(): argument has size != 1");
1008  )
1009 
1010  PCE<Storage> c(a.cijk(), 1);
1011  c.fastAccessCoeff(0) = std::cosh( a.fastAccessCoeff(0) );
1012 
1013  return c;
1014 }
1015 
1016 template <typename Storage>
1017 KOKKOS_INLINE_FUNCTION
1018 PCE<Storage>
1019 tanh(const PCE<Storage>& a)
1020 {
1021  KOKKOS_IF_ON_HOST(
1022  TEUCHOS_TEST_FOR_EXCEPTION(
1023  a.size() != 1, std::logic_error,
1024  "Sacado::UQ::PCE::tanh(): argument has size != 1");
1025  )
1026 
1027  PCE<Storage> c(a.cijk(), 1);
1028  c.fastAccessCoeff(0) = std::tanh( a.fastAccessCoeff(0) );
1029 
1030  return c;
1031 }
1032 
1033 template <typename Storage>
1034 KOKKOS_INLINE_FUNCTION
1035 PCE<Storage>
1036 acos(const PCE<Storage>& a)
1037 {
1038  KOKKOS_IF_ON_HOST(
1039  TEUCHOS_TEST_FOR_EXCEPTION(
1040  a.size() != 1, std::logic_error,
1041  "Sacado::UQ::PCE::acos(): argument has size != 1");
1042  )
1043 
1044  PCE<Storage> c(a.cijk(), 1);
1045  c.fastAccessCoeff(0) = std::acos( a.fastAccessCoeff(0) );
1046 
1047  return c;
1048 }
1049 
1050 template <typename Storage>
1051 KOKKOS_INLINE_FUNCTION
1052 PCE<Storage>
1053 asin(const PCE<Storage>& a)
1054 {
1055  KOKKOS_IF_ON_HOST(
1056  TEUCHOS_TEST_FOR_EXCEPTION(
1057  a.size() != 1, std::logic_error,
1058  "Sacado::UQ::PCE::asin(): argument has size != 1");
1059  )
1060 
1061  PCE<Storage> c(a.cijk(), 1);
1062  c.fastAccessCoeff(0) = std::asin( a.fastAccessCoeff(0) );
1063 
1064  return c;
1065 }
1066 
1067 template <typename Storage>
1068 KOKKOS_INLINE_FUNCTION
1069 PCE<Storage>
1070 atan(const PCE<Storage>& a)
1071 {
1072  KOKKOS_IF_ON_HOST(
1073  TEUCHOS_TEST_FOR_EXCEPTION(
1074  a.size() != 1, std::logic_error,
1075  "Sacado::UQ::PCE::atan(): argument has size != 1");
1076  )
1077 
1078  PCE<Storage> c(a.cijk(), 1);
1079  c.fastAccessCoeff(0) = std::atan( a.fastAccessCoeff(0) );
1080 
1081  return c;
1082 }
1083 
1084 /*
1085 template <typename Storage>
1086 KOKKOS_INLINE_FUNCTION
1087 PCE<Storage>
1088 acosh(const PCE<Storage>& a)
1089 {
1090  KOKKOS_IF_ON_HOST(
1091  TEUCHOS_TEST_FOR_EXCEPTION(
1092  a.size() != 1, std::logic_error,
1093  "Sacado::UQ::PCE::acosh(): argument has size != 1");
1094  )
1095 
1096  PCE<Storage> c(a.cijk(), 1);
1097  c.fastAccessCoeff(0) = std::acosh( a.fastAccessCoeff(0) );
1098 
1099  return c;
1100 }
1101 
1102 template <typename Storage>
1103 KOKKOS_INLINE_FUNCTION
1104 PCE<Storage>
1105 asinh(const PCE<Storage>& a)
1106 {
1107  KOKKOS_IF_ON_HOST(
1108  TEUCHOS_TEST_FOR_EXCEPTION(
1109  a.size() != 1, std::logic_error,
1110  "Sacado::UQ::PCE::asinh(): argument has size != 1");
1111  )
1112 
1113  PCE<Storage> c(a.cijk(), 1);
1114  c.fastAccessCoeff(0) = std::asinh( a.fastAccessCoeff(0) );
1115 
1116  return c;
1117 }
1118 
1119 template <typename Storage>
1120 KOKKOS_INLINE_FUNCTION
1121 PCE<Storage>
1122 atanh(const PCE<Storage>& a)
1123 {
1124  KOKKOS_IF_ON_HOST(
1125  TEUCHOS_TEST_FOR_EXCEPTION(
1126  a.size() != 1, std::logic_error,
1127  "Sacado::UQ::PCE::atanh(): argument has size != 1");
1128  )
1129 
1130  PCE<Storage> c(a.cijk(), 1);
1131  c.fastAccessCoeff(0) = std::atanh( a.fastAccessCoeff(0) );
1132 
1133  return c;
1134 }
1135 */
1136 
1137 template <typename Storage>
1138 KOKKOS_INLINE_FUNCTION
1139 PCE<Storage>
1140 fabs(const PCE<Storage>& a)
1141 {
1142  PCE<Storage> c(a.cijk(), 1);
1143  c.fastAccessCoeff(0) = a.two_norm();
1144  return c;
1145 }
1146 
1147 template <typename Storage>
1148 KOKKOS_INLINE_FUNCTION
1149 PCE<Storage>
1150 abs(const PCE<Storage>& a)
1151 {
1152  PCE<Storage> c(a.cijk(), 1);
1153  c.fastAccessCoeff(0) = a.two_norm();
1154  return c;
1155 }
1156 
1157 template <typename Storage>
1158 KOKKOS_INLINE_FUNCTION
1159 PCE<Storage>
1160 ceil(const PCE<Storage>& a)
1161 {
1162  PCE<Storage> c(a.cijk(), 1);
1163  c.fastAccessCoeff(0) = std::ceil( a.fastAccessCoeff(0) );
1164  return c;
1165 }
1166 
1167 // template <typename Storage>
1168 // KOKKOS_INLINE_FUNCTION
1169 // PCE<Storage>
1170 // max(const PCE<Storage>& a, const PCE<Storage>& b)
1171 // {
1172 // if (a.two_norm() >= b.two_norm())
1173 // return a;
1174 // return b;
1175 // }
1176 
1177 template <typename Storage>
1178 KOKKOS_INLINE_FUNCTION
1179 PCE<Storage>
1180 max(const typename PCE<Storage>::value_type& a,
1181  const PCE<Storage>& b)
1182 {
1183  if (a >= b.two_norm()) {
1184  PCE<Storage> c(b.cijk(), 1);
1185  c.fastAccessCoeff(0) = a;
1186  return c;
1187  }
1188  return b;
1189 }
1190 
1191 template <typename Storage>
1192 KOKKOS_INLINE_FUNCTION
1193 PCE<Storage>
1194 max(const PCE<Storage>& a,
1195  const typename PCE<Storage>::value_type& b)
1196 {
1197  if (a.two_norm() >= b)
1198  return a;
1199  PCE<Storage> c(a.cijk(), 1);
1200  c.fastAccessCoeff(0) = b;
1201  return c;
1202 }
1203 
1204 // template <typename Storage>
1205 // KOKKOS_INLINE_FUNCTION
1206 // PCE<Storage>
1207 // min(const PCE<Storage>& a, const PCE<Storage>& b)
1208 // {
1209 // if (a.two_norm() <= b.two_norm())
1210 // return a;
1211 // return b;
1212 // }
1213 
1214 template <typename Storage>
1215 KOKKOS_INLINE_FUNCTION
1216 PCE<Storage>
1217 min(const typename PCE<Storage>::value_type& a,
1218  const PCE<Storage>& b)
1219 {
1220  if (a <= b.two_norm()) {
1221  PCE<Storage> c(b.cijk(), 1);
1222  c.fastAccessCoeff(0) = a;
1223  return c;
1224  }
1225  return b;
1226 }
1227 
1228 template <typename Storage>
1229 KOKKOS_INLINE_FUNCTION
1230 PCE<Storage>
1231 min(const PCE<Storage>& a,
1232  const typename PCE<Storage>::value_type& b)
1233 {
1234  if (a.two_norm() <= b)
1235  return a;
1236  PCE<Storage> c(a.cijk(), 1);
1237  c.fastAccessCoeff(0) = b;
1238  return c;
1239 }
1240 
1241 template <typename Storage>
1242 KOKKOS_INLINE_FUNCTION
1243 bool
1244 operator==(const PCE<Storage>& a, const PCE<Storage>& b)
1245 {
1246  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1247  const ordinal_type asz = a.size();
1248  const ordinal_type bsz = b.size();
1249  const ordinal_type n = asz > bsz ? asz : bsz;
1250  for (ordinal_type i=0; i<n; i++)
1251  if (a.coeff(i) != b.coeff(i))
1252  return false;
1253  return true;
1254 }
1255 
1256 template <typename Storage>
1257 KOKKOS_INLINE_FUNCTION
1258 bool
1259 operator==(const typename PCE<Storage>::value_type& a,
1260  const PCE<Storage>& b)
1261 {
1262  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1263  const ordinal_type n = b.size();
1264  if (a != b.coeff(0))
1265  return false;
1266  for (ordinal_type i=1; i<n; i++)
1267  if (b.fastAccessCoeff(i) != typename PCE<Storage>::value_type(0.0))
1268  return false;
1269  return true;
1270 }
1271 
1272 template <typename Storage>
1273 KOKKOS_INLINE_FUNCTION
1274 bool
1275 operator==(const PCE<Storage>& a,
1276  const typename PCE<Storage>::value_type& b)
1277 {
1278  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1279  const ordinal_type n = a.size();
1280  if (a.coeff(0) != b)
1281  return false;
1282  for (ordinal_type i=1; i<n; i++)
1283  if (a.fastAccessCoeff(i) != typename PCE<Storage>::value_type(0.0))
1284  return false;
1285  return true;
1286 }
1287 
1288 template <typename Storage>
1289 KOKKOS_INLINE_FUNCTION
1290 bool
1291 operator!=(const PCE<Storage>& a, const PCE<Storage>& b)
1292 {
1293  return !(a == b);
1294 }
1295 
1296 template <typename Storage>
1297 KOKKOS_INLINE_FUNCTION
1298 bool
1299 operator!=(const typename PCE<Storage>::value_type& a,
1300  const PCE<Storage>& b)
1301 {
1302  return !(a == b);
1303 }
1304 
1305 template <typename Storage>
1306 KOKKOS_INLINE_FUNCTION
1307 bool
1308 operator!=(const PCE<Storage>& a,
1309  const typename PCE<Storage>::value_type& b)
1310 {
1311  return !(a == b);
1312 }
1313 
1314 template <typename Storage>
1315 KOKKOS_INLINE_FUNCTION
1316 bool
1317 operator<=(const PCE<Storage>& a, const PCE<Storage>& b)
1318 {
1319  return a.two_norm() <= b.two_norm();
1320 }
1321 
1322 template <typename Storage>
1323 KOKKOS_INLINE_FUNCTION
1324 bool
1325 operator<=(const typename PCE<Storage>::value_type& a,
1326  const PCE<Storage>& b)
1327 {
1328  return a <= b.two_norm();
1329 }
1330 
1331 template <typename Storage>
1332 KOKKOS_INLINE_FUNCTION
1333 bool
1334 operator<=(const PCE<Storage>& a,
1335  const typename PCE<Storage>::value_type& b)
1336 {
1337  return a.two_norm() <= b;
1338 }
1339 
1340 template <typename Storage>
1341 KOKKOS_INLINE_FUNCTION
1342 bool
1343 operator>=(const PCE<Storage>& a, const PCE<Storage>& b)
1344 {
1345  return a.two_norm() >= b.two_norm();
1346 }
1347 
1348 template <typename Storage>
1349 KOKKOS_INLINE_FUNCTION
1350 bool
1351 operator>=(const typename PCE<Storage>::value_type& a,
1352  const PCE<Storage>& b)
1353 {
1354  return a >= b.two_norm();
1355 }
1356 
1357 template <typename Storage>
1358 KOKKOS_INLINE_FUNCTION
1359 bool
1360 operator>=(const PCE<Storage>& a,
1361  const typename PCE<Storage>::value_type& b)
1362 {
1363  return a.two_norm() >= b;
1364 }
1365 
1366 template <typename Storage>
1367 KOKKOS_INLINE_FUNCTION
1368 bool
1369 operator<(const PCE<Storage>& a, const PCE<Storage>& b)
1370 {
1371  return a.two_norm() < b.two_norm();
1372 }
1373 
1374 template <typename Storage>
1375 KOKKOS_INLINE_FUNCTION
1376 bool
1377 operator<(const typename PCE<Storage>::value_type& a,
1378  const PCE<Storage>& b)
1379 {
1380  return a < b.two_norm();
1381 }
1382 
1383 template <typename Storage>
1384 KOKKOS_INLINE_FUNCTION
1385 bool
1386 operator<(const PCE<Storage>& a,
1387  const typename PCE<Storage>::value_type& b)
1388 {
1389  return a.two_norm() < b;
1390 }
1391 
1392 template <typename Storage>
1393 KOKKOS_INLINE_FUNCTION
1394 bool
1395 operator>(const PCE<Storage>& a, const PCE<Storage>& b)
1396 {
1397  return a.two_norm() > b.two_norm();
1398 }
1399 
1400 template <typename Storage>
1401 KOKKOS_INLINE_FUNCTION
1402 bool
1403 operator>(const typename PCE<Storage>::value_type& a,
1404  const PCE<Storage>& b)
1405 {
1406  return a > b.two_norm();
1407 }
1408 
1409 template <typename Storage>
1410 KOKKOS_INLINE_FUNCTION
1411 bool
1412 operator>(const PCE<Storage>& a,
1413  const typename PCE<Storage>::value_type& b)
1414 {
1415  return a.two_norm() > b;
1416 }
1417 
1418 template <typename Storage>
1419 KOKKOS_INLINE_FUNCTION
1420 bool toBool(const PCE<Storage>& x) {
1421  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1422  bool is_zero = true;
1423  const ordinal_type sz = x.size();
1424  for (ordinal_type i=0; i<sz; i++)
1425  is_zero = is_zero && (x.fastAccessCoeff(i) == 0.0);
1426  return !is_zero;
1427 }
1428 
1429 template <typename Storage>
1430 KOKKOS_INLINE_FUNCTION
1431 bool
1432 operator && (const PCE<Storage>& x1, const PCE<Storage>& x2)
1433 {
1434  return toBool(x1) && toBool(x2);
1435 }
1436 
1437 template <typename Storage>
1438 KOKKOS_INLINE_FUNCTION
1439 bool
1440 operator && (const typename PCE<Storage>::value_type& a,
1441  const PCE<Storage>& x2)
1442 {
1443  return a && toBool(x2);
1444 }
1445 
1446 template <typename Storage>
1447 KOKKOS_INLINE_FUNCTION
1448 bool
1449 operator && (const PCE<Storage>& x1,
1450  const typename PCE<Storage>::value_type& b)
1451 {
1452  return toBool(x1) && b;
1453 }
1454 
1455 template <typename Storage>
1456 KOKKOS_INLINE_FUNCTION
1457 bool
1458 operator || (const PCE<Storage>& x1, const PCE<Storage>& x2)
1459 {
1460  return toBool(x1) || toBool(x2);
1461 }
1462 
1463 template <typename Storage>
1464 KOKKOS_INLINE_FUNCTION
1465 bool
1466 operator || (const typename PCE<Storage>::value_type& a,
1467  const PCE<Storage>& x2)
1468 {
1469  return a || toBool(x2);
1470 }
1471 
1472 template <typename Storage>
1473 KOKKOS_INLINE_FUNCTION
1474 bool
1475 operator || (const PCE<Storage>& x1,
1476  const typename PCE<Storage>::value_type& b)
1477 {
1478  return toBool(x1) || b;
1479 }
1480 
1481 template <typename Storage>
1482 std::ostream&
1483 operator << (std::ostream& os, const PCE<Storage>& a)
1484 {
1485  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1486 
1487  os << "[ ";
1488 
1489  for (ordinal_type i=0; i<a.size(); i++) {
1490  os << a.coeff(i) << " ";
1491  }
1492 
1493  os << "]";
1494  return os;
1495 }
1496 
1497 template <typename Storage>
1498 std::istream&
1499 operator >> (std::istream& is, PCE<Storage>& a)
1500 {
1501  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1502 
1503  // Read in the opening "["
1504  char bracket;
1505  is >> bracket;
1506 
1507  for (ordinal_type i=0; i<a.size(); i++) {
1508  is >> a.fastAccessCoeff(i);
1509  }
1510 
1511  // Read in the closing "]"
1512 
1513  is >> bracket;
1514  return is;
1515 }
1516 
1517 template <typename Storage>
1518 KOKKOS_INLINE_FUNCTION
1519 void
1520 CG_divide(const PCE<Storage>& a, const PCE<Storage>& b, PCE<Storage>& c) {
1521  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1522  typedef typename PCE<Storage>::value_type value_type;
1523 
1524  const ordinal_type size = c.size();
1525 
1526  //Needed scalars
1527  value_type alpha, beta, rTz, rTz_old, resid;
1528 
1529  //Needed temporary PCEs
1530  PCE<Storage> r(a.cijk(),size);
1531  PCE<Storage> p(a.cijk(),size);
1532  PCE<Storage> bp(a.cijk(),size);
1533  PCE<Storage> z(a.cijk(),size);
1534 
1535  //compute residual = a - b*c
1536  r = a - b*c;
1537  z = r/b.coeff(0);
1538  p = z;
1539  resid = r.two_norm();
1540  //Compute <r,z>=rTz (L2 inner product)
1541  rTz = r.inner_product(z);
1542  ordinal_type k = 0;
1543  value_type tol = 1e-6;
1544  while ( resid > tol && k < 100){
1545  bp = b*p;
1546  //Compute alpha = <r,z>/<p,b*p>
1547  alpha = rTz/p.inner_product(bp);
1548  //Update the solution c = c + alpha*p
1549  c = c + alpha*p;
1550  rTz_old = rTz;
1551  //Compute the new residual r = r - alpha*b*p
1552  r = r - alpha*bp;
1553  resid = r.two_norm();
1554  //Compute beta = rTz_new/ rTz_old and new p
1555  z = r/b.coeff(0);
1556  rTz = r.inner_product(z);
1557  beta = rTz/rTz_old;
1558  p = z + beta*p;
1559  k++;
1560  }
1561  }
1562 
1563 } // namespace UQ
1564 } // namespace Sacado
Sacado::UQ::sqrt
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:847
Sacado::UQ::fabs
KOKKOS_INLINE_FUNCTION PCE< Storage > fabs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1140
Sacado::UQ::tan
KOKKOS_INLINE_FUNCTION PCE< Storage > tan(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:968
Sacado::UQ::operator||
KOKKOS_INLINE_FUNCTION bool operator||(const PCE< Storage > &x1, const PCE< Storage > &x2)
Definition: Sacado_UQ_PCE_Imp.hpp:1458
Sacado::UQ::sinh
KOKKOS_INLINE_FUNCTION PCE< Storage > sinh(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:985
Teuchos_SerialDenseMatrix.hpp
Sacado::PCE::operator-
OrthogPoly< T, Storage > operator-(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
Definition: Sacado_PCE_OrthogPolyImp.hpp:310
Sacado::UQ::pow
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:881
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Sacado::UQ::CG_divide
KOKKOS_INLINE_FUNCTION void CG_divide(const PCE< Storage > &a, const PCE< Storage > &b, PCE< Storage > &c)
Definition: Sacado_UQ_PCE_Imp.hpp:1520
Sacado::UQ::operator&&
KOKKOS_INLINE_FUNCTION bool operator&&(const PCE< Storage > &x1, const PCE< Storage > &x2)
Definition: Sacado_UQ_PCE_Imp.hpp:1432
Sacado::UQ::tanh
KOKKOS_INLINE_FUNCTION PCE< Storage > tanh(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1019
Sacado::UQ::cbrt
KOKKOS_INLINE_FUNCTION PCE< Storage > cbrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:864
Sacado::UQ::operator!=
KOKKOS_INLINE_FUNCTION bool operator!=(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1291
Sacado::UQ::acos
KOKKOS_INLINE_FUNCTION PCE< Storage > acos(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1036
Sacado::UQ::min
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1217
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Sacado::UQ::operator>
KOKKOS_INLINE_FUNCTION bool operator>(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1395
Sacado::UQ::operator>>
std::istream & operator>>(std::istream &is, PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1499
Sacado::UQ::operator/
KOKKOS_INLINE_FUNCTION PCE< Storage > operator/(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:717
Sacado::UQ::ceil
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1160
Sacado::UQ::operator-
KOKKOS_INLINE_FUNCTION PCE< Storage > operator-(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:541
Sacado::UQ::operator>=
KOKKOS_INLINE_FUNCTION bool operator>=(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1343
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:1180
Sacado::UQ::cosh
KOKKOS_INLINE_FUNCTION PCE< Storage > cosh(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1002
Sacado::is_constant
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
Definition: Stokhos_Is_Constant.hpp:21
Sacado::UQ::operator==
KOKKOS_INLINE_FUNCTION bool operator==(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1244
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
Sacado::UQ::atan
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1070
Sacado::UQ::exp
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:796
Sacado::UQ::sin
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:934
Sacado::UQ::operator+
KOKKOS_INLINE_FUNCTION PCE< Storage > operator+(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:466
Sacado::UQ::operator*
KOKKOS_INLINE_FUNCTION PCE< Storage > operator*(const PCE< Storage > &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:616
Sacado::UQ::log
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:813
Sacado::UQ::log10
KOKKOS_INLINE_FUNCTION PCE< Storage > log10(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:830
n
int n
Sacado::UQ::asin
KOKKOS_INLINE_FUNCTION PCE< Storage > asin(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1053
Sacado::UQ::cos
KOKKOS_INLINE_FUNCTION PCE< Storage > cos(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:951
Sacado::UQ::toBool
KOKKOS_INLINE_FUNCTION bool toBool(const PCE< Storage > &x)
Definition: Sacado_UQ_PCE_Imp.hpp:1420
Generated on Fri Nov 22 2024 09:18:09 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5