Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
45namespace Sacado {
46namespace UQ {
47
48/*
49template <typename Storage>
50KOKKOS_INLINE_FUNCTION
51typename PCE<Storage>::value_type
52PCE<Storage>::
53evaluate(const Teuchos::Array<typename PCE<Storage>::value_type>& point) const
54{
55 return s_.evaluate(point);
56}
57
58template <typename Storage>
59KOKKOS_INLINE_FUNCTION
60typename PCE<Storage>::value_type
61PCE<Storage>::
62evaluate(
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
70template <typename Storage>
71KOKKOS_INLINE_FUNCTION
72typename PCE<Storage>::value_type
73PCE<Storage>::
74standard_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
83template <typename Storage>
84KOKKOS_INLINE_FUNCTION
85typename PCE<Storage>::value_type
86PCE<Storage>::
87two_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
96template <typename Storage>
97KOKKOS_INLINE_FUNCTION
98typename PCE<Storage>::value_type
99PCE<Storage>::
100inner_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
112template <typename Storage>
113KOKKOS_INLINE_FUNCTION
114bool
115PCE<Storage>::
116isEqualTo(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
126template <typename Storage>
127KOKKOS_INLINE_FUNCTION
128bool
129PCE<Storage>::
130isEqualTo(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
140template <typename Storage>
141KOKKOS_INLINE_FUNCTION
142PCE<Storage>&
143PCE<Storage>::
144operator=(const typename PCE<Storage>::value_type v)
145{
146 s_.init(value_type(0));
147 s_[0] = v;
148 return *this;
149}
150
151template <typename Storage>
152KOKKOS_INLINE_FUNCTION
153/*volatile*/ PCE<Storage>&
154PCE<Storage>::
155operator=(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
162template <typename Storage>
163KOKKOS_INLINE_FUNCTION
164PCE<Storage>&
165PCE<Storage>::
166operator=(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
188template <typename Storage>
189KOKKOS_INLINE_FUNCTION
190PCE<Storage>&
191PCE<Storage>::
192operator=(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
214template <typename Storage>
215KOKKOS_INLINE_FUNCTION
216/*volatile*/ PCE<Storage>&
217PCE<Storage>::
218operator=(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
240template <typename Storage>
241KOKKOS_INLINE_FUNCTION
242/*volatile*/ PCE<Storage>&
243PCE<Storage>::
244operator=(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
267template <typename Storage>
268KOKKOS_INLINE_FUNCTION
269PCE<Storage>
271operator-() 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
282template <typename Storage>
283KOKKOS_INLINE_FUNCTION
284PCE<Storage>
286operator-() 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
297template <typename Storage>
298KOKKOS_INLINE_FUNCTION
299PCE<Storage>&
300PCE<Storage>::
301operator*=(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
310template <typename Storage>
311KOKKOS_INLINE_FUNCTION
312/*volatile*/ PCE<Storage>&
313PCE<Storage>::
314operator*=(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
323template <typename Storage>
324KOKKOS_INLINE_FUNCTION
325PCE<Storage>&
326PCE<Storage>::
327operator/=(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
336template <typename Storage>
337KOKKOS_INLINE_FUNCTION
338/*volatile*/ PCE<Storage>&
339PCE<Storage>::
340operator/=(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
349template <typename Storage>
350KOKKOS_INLINE_FUNCTION
351PCE<Storage>&
352PCE<Storage>::
353operator+=(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
367template <typename Storage>
368KOKKOS_INLINE_FUNCTION
369PCE<Storage>&
370PCE<Storage>::
371operator-=(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
385template <typename Storage>
386KOKKOS_INLINE_FUNCTION
387PCE<Storage>&
388PCE<Storage>::
389operator*=(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
445template <typename Storage>
446KOKKOS_INLINE_FUNCTION
447PCE<Storage>&
448PCE<Storage>::
449operator/=(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
494template <typename Storage>
495KOKKOS_INLINE_FUNCTION
496PCE<Storage>
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
529template <typename Storage>
530KOKKOS_INLINE_FUNCTION
531PCE<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
549template <typename Storage>
550KOKKOS_INLINE_FUNCTION
551PCE<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
569template <typename Storage>
570KOKKOS_INLINE_FUNCTION
571PCE<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
604template <typename Storage>
605KOKKOS_INLINE_FUNCTION
606PCE<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
624template <typename Storage>
625KOKKOS_INLINE_FUNCTION
626PCE<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
644template <typename Storage>
645KOKKOS_INLINE_FUNCTION
646PCE<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
707template <typename Storage>
708KOKKOS_INLINE_FUNCTION
709PCE<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
726template <typename Storage>
727KOKKOS_INLINE_FUNCTION
728PCE<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
745template <typename Storage>
746KOKKOS_INLINE_FUNCTION
747PCE<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__)
761TEUCHOS_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
794template <typename Storage>
795KOKKOS_INLINE_FUNCTION
796PCE<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
805template <typename Storage>
806KOKKOS_INLINE_FUNCTION
807PCE<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
824template <typename Storage>
825KOKKOS_INLINE_FUNCTION
826PCE<Storage>
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
841template <typename Storage>
842KOKKOS_INLINE_FUNCTION
843PCE<Storage>
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
858template <typename Storage>
859KOKKOS_INLINE_FUNCTION
860PCE<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
875template <typename Storage>
876KOKKOS_INLINE_FUNCTION
877PCE<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
892template <typename Storage>
893KOKKOS_INLINE_FUNCTION
894PCE<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
909template <typename Storage>
910KOKKOS_INLINE_FUNCTION
911PCE<Storage>
912pow(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
926template <typename Storage>
927KOKKOS_INLINE_FUNCTION
928PCE<Storage>
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
944template <typename Storage>
945KOKKOS_INLINE_FUNCTION
946PCE<Storage>
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
962template <typename Storage>
963KOKKOS_INLINE_FUNCTION
964PCE<Storage>
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
979template <typename Storage>
980KOKKOS_INLINE_FUNCTION
981PCE<Storage>
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
996template <typename Storage>
997KOKKOS_INLINE_FUNCTION
998PCE<Storage>
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
1013template <typename Storage>
1014KOKKOS_INLINE_FUNCTION
1015PCE<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
1030template <typename Storage>
1031KOKKOS_INLINE_FUNCTION
1032PCE<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
1047template <typename Storage>
1048KOKKOS_INLINE_FUNCTION
1049PCE<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
1064template <typename Storage>
1065KOKKOS_INLINE_FUNCTION
1066PCE<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
1081template <typename Storage>
1082KOKKOS_INLINE_FUNCTION
1083PCE<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
1098template <typename Storage>
1099KOKKOS_INLINE_FUNCTION
1100PCE<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/*
1116template <typename Storage>
1117KOKKOS_INLINE_FUNCTION
1118PCE<Storage>
1119acosh(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
1133template <typename Storage>
1134KOKKOS_INLINE_FUNCTION
1135PCE<Storage>
1136asinh(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
1150template <typename Storage>
1151KOKKOS_INLINE_FUNCTION
1152PCE<Storage>
1153atanh(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
1168template <typename Storage>
1169KOKKOS_INLINE_FUNCTION
1170PCE<Storage>
1172{
1173 PCE<Storage> c(a.cijk(), 1);
1174 c.fastAccessCoeff(0) = a.two_norm();
1175 return c;
1176}
1177
1178template <typename Storage>
1179KOKKOS_INLINE_FUNCTION
1180PCE<Storage>
1182{
1183 PCE<Storage> c(a.cijk(), 1);
1184 c.fastAccessCoeff(0) = a.two_norm();
1185 return c;
1186}
1187
1188template <typename Storage>
1189KOKKOS_INLINE_FUNCTION
1190PCE<Storage>
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
1208template <typename Storage>
1209KOKKOS_INLINE_FUNCTION
1210PCE<Storage>
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
1222template <typename Storage>
1223PCE<Storage>
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
1244template <typename Storage>
1245KOKKOS_INLINE_FUNCTION
1246PCE<Storage>
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
1258template <typename Storage>
1259KOKKOS_INLINE_FUNCTION
1260PCE<Storage>
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
1271template <typename Storage>
1272KOKKOS_INLINE_FUNCTION
1273bool
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
1286template <typename Storage>
1287KOKKOS_INLINE_FUNCTION
1288bool
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
1302template <typename Storage>
1303KOKKOS_INLINE_FUNCTION
1304bool
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
1318template <typename Storage>
1319KOKKOS_INLINE_FUNCTION
1320bool
1322{
1323 return !(a == b);
1324}
1325
1326template <typename Storage>
1327KOKKOS_INLINE_FUNCTION
1328bool
1330 const PCE<Storage>& b)
1331{
1332 return !(a == b);
1333}
1334
1335template <typename Storage>
1336KOKKOS_INLINE_FUNCTION
1337bool
1339 const typename PCE<Storage>::value_type& b)
1340{
1341 return !(a == b);
1342}
1343
1344template <typename Storage>
1345KOKKOS_INLINE_FUNCTION
1346bool
1347operator<=(const PCE<Storage>& a, const PCE<Storage>& b)
1348{
1349 return a.two_norm() <= b.two_norm();
1350}
1351
1352template <typename Storage>
1353KOKKOS_INLINE_FUNCTION
1354bool
1355operator<=(const typename PCE<Storage>::value_type& a,
1356 const PCE<Storage>& b)
1357{
1358 return a <= b.two_norm();
1359}
1360
1361template <typename Storage>
1362KOKKOS_INLINE_FUNCTION
1363bool
1364operator<=(const PCE<Storage>& a,
1365 const typename PCE<Storage>::value_type& b)
1366{
1367 return a.two_norm() <= b;
1368}
1369
1370template <typename Storage>
1371KOKKOS_INLINE_FUNCTION
1372bool
1374{
1375 return a.two_norm() >= b.two_norm();
1376}
1377
1378template <typename Storage>
1379KOKKOS_INLINE_FUNCTION
1380bool
1382 const PCE<Storage>& b)
1383{
1384 return a >= b.two_norm();
1385}
1386
1387template <typename Storage>
1388KOKKOS_INLINE_FUNCTION
1389bool
1391 const typename PCE<Storage>::value_type& b)
1392{
1393 return a.two_norm() >= b;
1394}
1395
1396template <typename Storage>
1397KOKKOS_INLINE_FUNCTION
1398bool
1399operator<(const PCE<Storage>& a, const PCE<Storage>& b)
1400{
1401 return a.two_norm() < b.two_norm();
1402}
1403
1404template <typename Storage>
1405KOKKOS_INLINE_FUNCTION
1406bool
1407operator<(const typename PCE<Storage>::value_type& a,
1408 const PCE<Storage>& b)
1409{
1410 return a < b.two_norm();
1411}
1412
1413template <typename Storage>
1414KOKKOS_INLINE_FUNCTION
1415bool
1416operator<(const PCE<Storage>& a,
1417 const typename PCE<Storage>::value_type& b)
1418{
1419 return a.two_norm() < b;
1420}
1421
1422template <typename Storage>
1423KOKKOS_INLINE_FUNCTION
1424bool
1426{
1427 return a.two_norm() > b.two_norm();
1428}
1429
1430template <typename Storage>
1431KOKKOS_INLINE_FUNCTION
1432bool
1434 const PCE<Storage>& b)
1435{
1436 return a > b.two_norm();
1437}
1438
1439template <typename Storage>
1440KOKKOS_INLINE_FUNCTION
1441bool
1443 const typename PCE<Storage>::value_type& b)
1444{
1445 return a.two_norm() > b;
1446}
1447
1448template <typename Storage>
1449KOKKOS_INLINE_FUNCTION
1450bool 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
1459template <typename Storage>
1460KOKKOS_INLINE_FUNCTION
1461bool
1463{
1464 return toBool(x1) && toBool(x2);
1465}
1466
1467template <typename Storage>
1468KOKKOS_INLINE_FUNCTION
1469bool
1471 const PCE<Storage>& x2)
1472{
1473 return a && toBool(x2);
1474}
1475
1476template <typename Storage>
1477KOKKOS_INLINE_FUNCTION
1478bool
1480 const typename PCE<Storage>::value_type& b)
1481{
1482 return toBool(x1) && b;
1483}
1484
1485template <typename Storage>
1486KOKKOS_INLINE_FUNCTION
1487bool
1489{
1490 return toBool(x1) || toBool(x2);
1491}
1492
1493template <typename Storage>
1494KOKKOS_INLINE_FUNCTION
1495bool
1497 const PCE<Storage>& x2)
1498{
1499 return a || toBool(x2);
1500}
1501
1502template <typename Storage>
1503KOKKOS_INLINE_FUNCTION
1504bool
1506 const typename PCE<Storage>::value_type& b)
1507{
1508 return toBool(x1) || b;
1509}
1510
1511template <typename Storage>
1512std::ostream&
1513operator << (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
1527template <typename Storage>
1528std::istream&
1529operator >> (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
1547template <typename Storage>
1548void
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
OrthogPoly< T, Storage > operator-(const OrthogPoly< T, Storage > &a, const OrthogPoly< T, Storage > &b)
KOKKOS_INLINE_FUNCTION bool operator||(const PCE< Storage > &x1, const PCE< Storage > &x2)
KOKKOS_INLINE_FUNCTION bool toBool(const PCE< Storage > &x)
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 > tan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > tanh(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cos(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator<(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
std::ostream & operator<<(std::ostream &os, const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cosh(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool 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 > abs(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > sinh(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 > acos(const PCE< Storage > &a)
std::istream & operator>>(std::istream &is, PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator!=(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(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 > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cbrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > fabs(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 bool operator>=(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator<=(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION bool operator&&(const PCE< Storage > &x1, const PCE< Storage > &x2)
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > operator-(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267