42#ifndef STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
43#define STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
48#include "Teuchos_SerialDenseMatrix.hpp"
49#include "Teuchos_BLAS.hpp"
57 template <
typename ordinal_t,
59 typename point_compare_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);
101 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > bases_ = product_basis.
getCoordinateBases();
105 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(bases_);
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;
117 Teuchos::Array< Teuchos::Array<value_type> > gp(
dim);
118 Teuchos::Array< Teuchos::Array<value_type> > gw(
dim);
119 Teuchos::Array< Teuchos::Array< Teuchos::Array<value_type> > > gv(
dim);
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);
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 it ==
points.end(), std::logic_error,
"Invalid term " <<
point);
254 return it->second.second;
270 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
271 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
273 bool trans =
false)
const {
290 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
291 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
293 bool trans =
false)
const {
299 can_use_pst = can_use_pst && (input.numCols() ==
pce2qp.numRows());
301 can_use_pst = can_use_pst && (input.numRows() ==
pce2qp.numRows());
313 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
315 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
316 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
320 TEUCHOS_ASSERT(input.numCols() <= A.numCols());
321 TEUCHOS_ASSERT(result.numCols() == A.numRows());
322 TEUCHOS_ASSERT(result.numRows() == input.numRows());
323 blas.GEMM(Teuchos::NO_TRANS, Teuchos::TRANS, input.numRows(),
324 A.numRows(), input.numCols(), alpha, input.values(),
325 input.stride(), A.values(), A.stride(), beta,
326 result.values(), result.stride());
329 TEUCHOS_ASSERT(input.numRows() <= A.numCols());
330 TEUCHOS_ASSERT(result.numRows() == A.numRows());
331 TEUCHOS_ASSERT(result.numCols() == input.numCols());
332 blas.GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, A.numRows(),
333 input.numCols(), input.numRows(), alpha, A.values(),
334 A.stride(), input.values(), input.stride(), beta,
335 result.values(), result.stride());
367 const Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >& Ak,
369 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
370 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
374 bool reorder_result)
const {
378 TEUCHOS_ASSERT(input.numRows() == result.numRows());
381 m = result.numCols();
384 TEUCHOS_ASSERT(input.numCols() == result.numCols());
387 m = result.numRows();
392 M *= Ak[k].numRows();
393 N *= Ak[k].numCols();
395 TEUCHOS_ASSERT(n == N);
396 TEUCHOS_ASSERT(m == M);
400 Teuchos::Array<value_type> tmp1(sz*p), tmp2(sz*p);
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);
441 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x(
442 Teuchos::View, &tmp1[0], n, n, nk*p);
445 Teuchos::SerialDenseMatrix<ordinal_type,value_type> y(
446 Teuchos::View, &tmp2[0], nk, nk, n*p);
450 y(i,
j+l*n) = x(
j,i+l*nk);
453 Teuchos::SerialDenseMatrix<ordinal_type,value_type> z(
454 Teuchos::View, &tmp1[0], mk, mk, n*p);
456 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Ak[k], y, 0.0);
457 TEUCHOS_ASSERT(ret == 0);
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];
516 Teuchos::Array<ordinal_type>
perm;
519 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
qp2pce;
522 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
pce2qp;
525 Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >
qp2pce_k;
528 Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >
pce2qp_k;
531 Teuchos::BLAS<ordinal_type,value_type>
blas;
A comparison functor implementing a strict weak ordering based lexographic ordering.
A multidimensional index.
ordinal_type size() const
Size.
virtual ordinal_type size() const =0
Return total size of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
virtual MultiIndex< ordinal_type > getMaxOrders() const =0
Return maximum order allowable for each coordinate basis.
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const =0
Get orders of each coordinate polynomial given an index i.
An operator interface for building pseudo-spectral approximations.
point_map_type::const_iterator const_iterator
point_set_type::iterator set_iterator
point_set_type::const_iterator const_set_iterator
Teuchos::Array< point_type > point_map_type
std::map< point_type, std::pair< value_type, ordinal_type >, point_compare_type > point_set_type
point_map_type::iterator iterator
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials.
Container storing a term in a generalized tensor product.
Iterator class for iterating over elements of the index set.
A tensor product index set.
const_iterator begin() const
Return iterator for first element in the set.
const_iterator end() const
Return iterator for end of the index set.
An operator for building pseudo-spectral coefficients using tensor-product quadrature.
iterator begin()
Iterator to begining of point set.
ordinal_type coeff_size() const
Number of coefficients.
const_iterator end() const
Iterator to end of point set.
set_iterator set_begin()
Iterator to begining of point set.
base_type::iterator iterator
base_type::point_map_type point_map_type
base_type::const_set_iterator const_set_iterator
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.
point_set_type points
Quadrature points.
virtual ~TensorProductPseudoSpectralOperator()
Destructor.
base_type::set_iterator set_iterator
const_set_iterator set_end() const
Iterator to end of point set.
point_map_type point_map
Map index to point term.
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.
const_iterator begin() const
Iterator to begining of point set.
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.
base_type::point_type point_type
set_iterator set_end()
Iterator to end of point set.
ordinal_type dim
Dimension.
base_type::const_iterator const_iterator
const point_type & point(ordinal_type n) const
Get point for given index.
PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > base_type
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.
bool reorder
Do we need to reorder coefficients for PST.
ordinal_type coeff_sz
Number of coefficients.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > qp2pce
Matrix mapping points to coefficients.
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > pce2qp_k
Matrix mapping coefficients to points for each dimension for PST.
Teuchos::BLAS< ordinal_type, value_type > blas
BLAS wrappers.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > pce2qp
Matrix mapping coefficients to points.
MultiIndex< ordinal_type > multiindex_type
ordinal_type index(const point_type &point) const
Get point index for given point.
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > qp2pce_k
Matrix mapping points to coefficients for each dimension for PST.
const_set_iterator set_begin() const
Iterator to begining of point set.
bool use_pst
Use partial-summation-transformation.
iterator end()
Iterator to end of point set.
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.
ordinal_type point_size() const
Number of points.
Teuchos::Array< ordinal_type > perm
Permutation array when reordering for PST.
base_type::point_set_type point_set_type
Top-level namespace for Stokhos classes and functions.
LexographicLess< TensorProductElement< ordinal_type, value_type >, FloatingPointLess< value_type > > type