Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_ETPCE_OrthogPolyImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Sacado_DynamicArrayTraits.hpp"
43#include "Teuchos_TimeMonitor.hpp"
45
46namespace Sacado {
47namespace ETPCE {
48
49template <int NN, typename ExprT>
51 static const int N = NN;
52 typedef typename ExprT::value_type value_type;
53 const ExprT& expr_;
54 ExprQuadFuncWrapper(const ExprT& expr) : expr_(expr) {}
55 template <typename tuple_type>
56 KERNEL_PREFIX value_type operator() (tuple_type x) const {
57 return expr_.template eval_sample<0>(x);
58 }
59};
60
61template <typename T, typename Storage>
62template <typename S>
63void
64OrthogPolyImpl<T,Storage>::
65expressionCopy(const Expr<S>& x)
66{
67#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
68 TEUCHOS_FUNC_TIME_MONITOR("ETPCE ExpressionCopy(" << x.name() << ")");
69#endif
70 int p = x.order();
71 if (p == 0) {
72 (*th_)[0] = x.val();
73 }
74 else if (p <= 2) {
75 int sz = th_->size();
76 approx_type* tc = th_.get();
77 bool on_rhs = false;
78
79 // Check if *this is on the right-hand-side of the expression
80 if (p == 2) {
81 const int N = Expr<S>::num_args;
82 for (int i=0; i<N; i++) {
83 const approx_type* opa = &(x.getArg(i));
84 if (opa == &(*th_)) {
85 on_rhs = true;
86 break;
87 }
88 }
89 }
90
91 // If we are on the RHS, we have to put the results in a temporary
92 if (on_rhs)
93 tc = new approx_type(expansion_->getBasis(), sz);
94
95 (*tc)[0] = x.val();
96 if (x.has_fast_access(sz)) {
97 for (int i=1; i<sz; i++)
98 (*tc)[i] = x.fast_higher_order_coeff(i);
99 }
100 else {
101 for (int i=1; i<sz; i++)
102 (*tc)[i] = x.higher_order_coeff(i);
103 }
104
105 // Set underlying OPA if we created a temporary
106 if (on_rhs) {
107 th_ = Sacado::Handle<approx_type>(tc);
108 }
109 }
110 else {
111 const int N = Expr<S>::num_args;
112 const approx_type* opas[N];
113 for (int i=0; i<N; i++)
114 opas[i] = &(x.getArg(i));
115 ExprQuadFuncWrapper< N, Expr<S> > func(x);
116 quad_expansion_->nary_op(func, *th_, opas);
117 }
118}
119
120template <typename T, typename Storage>
121OrthogPolyImpl<T,Storage>::
122OrthogPolyImpl() :
123 expansion_(),
124 quad_expansion_(),
125 th_(new Stokhos::OrthogPolyApprox<int,value_type,Storage>)
126{
127 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
128 expansion_ = const_expansion_;
129}
130
131template <typename T, typename Storage>
132OrthogPolyImpl<T,Storage>::
133OrthogPolyImpl(const typename OrthogPolyImpl<T,Storage>::value_type& x) :
134 expansion_(),
135 quad_expansion_(),
136 th_(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(Teuchos::null, 1, &x))
137{
138 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
139 expansion_ = const_expansion_;
140}
141
142template <typename T, typename Storage>
143OrthogPolyImpl<T,Storage>::
144OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion) :
145 expansion_(expansion),
146 quad_expansion_(Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_)),
147 th_(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(expansion_->getBasis()))
148{
149 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
150}
151
152template <typename T, typename Storage>
153OrthogPolyImpl<T,Storage>::
154OrthogPolyImpl(const Teuchos::RCP<expansion_type>& expansion,
155 ordinal_type sz) :
156 expansion_(expansion),
157 quad_expansion_(Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_)),
158 th_(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(expansion_->getBasis(), sz))
159{
160 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
161}
162
163template <typename T, typename Storage>
164OrthogPolyImpl<T,Storage>::
165OrthogPolyImpl(const OrthogPolyImpl<T,Storage>& x) :
166 expansion_(x.expansion_),
167 quad_expansion_(x.quad_expansion_),
168 th_(x.th_)
169{
170 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
171}
172
173template <typename T, typename Storage>
174template <typename S>
175OrthogPolyImpl<T,Storage>::
176OrthogPolyImpl(const Expr<S>& x) :
177 expansion_(x.expansion()),
178 quad_expansion_(Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_)),
179 th_(new Stokhos::OrthogPolyApprox<int,value_type,Storage>(expansion_->getBasis(), x.size()))
180{
181 const_expansion_ = Teuchos::rcp(new Stokhos::ConstantOrthogPolyExpansion<int,T>);
182 expressionCopy(x);
183}
184
185template <typename T, typename Storage>
186void
187OrthogPolyImpl<T,Storage>::
188reset(const Teuchos::RCP<expansion_type>& expansion)
189{
190 expansion_ = expansion;
191 quad_expansion_ = Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_);
192 th_->reset(expansion_->getBasis());
193}
194
195template <typename T, typename Storage>
196void
197OrthogPolyImpl<T,Storage>::
198reset(const Teuchos::RCP<expansion_type>& expansion, ordinal_type sz)
199{
200 expansion_ = expansion;
201 quad_expansion_ = Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_);
202 th_->reset(expansion_->getBasis(), sz);
203}
204
205template <typename T, typename Storage>
206typename OrthogPolyImpl<T,Storage>::value_type
207OrthogPolyImpl<T,Storage>::
208evaluate(const Teuchos::Array<typename OrthogPolyImpl<T,Storage>::value_type>& point) const
209{
210 return th_->evaluate(point);
211}
212
213template <typename T, typename Storage>
214typename OrthogPolyImpl<T,Storage>::value_type
215OrthogPolyImpl<T,Storage>::
216evaluate(
217 const Teuchos::Array<typename OrthogPolyImpl<T,Storage>::value_type>& point,
218 const Teuchos::Array<typename OrthogPolyImpl<T,Storage>::value_type>& bvals) const
219{
220 return th_->evaluate(point, bvals);
221}
222
223template <typename T, typename Storage>
224template <typename S>
225bool
226OrthogPolyImpl<T,Storage>::
227isEqualTo(const Expr<S>& x) const {
228 typedef IsEqual<value_type> IE;
229 if (x.size() != this->size()) return false;
230 // Allow expansions to be different if their size is 1 and one
231 // of them is a constant expansion
232 if (expansion_ != x.expansion_) {
233 if (x.size() != 1)
234 return false;
235 if ((expansion_ != const_expansion_) &&
236 (x.expansion_ != x.const_expansion_))
237 return false;
238 }
239 bool eq = true;
240 for (int i=0; i<this->size(); i++)
241 eq = eq && IE::eval(x.coeff(i), this->coeff(i));
242 return eq;
243}
244
245template <typename T, typename Storage>
246OrthogPolyImpl<T,Storage>&
247OrthogPolyImpl<T,Storage>::
248operator=(const typename OrthogPolyImpl<T,Storage>::value_type& v)
249{
250 th_.makeOwnCopy();
251 *th_ = v;
252 return *this;
253}
254
255template <typename T, typename Storage>
256OrthogPolyImpl<T,Storage>&
257OrthogPolyImpl<T,Storage>::
258operator=(const OrthogPolyImpl<T,Storage>& x)
259{
260 expansion_ = x.expansion_;
261 quad_expansion_ = x.quad_expansion_;
262 th_ = x.th_;
263 return *this;
264}
265
266template <typename T, typename Storage>
267template <typename S>
268OrthogPolyImpl<T,Storage>&
269OrthogPolyImpl<T,Storage>::
270operator=(const Expr<S>& x)
271{
272 th_.makeOwnCopy();
273 expansion_ = x.expansion();
274 quad_expansion_ = Teuchos::rcp_dynamic_cast<quad_expansion_type>(expansion_);
275 th_->reset(expansion_->getBasis(), x.size());
276 expressionCopy(x);
277 return *this;
278}
279
280template <typename T, typename Storage>
281OrthogPolyImpl<T,Storage>&
282OrthogPolyImpl<T,Storage>::
283operator+=(const typename OrthogPolyImpl<T,Storage>::value_type& v)
284{
285 th_.makeOwnCopy();
286 expansion_->plusEqual(*th_, v);
287 return *this;
288}
289
290template <typename T, typename Storage>
291OrthogPolyImpl<T,Storage>&
292OrthogPolyImpl<T,Storage>::
293operator-=(const typename OrthogPolyImpl<T,Storage>::value_type& v)
294{
295 th_.makeOwnCopy();
296 expansion_->minusEqual(*th_, v);
297 return *this;
298}
299
300template <typename T, typename Storage>
301OrthogPolyImpl<T,Storage>&
302OrthogPolyImpl<T,Storage>::
303operator*=(const typename OrthogPolyImpl<T,Storage>::value_type& v)
304{
305 th_.makeOwnCopy();
306 expansion_->timesEqual(*th_, v);
307 return *this;
308}
309
310template <typename T, typename Storage>
311OrthogPolyImpl<T,Storage>&
312OrthogPolyImpl<T,Storage>::
313operator/=(const typename OrthogPolyImpl<T,Storage>::value_type& v)
314{
315 th_.makeOwnCopy();
316 expansion_->divideEqual(*th_, v);
317 return *this;
318}
319
320template <typename T, typename Storage>
321std::ostream&
322operator << (std::ostream& os, const OrthogPoly<T,Storage>& a)
323{
324 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
325
326 os << "[ ";
327
328 for (ordinal_type i=0; i<a.size(); i++) {
329 os << a.coeff(i) << " ";
330 }
331
332 os << "]\n";
333 return os;
334}
335
336template <typename T, typename Storage>
337std::istream&
338operator >> (std::istream& is, OrthogPoly<T,Storage>& a)
339{
340 typedef typename OrthogPoly<T,Storage>::ordinal_type ordinal_type;
341
342 // Read in the opening "["
343 char bracket;
344 is >> bracket;
345
346 for (ordinal_type i=0; i<a.size(); i++) {
347 is >> a.fastAccessCoeff(i);
348 }
349
350 // Read in the closing "]"
351
352 is >> bracket;
353 return is;
354}
355
356} // namespace PCE
357} // namespace Sacado
Orthogonal polynomial expansion class for constant (size 1) expansions.
std::istream & operator>>(std::istream &is, OrthogPoly< T, Storage > &a)
std::ostream & operator<<(std::ostream &os, const OrthogPoly< T, Storage > &a)
Top-level namespace for Stokhos classes and functions.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition csr_vector.h:265
KERNEL_PREFIX value_type operator()(tuple_type x) const