Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_OrthogPolyExpansionBaseImp.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 "Stokhos_DenseDirectDivisionExpansionStrategy.hpp"
43 #include "Stokhos_SPDDenseDirectDivisionExpansionStrategy.hpp"
44 #include "Stokhos_MeanBasedDivisionExpansionStrategy.hpp"
45 #include "Stokhos_GMRESDivisionExpansionStrategy.hpp"
46 #include "Stokhos_CGDivisionExpansionStrategy.hpp"
47 #include "Teuchos_Assert.hpp"
48 #include "Teuchos_TimeMonitor.hpp"
49 #include "Stokhos_DynamicArrayTraits.hpp"
50 
51 template <typename ordinal_type, typename value_type, typename node_type>
52 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
53 OrthogPolyExpansionBase(
54  const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
55  const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
56  const Teuchos::RCP<Teuchos::ParameterList>& params_) :
57  basis(basis_),
58  Cijk(Cijk_),
59  params(params_),
60  sz(basis->size())
61 {
62  if (params == Teuchos::null)
63  params = Teuchos::rcp(new Teuchos::ParameterList);
64 
65  // Create division strategy
66  std::string name = params->get("Division Strategy","Dense Direct");
67  double tol = params->get("Division Tolerance", 1e-6);
68  int prec_iter = params->get("prec_iter", 1);
69  int max_it = params->get("max_it_div", 50);
70  std::string prec = params->get("Prec Strategy","None");
71  int PrecNum;
72  if (prec == "None")
73  PrecNum=0;
74  else if (prec == "Diag")
75  PrecNum=1;
76  else if (prec == "Jacobi")
77  PrecNum=2;
78  else if (prec == "GS")
79  PrecNum=3;
80  else if (prec == "Schur")
81  PrecNum=4;
82  else
83  PrecNum=-1;
84  std::string linopt = params->get("Prec option", "whole");
85  int linear = 0;
86  if (linopt == "linear")
87  linear = 1;
88 
89  std::string schuropt = params->get("Schur option", "diag");
90  int diag = 0;
91  if (schuropt == "full")
92  diag = 1;
93 
94  int equil = params->get("Equilibrate", 0);
95 
96 
97  if (name == "Dense Direct")
98  division_strategy =
99  Teuchos::rcp(new DenseDirectDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk));
100  else if (name == "SPD Dense Direct")
101  division_strategy =
102  Teuchos::rcp(new SPDDenseDirectDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk));
103  else if (name == "Mean-Based")
104  division_strategy =
105  Teuchos::rcp(new MeanBasedDivisionExpansionStrategy<ordinal_type,value_type,node_type>());
106  else if (name == "GMRES")
107  division_strategy =
108  Teuchos::rcp(new GMRESDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk, prec_iter, tol, PrecNum, max_it, linear, diag, equil));
109  else if (name == "CG")
110  division_strategy =
111  Teuchos::rcp(new CGDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk, prec_iter, tol, PrecNum, max_it, linear, diag, equil));
112 
113 // TEUCHOS_TEST_FOR_EXCEPTION(
114 // true, std::logic_error,
115 // "Invalid division strategy name" << name);
116 }
117 
118 template <typename ordinal_type, typename value_type, typename node_type>
119 void
120 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
121 unaryMinus(
122  Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
123  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
124 {
125 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
126  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::unaryMinus(OPA)");
127 #endif
128  ordinal_type pc = a.size();
129  if (c.size() != pc)
130  c.resize(pc);
131 
132  value_type* cc = c.coeff();
133  const value_type* ca = a.coeff();
134 
135  for (ordinal_type i=0; i<pc; i++)
136  {
137  cc[i] = -ca[i];
138  }
139 }
140 template <typename ordinal_type, typename value_type, typename node_type>
141 void
142 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
143 plusEqual(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
144  const value_type& val)
145 {
146  c[0] += val;
147 }
148 
149 template <typename ordinal_type, typename value_type, typename node_type>
150 void
151 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
152 minusEqual(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
153  const value_type& val)
154 {
155  c[0] -= val;
156 }
157 
158 template <typename ordinal_type, typename value_type, typename node_type>
159 void
160 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
161 timesEqual(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
162  const value_type& val)
163 {
164 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
165  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::timesEqual(const)");
166 #endif
167  ordinal_type pc = c.size();
168  value_type* cc = c.coeff();
169  for (ordinal_type i=0; i<pc; i++)
170  cc[i] *= val;
171 }
172 
173 template <typename ordinal_type, typename value_type, typename node_type>
174 void
175 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
176 divideEqual(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
177  const value_type& val)
178 {
179 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
180  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::divideEqual(const)");
181 #endif
182  ordinal_type pc = c.size();
183  value_type* cc=c.coeff();
184  for (ordinal_type i=0; i<pc; i++){
185  cc[i] /= val;
186 }
187 }
188 template <typename ordinal_type, typename value_type, typename node_type>
189 void
190 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
191 plusEqual(
192  Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
193  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
194 {
195 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
196  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::plusEqual(OPA)");
197 #endif
198  ordinal_type xp = x.size();
199  if (c.size() < xp)
200  c.resize(xp);
201 
202  value_type* cc = c.coeff();
203  const value_type* xc = x.coeff();
204  for (ordinal_type i=0; i<xp; i++)
205  cc[i] += xc[i];
206 }
207 
208 template <typename ordinal_type, typename value_type, typename node_type>
209 void
210 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
211 minusEqual(
212  Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
213  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
214 {
215 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
216  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::minusEqual(OPA)");
217 #endif
218  ordinal_type xp = x.size();
219  if (c.size() < xp)
220  c.resize(xp);
221 
222  value_type* cc = c.coeff();
223  const value_type* xc = x.coeff();
224  for (ordinal_type i=0; i<xp; i++)
225  cc[i] -= xc[i];
226 }
227 
228 template <typename ordinal_type, typename value_type, typename node_type>
229 void
230 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
231 timesEqual(
232  Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
233  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
234 {
235 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
236  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::timesEqual(OPA)");
237 #endif
238  ordinal_type p = c.size();
239  ordinal_type xp = x.size();
240  ordinal_type pc;
241  if (p > 1 && xp > 1)
242  pc = sz;
243  else
244  pc = p*xp;
245  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
246  "Stokhos::OrthogPolyExpansionBase::timesEqual()" <<
247  ": Expansion size (" << sz <<
248  ") is too small for computation.");
249  if (c.size() != pc)
250  c.resize(pc);
251 
252  value_type* cc = c.coeff();
253  const value_type* xc = x.coeff();
254 
255  if (p > 1 && xp > 1) {
256  // Copy c coefficients into temporary array
257  value_type* tc = Stokhos::ds_array<value_type>::get_and_fill(cc,p);
258 
259  typename Cijk_type::i_iterator i_begin = Cijk->i_begin();
260  typename Cijk_type::i_iterator i_end = Cijk->i_end();
261  if (pc < Cijk->num_i())
262  i_end = Cijk->find_i(pc);
263  ordinal_type k_lim = p;
264  ordinal_type j_lim = xp;
265  const value_type* kc = tc;
266  const value_type* jc = xc;
267  if (xp < p) {
268  k_lim = xp;
269  j_lim = p;
270  kc = xc;
271  jc = tc;
272  }
273 
274  value_type tmp, cijk;
275  ordinal_type i,j,k;
276  for (typename Cijk_type::i_iterator i_it=i_begin; i_it!=i_end; ++i_it) {
277  i = index(i_it);
278  tmp = value_type(0.0);
279  for (typename Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
280  k_it != Cijk->k_end(i_it); ++k_it) {
281  k = index(k_it);
282  if (k < k_lim) {
283  for (typename Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
284  j_it != Cijk->j_end(k_it); ++j_it) {
285  j = index(j_it);
286  cijk = value(j_it);
287  if (j < j_lim)
288  tmp += cijk*kc[k]*jc[j];
289  }
290  }
291  }
292  cc[i] = tmp / basis->norm_squared(i);
293  }
294 
295  Stokhos::ds_array<value_type>::destroy_and_release(tc, p);
296  }
297  else if (p > 1) {
298  for (ordinal_type i=0; i<p; i++)
299  cc[i] *= xc[0];
300  }
301  else if (xp > 1) {
302  for (ordinal_type i=1; i<xp; i++)
303  cc[i] = cc[0]*xc[i];
304  cc[0] *= xc[0];
305  }
306  else {
307  cc[0] *= xc[0];
308  }
309 }
310 
311 template <typename ordinal_type, typename value_type, typename node_type>
312 void
313 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
314 divideEqual(
315  Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
316  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type >& x)
317 {
318  division_strategy->divide(c, 1.0, c, x, 0.0);
319 }
320 
321 template <typename ordinal_type, typename value_type, typename node_type>
322 void
323 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
324 plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
325  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
326  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
327 {
328 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
329  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::plus(OPA,OPA)");
330 #endif
331  ordinal_type pa = a.size();
332  ordinal_type pb = b.size();
333  ordinal_type pc = pa > pb ? pa : pb;
334  if (c.size() != pc)
335  c.resize(pc);
336 
337  const value_type* ca = a.coeff();
338  const value_type* cb = b.coeff();
339  value_type* cc = c.coeff();
340 
341  if (pa > pb) {
342  for (ordinal_type i=0; i<pb; i++)
343  cc[i] = ca[i] + cb[i];
344  for (ordinal_type i=pb; i<pc; i++)
345  cc[i] = ca[i];
346  }
347  else {
348  for (ordinal_type i=0; i<pa; i++)
349  cc[i] = ca[i] + cb[i];
350  for (ordinal_type i=pa; i<pc; i++)
351  cc[i] = cb[i];
352  }
353 }
354 
355 template <typename ordinal_type, typename value_type, typename node_type>
356 void
357 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
358 plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
359  const value_type& a,
360  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
361 {
362 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
363  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::plus(const,OPA)");
364 #endif
365  ordinal_type pc = b.size();
366  if (c.size() != pc)
367  c.resize(pc);
368 
369  const value_type* cb = b.coeff();
370  value_type* cc = c.coeff();
371 
372  cc[0] = a + cb[0];
373  for (ordinal_type i=1; i<pc; i++)
374  cc[i] = cb[i];
375 }
376 
377 template <typename ordinal_type, typename value_type, typename node_type>
378 void
379 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
380 plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
381  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
382  const value_type& b)
383 {
384 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
385  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::plus(OPA,const)");
386 #endif
387  ordinal_type pc = a.size();
388  if (c.size() != pc)
389  c.resize(pc);
390 
391  const value_type* ca = a.coeff();
392  value_type* cc = c.coeff();
393 
394  cc[0] = ca[0] + b;
395  for (ordinal_type i=1; i<pc; i++)
396  cc[i] = ca[i];
397 }
398 
399 template <typename ordinal_type, typename value_type, typename node_type>
400 void
401 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
402 minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
403  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
404  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
405 {
406 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
407  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::minus(OPA,OPA)");
408 #endif
409  ordinal_type pa = a.size();
410  ordinal_type pb = b.size();
411  ordinal_type pc = pa > pb ? pa : pb;
412  if (c.size() != pc)
413  c.resize(pc);
414 
415  const value_type* ca = a.coeff();
416  const value_type* cb = b.coeff();
417  value_type* cc = c.coeff();
418 
419  if (pa > pb) {
420  for (ordinal_type i=0; i<pb; i++)
421  cc[i] = ca[i] - cb[i];
422  for (ordinal_type i=pb; i<pc; i++)
423  cc[i] = ca[i];
424  }
425  else {
426  for (ordinal_type i=0; i<pa; i++)
427  cc[i] = ca[i] - cb[i];
428  for (ordinal_type i=pa; i<pc; i++)
429  cc[i] = -cb[i];
430  }
431 }
432 
433 template <typename ordinal_type, typename value_type, typename node_type>
434 void
435 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
436 minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
437  const value_type& a,
438  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
439 {
440 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
441  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::minus(const,OPA)");
442 #endif
443  ordinal_type pc = b.size();
444  if (c.size() != pc)
445  c.resize(pc);
446 
447  const value_type* cb = b.coeff();
448  value_type* cc = c.coeff();
449 
450  cc[0] = a - cb[0];
451  for (ordinal_type i=1; i<pc; i++)
452  cc[i] = -cb[i];
453 }
454 
455 template <typename ordinal_type, typename value_type, typename node_type>
456 void
457 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
458 minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
459  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
460  const value_type& b)
461 {
462 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
463  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::minus(OPA,const)");
464 #endif
465  ordinal_type pc = a.size();
466  if (c.size() != pc)
467  c.resize(pc);
468 
469  const value_type* ca = a.coeff();
470  value_type* cc = c.coeff();
471 
472  cc[0] = ca[0] - b;
473  for (ordinal_type i=1; i<pc; i++)
474  cc[i] = ca[i];
475 }
476 
477 template <typename ordinal_type, typename value_type, typename node_type>
478 void
479 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
480 times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
481  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
482  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
483 {
484 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
485  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::times(OPA,OPA)");
486 #endif
487  ordinal_type pa = a.size();
488  ordinal_type pb = b.size();
489  ordinal_type pc;
490  if (pa > 1 && pb > 1)
491  pc = sz;
492  else
493  pc = pa*pb;
494  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
495  "Stokhos::OrthogPolyExpansionBase::times()" <<
496  ": Expansion size (" << sz <<
497  ") is too small for computation.");
498  if (c.size() != pc)
499  c.resize(pc);
500 
501  const value_type* ca = a.coeff();
502  const value_type* cb = b.coeff();
503  value_type* cc = c.coeff();
504 
505  if (pa > 1 && pb > 1) {
506  typename Cijk_type::i_iterator i_begin = Cijk->i_begin();
507  typename Cijk_type::i_iterator i_end = Cijk->i_end();
508  if (pc < Cijk->num_i())
509  i_end = Cijk->find_i(pc);
510  ordinal_type k_lim = pa;
511  ordinal_type j_lim = pb;
512  const value_type* kc = ca;
513  const value_type* jc = cb;
514  if (pb < pa) {
515  k_lim = pb;
516  j_lim = pa;
517  kc = cb;
518  jc = ca;
519  }
520  value_type tmp, cijk;
521  ordinal_type i,j,k;
522  for (typename Cijk_type::i_iterator i_it=i_begin; i_it!=i_end; ++i_it) {
523  i = index(i_it);
524  tmp = value_type(0.0);
525  for (typename Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
526  k_it != Cijk->k_end(i_it); ++k_it) {
527  k = index(k_it);
528  if (k < k_lim) {
529  for (typename Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
530  j_it != Cijk->j_end(k_it); ++j_it) {
531  j = index(j_it);
532  cijk = value(j_it);
533  if (j < j_lim)
534  tmp += cijk*kc[k]*jc[j];
535  }
536  }
537  }
538  cc[i] = tmp / basis->norm_squared(i);
539  }
540  }
541  else if (pa > 1) {
542  for (ordinal_type i=0; i<pc; i++)
543  cc[i] = ca[i]*cb[0];
544  }
545  else if (pb > 1) {
546  for (ordinal_type i=0; i<pc; i++)
547  cc[i] = ca[0]*cb[i];
548  }
549  else {
550  cc[0] = ca[0]*cb[0];
551  }
552 }
553 
554 template <typename ordinal_type, typename value_type, typename node_type>
555 void
556 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
557 times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
558  const value_type& a,
559  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
560 {
561 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
562  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::times(const,OPA)");
563 #endif
564  ordinal_type pc = b.size();
565  if (c.size() != pc)
566  c.resize(pc);
567 
568  const value_type* cb = b.coeff();
569  value_type* cc = c.coeff();
570 
571  for (ordinal_type i=0; i<pc; i++)
572  cc[i] = a*cb[i];
573 }
574 
575 template <typename ordinal_type, typename value_type, typename node_type>
576 void
577 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
578 times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
579  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
580  const value_type& b)
581 {
582 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
583  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::times(OPA,const)");
584 #endif
585  ordinal_type pc = a.size();
586  if (c.size() != pc)
587  c.resize(pc);
588 
589  const value_type* ca = a.coeff();
590  value_type* cc = c.coeff();
591 
592  for (ordinal_type i=0; i<pc; i++)
593  cc[i] = ca[i]*b;
594 }
595 
596 template <typename ordinal_type, typename value_type, typename node_type>
597 void
598 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
599 divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
600  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
601  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
602 {
603  division_strategy->divide(c, 1.0, a, b, 0.0);
604 }
605 
606 template <typename ordinal_type, typename value_type, typename node_type>
607 void
608 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
609 divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
610  const value_type& a,
611  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
612 {
613  Stokhos::OrthogPolyApprox<ordinal_type,value_type> aa(b.basis(), 1, &a);
614  division_strategy->divide(c, 1.0, aa, b, 0.0);
615 }
616 
617 template <typename ordinal_type, typename value_type, typename node_type>
618 void
619 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
620 divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
621  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
622  const value_type& b)
623 {
624 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
625  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::divide(OPA,const)");
626 #endif
627  ordinal_type pc = a.size();
628  if (c.size() != pc)
629  c.resize(pc);
630 
631  const value_type* ca = a.coeff();
632  value_type* cc = c.coeff();
633 
634  for (ordinal_type i=0; i<pc; i++)
635  cc[i] = ca[i]/b;
636 }
637 
638 template <typename ordinal_type, typename value_type, typename node_type>
639 void
640 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
641 fabs(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
642  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
643 {
644 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
645  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::fabs(OPA)");
646 #endif
647  c.init(0.0);
648  c[0] = a.two_norm();
649 }
650 
651 template <typename ordinal_type, typename value_type, typename node_type>
652 void
653 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
654 abs(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
655  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
656 {
657 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
658  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::abs(OPA)");
659 #endif
660  c.init(0.0);
661  c[0] = a.two_norm();
662 }
663 
664 template <typename ordinal_type, typename value_type, typename node_type>
665 void
666 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
667 max(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
668  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
669  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
670 {
671 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
672  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::max(OPA,OPA)");
673 #endif
674  if (a.two_norm() >= b.two_norm())
675  c = a;
676  else
677  c = b;
678 }
679 
680 template <typename ordinal_type, typename value_type, typename node_type>
681 void
682 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
683 max(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
684  const value_type& a,
685  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
686 {
687 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
688  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::max(const,OPA)");
689 #endif
690  if (a >= b.two_norm()) {
691  c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
692  c[0] = a;
693  }
694  else
695  c = b;
696 }
697 
698 template <typename ordinal_type, typename value_type, typename node_type>
699 void
700 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
701 max(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
702  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
703  const value_type& b)
704 {
705 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
706  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::max(OPA,const)");
707 #endif
708  if (a.two_norm() >= b)
709  c = a;
710  else {
711  c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
712  c[0] = b;
713  }
714 }
715 
716 template <typename ordinal_type, typename value_type, typename node_type>
717 void
718 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
719 min(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
720  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
721  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
722 {
723 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
724  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::min(OPA,OPA)");
725 #endif
726  if (a.two_norm() <= b.two_norm())
727  c = a;
728  else
729  c = b;
730 }
731 
732 template <typename ordinal_type, typename value_type, typename node_type>
733 void
734 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
735 min(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
736  const value_type& a,
737  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
738 {
739 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
740  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::min(const,OPA)");
741 #endif
742  if (a <= b.two_norm()) {
743  c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
744  c[0] = a;
745  }
746  else
747  c = b;
748 }
749 
750 template <typename ordinal_type, typename value_type, typename node_type>
751 void
752 Stokhos::OrthogPolyExpansionBase<ordinal_type, value_type, node_type>::
753 min(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
754  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
755  const value_type& b)
756 {
757 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
758  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::OrthogPolyExpansionBase::min(OPA,const)");
759 #endif
760  if (a.two_norm() <= b)
761  c = a;
762  else {
763  c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
764  c[0] = b;
765  }
766 }
Stokhos::OrthogPolyExpansionBase::divideEqual
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:176
Stokhos::OrthogPolyExpansionBase::basis
Teuchos::RCP< const OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
Definition: Stokhos_OrthogPolyExpansionBase.hpp:186
TEUCHOS_FUNC_TIME_MONITOR
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Stokhos::OrthogPolyExpansionBase::params
Teuchos::RCP< Teuchos::ParameterList > params
Parameter list.
Definition: Stokhos_OrthogPolyExpansionBase.hpp:195
Stokhos::OrthogPolyApprox::resize
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Definition: Stokhos_OrthogPolyApproxImp.hpp:152
Stokhos::OrthogPolyExpansionBase::minusEqual
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:152
Stokhos::CGDivisionExpansionStrategy
Strategy interface for computing PCE of a/b using only b[0].
Definition: Stokhos_CGDivisionExpansionStrategy.hpp:72
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition: Stokhos_Sparse3Tensor.hpp:56
Teuchos::ParameterList::get
T & get(ParameterList &l, const std::string &name)
Stokhos_MeanBasedDivisionExpansionStrategy.hpp
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Stokhos::OrthogPolyApprox::init
void init(const value_type &v)
Initialize coefficients to value.
Definition: Stokhos_OrthogPolyApproxImp.hpp:106
Stokhos::OrthogPolyExpansionBase::divide
void divide(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:599
Stokhos::OrthogPolyExpansionBase::division_strategy
Teuchos::RCP< Stokhos::DivisionExpansionStrategy< ordinal_type, value_type, node_type > > division_strategy
Division expansion strategy.
Definition: Stokhos_OrthogPolyExpansionBase.hpp:198
Stokhos::ds_array::destroy_and_release
static void destroy_and_release(T *m, int sz)
Destroy array elements and release memory.
Definition: Stokhos_DynamicArrayTraits.hpp:115
Stokhos_GMRESDivisionExpansionStrategy.hpp
Stokhos::OrthogPolyExpansionBase::timesEqual
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:161
Stokhos::OrthogPolyExpansionBase::max
void max(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:667
Stokhos::OrthogPolyApprox::coeff
pointer coeff()
Return coefficient array.
Definition: Stokhos_OrthogPolyApproxImp.hpp:168
Stokhos::OrthogPolyApprox::two_norm
value_type two_norm() const
Compute the two-norm of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:292
Stokhos::OrthogPolyExpansionBase::plusEqual
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:143
linear
Definition: experimental/linear2d_diffusion_pce.cpp:146
Stokhos::ds_array::get_and_fill
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.
Definition: Stokhos_DynamicArrayTraits.hpp:82
Stokhos::SparseArrayIterator
Bi-directional iterator for traversing a sparse array.
Definition: Stokhos_SparseArray.hpp:56
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Stokhos::OrthogPolyBasis
Abstract base class for multivariate orthogonal polynomials.
Definition: Stokhos_OrthogPolyBasis.hpp:74
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Stokhos::OrthogPolyExpansionBase::min
void min(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:719
Stokhos::OrthogPolyApprox::basis
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Definition: Stokhos_OrthogPolyApproxImp.hpp:130
Stokhos::OrthogPolyExpansionBase::plus
void plus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:324
Stokhos::OrthogPolyExpansionBase::minus
void minus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:402
Stokhos::OrthogPolyExpansionBase::times
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:480
Stokhos_DenseDirectDivisionExpansionStrategy.hpp
Kokkos::cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:217
Stokhos::GMRESDivisionExpansionStrategy
Strategy interface for computing PCE of a/b using only b[0].
Definition: Stokhos_GMRESDivisionExpansionStrategy.hpp:69
Stokhos::OrthogPolyApprox
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Definition: Stokhos_OrthogPolyApprox.hpp:63
val
expr val()
Stokhos_DynamicArrayTraits.hpp
Stokhos::OrthogPolyExpansionBase::abs
void abs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:654
Stokhos::OrthogPolyApprox::size
ordinal_type size() const
Return size.
Definition: Stokhos_OrthogPolyApproxImp.hpp:160
diag
Definition: experimental/linear2d_diffusion_pce.cpp:140
Teuchos_TimeMonitor.hpp
Stokhos::SPDDenseDirectDivisionExpansionStrategy
Strategy interface for computing PCE of a/b using only b[0].
Definition: Stokhos_SPDDenseDirectDivisionExpansionStrategy.hpp:63
Stokhos_CGDivisionExpansionStrategy.hpp
Stokhos_SPDDenseDirectDivisionExpansionStrategy.hpp
Stokhos::OrthogPolyExpansionBase::Cijk
Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > Cijk
Triple-product tensor.
Definition: Stokhos_OrthogPolyExpansionBase.hpp:192
Stokhos::MeanBasedDivisionExpansionStrategy
Strategy interface for computing PCE of a/b using only b[0].
Definition: Stokhos_MeanBasedDivisionExpansionStrategy.hpp:57
Teuchos_Assert.hpp
Stokhos::OrthogPolyExpansionBase::unaryMinus
void unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:121
Stokhos::OrthogPolyExpansionBase::fabs
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:641
Stokhos::DenseDirectDivisionExpansionStrategy
Strategy interface for computing PCE of a/b using only b[0].
Definition: Stokhos_DenseDirectDivisionExpansionStrategy.hpp:63
Stokhos::OrthogPolyExpansionBase::OrthogPolyExpansionBase
OrthogPolyExpansionBase(const Teuchos::RCP< const OrthogPolyBasis< ordinal_type, value_type > > &basis, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Constructor.
Definition: Stokhos_OrthogPolyExpansionBaseImp.hpp:53
Generated by   doxygen 1.8.5