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