42 #ifndef STOKHOS_QUADRATURE_PSEUDO_SPECTRAL_OPERATOR_HPP 43 #define STOKHOS_QUADRATURE_PSEUDO_SPECTRAL_OPERATOR_HPP 47 #include "Teuchos_SerialDenseMatrix.hpp" 48 #include "Teuchos_BLAS.hpp" 56 template <
typename ordinal_t,
58 typename point_compare_type =
81 const point_compare_type& point_compare = point_compare_type()) :
85 const Teuchos::Array<value_type>& quad_weights = quad.
getQuadWeights();
86 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
88 const Teuchos::Array< Teuchos::Array<value_type> > & quad_vals =
99 thePoint[k] = quad_points[i][k];
101 points[thePoint] = std::make_pair(quad_weights[i],i);
109 typename point_set_type::iterator di =
points.begin();
110 typename point_set_type::iterator di_end =
points.end();
112 for (; di != di_end; ++di) {
115 qp2pce(i,jdx) = quad_weights[
j]*quad_vals[
j][i] /
117 pce2qp(jdx,i) = quad_vals[
j][i];
126 while (di != di_end) {
127 di->second.second = idx;
172 TEUCHOS_TEST_FOR_EXCEPTION(
173 it ==
points.end(), std::logic_error,
"Invalid term " <<
point);
174 return it->second.second;
190 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
191 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
193 bool trans =
false)
const {
196 ret = result.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, alpha,
199 ret = result.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha,
201 TEUCHOS_ASSERT(ret == 0);
214 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
215 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
217 bool trans =
false)
const {
219 TEUCHOS_ASSERT(input.numCols() <=
pce2qp.numCols());
220 TEUCHOS_ASSERT(result.numCols() ==
pce2qp.numRows());
221 TEUCHOS_ASSERT(result.numRows() == input.numRows());
222 blas.GEMM(Teuchos::NO_TRANS, Teuchos::TRANS, input.numRows(),
223 pce2qp.numRows(), input.numCols(), alpha, input.values(),
225 result.values(), result.stride());
228 TEUCHOS_ASSERT(input.numRows() <=
pce2qp.numCols());
229 TEUCHOS_ASSERT(result.numRows() ==
pce2qp.numRows());
230 TEUCHOS_ASSERT(result.numCols() == input.numCols());
231 blas.GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS,
pce2qp.numRows(),
232 input.numCols(), input.numRows(), alpha,
pce2qp.values(),
233 pce2qp.stride(), input.values(), input.stride(), beta,
234 result.values(), result.stride());
250 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
qp2pce;
253 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
pce2qp;
256 Teuchos::BLAS<ordinal_type,value_type>
blas;
set_iterator set_begin()
Iterator to begining of point set.
Teuchos::BLAS< ordinal_type, value_type > blas
BLAS wrappers.
point_set_type::iterator set_iterator
std::map< point_type, std::pair< value_type, ordinal_type >, point_compare_type > point_set_type
Teuchos::SerialDenseMatrix< ordinal_type, value_type > pce2qp
Matrix mapping coefficients to points.
ordinal_type coeff_sz
Number of coefficients.
Container storing a term in a generalized tensor product.
base_type::point_set_type point_set_type
Teuchos::Array< point_type > point_map_type
base_type::iterator iterator
An operator interface for building pseudo-spectral approximations.
ordinal_type point_size() const
Number of points.
virtual ordinal_type size() const =0
Get number of quadrature points.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
base_type::const_iterator const_iterator
iterator begin()
Iterator to begining of point set.
MultiIndex< ordinal_type > multiindex_type
const_iterator end() const
Iterator to end of point set.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > qp2pce
Matrix mapping points to coefficients.
const point_type & point(ordinal_type n) const
Get point for given index.
virtual ordinal_type dimension() const =0
Return dimension of basis.
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.
point_map_type point_map
Map index to point term.
virtual ~QuadraturePseudoSpectralOperator()
Destructor.
virtual const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
base_type::const_set_iterator const_set_iterator
Abstract base class for multivariate orthogonal polynomials.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const =0
Get quadrature points.
Abstract base class for quadrature methods.
const_iterator begin() const
Iterator to begining 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.
Top-level namespace for Stokhos classes and functions.
PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > base_type
iterator end()
Iterator to end of point set.
base_type::point_map_type point_map_type
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
base_type::point_type point_type
point_map_type::iterator iterator
const_set_iterator set_end() const
Iterator to end of point set.
base_type::set_iterator set_iterator
const_set_iterator set_begin() const
Iterator to begining of point set.
ordinal_type index(const point_type &point) const
Get point index for given point.
point_map_type::const_iterator const_iterator
ordinal_type coeff_size() const
Number of coefficients.
LexographicLess< TensorProductElement< ordinal_type, value_type >, FloatingPointLess< value_type > > type
An operator for building pseudo-spectral coefficients using an arbitrary quadrature rule...
point_set_type points
Quadrature points.
virtual ordinal_type size() const =0
Return total size of basis.
set_iterator set_end()
Iterator to end of point set.
point_set_type::const_iterator const_set_iterator
QuadraturePseudoSpectralOperator(const OrthogPolyBasis< ordinal_type, value_type > &basis, const Quadrature< ordinal_type, value_type > &quad, const point_compare_type &point_compare=point_compare_type())
Constructor.