42 #ifndef STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
43 #define STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
57 template <
typename ordinal_t,
59 typename point_compare_type =
60 typename DefaultPointCompare<ordinal_t,value_t>::type>
81 bool use_pst_ =
false,
83 const point_compare_type& point_compare = point_compare_type()) :
86 dim(product_basis.dimension()),
96 const tb_type *tb =
dynamic_cast<const tb_type*
>(&product_basis);
108 if (bases[i]->order() < max_orders[i])
109 bases[i] = bases[i]->cloneWithOrder(max_orders[i]);
113 if (multiindex.dimension() == 0)
114 multiindex = max_orders;
126 bases[k]->getQuadPoints(2*multiindex[k], gp[k], gw[k], gv[k]);
127 n[k] = gp[k].
size()-1;
140 bases[k]->norm_squared(i);
150 index_iterator point_iterator = pointIndexSet.
begin();
151 index_iterator point_end = pointIndexSet.
end();
153 while (point_iterator != point_end) {
156 point[k] = gp[k][(*point_iterator)[k]];
157 w *= gw[k][(*point_iterator)[k]];
167 typename point_set_type::iterator di =
points.begin();
168 typename point_set_type::iterator di_end =
points.end();
169 while (di != di_end) {
170 di->second.second = idx;
179 typedef std::map<multiindex_type,ordinal_type,lexo_type> reorder_type;
180 reorder_type reorder_map;
182 reorder_map[product_basis.
term(i)] = i;
185 for (
typename reorder_type::iterator it = reorder_map.begin();
186 it != reorder_map.end(); ++it)
187 perm[idx++] = it->second;
198 for (point_iterator = pointIndexSet.
begin();
199 point_iterator != point_end;
202 point[k] = gp[k][(*point_iterator)[k]];
208 qp2pce(i,j) *= gw[k][l]*gv[k][l][m] / bases[k]->norm_squared(m);
209 pce2qp(j,i) *= gv[k][l][m];
253 it ==
points.end(), std::logic_error,
"Invalid term " <<
point);
254 return it->second.second;
273 bool trans =
false)
const {
293 bool trans =
false)
const {
374 bool reorder_result)
const {
392 M *= Ak[k].numRows();
393 N *= Ak[k].numCols();
409 tmp1[i+
j*n] = input(
j,idx);
415 tmp1[i+
j*n] = input(
j,i);
424 tmp1[i+
j*n] = input(idx,
j);
430 tmp1[i+
j*n] = input(i,
j);
450 y(i,
j+l*n) = x(
j,i+l*nk);
464 if (reorder_result) {
469 result(
j,idx) = beta*result(
j,idx) + alpha*tmp1[i+
j*m];
475 result(
j,i) = beta*result(
j,i) + alpha*tmp1[i+
j*m];
479 if (reorder_result) {
484 result(idx,
j) = beta*result(idx,
j) + alpha*tmp1[i+
j*m];
490 result(i,
j) = beta*result(i,
j) + alpha*tmp1[i+
j*m];
ScalarType * values() const
ordinal_type dim
Dimension.
point_set_type::iterator set_iterator
base_type::const_iterator const_iterator
ordinal_type coeff_size() const
Number of coefficients.
std::map< point_type, std::pair< value_type, ordinal_type >, point_compare_type > point_set_type
const point_type & point(ordinal_type n) const
Get point for given index.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > qp2pce
Matrix mapping points to coefficients.
bool use_pst
Use partial-summation-transformation.
Container storing a term in a generalized tensor product.
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > qp2pce_k
Matrix mapping points to coefficients for each dimension for PST.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
An operator interface for building pseudo-spectral approximations.
set_iterator set_begin()
Iterator to begining of point set.
base_type::const_set_iterator const_set_iterator
point_set_type points
Quadrature points.
bool reorder
Do we need to reorder coefficients for PST.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
iterator end()
Iterator to end of point set.
base_type::point_type point_type
const_iterator end() const
Iterator to end of point set.
set_iterator set_end()
Iterator to end of point set.
Teuchos::BLAS< ordinal_type, value_type > blas
BLAS wrappers.
PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > base_type
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const =0
Get orders of each coordinate polynomial given an index i.
point_map_type point_map
Map index to point term.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
TensorProductPseudoSpectralOperator(const ProductBasis< ordinal_type, value_type > &product_basis, bool use_pst_=false, multiindex_type multiindex=multiindex_type(), const point_compare_type &point_compare=point_compare_type())
Constructor.
base_type::iterator iterator
MultiIndex< ordinal_type > multiindex_type
base_type::point_map_type point_map_type
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
base_type::point_set_type point_set_type
void apply_direct(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const value_type &alpha, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result, const value_type &beta, bool trans) const
Apply transformation operator using direct method.
void resize(size_type new_size, const value_type &x=value_type())
const_iterator begin() const
Return iterator for first element in the set.
const_set_iterator set_begin() const
Iterator to begining of point set.
const_iterator begin() const
Iterator to begining of point set.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > pce2qp
Matrix mapping coefficients to points.
virtual void transformQP2PCE(const value_type &alpha, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result, const value_type &beta, bool trans=false) const
Transform values at quadrature points to PCE coefficients.
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
OrdinalType numCols() const
void apply_pst(const Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > &Ak, const value_type &alpha, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result, const value_type &beta, bool trans, bool reorder_input, bool reorder_result) const
Apply tranformation operator using PST method.
virtual void transformPCE2QP(const value_type &alpha, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result, const value_type &beta, bool trans=false) const
Transform PCE coefficients to quadrature values.
base_type::set_iterator set_iterator
ordinal_type coeff_sz
Number of coefficients.
int reshape(OrdinalType numRows, OrdinalType numCols)
point_map_type::iterator iterator
ordinal_type point_size() const
Number of points.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Iterator class for iterating over elements of the index set.
An operator for building pseudo-spectral coefficients using tensor-product quadrature.
virtual ~TensorProductPseudoSpectralOperator()
Destructor.
iterator begin()
Iterator to begining of point set.
#define TEUCHOS_ASSERT(assertion_test)
const_set_iterator set_end() const
Iterator to end of point set.
point_map_type::const_iterator const_iterator
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > pce2qp_k
Matrix mapping coefficients to points for each dimension for PST.
A tensor product index set.
const_iterator end() const
Return iterator for end of the index set.
ordinal_type index(const point_type &point) const
Get point index for given point.
Teuchos::Array< ordinal_type > perm
Permutation array when reordering for PST.
virtual ordinal_type size() const =0
Return total size of basis.
OrdinalType stride() const
virtual MultiIndex< ordinal_type > getMaxOrders() const =0
Return maximum order allowable for each coordinate basis.
OrdinalType numRows() const
point_set_type::const_iterator const_set_iterator