Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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
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
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
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
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
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
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>&
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
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
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>&
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>&
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
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
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
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>&
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
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>&
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
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
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
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
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
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
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 // max(const PCE<Storage>& a, const PCE<Storage>& b)
1192 // {
1193 // if (a.two_norm() >= b.two_norm())
1194 // return a;
1195 // return b;
1196 // }
1197 
1198 template <typename Storage>
1199 KOKKOS_INLINE_FUNCTION
1200 PCE<Storage>
1201 max(const typename PCE<Storage>::value_type& a,
1202  const PCE<Storage>& b)
1203 {
1204  if (a >= b.two_norm()) {
1205  PCE<Storage> c(b.cijk(), 1);
1206  c.fastAccessCoeff(0) = a;
1207  return c;
1208  }
1209  return b;
1210 }
1211 
1212 template <typename Storage>
1213 PCE<Storage>
1215  const typename PCE<Storage>::value_type& b)
1216 {
1217  if (a.two_norm() >= b)
1218  return a;
1219  PCE<Storage> c(a.cijk(), 1);
1220  c.fastAccessCoeff(0) = b;
1221  return c;
1222 }
1223 
1224 // template <typename Storage>
1225 // KOKKOS_INLINE_FUNCTION
1226 // PCE<Storage>
1227 // min(const PCE<Storage>& a, const PCE<Storage>& b)
1228 // {
1229 // if (a.two_norm() <= b.two_norm())
1230 // return a;
1231 // return b;
1232 // }
1233 
1234 template <typename Storage>
1235 KOKKOS_INLINE_FUNCTION
1236 PCE<Storage>
1237 min(const typename PCE<Storage>::value_type& a,
1238  const PCE<Storage>& b)
1239 {
1240  if (a <= b.two_norm()) {
1241  PCE<Storage> c(b.cijk(), 1);
1242  c.fastAccessCoeff(0) = a;
1243  return c;
1244  }
1245  return b;
1246 }
1247 
1248 template <typename Storage>
1249 KOKKOS_INLINE_FUNCTION
1250 PCE<Storage>
1252  const typename PCE<Storage>::value_type& b)
1253 {
1254  if (a.two_norm() <= b)
1255  return a;
1256  PCE<Storage> c(a.cijk(), 1);
1257  c.fastAccessCoeff(0) = b;
1258  return c;
1259 }
1260 
1261 template <typename Storage>
1262 KOKKOS_INLINE_FUNCTION
1263 bool
1265 {
1266  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1267  const ordinal_type asz = a.size();
1268  const ordinal_type bsz = b.size();
1269  const ordinal_type n = asz > bsz ? asz : bsz;
1270  for (ordinal_type i=0; i<n; i++)
1271  if (a.coeff(i) != b.coeff(i))
1272  return false;
1273  return true;
1274 }
1275 
1276 template <typename Storage>
1277 KOKKOS_INLINE_FUNCTION
1278 bool
1280  const PCE<Storage>& b)
1281 {
1282  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1283  const ordinal_type n = b.size();
1284  if (a != b.coeff(0))
1285  return false;
1286  for (ordinal_type i=1; i<n; i++)
1287  if (b.fastAccessCoeff(i) != typename PCE<Storage>::value_type(0.0))
1288  return false;
1289  return true;
1290 }
1291 
1292 template <typename Storage>
1293 KOKKOS_INLINE_FUNCTION
1294 bool
1296  const typename PCE<Storage>::value_type& b)
1297 {
1298  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1299  const ordinal_type n = a.size();
1300  if (a.coeff(0) != b)
1301  return false;
1302  for (ordinal_type i=1; i<n; i++)
1303  if (a.fastAccessCoeff(i) != typename PCE<Storage>::value_type(0.0))
1304  return false;
1305  return true;
1306 }
1307 
1308 template <typename Storage>
1309 KOKKOS_INLINE_FUNCTION
1310 bool
1312 {
1313  return !(a == b);
1314 }
1315 
1316 template <typename Storage>
1317 KOKKOS_INLINE_FUNCTION
1318 bool
1320  const PCE<Storage>& b)
1321 {
1322  return !(a == b);
1323 }
1324 
1325 template <typename Storage>
1326 KOKKOS_INLINE_FUNCTION
1327 bool
1329  const typename PCE<Storage>::value_type& b)
1330 {
1331  return !(a == b);
1332 }
1333 
1334 template <typename Storage>
1335 KOKKOS_INLINE_FUNCTION
1336 bool
1337 operator<=(const PCE<Storage>& a, const PCE<Storage>& b)
1338 {
1339  return a.two_norm() <= b.two_norm();
1340 }
1341 
1342 template <typename Storage>
1343 KOKKOS_INLINE_FUNCTION
1344 bool
1346  const PCE<Storage>& b)
1347 {
1348  return a <= b.two_norm();
1349 }
1350 
1351 template <typename Storage>
1352 KOKKOS_INLINE_FUNCTION
1353 bool
1354 operator<=(const PCE<Storage>& a,
1355  const typename PCE<Storage>::value_type& b)
1356 {
1357  return a.two_norm() <= b;
1358 }
1359 
1360 template <typename Storage>
1361 KOKKOS_INLINE_FUNCTION
1362 bool
1364 {
1365  return a.two_norm() >= b.two_norm();
1366 }
1367 
1368 template <typename Storage>
1369 KOKKOS_INLINE_FUNCTION
1370 bool
1372  const PCE<Storage>& b)
1373 {
1374  return a >= b.two_norm();
1375 }
1376 
1377 template <typename Storage>
1378 KOKKOS_INLINE_FUNCTION
1379 bool
1381  const typename PCE<Storage>::value_type& b)
1382 {
1383  return a.two_norm() >= b;
1384 }
1385 
1386 template <typename Storage>
1387 KOKKOS_INLINE_FUNCTION
1388 bool
1389 operator<(const PCE<Storage>& a, const PCE<Storage>& b)
1390 {
1391  return a.two_norm() < b.two_norm();
1392 }
1393 
1394 template <typename Storage>
1395 KOKKOS_INLINE_FUNCTION
1396 bool
1398  const PCE<Storage>& b)
1399 {
1400  return a < b.two_norm();
1401 }
1402 
1403 template <typename Storage>
1404 KOKKOS_INLINE_FUNCTION
1405 bool
1406 operator<(const PCE<Storage>& a,
1407  const typename PCE<Storage>::value_type& b)
1408 {
1409  return a.two_norm() < b;
1410 }
1411 
1412 template <typename Storage>
1413 KOKKOS_INLINE_FUNCTION
1414 bool
1416 {
1417  return a.two_norm() > b.two_norm();
1418 }
1419 
1420 template <typename Storage>
1421 KOKKOS_INLINE_FUNCTION
1422 bool
1424  const PCE<Storage>& b)
1425 {
1426  return a > b.two_norm();
1427 }
1428 
1429 template <typename Storage>
1430 KOKKOS_INLINE_FUNCTION
1431 bool
1433  const typename PCE<Storage>::value_type& b)
1434 {
1435  return a.two_norm() > b;
1436 }
1437 
1438 template <typename Storage>
1439 KOKKOS_INLINE_FUNCTION
1440 bool toBool(const PCE<Storage>& x) {
1441  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1442  bool is_zero = true;
1443  const ordinal_type sz = x.size();
1444  for (ordinal_type i=0; i<sz; i++)
1445  is_zero = is_zero && (x.fastAccessCoeff(i) == 0.0);
1446  return !is_zero;
1447 }
1448 
1449 template <typename Storage>
1450 KOKKOS_INLINE_FUNCTION
1451 bool
1453 {
1454  return toBool(x1) && toBool(x2);
1455 }
1456 
1457 template <typename Storage>
1458 KOKKOS_INLINE_FUNCTION
1459 bool
1461  const PCE<Storage>& x2)
1462 {
1463  return a && toBool(x2);
1464 }
1465 
1466 template <typename Storage>
1467 KOKKOS_INLINE_FUNCTION
1468 bool
1470  const typename PCE<Storage>::value_type& b)
1471 {
1472  return toBool(x1) && b;
1473 }
1474 
1475 template <typename Storage>
1476 KOKKOS_INLINE_FUNCTION
1477 bool
1479 {
1480  return toBool(x1) || toBool(x2);
1481 }
1482 
1483 template <typename Storage>
1484 KOKKOS_INLINE_FUNCTION
1485 bool
1487  const PCE<Storage>& x2)
1488 {
1489  return a || toBool(x2);
1490 }
1491 
1492 template <typename Storage>
1493 KOKKOS_INLINE_FUNCTION
1494 bool
1496  const typename PCE<Storage>::value_type& b)
1497 {
1498  return toBool(x1) || b;
1499 }
1500 
1501 template <typename Storage>
1502 std::ostream&
1503 operator << (std::ostream& os, const PCE<Storage>& a)
1504 {
1505  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1506 
1507  os << "[ ";
1508 
1509  for (ordinal_type i=0; i<a.size(); i++) {
1510  os << a.coeff(i) << " ";
1511  }
1512 
1513  os << "]";
1514  return os;
1515 }
1516 
1517 template <typename Storage>
1518 std::istream&
1519 operator >> (std::istream& is, PCE<Storage>& a)
1520 {
1521  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1522 
1523  // Read in the opening "["
1524  char bracket;
1525  is >> bracket;
1526 
1527  for (ordinal_type i=0; i<a.size(); i++) {
1528  is >> a.fastAccessCoeff(i);
1529  }
1530 
1531  // Read in the closing "]"
1532 
1533  is >> bracket;
1534  return is;
1535 }
1536 
1537 template <typename Storage>
1538 void
1540  typedef typename PCE<Storage>::ordinal_type ordinal_type;
1541  typedef typename PCE<Storage>::value_type value_type;
1542 
1543  const ordinal_type size = c.size();
1544 
1545  //Needed scalars
1546  value_type alpha, beta, rTz, rTz_old, resid;
1547 
1548  //Needed temporary PCEs
1549  PCE<Storage> r(a.cijk(),size);
1550  PCE<Storage> p(a.cijk(),size);
1551  PCE<Storage> bp(a.cijk(),size);
1552  PCE<Storage> z(a.cijk(),size);
1553 
1554  //compute residual = a - b*c
1555  r = a - b*c;
1556  z = r/b.coeff(0);
1557  p = z;
1558  resid = r.two_norm();
1559  //Compute <r,z>=rTz (L2 inner product)
1560  rTz = r.inner_product(z);
1561  ordinal_type k = 0;
1562  value_type tol = 1e-6;
1563  while ( resid > tol && k < 100){
1564  bp = b*p;
1565  //Compute alpha = <r,z>/<p,b*p>
1566  alpha = rTz/p.inner_product(bp);
1567  //Update the solution c = c + alpha*p
1568  c = c + alpha*p;
1569  rTz_old = rTz;
1570  //Compute the new residual r = r - alpha*b*p
1571  r = r - alpha*bp;
1572  resid = r.two_norm();
1573  //Compute beta = rTz_new/ rTz_old and new p
1574  z = r/b.coeff(0);
1575  rTz = r.inner_product(z);
1576  beta = rTz/rTz_old;
1577  p = z + beta*p;
1578  k++;
1579  }
1580  }
1581 
1582 } // namespace UQ
1583 } // namespace Sacado
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > fabs(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > tan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator||(const PCE< Storage > &x1, const PCE< Storage > &x2)
KOKKOS_INLINE_FUNCTION PCE< Storage > sinh(const PCE< Storage > &a)
OrthogPoly< T, Storage > operator-(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const typename PCE< Storage >::value_type &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > tanh(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cbrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator!=(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > acos(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator&&(const PCE< Storage > &x1, const PCE< Storage > &x2)
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION bool operator>(const PCE< Storage > &a, const PCE< Storage > &b)
std::istream & operator>>(std::istream &is, PCE< Storage > &a)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
KOKKOS_INLINE_FUNCTION PCE< Storage > operator/(const PCE< Storage > &a, const PCE< Storage > &b)
void CG_divide(const PCE< Storage > &a, const PCE< Storage > &b, PCE< Storage > &c)
KOKKOS_INLINE_FUNCTION PCE< Storage > operator-(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION bool operator>=(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > cosh(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
KOKKOS_INLINE_FUNCTION bool operator==(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > operator+(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > operator*(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > log10(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > asin(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cos(const PCE< Storage > &a)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
KOKKOS_INLINE_FUNCTION bool toBool(const PCE< Storage > &x)