Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SmolyakPseudoSpectralOperatorImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41template <typename ordinal_type, typename value_type,
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,
48 bool use_pst,
49 const point_compare_type& point_compare) :
50 use_smolyak(use_smolyak_apply),
51 points(point_compare) {
52
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;
55
56 // Generate sparse grid and tensor operators
57 coeff_sz = smolyak_basis.size();
58 ordinal_type dim = smolyak_basis.dimension();
59 ordinal_type num_terms = smolyak_basis.getNumSmolyakTerms();
60 for (ordinal_type i=0; i<num_terms; ++i) {
61
62 // Get tensor product basis for given term
63 Teuchos::RCP<const tensor_product_basis_type> tp_basis =
64 smolyak_basis.getTensorProductBasis(i);
65
66 // Get coefficient multi-index defining basis orders
67 multiindex_type coeff_index = tp_basis->getMaxOrders();
68
69 // Apply growth rule to cofficient multi-index
70 multiindex_type point_growth_index(dim);
71 for (ordinal_type j=0; j<dim; ++j) {
72 point_growth_index[j] =
73 smolyak_basis.getCoordinateBases()[j]->pointGrowth(coeff_index[j]);
74 }
75
76 // Build tensor product operator for given index
77 Teuchos::RCP<operator_type> op =
78 Teuchos::rcp(new operator_type(*tp_basis, use_pst,
79 point_growth_index));
80 if (use_smolyak)
81 operators.push_back(op);
82
83 // Get smolyak cofficient for given index
84 value_type c = smolyak_basis.getSmolyakCoefficient(i);
85 if (use_smolyak)
86 smolyak_coeffs.push_back(c);
87
88 // Include points in union over all sets
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) {
92 const point_type& point = op_set_iterator->first;
93 value_type w = op_set_iterator->second.first;
94 set_iterator si = points.find(point);
95 if (si == points.end())
96 points[point] = std::make_pair(c*w,ordinal_type(0));
97 else
98 si->second.first += c*w;
99 }
100
101 }
102
103 // Generate linear ordering of points
104 ordinal_type idx = 0;
105 point_map.resize(points.size());
106 for (set_iterator si = points.begin(); si != points.end(); ++si) {
107 si->second.second = idx;
108 point_map[idx] = si->first;
109 ++idx;
110 }
111
112 if (use_smolyak) {
113
114 // Build gather/scatter maps into global domain/range for each operator
115 gather_maps.resize(operators.size());
116 scatter_maps.resize(operators.size());
117 for (ordinal_type i=0; i<operators.size(); i++) {
118 Teuchos::RCP<operator_type> op = operators[i];
119
120 gather_maps[i].reserve(op->point_size());
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) {
124 gather_maps[i].push_back(points[*op_iterator].second);
125 }
126
127 Teuchos::RCP<const tensor_product_basis_type> tp_basis =
128 smolyak_basis.getTensorProductBasis(i);
129 ordinal_type op_coeff_sz = tp_basis->size();
130 scatter_maps[i].reserve(op_coeff_sz);
131 for (ordinal_type j=0; j<op_coeff_sz; ++j) {
132 scatter_maps[i].push_back(smolyak_basis.index(tp_basis->term(j)));
133 }
134 }
135 }
136
137 //else {
138
139 // Generate quadrature operator
140 ordinal_type nqp = points.size();
142 qp2pce.reshape(npc,nqp);
143 pce2qp.reshape(nqp,npc);
144 qp2pce.putScalar(1.0);
145 pce2qp.putScalar(1.0);
146 Teuchos::Array<value_type> vals(npc);
147 for (set_iterator si = points.begin(); si != points.end(); ++si) {
148 ordinal_type j = si->second.second;
149 value_type w = si->second.first;
150 point_type point = si->first;
151 smolyak_basis.evaluateBases(point, vals);
152 for (ordinal_type i=0; i<npc; ++i) {
153 qp2pce(i,j) = w*vals[i] / smolyak_basis.norm_squared(i);
154 pce2qp(j,i) = vals[i];
155 }
156 }
157 //}
158
159}
160
161template <typename ordinal_type, typename value_type,
162 typename point_compare_type>
163void
166 const value_type& alpha,
167 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
168 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
169 const value_type& beta,
170 bool trans) const {
171
172 if (use_smolyak)
173 transformQP2PCE_smolyak(alpha, input, result, beta, trans);
174 else
175 apply_direct(qp2pce, alpha, input, result, beta, trans);
176}
177
178template <typename ordinal_type, typename value_type,
179 typename point_compare_type>
180void
183 const value_type& alpha,
184 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
185 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
186 const value_type& beta,
187 bool trans) const {
188
189 // Currently we use the direct method for mapping PCE->QP because the
190 // current implementation doesn't work. Need to evaluate tensor bases
191 // on all quad points, not just the quad points associated with that
192 // basis.
193
194 // if (use_smolyak)
195 // transformPCE2QP_smolyak(alpha, input, result, beta, trans);
196 // else
197 apply_direct(pce2qp, alpha, input, result, beta, trans);
198}
199
200template <typename ordinal_type, typename value_type,
201 typename point_compare_type>
202void
205 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
206 const value_type& alpha,
207 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
208 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
209 const value_type& beta,
210 bool trans) const {
211 if (trans) {
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());
219 }
220 else {
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());
228 }
229}
230
231template <typename ordinal_type, typename value_type,
232 typename point_compare_type>
233void
236 const value_type& alpha,
237 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
238 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
239 const value_type& beta,
240 bool trans) const {
241 Teuchos::SerialDenseMatrix<ordinal_type,value_type> op_input, op_result;
242 result.scale(beta);
243 for (ordinal_type i=0; i<operators.size(); i++) {
244 Teuchos::RCP<operator_type> op = operators[i];
245 if (trans) {
246 op_input.reshape(input.numRows(), op->point_size());
247 op_result.reshape(result.numRows(), op->coeff_size());
248 }
249 else {
250 op_input.reshape(op->point_size(), input.numCols());
251 op_result.reshape(op->coeff_size(), result.numCols());
252 }
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);
256 }
257}
258
259template <typename ordinal_type, typename value_type,
260 typename point_compare_type>
261void
264 const value_type& alpha,
265 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
266 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
267 const value_type& beta,
268 bool trans) const {
269 Teuchos::SerialDenseMatrix<ordinal_type,value_type> op_input, op_result;
270 result.scale(beta);
271
272 for (ordinal_type i=0; i<operators.size(); i++) {
273 Teuchos::RCP<operator_type> op = operators[i];
274 if (trans) {
275 op_input.reshape(input.numRows(), op->coeff_size());
276 op_result.reshape(result.numRows(), op->point_size());
277 }
278 else {
279 op_input.reshape(op->coeff_size(), input.numCols());
280 op_result.reshape(op->point_size(), result.numCols());
281 }
282
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);
286 }
287}
288
289template <typename ordinal_type, typename value_type,
290 typename point_compare_type>
291void
293gather(
294 const Teuchos::Array<ordinal_type>& map,
295 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
296 bool trans,
297 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result) const {
298 if (trans) {
299 for (ordinal_type j=0; j<map.size(); j++)
300 for (ordinal_type i=0; i<input.numRows(); i++)
301 result(i,j) = input(i,map[j]);
302 }
303 else {
304 for (ordinal_type j=0; j<input.numCols(); j++)
305 for (ordinal_type i=0; i<map.size(); i++)
306 result(i,j) = input(map[i],j);
307 }
308}
309
310template <typename ordinal_type, typename value_type,
311 typename point_compare_type>
312void
314scatter(
315 const Teuchos::Array<ordinal_type>& map,
316 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
317 bool trans,
318 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result) const {
319 if (trans) {
320 for (ordinal_type j=0; j<map.size(); j++)
321 for (ordinal_type i=0; i<input.numRows(); i++)
322 result(i,map[j]) += input(i,j);
323 }
324 else {
325 for (ordinal_type j=0; j<input.numCols(); j++)
326 for (ordinal_type i=0; i<map.size(); i++)
327 result(map[i],j) += input(i,j);
328 }
329}
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.
Teuchos::Array< value_type > smolyak_coeffs
Smolyak coefficients.
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.
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
Container storing a term in a generalized tensor product.