42 typename point_compare_type>
43template <
typename coeff_compare_type>
46 const SmolyakBasis<ordinal_type,value_type,coeff_compare_type>& smolyak_basis,
47 bool use_smolyak_apply,
49 const point_compare_type& point_compare) :
50 use_smolyak(use_smolyak_apply),
51 points(point_compare) {
53 typedef SmolyakBasis<ordinal_type,value_type,coeff_compare_type> smolyak_basis_type;
54 typedef typename smolyak_basis_type::tensor_product_basis_type tensor_product_basis_type;
59 ordinal_type num_terms = smolyak_basis.getNumSmolyakTerms();
63 Teuchos::RCP<const tensor_product_basis_type> tp_basis =
64 smolyak_basis.getTensorProductBasis(i);
72 point_growth_index[
j] =
73 smolyak_basis.getCoordinateBases()[
j]->pointGrowth(coeff_index[
j]);
77 Teuchos::RCP<operator_type> op =
84 value_type c = smolyak_basis.getSmolyakCoefficient(i);
89 typename operator_type::set_iterator op_set_iterator = op->set_begin();
90 typename operator_type::set_iterator op_set_end = op->set_end();
91 for (; op_set_iterator != op_set_end; ++op_set_iterator) {
98 si->second.first += c*w;
107 si->second.second = idx;
118 Teuchos::RCP<operator_type> op =
operators[i];
121 typename operator_type::iterator op_iterator = op->begin();
122 typename operator_type::iterator op_end = op->end();
123 for (; op_iterator != op_end; ++op_iterator) {
127 Teuchos::RCP<const tensor_product_basis_type> tp_basis =
128 smolyak_basis.getTensorProductBasis(i);
132 scatter_maps[i].push_back(smolyak_basis.index(tp_basis->term(
j)));
146 Teuchos::Array<value_type> vals(npc);
151 smolyak_basis.evaluateBases(
point, vals);
153 qp2pce(i,
j) = w*vals[i] / smolyak_basis.norm_squared(i);
161template <
typename ordinal_type,
typename value_type,
162 typename point_compare_type>
167 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
168 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
173 transformQP2PCE_smolyak(alpha, input, result, beta, trans);
175 apply_direct(qp2pce, alpha, input, result, beta, trans);
178template <
typename ordinal_type,
typename value_type,
179 typename point_compare_type>
184 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
185 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
197 apply_direct(pce2qp, alpha, input, result, beta, trans);
200template <
typename ordinal_type,
typename value_type,
201 typename point_compare_type>
205 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
207 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
208 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
212 TEUCHOS_ASSERT(input.numCols() <= A.numCols());
213 TEUCHOS_ASSERT(result.numCols() == A.numRows());
214 TEUCHOS_ASSERT(result.numRows() == input.numRows());
215 blas.GEMM(Teuchos::NO_TRANS, Teuchos::TRANS, input.numRows(),
216 A.numRows(), input.numCols(), alpha, input.values(),
217 input.stride(), A.values(), A.stride(), beta,
218 result.values(), result.stride());
221 TEUCHOS_ASSERT(input.numRows() <= A.numCols());
222 TEUCHOS_ASSERT(result.numRows() == A.numRows());
223 TEUCHOS_ASSERT(result.numCols() == input.numCols());
224 blas.GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, A.numRows(),
225 input.numCols(), input.numRows(), alpha, A.values(),
226 A.stride(), input.values(), input.stride(), beta,
227 result.values(), result.stride());
231template <
typename ordinal_type,
typename value_type,
232 typename point_compare_type>
237 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
238 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
241 Teuchos::SerialDenseMatrix<ordinal_type,value_type> op_input, op_result;
244 Teuchos::RCP<operator_type> op = operators[i];
246 op_input.reshape(input.numRows(), op->point_size());
247 op_result.reshape(result.numRows(), op->coeff_size());
250 op_input.reshape(op->point_size(), input.numCols());
251 op_result.reshape(op->coeff_size(), result.numCols());
253 gather(gather_maps[i], input, trans, op_input);
254 op->transformQP2PCE(smolyak_coeffs[i], op_input, op_result, 0.0, trans);
255 scatter(scatter_maps[i], op_result, trans, result);
259template <
typename ordinal_type,
typename value_type,
260 typename point_compare_type>
265 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
266 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
269 Teuchos::SerialDenseMatrix<ordinal_type,value_type> op_input, op_result;
273 Teuchos::RCP<operator_type> op = operators[i];
275 op_input.reshape(input.numRows(), op->coeff_size());
276 op_result.reshape(result.numRows(), op->point_size());
279 op_input.reshape(op->coeff_size(), input.numCols());
280 op_result.reshape(op->point_size(), result.numCols());
283 gather(scatter_maps[i], input, trans, op_input);
284 op->transformPCE2QP(smolyak_coeffs[i], op_input, op_result, 0.0, trans);
285 scatter(gather_maps[i], op_result, trans, result);
289template <
typename ordinal_type,
typename value_type,
290 typename point_compare_type>
294 const Teuchos::Array<ordinal_type>& map,
295 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
297 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result)
const {
301 result(i,
j) = input(i,map[
j]);
306 result(i,
j) = input(map[i],
j);
310template <
typename ordinal_type,
typename value_type,
311 typename point_compare_type>
315 const Teuchos::Array<ordinal_type>& map,
316 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
318 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result)
const {
322 result(i,map[
j]) += input(i,
j);
327 result(map[i],
j) += input(i,
j);
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.
operator_set_type operators
Tensor pseudospectral operators.
void transformPCE2QP_smolyak(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
Transform PCE coefficients to values at quadrature points using Smolyak formula.
point_map_type point_map
Map index to sparse grid term.
void gather(const Teuchos::Array< ordinal_type > &map, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, bool trans, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result) const
Teuchos::SerialDenseMatrix< ordinal_type, value_type > pce2qp
Matrix mapping coefficients to points.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > qp2pce
Matrix mapping points to coefficients.
Teuchos::Array< Teuchos::Array< ordinal_type > > gather_maps
Gather maps for each operator for Smolyak apply.
point_set_type points
Smolyak sparse grid.
Teuchos::Array< value_type > smolyak_coeffs
Smolyak coefficients.
bool use_smolyak
Use Smolyak apply method.
SmolyakPseudoSpectralOperator(const SmolyakBasis< ordinal_type, value_type, coeff_compare_type > &smolyak_basis, bool use_smolyak_apply=true, bool use_pst=true, const point_compare_type &point_compare=point_compare_type())
Constructor.
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.
Teuchos::Array< Teuchos::Array< ordinal_type > > scatter_maps
Scatter maps for each operator for Smolyak apply.
base_type::set_iterator set_iterator
const point_type & point(ordinal_type n) const
Get point for given index.
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 scatter(const Teuchos::Array< ordinal_type > &map, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &input, bool trans, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &result) const
void transformQP2PCE_smolyak(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
Transform values at quadrature points to PCE coefficients using Smolyak formula.
TensorProductPseudoSpectralOperator< ordinal_type, value_type > operator_type
ordinal_type coeff_sz
Number of coefficients.
Container storing a term in a generalized tensor product.