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