47 #include "Teuchos_Assert.hpp" 48 #include "Teuchos_TimeMonitor.hpp" 51 template <
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));
118 template <
typename ordinal_type,
typename value_type,
typename node_type>
125 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 126 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::unaryMinus(OPA)");
140 template <
typename ordinal_type,
typename value_type,
typename node_type>
149 template <
typename ordinal_type,
typename value_type,
typename node_type>
158 template <
typename ordinal_type,
typename value_type,
typename node_type>
164 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 165 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::timesEqual(const)");
173 template <
typename ordinal_type,
typename value_type,
typename node_type>
179 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 180 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divideEqual(const)");
188 template <
typename ordinal_type,
typename value_type,
typename node_type>
195 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 196 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::plusEqual(OPA)");
208 template <
typename ordinal_type,
typename value_type,
typename node_type>
215 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 216 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::minusEqual(OPA)");
228 template <
typename ordinal_type,
typename value_type,
typename node_type>
235 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 236 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::timesEqual(OPA)");
245 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
246 "Stokhos::OrthogPolyExpansionBase::timesEqual()" <<
247 ": Expansion size (" << sz <<
248 ") is too small for computation.");
255 if (p > 1 && xp > 1) {
261 if (pc < Cijk->num_i())
262 i_end = Cijk->find_i(pc);
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);
311 template <
typename ordinal_type,
typename value_type,
typename node_type>
318 division_strategy->divide(c, 1.0, c,
x, 0.0);
321 template <
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)");
343 cc[i] = ca[i] + cb[i];
349 cc[i] = ca[i] + cb[i];
355 template <
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)");
377 template <
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)");
399 template <
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)");
421 cc[i] = ca[i] - cb[i];
427 cc[i] = ca[i] - cb[i];
433 template <
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)");
455 template <
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)");
477 template <
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)");
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.");
505 if (pa > 1 && pb > 1) {
508 if (pc < Cijk->num_i())
509 i_end = Cijk->find_i(pc);
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);
554 template <
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)");
575 template <
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)");
596 template <
typename ordinal_type,
typename value_type,
typename node_type>
603 division_strategy->divide(c, 1.0, a, b, 0.0);
606 template <
typename ordinal_type,
typename value_type,
typename node_type>
614 division_strategy->divide(c, 1.0, aa, b, 0.0);
617 template <
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)");
638 template <
typename ordinal_type,
typename value_type,
typename node_type>
644 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 645 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::fabs(OPA)");
651 template <
typename ordinal_type,
typename value_type,
typename node_type>
657 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 658 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::abs(OPA)");
664 template <
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)");
680 template <
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)");
698 template <
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)");
716 template <
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)");
732 template <
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)");
750 template <
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)");
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
void init(const value_type &v)
Initialize coefficients to value.
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)
static void destroy_and_release(T *m, int sz)
Destroy array elements and release memory.
void timesEqual(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)
pointer coeff()
Return coefficient array.
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
ordinal_type size() const
Return size.
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.
Bi-directional iterator for traversing a sparse array.
Abstract base class for multivariate orthogonal polynomials.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
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 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 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 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)
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)
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type two_norm() const
Compute the two-norm of expansion.
SparseArrayIterator< index_iterator, value_iterator >::value_reference value(const SparseArrayIterator< index_iterator, value_iterator > &it)
void abs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
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.