47#include "Teuchos_Assert.hpp"
48#include "Teuchos_TimeMonitor.hpp"
51template <
typename ordinal_type,
typename value_type,
typename node_type>
56 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
62 if (params == Teuchos::null)
63 params = Teuchos::rcp(
new Teuchos::ParameterList);
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");
74 else if (prec ==
"Diag")
76 else if (prec ==
"Jacobi")
78 else if (prec ==
"GS")
80 else if (prec ==
"Schur")
84 std::string linopt = params->get(
"Prec option",
"whole");
86 if (linopt ==
"linear")
89 std::string schuropt = params->get(
"Schur option",
"diag");
91 if (schuropt ==
"full")
94 int equil = params->get(
"Equilibrate", 0);
97 if (name ==
"Dense Direct")
99 Teuchos::rcp(
new DenseDirectDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk));
100 else if (name ==
"SPD Dense Direct")
102 Teuchos::rcp(
new SPDDenseDirectDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk));
103 else if (name ==
"Mean-Based")
105 Teuchos::rcp(
new MeanBasedDivisionExpansionStrategy<ordinal_type,value_type,node_type>());
106 else if (name ==
"GMRES")
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")
111 Teuchos::rcp(
new CGDivisionExpansionStrategy<ordinal_type,value_type,node_type>(this->basis, this->Cijk, prec_iter, tol, PrecNum, max_it,
linear,
diag, equil));
118template <
typename ordinal_type,
typename value_type,
typename node_type>
125#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
126 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::unaryMinus(OPA)");
128 ordinal_type pc = a.
size();
133 const value_type* ca = a.
coeff();
135 for (ordinal_type i=0; i<pc; i++)
140template <
typename ordinal_type,
typename value_type,
typename node_type>
149template <
typename ordinal_type,
typename value_type,
typename node_type>
153 const value_type&
val)
158template <
typename ordinal_type,
typename value_type,
typename node_type>
164#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
165 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::timesEqual(const)");
167 ordinal_type pc = c.
size();
169 for (ordinal_type i=0; i<pc; i++)
173template <
typename ordinal_type,
typename value_type,
typename node_type>
177 const value_type&
val)
179#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
180 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divideEqual(const)");
182 ordinal_type pc = c.
size();
183 value_type* cc=c.
coeff();
184 for (ordinal_type i=0; i<pc; i++){
188template <
typename ordinal_type,
typename value_type,
typename node_type>
195#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
196 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::plusEqual(OPA)");
198 ordinal_type xp = x.size();
202 value_type* cc = c.
coeff();
203 const value_type* xc = x.coeff();
204 for (ordinal_type i=0; i<xp; i++)
208template <
typename ordinal_type,
typename value_type,
typename node_type>
215#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
216 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::minusEqual(OPA)");
218 ordinal_type xp = x.size();
222 value_type* cc = c.
coeff();
223 const value_type* xc = x.coeff();
224 for (ordinal_type i=0; i<xp; i++)
228template <
typename ordinal_type,
typename value_type,
typename node_type>
235#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
236 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::timesEqual(OPA)");
238 ordinal_type p = c.
size();
239 ordinal_type xp = x.size();
245 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
246 "Stokhos::OrthogPolyExpansionBase::timesEqual()" <<
247 ": Expansion size (" << sz <<
248 ") is too small for computation.");
252 value_type* cc = c.
coeff();
253 const value_type* xc = x.coeff();
255 if (p > 1 && xp > 1) {
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;
274 value_type tmp, cijk;
278 tmp = value_type(0.0);
280 k_it != Cijk->k_end(i_it); ++k_it) {
284 j_it != Cijk->j_end(k_it); ++j_it) {
288 tmp += cijk*kc[k]*jc[
j];
292 cc[i] = tmp / basis->norm_squared(i);
298 for (ordinal_type i=0; i<p; i++)
302 for (ordinal_type i=1; i<xp; i++)
311template <
typename ordinal_type,
typename value_type,
typename node_type>
318 division_strategy->divide(c, 1.0, c, x, 0.0);
321template <
typename ordinal_type,
typename value_type,
typename node_type>
328#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
329 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::plus(OPA,OPA)");
331 ordinal_type pa = a.
size();
332 ordinal_type pb = b.
size();
333 ordinal_type pc = pa > pb ? pa : pb;
337 const value_type* ca = a.
coeff();
338 const value_type* cb = b.
coeff();
339 value_type* cc = c.
coeff();
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++)
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++)
355template <
typename ordinal_type,
typename value_type,
typename node_type>
362#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
363 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::plus(const,OPA)");
365 ordinal_type pc = b.
size();
369 const value_type* cb = b.
coeff();
370 value_type* cc = c.
coeff();
373 for (ordinal_type i=1; i<pc; i++)
377template <
typename ordinal_type,
typename value_type,
typename node_type>
384#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
385 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::plus(OPA,const)");
387 ordinal_type pc = a.
size();
391 const value_type* ca = a.
coeff();
392 value_type* cc = c.
coeff();
395 for (ordinal_type i=1; i<pc; i++)
399template <
typename ordinal_type,
typename value_type,
typename node_type>
406#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
407 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::minus(OPA,OPA)");
409 ordinal_type pa = a.
size();
410 ordinal_type pb = b.
size();
411 ordinal_type pc = pa > pb ? pa : pb;
415 const value_type* ca = a.
coeff();
416 const value_type* cb = b.
coeff();
417 value_type* cc = c.
coeff();
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++)
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++)
433template <
typename ordinal_type,
typename value_type,
typename node_type>
440#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
441 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::minus(const,OPA)");
443 ordinal_type pc = b.
size();
447 const value_type* cb = b.
coeff();
448 value_type* cc = c.
coeff();
451 for (ordinal_type i=1; i<pc; i++)
455template <
typename ordinal_type,
typename value_type,
typename node_type>
462#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
463 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::minus(OPA,const)");
465 ordinal_type pc = a.
size();
469 const value_type* ca = a.
coeff();
470 value_type* cc = c.
coeff();
473 for (ordinal_type i=1; i<pc; i++)
477template <
typename ordinal_type,
typename value_type,
typename node_type>
484#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
485 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::times(OPA,OPA)");
487 ordinal_type pa = a.
size();
488 ordinal_type pb = b.
size();
490 if (pa > 1 && pb > 1)
494 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
495 "Stokhos::OrthogPolyExpansionBase::times()" <<
496 ": Expansion size (" << sz <<
497 ") is too small for computation.");
501 const value_type* ca = a.
coeff();
502 const value_type* cb = b.
coeff();
503 value_type* cc = c.
coeff();
505 if (pa > 1 && pb > 1) {
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;
520 value_type tmp, cijk;
524 tmp = value_type(0.0);
526 k_it != Cijk->k_end(i_it); ++k_it) {
530 j_it != Cijk->j_end(k_it); ++j_it) {
534 tmp += cijk*kc[k]*jc[
j];
538 cc[i] = tmp / basis->norm_squared(i);
542 for (ordinal_type i=0; i<pc; i++)
546 for (ordinal_type i=0; i<pc; i++)
554template <
typename ordinal_type,
typename value_type,
typename node_type>
561#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
562 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::times(const,OPA)");
564 ordinal_type pc = b.
size();
568 const value_type* cb = b.
coeff();
569 value_type* cc = c.
coeff();
571 for (ordinal_type i=0; i<pc; i++)
575template <
typename ordinal_type,
typename value_type,
typename node_type>
582#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
583 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::times(OPA,const)");
585 ordinal_type pc = a.
size();
589 const value_type* ca = a.
coeff();
590 value_type* cc = c.
coeff();
592 for (ordinal_type i=0; i<pc; i++)
596template <
typename ordinal_type,
typename value_type,
typename node_type>
603 division_strategy->divide(c, 1.0, a, b, 0.0);
606template <
typename ordinal_type,
typename value_type,
typename node_type>
614 division_strategy->divide(c, 1.0, aa, b, 0.0);
617template <
typename ordinal_type,
typename value_type,
typename node_type>
624#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
625 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divide(OPA,const)");
627 ordinal_type pc = a.
size();
631 const value_type* ca = a.
coeff();
632 value_type* cc = c.
coeff();
634 for (ordinal_type i=0; i<pc; i++)
638template <
typename ordinal_type,
typename value_type,
typename node_type>
644#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
645 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::fabs(OPA)");
651template <
typename ordinal_type,
typename value_type,
typename node_type>
657#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
658 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::abs(OPA)");
664template <
typename ordinal_type,
typename value_type,
typename node_type>
671#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
672 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::max(OPA,OPA)");
680template <
typename ordinal_type,
typename value_type,
typename node_type>
687#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
688 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::max(const,OPA)");
691 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
698template <
typename ordinal_type,
typename value_type,
typename node_type>
705#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
706 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::max(OPA,const)");
711 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
716template <
typename ordinal_type,
typename value_type,
typename node_type>
723#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
724 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::min(OPA,OPA)");
732template <
typename ordinal_type,
typename value_type,
typename node_type>
739#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
740 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::min(const,OPA)");
743 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
750template <
typename ordinal_type,
typename value_type,
typename node_type>
757#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
758 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::min(OPA,const)");
763 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void init(const value_type &v)
Initialize coefficients to value.
value_type two_norm() const
Compute the two-norm of expansion.
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
pointer coeff()
Return coefficient array.
ordinal_type size() const
Return size.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Abstract base class for multivariate orthogonal polynomials.
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
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)
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 > ¶ms=Teuchos::null)
Constructor.
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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)
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
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)
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)
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)
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)
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void abs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
j_sparse_array::const_iterator ikj_iterator
Iterator for looping over j entries given i and k.
ikj_sparse_array::const_iterator i_iterator
Iterator for looping over i entries.
Bi-directional iterator for traversing a sparse array.
static void destroy_and_release(T *m, int sz)
Destroy array elements and release memory.
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.