41#include "Teuchos_TimeMonitor.hpp"
43template <
typename ordinal_type,
typename value_type>
47 const value_type& sparse_tol_,
48 bool use_old_cijk_alg_,
49 const Teuchos::RCP< Teuchos::Array<value_type> >& deriv_coeffs_) :
55 sparse_tol(sparse_tol_),
56 use_old_cijk_alg(use_old_cijk_alg_),
57 deriv_coeffs(deriv_coeffs_),
63 for (ordinal_type i=0; i<
d; i++) {
75 for (ordinal_type k=0; k<
sz; k++) {
76 nrm = value_type(1.0);
77 for (ordinal_type i=0; i<
d; i++)
83 name =
"Complete polynomial basis (";
84 for (ordinal_type i=0; i<
d-1; i++)
90 for (ordinal_type
j=0;
j<
d;
j++)
95 deriv_coeffs = Teuchos::rcp(
new Teuchos::Array<value_type>(
d));
96 for (ordinal_type
j=0;
j<
d;
j++)
101template <
typename ordinal_type,
typename value_type>
107template <
typename ordinal_type,
typename value_type>
115template <
typename ordinal_type,
typename value_type>
123template <
typename ordinal_type,
typename value_type>
131template <
typename ordinal_type,
typename value_type>
132const Teuchos::Array<value_type>&
139template <
typename ordinal_type,
typename value_type>
147template <
typename ordinal_type,
typename value_type>
148Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
152#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
153 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Fill Time");
155 if (use_old_cijk_alg)
156 return computeTripleProductTensorOld(sz);
158 return computeTripleProductTensorNew(sz);
161template <
typename ordinal_type,
typename value_type>
162Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
166#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
167 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Fill Time");
169 if (use_old_cijk_alg)
170 return computeTripleProductTensorOld(d+1);
172 return computeTripleProductTensorNew(d+1);
175template <
typename ordinal_type,
typename value_type>
176Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
181 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
185 Teuchos::Array< Teuchos::RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
186 for (ordinal_type i=0; i<d; i++)
187 Cijk_1d[i] = bases[i]->computeTripleProductTensor();
189 for (ordinal_type
j=0;
j<sz;
j++) {
190 for (ordinal_type i=0; i<sz; i++) {
191 for (ordinal_type k=0; k<order; k++) {
192 value_type c = value_type(1.0);
193 for (ordinal_type l=0; l<d; l++)
194 c *= (*Cijk_1d[l])(terms[i][l],terms[
j][l],terms[k][l]);
195 if (std::abs(c/norms[i]) > sparse_tol)
196 Cijk->add_term(i,
j,k,c);
201 Cijk->fillComplete();
206template <
typename ordinal_type,
typename value_type>
207Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
221 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
228 ordinal_type k_lim = 0;
229 for (ordinal_type i=0; i<d; i++)
230 k_lim = k_lim + term[i];
234 Teuchos::Array< Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
235 for (ordinal_type i=0; i<d; i++) {
236 if (k_lim <= basis_orders[i]+1)
237 Cijk_1d[i] = bases[i]->computeSparseTripleProductTensor(k_lim);
239 Cijk_1d[i] = bases[i]->computeSparseTripleProductTensor(basis_orders[i]+1);
249 Teuchos::Array<k_iterator> k_iterators(d, Cijk_1d[0]->k_begin());
250 Teuchos::Array<kj_iterator > j_iterators(d, Cijk_1d[0]->j_begin(k_iterators[0]));
251 Teuchos::Array<kji_iterator > i_iterators(d, Cijk_1d[0]->i_begin(j_iterators[0]));
253 ordinal_type sum_i = 0;
254 ordinal_type sum_j = 0;
255 ordinal_type sum_k = 0;
256 for (ordinal_type dim=0; dim<d; dim++) {
257 k_iterators[dim] = Cijk_1d[dim]->k_begin();
258 j_iterators[dim] = Cijk_1d[dim]->j_begin(k_iterators[dim]);
259 i_iterators[dim] = Cijk_1d[dim]->i_begin(j_iterators[dim]);
260 terms_i[dim] = Stokhos::index(i_iterators[dim]);
261 terms_j[dim] = Stokhos::index(j_iterators[dim]);
262 terms_k[dim] = Stokhos::index(k_iterators[dim]);
263 sum_i += terms_i[dim];
264 sum_j += terms_j[dim];
265 sum_k += terms_k[dim];
275 ordinal_type cnt = 0;
279 if (sum_i <= p && sum_j <= p && sum_k <= p) {
281 K = CPBUtils::compute_index(terms_k, terms, num_terms, p);
284 I = CPBUtils::compute_index(terms_i, terms, num_terms, p);
286 J = CPBUtils::compute_index(terms_j, terms, num_terms, p);
287 value_type c = value_type(1.0);
288 for (ordinal_type dim=0; dim<d; dim++)
289 c *= value(i_iterators[dim]);
290 if (std::abs(c/norms[I]) > sparse_tol)
291 Cijk->add_term(I,J,K,c);
296 ordinal_type cdim = 0;
301 while (inc && cdim < d) {
302 ordinal_type cur_dim = cdim;
305 if (i_iterators[cdim] != Cijk_1d[cdim]->i_end(j_iterators[cdim])) {
306 terms_i[cdim] = Stokhos::index(i_iterators[cdim]);
308 for (
int dim=0; dim<d; dim++)
309 sum_i += terms_i[dim];
311 if (i_iterators[cdim] == Cijk_1d[cdim]->i_end(j_iterators[cdim]) ||
315 if (j_iterators[cdim] != Cijk_1d[cdim]->j_end(k_iterators[cdim])) {
316 terms_j[cdim] = Stokhos::index(j_iterators[cdim]);
318 for (
int dim=0; dim<d; dim++)
319 sum_j += terms_j[dim];
321 if (j_iterators[cdim] == Cijk_1d[cdim]->j_end(k_iterators[cdim]) ||
325 if (k_iterators[cdim] != Cijk_1d[cdim]->k_end()) {
326 terms_k[cdim] = Stokhos::index(k_iterators[cdim]);
328 for (
int dim=0; dim<d; dim++)
329 sum_k += terms_k[dim];
331 if (k_iterators[cdim] == Cijk_1d[cdim]->k_end() || sum_k > p) {
332 k_iterators[cdim] = Cijk_1d[cdim]->k_begin();
334 terms_k[cur_dim] = Stokhos::index(k_iterators[cur_dim]);
336 for (
int dim=0; dim<d; dim++)
337 sum_k += terms_k[dim];
341 j_iterators[cur_dim] =
342 Cijk_1d[cur_dim]->j_begin(k_iterators[cur_dim]);
343 terms_j[cur_dim] = Stokhos::index(j_iterators[cur_dim]);
345 for (
int dim=0; dim<d; dim++)
346 sum_j += terms_j[dim];
350 i_iterators[cur_dim] = Cijk_1d[cur_dim]->i_begin(j_iterators[cur_dim]);
351 terms_i[cur_dim] = Stokhos::index(i_iterators[cur_dim]);
353 for (
int dim=0; dim<d; dim++)
354 sum_i += terms_i[dim];
359 if (sum_i > p || sum_j > p || sum_k > p)
369 Cijk->fillComplete();
374template <
typename ordinal_type,
typename value_type>
375Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
378 const Teuchos::RCP<
const Teuchos::SerialDenseMatrix<ordinal_type, value_type> >& Bij,
382 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Dijk =
384 for (ordinal_type i=0; i<sz; i++)
385 for (ordinal_type
j=0;
j<sz;
j++)
386 for (ordinal_type k=0; k<sz; k++)
387 (*Dijk)(i,
j,k) = value_type(0.0);
391 for (ordinal_type k=0; k<sz; k++) {
393 m_it!=Cijk->k_end(); ++m_it) {
394 m = Stokhos::index(m_it);
396 j_it != Cijk->j_end(m_it); ++j_it) {
397 j = Stokhos::index(j_it);
399 i_it != Cijk->i_end(j_it); ++i_it) {
400 i = Stokhos::index(i_it);
401 c = Stokhos::value(i_it);
402 (*Dijk)(i,
j,k) += (*Bij)(m,k)*c/norms[m];
411template <
typename ordinal_type,
typename value_type>
412Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
417 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > Bij =
418 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(sz,
422 Teuchos::Array< Teuchos::RCP<const Teuchos::SerialDenseMatrix<ordinal_type,value_type> > > Bij_1d(d);
423 for (ordinal_type i=0; i<d; i++)
424 Bij_1d[i] = bases[i]->computeDerivDoubleProductTensor();
426 for (ordinal_type i=0; i<sz; i++) {
427 for (ordinal_type k=0; k<sz; k++) {
428 value_type t = value_type(1.0);
429 value_type c = value_type(0.0);
430 for (ordinal_type
j=0;
j<d;
j++) {
431 bool is_zero =
false;
432 for (ordinal_type l=0; l<d; l++) {
433 if (l !=
j && terms[i][l] != terms[k][l])
436 t *= bases[l]->norm_squared(terms[k][l]);
439 c += t*(*deriv_coeffs)[
j]*(*Bij_1d[
j])(terms[k][
j],terms[i][
j]);
448template <
typename ordinal_type,
typename value_type>
455 value_type z = value_type(1.0);
456 for (ordinal_type
j=0;
j<d;
j++)
457 z = z * bases[
j]->evaluate(value_type(0.0), terms[i][
j]);
462template <
typename ordinal_type,
typename value_type>
465evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
466 Teuchos::Array<value_type>& basis_vals)
const
468 for (ordinal_type
j=0;
j<d;
j++)
469 bases[
j]->evaluateBases(point[
j], basis_eval_tmp[
j]);
472 for (ordinal_type i=0; i<sz; i++) {
473 value_type t = value_type(1.0);
474 for (ordinal_type
j=0;
j<d;
j++)
475 t *= basis_eval_tmp[
j][terms[i][
j]];
480template <
typename ordinal_type,
typename value_type>
483print(std::ostream& os)
const
485 os <<
"Complete basis of order " << p <<
", dimension " << d
486 <<
", and size " << sz <<
". Component bases:\n";
487 for (ordinal_type i=0; i<d; i++)
489 os <<
"Basis vector norms (squared):\n\t";
490 for (ordinal_type i=0; i<static_cast<ordinal_type>(norms.size()); i++)
491 os << norms[i] <<
" ";
495template <
typename ordinal_type,
typename value_type>
498term(ordinal_type i)
const
503template <
typename ordinal_type,
typename value_type>
508 return CPBUtils::compute_index(term, terms, num_terms, p);
511template <
typename ordinal_type,
typename value_type>
519template <
typename ordinal_type,
typename value_type>
520Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
527template <
typename ordinal_type,
typename value_type>
533 for (ordinal_type i=0; i<d; ++i)
534 max_orders[i] = basis_orders[i];
static ordinal_type compute_terms(ordinal_type p, ordinal_type d, ordinal_type &sz, Teuchos::Array< MultiIndex< ordinal_type > > &terms, Teuchos::Array< ordinal_type > &num_terms)
Compute the 2-D array of basis terms which maps a basis index into the orders for each basis dimensio...
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
virtual ~CompletePolynomialBasis()
Destructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
std::string name
Name of basis.
ordinal_type sz
Total size of basis.
virtual void print(std::ostream &os) const
Print basis to stream os.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > bases
Array of bases.
Teuchos::Array< MultiIndex< ordinal_type > > terms
2-D array of basis terms
CompletePolynomialBasis(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const value_type &sparse_tol=1.0e-12, bool use_old_cijk_alg=false, const Teuchos::RCP< Teuchos::Array< value_type > > &deriv_coeffs=Teuchos::null)
Constructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual ordinal_type size() const
Return total size of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorNew(ordinal_type order) const
Compute triple product tensor using new algorithm.
virtual Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > computeDerivDoubleProductTensor() const
Compute double product tensor where represents the derivative of in the direction .
ordinal_type order() const
Return order of basis.
Teuchos::RCP< Teuchos::Array< value_type > > deriv_coeffs
Coefficients for derivative.
virtual const std::string & getName() const
Return string name of basis.
Teuchos::Array< value_type > norms
Norms.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorOld(ordinal_type order) const
Compute triple product tensor using old algorithm.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeDerivTripleProductTensor(const Teuchos::RCP< const Teuchos::SerialDenseMatrix< ordinal_type, value_type > > &Bij, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk) const
Compute triple product tensor where represents the derivative of in the direction .
ordinal_type dimension() const
Return dimension of basis.
Teuchos::Array< Teuchos::Array< value_type > > basis_eval_tmp
Temporary array used in basis evaluation.
Teuchos::Array< ordinal_type > basis_orders
Array storing order of each basis.
ordinal_type p
Total order of basis.
ordinal_type d
Total dimension of basis.
Data structure storing a dense 3-tensor C(i,j,k).
A multidimensional index.
Abstract base class for 1-D orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
Bi-directional iterator for traversing a sparse array.