44#include "Teuchos_Assert.hpp"
46#include "Teuchos_TimeMonitor.hpp"
47#include "Teuchos_Tuple.hpp"
49template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
54 const Teuchos::RCP<
const PseudoSpectralOperator<ordinal_type, value_type, point_compare_type> >& ps_op_,
55 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
58 sz(ps_op->coeff_size()),
59 nqp(ps_op->point_size()),
64 Teuchos::RCP<Teuchos::ParameterList> params = params_;
65 if (params == Teuchos::null)
66 params = Teuchos::rcp(
new Teuchos::ParameterList);
67 use_quad_for_times = params->get(
"Use Quadrature for Times",
false);
68 use_quad_for_division = params->get(
"Use Quadrature for Division",
true);
71template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
72template <
typename FuncT>
79 ordinal_type pa = a.
size();
94#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
95 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- Unary Polynomial Evaluation");
99 SDV a_sdv(Teuchos::View,
const_cast<value_type*
>(a.
coeff()), pa);
100 ps_op->transformPCE2QP(1.0, a_sdv, avals, 0.0,
false);
105#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
106 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- Unary Function Evaluation");
110 for (ordinal_type qp=0; qp<nqp; qp++) {
111 fvals[qp] = func(avals[qp]);
117#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
118 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- Unary Polynomial Integration");
122 SDV c_sdv(Teuchos::View, c.
coeff(), pc);
123 ps_op->transformQP2PCE(1.0, fvals, c_sdv, 0.0,
false);
128template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
129template <
typename FuncT>
137 ordinal_type pa = a.
size();
138 ordinal_type pb = b.
size();
140 if (pa == 1 && pb == 1)
148 c[0] = func(a[0], b[0]);
153#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
154 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PP Binary Polynomial Evaluation");
158 SDV a_sdv(Teuchos::View,
const_cast<value_type*
>(a.
coeff()), pa);
159 SDV b_sdv(Teuchos::View,
const_cast<value_type*
>(b.
coeff()), pb);
160 ps_op->transformPCE2QP(1.0, a_sdv, avals, 0.0,
false);
161 ps_op->transformPCE2QP(1.0, b_sdv, bvals, 0.0,
false);
166#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
167 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PP Binary Function Evaluation");
171 for (ordinal_type qp=0; qp<nqp; qp++)
172 fvals[qp] = func(avals[qp], bvals[qp]);
177#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
178 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PP Binary Polynomial Integration");
182 SDV c_sdv(Teuchos::View, c.
coeff(), pc);
183 ps_op->transformQP2PCE(1.0, fvals, c_sdv, 0.0,
false);
188template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
189template <
typename FuncT>
197 ordinal_type pb = b.
size();
207 c[0] = func(a, b[0]);
212#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
213 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- CP Binary Polynomial Evaluation");
217 SDV b_sdv(Teuchos::View,
const_cast<value_type*
>(b.
coeff()), pb);
218 ps_op->transformPCE2QP(1.0, b_sdv, bvals, 0.0,
false);
223#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
224 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- CP Binary Function Evaluation");
228 for (ordinal_type qp=0; qp<nqp; qp++)
229 fvals[qp] = func(a, bvals[qp]);
234#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
235 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- CP Binary Polynomial Integration");
239 SDV c_sdv(Teuchos::View, c.
coeff(), pc);
240 ps_op->transformQP2PCE(1.0, fvals, c_sdv, 0.0,
false);
245template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
246template <
typename FuncT>
254 ordinal_type pa = a.
size();
264 c[0] = func(a[0], b);
269#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
270 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PC Binary Polynomial Evaluation");
274 SDV a_sdv(Teuchos::View,
const_cast<value_type*
>(a.
coeff()), pa);
275 ps_op->transformPCE2QP(1.0, a_sdv, avals, 0.0,
false);
280#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
281 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PC Binary Function Evaluation");
285 for (ordinal_type qp=0; qp<nqp; qp++)
286 fvals[qp] = func(avals[qp], b);
291#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
292 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- PC Binary Polynomial Integration");
296 SDV c_sdv(Teuchos::View, c.
coeff(), pc);
297 ps_op->transformQP2PCE(1.0, fvals, c_sdv, 0.0,
false);
302template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
303template <
typename FuncT>
310 const int N = FuncT::N;
311 bool is_constant =
true;
312 for (
int i=0; i<N; i++) {
313 if (na[i]->size() > 1) {
328 for (
int i=0; i<N; i++)
329 val[i] = (*na[i])[0];
334 if (N >= navals.size())
336 if (navals[N].size() != N) {
338 for (
int i=0; i<N; i++)
339 navals[N][i].resize(nqp);
343#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
344 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- N(" << N <<
")-ary Polynomial Evaluation");
348 for (
int i=0; i<N; i++) {
349 SDV sdv(Teuchos::View,
const_cast<value_type*
>(na[i]->coeff()),
351 ps_op->transformPCE2QP(1.0, sdv, navals[N][i], 0.0,
false);
357#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
358 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- N(" << N <<
")-ary Function Evaluation");
363 for (ordinal_type qp=0; qp<nqp; qp++) {
364 for (
int i=0; i<N; i++)
365 val[i] = navals[N][i][qp];
366 fvals[qp] = func(
val);
372#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
373 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::PSExp -- N(" << N <<
")-ary Polynomial Integration");
377 SDV c_sdv(Teuchos::View, c.
coeff(), pc);
378 ps_op->transformQP2PCE(1.0, fvals, c_sdv, 0.0,
false);
383template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
393template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
403template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
410 if (use_quad_for_times)
416template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
423#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
424 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divideEqual(OPA)");
427 ordinal_type p = c.
size();
428 value_type* cc = c.
coeff();
429 const value_type* xc = x.coeff();
430 for (ordinal_type i=0; i<p; i++)
434 if (use_quad_for_division)
441template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
448 if (use_quad_for_times)
454template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
464template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
474template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
481#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
482 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divide(OPA,OPA)");
485 ordinal_type pc = a.
size();
489 const value_type* ca = a.
coeff();
490 const value_type* cb = b.
coeff();
491 value_type* cc = c.
coeff();
493 for (ordinal_type i=0; i<pc; i++)
497 if (use_quad_for_division)
504template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
511 if (use_quad_for_division)
517template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
527template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
536template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
545template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
554template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
563template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
572template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
582template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
592template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
602template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
611template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
620template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
629template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
638template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
647template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
656template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
665template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
674template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
683template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
693template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
703template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
713template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
722template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
731template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
740template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
741template <
typename ExprT1,
typename ExprT2>
746 ordinal_type pa = a.size();
747 ordinal_type pb = b.size();
749 if (pa > 1 && pb > 1) {
750 ordinal_type k_lim = pa;
751 ordinal_type j_lim = pb;
758 TEUCHOS_TEST_FOR_EXCEPTION(i_it == this->Cijk->i_end(), std::logic_error,
759 "Stokhos::PseudoSpectralOrthogPolyExpansion::compute_times_coeff()"
760 <<
": Index " << i <<
" is out of range [0,"
761 << this->Cijk->num_i() <<
")!");
763 value_type cc = value_type(0);
764 value_type aa, bb, cijk;
767 k_it != this->Cijk->k_end(i_it); ++k_it) {
774 aa = a.higher_order_coeff(k);
780 aa = b.higher_order_coeff(k);
783 j_it != this->Cijk->j_end(k_it); ++j_it) {
791 bb = b.higher_order_coeff(
j);
797 bb = a.higher_order_coeff(
j);
804 return cc / this->basis->norm_squared(i);
807 return a.val() * b.val();
809 return a.higher_order_coeff(i)*b.val();
812 return a.val()*b.higher_order_coeff(i);
816template <
typename ordinal_type,
typename value_type,
typename po
int_compare_type,
typename node_type>
817template <
typename ExprT1,
typename ExprT2>
824 TEUCHOS_TEST_FOR_EXCEPTION(i_it == this->Cijk->i_end(), std::logic_error,
825 "Stokhos::PseudoSpectralOrthogPolyExpansion::fast_ompute_times_coeff()"
826 <<
": Index " << i <<
" is out of range [0,"
827 << this->Cijk->num_i() <<
")!");
829 value_type cc = value_type(0);
830 value_type aa, bb, cijk;
833 k_it != this->Cijk->k_end(i_it); ++k_it) {
838 aa = a.higher_order_coeff(k);
840 j_it != this->Cijk->j_end(k_it); ++j_it) {
846 bb = b.higher_order_coeff(
j);
851 return cc / this->basis->norm_squared(i);
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
pointer coeff()
Return coefficient array.
ordinal_type size() const
Return size.
Abstract base class for multivariate orthogonal polynomials.
Base class for consolidating common expansion implementations.
void binary_op(const FuncT &func, 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)
Nonlinear binary function.
void asinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void asin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
value_type fast_compute_times_coeff(ordinal_type k, const ExprT1 &a, const ExprT2 &b) const
value_type compute_times_coeff(ordinal_type k, const ExprT1 &a, const ExprT2 &b) const
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void acosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void tanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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 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 unary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Nonlinear unary function.
void pow(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 nary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > **a)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void cbrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void acos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void log(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void atan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Teuchos::SerialDenseVector< ordinal_type, value_type > SDV
Short-hand for SerialDenseVector.
void sqrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void cosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void tan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void atan2(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 atanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
PseudoSpectralOrthogPolyExpansion(const Teuchos::RCP< const OrthogPolyBasis< ordinal_type, value_type > > &basis, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::RCP< const PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > > &ps_op, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Constructor.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void log10(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void sinh(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.
ikj_sparse_array::const_iterator i_iterator
Iterator for looping over i entries.
Bi-directional iterator for traversing a sparse array.