44#include "Teuchos_Assert.hpp"
46#include "Teuchos_TimeMonitor.hpp"
47#include "Teuchos_Tuple.hpp"
49template <
typename ordinal_type,
typename value_type,
typename node_type>
54 const Teuchos::RCP<
const Quadrature<ordinal_type, value_type> >& quad_,
55 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
58 sz(this->basis->size()),
60 quad_points(quad->getQuadPoints()),
61 quad_weights(quad->getQuadWeights()),
62 quad_values(quad->getBasisAtQuadPoints()),
63 norms(this->basis->norm_squared()),
64 nqp(quad_points.size()),
71 for (ordinal_type qp=0; qp<nqp; qp++)
72 for (ordinal_type i=0; i<sz; i++) {
73 qv[qp*sz+i] = quad_values[qp][i];
74 sqv[qp*sz+i] = quad_values[qp][i]/norms[i];
77 Teuchos::RCP<Teuchos::ParameterList> params = params_;
78 if (params == Teuchos::null)
79 params = Teuchos::rcp(
new Teuchos::ParameterList);
80 use_quad_for_times = params->get(
"Use Quadrature for Times",
false);
81 use_quad_for_division = params->get(
"Use Quadrature for Division",
true);
84template <
typename ordinal_type,
typename value_type,
typename node_type>
85template <
typename FuncT>
92 ordinal_type pa = a.
size();
107#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
108 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- Unary Polynomial Evaluation");
112 blas.GEMV(Teuchos::TRANS, pa, nqp, 1.0, &qv[0], sz, a.
coeff(), 1, 0.0,
118#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
119 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- Unary Function Evaluation");
123 for (ordinal_type qp=0; qp<nqp; qp++) {
124 if (quad_weights[qp] != value_type(0))
125 fvals[qp] = func(avals[qp])*quad_weights[qp];
127 fvals[qp] = value_type(0);
133#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
134 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- Unary Polynomial Integration");
138 blas.GEMV(Teuchos::NO_TRANS, pc, nqp, 1.0, &sqv[0], sz, &fvals[0], 1, 0.0,
144template <
typename ordinal_type,
typename value_type,
typename node_type>
145template <
typename FuncT>
153 ordinal_type pa = a.
size();
154 ordinal_type pb = b.
size();
156 if (pa == 1 && pb == 1)
164 c[0] = func(a[0], b[0]);
169#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
170 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PP Binary Polynomial Evaluation");
174 blas.GEMV(Teuchos::TRANS, pa, nqp, 1.0, &qv[0], sz, a.
coeff(), 1, 0.0,
176 blas.GEMV(Teuchos::TRANS, pb, nqp, 1.0, &qv[0], sz, b.
coeff(), 1, 0.0,
182#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
183 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PP Binary Function Evaluation");
187 for (ordinal_type qp=0; qp<nqp; qp++)
188 if (quad_weights[qp] != value_type(0))
189 fvals[qp] = func(avals[qp], bvals[qp])*quad_weights[qp];
191 fvals[qp] = value_type(0);
196#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
197 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PP Binary Polynomial Integration");
201 blas.GEMV(Teuchos::NO_TRANS, pc, nqp, 1.0, &sqv[0], sz, &fvals[0], 1, 0.0,
207template <
typename ordinal_type,
typename value_type,
typename node_type>
208template <
typename FuncT>
216 ordinal_type pb = b.
size();
226 c[0] = func(a, b[0]);
231#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
232 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- CP Binary Polynomial Evaluation");
236 blas.GEMV(Teuchos::TRANS, pb, nqp, 1.0, &qv[0], sz, b.
coeff(), 1, 0.0,
242#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
243 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- CP Binary Function Evaluation");
247 for (ordinal_type qp=0; qp<nqp; qp++)
248 if (quad_weights[qp] != value_type(0))
249 fvals[qp] = func(a, bvals[qp])*quad_weights[qp];
251 fvals[qp] = value_type(0);
256#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
257 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- CP Binary Polynomial Integration");
261 blas.GEMV(Teuchos::NO_TRANS, pc, nqp, 1.0, &sqv[0], sz, &fvals[0], 1, 0.0,
267template <
typename ordinal_type,
typename value_type,
typename node_type>
268template <
typename FuncT>
276 ordinal_type pa = a.
size();
286 c[0] = func(a[0], b);
291#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
292 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PC Binary Polynomial Evaluation");
296 blas.GEMV(Teuchos::TRANS, pa, nqp, 1.0, &qv[0], sz, a.
coeff(), 1, 0.0,
302#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
303 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PC Binary Function Evaluation");
307 for (ordinal_type qp=0; qp<nqp; qp++)
308 if (quad_weights[qp] != value_type(0))
309 fvals[qp] = func(avals[qp], b)*quad_weights[qp];
311 fvals[qp] = value_type(0);
316#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
317 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- PC Binary Polynomial Integration");
321 blas.GEMV(Teuchos::NO_TRANS, pc, nqp, 1.0, &sqv[0], sz, &fvals[0], 1, 0.0,
327template <
typename ordinal_type,
typename value_type,
typename node_type>
328template <
typename FuncT>
335 const int N = FuncT::N;
336 bool is_constant =
true;
337 for (
int i=0; i<N; i++) {
338 if (na[i]->size() > 1) {
353 for (
int i=0; i<N; i++)
354 val[i] = (*na[i])[0];
359 if (N >= navals.size())
361 if (navals[N].size() != N) {
363 for (
int i=0; i<N; i++)
364 navals[N][i].resize(nqp);
368#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
369 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- N(" << N <<
")-ary Polynomial Evaluation");
373 for (
int i=0; i<N; i++) {
375 ordinal_type pa = na[i]->
size();
376 blas.GEMV(Teuchos::TRANS, pa, nqp, 1.0, &qv[0], sz, na[i]->coeff(), 1, 0.0,
377 &navals[N][i][0], 1);
383#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
384 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- N(" << N <<
")-ary Function Evaluation");
389 for (ordinal_type qp=0; qp<nqp; qp++) {
390 if (quad_weights[qp] != value_type(0)) {
391 for (
int i=0; i<N; i++)
392 val[i] = navals[N][i][qp];
393 fvals[qp] = func(
val)*quad_weights[qp];
396 fvals[qp] = value_type(0);
402#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
403 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::QuadExp -- N(" << N <<
")-ary Polynomial Integration");
407 blas.GEMV(Teuchos::NO_TRANS, pc, nqp, 1.0, &sqv[0], sz, &fvals[0], 1, 0.0,
413template <
typename ordinal_type,
typename value_type,
typename node_type>
423template <
typename ordinal_type,
typename value_type,
typename node_type>
433template <
typename ordinal_type,
typename value_type,
typename node_type>
440 if (use_quad_for_times)
446template <
typename ordinal_type,
typename value_type,
typename node_type>
453#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
454 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divideEqual(OPA)");
457 ordinal_type p = c.
size();
458 value_type* cc = c.
coeff();
459 const value_type* xc = x.coeff();
460 for (ordinal_type i=0; i<p; i++)
464 if (use_quad_for_division)
471template <
typename ordinal_type,
typename value_type,
typename node_type>
478 if (use_quad_for_times)
484template <
typename ordinal_type,
typename value_type,
typename node_type>
494template <
typename ordinal_type,
typename value_type,
typename node_type>
504template <
typename ordinal_type,
typename value_type,
typename node_type>
511#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
512 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::OrthogPolyExpansionBase::divide(OPA,OPA)");
515 ordinal_type pc = a.
size();
519 const value_type* ca = a.
coeff();
520 const value_type* cb = b.
coeff();
521 value_type* cc = c.
coeff();
523 for (ordinal_type i=0; i<pc; i++)
527 if (use_quad_for_division)
534template <
typename ordinal_type,
typename value_type,
typename node_type>
541 if (use_quad_for_division)
547template <
typename ordinal_type,
typename value_type,
typename node_type>
557template <
typename ordinal_type,
typename value_type,
typename node_type>
566template <
typename ordinal_type,
typename value_type,
typename node_type>
575template <
typename ordinal_type,
typename value_type,
typename node_type>
584template <
typename ordinal_type,
typename value_type,
typename node_type>
593template <
typename ordinal_type,
typename value_type,
typename node_type>
602template <
typename ordinal_type,
typename value_type,
typename node_type>
612template <
typename ordinal_type,
typename value_type,
typename node_type>
622template <
typename ordinal_type,
typename value_type,
typename node_type>
632template <
typename ordinal_type,
typename value_type,
typename node_type>
641template <
typename ordinal_type,
typename value_type,
typename node_type>
650template <
typename ordinal_type,
typename value_type,
typename node_type>
659template <
typename ordinal_type,
typename value_type,
typename node_type>
668template <
typename ordinal_type,
typename value_type,
typename node_type>
677template <
typename ordinal_type,
typename value_type,
typename node_type>
686template <
typename ordinal_type,
typename value_type,
typename node_type>
695template <
typename ordinal_type,
typename value_type,
typename node_type>
704template <
typename ordinal_type,
typename value_type,
typename node_type>
713template <
typename ordinal_type,
typename value_type,
typename node_type>
723template <
typename ordinal_type,
typename value_type,
typename node_type>
733template <
typename ordinal_type,
typename value_type,
typename node_type>
743template <
typename ordinal_type,
typename value_type,
typename node_type>
752template <
typename ordinal_type,
typename value_type,
typename node_type>
761template <
typename ordinal_type,
typename value_type,
typename node_type>
770template <
typename ordinal_type,
typename value_type,
typename node_type>
771template <
typename ExprT1,
typename ExprT2>
776 ordinal_type pa = a.size();
777 ordinal_type pb = b.size();
779 if (pa > 1 && pb > 1) {
780 ordinal_type k_lim = pa;
781 ordinal_type j_lim = pb;
788 TEUCHOS_TEST_FOR_EXCEPTION(i_it == this->Cijk->i_end(), std::logic_error,
789 "Stokhos::QuadOrthogPolyExpansion::compute_times_coeff()"
790 <<
": Index " << i <<
" is out of range [0,"
791 << this->Cijk->num_i() <<
")!");
793 value_type cc = value_type(0);
794 value_type aa, bb, cijk;
797 k_it != this->Cijk->k_end(i_it); ++k_it) {
804 aa = a.higher_order_coeff(k);
810 aa = b.higher_order_coeff(k);
813 j_it != this->Cijk->j_end(k_it); ++j_it) {
821 bb = b.higher_order_coeff(
j);
827 bb = a.higher_order_coeff(
j);
834 return cc / norms[i];
837 return a.val() * b.val();
839 return a.higher_order_coeff(i)*b.val();
842 return a.val()*b.higher_order_coeff(i);
846template <
typename ordinal_type,
typename value_type,
typename node_type>
847template <
typename ExprT1,
typename ExprT2>
854 TEUCHOS_TEST_FOR_EXCEPTION(i_it == this->Cijk->i_end(), std::logic_error,
855 "Stokhos::QuadOrthogPolyExpansion::fast_ompute_times_coeff()"
856 <<
": Index " << i <<
" is out of range [0,"
857 << this->Cijk->num_i() <<
")!");
859 value_type cc = value_type(0);
860 value_type aa, bb, cijk;
863 k_it != this->Cijk->k_end(i_it); ++k_it) {
868 aa = a.higher_order_coeff(k);
870 j_it != this->Cijk->j_end(k_it); ++j_it) {
876 bb = b.higher_order_coeff(
j);
881 return cc / norms[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 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 tan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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 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 tanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void sqrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void acosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
QuadOrthogPolyExpansion(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 Quadrature< ordinal_type, value_type > > &quad, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Constructor.
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 nary_op(const FuncT &func, 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 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 sinh(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 asinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void cbrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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 atan(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)
value_type fast_compute_times_coeff(ordinal_type k, const ExprT1 &a, const ExprT2 &b) const
void acos(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)
value_type compute_times_coeff(ordinal_type k, const ExprT1 &a, const ExprT2 &b) const
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 asin(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 sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void atanh(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.