Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_TensorProductPseudoSpectralOperator.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
41
42#ifndef STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
43#define STOKHOS_TENSOR_PRODUCT_PSEUDO_SPECTRAL_OPERATOR_HPP
44
48#include "Teuchos_SerialDenseMatrix.hpp"
49#include "Teuchos_BLAS.hpp"
50
51namespace Stokhos {
52
57 template <typename ordinal_t,
58 typename value_t,
59 typename point_compare_type =
62 public PseudoSpectralOperator<ordinal_t,value_t,point_compare_type> {
63 public:
64
65 typedef ordinal_t ordinal_type;
66 typedef value_t value_type;
71 typedef typename base_type::iterator iterator;
75
77
80 const ProductBasis<ordinal_type,value_type>& product_basis,
81 bool use_pst_ = false,
82 multiindex_type multiindex = multiindex_type(),
83 const point_compare_type& point_compare = point_compare_type()) :
84 use_pst(use_pst_),
85 coeff_sz(product_basis.size()),
86 dim(product_basis.dimension()),
87 points(point_compare),
88 reorder(false) {
89
90 // Check if we need to reorder for PST
91 // This isn't a perfect test since it doesn't catch not tensor-product
92 // bases
94 if (use_pst) {
96 const tb_type *tb = dynamic_cast<const tb_type*>(&product_basis);
97 if (tb == NULL)
98 reorder = true;
99 }
100
101 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > bases_ = product_basis.getCoordinateBases();
102
103 // Make sure order of each 1-D basis is large enough for
104 // supplied set of coefficients
105 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(bases_);
106 multiindex_type max_orders = product_basis.getMaxOrders();
107 for (ordinal_type i=0; i<dim; ++i) {
108 if (bases[i]->order() < max_orders[i])
109 bases[i] = bases[i]->cloneWithOrder(max_orders[i]);
110 }
111
112 // Check for default quadrature order
113 if (multiindex.dimension() == 0)
114 multiindex = max_orders;
115
116 // Compute quad points, weights, values
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);
121 if (use_pst) {
122 qp2pce_k.resize(dim);
123 pce2qp_k.resize(dim);
124 }
125 for (ordinal_type k=0; k<dim; k++) {
126 bases[k]->getQuadPoints(2*multiindex[k], gp[k], gw[k], gv[k]);
127 n[k] = gp[k].size()-1;
128
129 // Generate quadrature operators for PST
130 if (use_pst) {
131 ordinal_type npc = bases[k]->size();
132 ordinal_type nqp = gp[k].size();
133 qp2pce_k[k].reshape(npc,nqp);
134 pce2qp_k[k].reshape(nqp,npc);
135 qp2pce_k[k].putScalar(1.0);
136 pce2qp_k[k].putScalar(1.0);
137 for (ordinal_type j=0; j<nqp; j++) {
138 for (ordinal_type i=0; i<npc; i++) {
139 qp2pce_k[k](i,j) *= gw[k][j]*gv[k][j][i] /
140 bases[k]->norm_squared(i);
141 pce2qp_k[k](j,i) *= gv[k][j][i];
142 }
143 }
144 }
145 }
146
147 // Generate points and weights
149 typedef typename TensorProductIndexSet<ordinal_type>::iterator index_iterator;
150 index_iterator point_iterator = pointIndexSet.begin();
151 index_iterator point_end = pointIndexSet.end();
153 while (point_iterator != point_end) {
154 value_type w = 1.0;
155 for (ordinal_type k=0; k<dim; k++) {
156 point[k] = gp[k][(*point_iterator)[k]];
157 w *= gw[k][(*point_iterator)[k]];
158 }
159 points[point] = std::make_pair(w,ordinal_type(0));
160 ++point_iterator;
161 }
162
163 // Generate linear ordering of points
164 ordinal_type nqp = points.size();
165 point_map.resize(nqp);
166 ordinal_type idx=0;
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;
171 point_map[idx] = di->first;
172 ++idx;
173 ++di;
174 }
175
176 // Build permutation array for reordering if coefficients aren't
177 // lexographically ordered
178 if (use_pst && reorder) {
179 typedef std::map<multiindex_type,ordinal_type,lexo_type> reorder_type;
180 reorder_type reorder_map;
181 for (ordinal_type i=0; i<coeff_sz; ++i)
182 reorder_map[product_basis.term(i)] = i;
183 perm.resize(coeff_sz);
184 ordinal_type idx = 0;
185 for (typename reorder_type::iterator it = reorder_map.begin();
186 it != reorder_map.end(); ++it)
187 perm[idx++] = it->second;
188 TEUCHOS_ASSERT(idx == coeff_sz);
189 }
190
191 // Generate quadrature operator for direct apply
192 // Always do this because we may not use PST in some situations
193 ordinal_type npc = product_basis.size();
194 qp2pce.reshape(npc,nqp);
195 pce2qp.reshape(nqp,npc);
196 qp2pce.putScalar(1.0);
197 pce2qp.putScalar(1.0);
198 for (point_iterator = pointIndexSet.begin();
199 point_iterator != point_end;
200 ++point_iterator) {
201 for (ordinal_type k=0; k<dim; k++)
202 point[k] = gp[k][(*point_iterator)[k]];
203 ordinal_type j = points[point].second;
204 for (ordinal_type i=0; i<npc; i++) {
205 for (ordinal_type k=0; k<dim; k++) {
206 ordinal_type l = (*point_iterator)[k];
207 ordinal_type m = product_basis.term(i)[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];
210 }
211 }
212 }
213
214 }
215
218
220 ordinal_type point_size() const { return points.size(); }
221
223 ordinal_type coeff_size() const { return coeff_sz; }
224
226 iterator begin() { return point_map.begin(); }
227
229 iterator end() { return point_map.end(); }
230
232 const_iterator begin() const { return point_map.begin(); }
233
235 const_iterator end() const { return point_map.end(); }
236
238 set_iterator set_begin() { return points.begin(); }
239
241 set_iterator set_end() { return points.end(); }
242
244 const_set_iterator set_begin() const { return points.begin(); }
245
247 const_set_iterator set_end() const { return points.end(); }
248
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 it == points.end(), std::logic_error, "Invalid term " << point);
254 return it->second.second;
255 }
256
258 const point_type& point(ordinal_type n) const { return point_map[n]; }
259
261
268 virtual void transformQP2PCE(
269 const value_type& alpha,
270 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
271 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
272 const value_type& beta,
273 bool trans = false) const {
274 if (use_pst)
275 apply_pst(qp2pce_k, alpha, input, result, beta, trans, false, reorder);
276 else
277 apply_direct(qp2pce, alpha, input, result, beta, trans);
278 }
279
281
288 virtual void transformPCE2QP(
289 const value_type& alpha,
290 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
291 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
292 const value_type& beta,
293 bool trans = false) const {
294
295 // Only use PST if it was requested and number of rows/columns match
296 // (since we may allow fewer rows in input than columns of pce2qp)
297 bool can_use_pst = use_pst;
298 if (trans)
299 can_use_pst = can_use_pst && (input.numCols() == pce2qp.numRows());
300 else
301 can_use_pst = can_use_pst && (input.numRows() == pce2qp.numRows());
302
303 if (can_use_pst)
304 apply_pst(pce2qp_k, alpha, input, result, beta, trans, reorder, false);
305 else
306 apply_direct(pce2qp, alpha, input, result, beta, trans);
307 }
308
309 protected:
310
313 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
314 const value_type& alpha,
315 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
316 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
317 const value_type& beta,
318 bool trans) const {
319 if (trans) {
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());
327 }
328 else {
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());
336 }
337 }
338
340
367 const Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >& Ak,
368 const value_type& alpha,
369 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
370 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
371 const value_type& beta,
372 bool trans,
373 bool reorder_input,
374 bool reorder_result) const {
375
376 ordinal_type n, m, p;
377 if (trans) {
378 TEUCHOS_ASSERT(input.numRows() == result.numRows());
379 n = input.numCols();
380 p = input.numRows();
381 m = result.numCols();
382 }
383 else {
384 TEUCHOS_ASSERT(input.numCols() == result.numCols());
385 n = input.numRows();
386 p = input.numCols();
387 m = result.numRows();
388 }
389 ordinal_type M = 1;
390 ordinal_type N = 1;
391 for (ordinal_type k=0; k<dim; ++k) {
392 M *= Ak[k].numRows();
393 N *= Ak[k].numCols();
394 }
395 TEUCHOS_ASSERT(n == N);
396 TEUCHOS_ASSERT(m == M);
397 ordinal_type sz = std::max(n,m);
398
399 // Temporary buffers used in algorithm
400 Teuchos::Array<value_type> tmp1(sz*p), tmp2(sz*p);
401
402 // Copy input to tmp1 (transpose if necessary)
403 if (trans) {
404 if (reorder_input) {
405 ordinal_type idx;
406 for (ordinal_type j=0; j<p; j++)
407 for (ordinal_type i=0; i<n; i++) {
408 idx = perm[i];
409 tmp1[i+j*n] = input(j,idx);
410 }
411 }
412 else {
413 for (ordinal_type j=0; j<p; j++)
414 for (ordinal_type i=0; i<n; i++)
415 tmp1[i+j*n] = input(j,i);
416 }
417 }
418 else {
419 if (reorder_input) {
420 ordinal_type idx;
421 for (ordinal_type j=0; j<p; j++)
422 for (ordinal_type i=0; i<n; i++) {
423 idx = perm[i];
424 tmp1[i+j*n] = input(idx,j);
425 }
426 }
427 else {
428 for (ordinal_type j=0; j<p; j++)
429 for (ordinal_type i=0; i<n; i++)
430 tmp1[i+j*n] = input(i,j);
431 }
432 }
433
434 // Loop over each term in Kronecker product
435 for (ordinal_type k=0; k<dim; k++) {
436 ordinal_type mk = Ak[k].numRows();
437 ordinal_type nk = Ak[k].numCols();
438 n = n / nk;
439
440 // x = reshape(x, n, n_k)
441 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x(
442 Teuchos::View, &tmp1[0], n, n, nk*p);
443
444 // y = x'
445 Teuchos::SerialDenseMatrix<ordinal_type,value_type> y(
446 Teuchos::View, &tmp2[0], nk, nk, n*p);
447 for (ordinal_type l=0; l<p; l++)
448 for (ordinal_type j=0; j<n; j++)
449 for (ordinal_type i=0; i<nk; i++)
450 y(i,j+l*n) = x(j,i+l*nk);
451
452 // x = A_k*x
453 Teuchos::SerialDenseMatrix<ordinal_type,value_type> z(
454 Teuchos::View, &tmp1[0], mk, mk, n*p);
455 ordinal_type ret =
456 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Ak[k], y, 0.0);
457 TEUCHOS_ASSERT(ret == 0);
458
459 n = n * mk;
460 }
461
462 // Sum answer into result (transposing and/or reordering if necessary)
463 if (trans) {
464 if (reorder_result) {
465 ordinal_type idx;
466 for (ordinal_type j=0; j<p; j++)
467 for (ordinal_type i=0; i<m; i++) {
468 idx = perm[i];
469 result(j,idx) = beta*result(j,idx) + alpha*tmp1[i+j*m];
470 }
471 }
472 else {
473 for (ordinal_type j=0; j<p; j++)
474 for (ordinal_type i=0; i<m; i++)
475 result(j,i) = beta*result(j,i) + alpha*tmp1[i+j*m];
476 }
477 }
478 else {
479 if (reorder_result) {
480 ordinal_type idx;
481 for (ordinal_type j=0; j<p; j++)
482 for (ordinal_type i=0; i<m; i++) {
483 idx = perm[i];
484 result(idx,j) = beta*result(idx,j) + alpha*tmp1[i+j*m];
485 }
486 }
487 else {
488 for (ordinal_type j=0; j<p; j++)
489 for (ordinal_type i=0; i<m; i++)
490 result(i,j) = beta*result(i,j) + alpha*tmp1[i+j*m];
491 }
492 }
493 }
494
495 protected:
496
499
502
505
508
511
514
516 Teuchos::Array<ordinal_type> perm;
517
519 Teuchos::SerialDenseMatrix<ordinal_type,value_type> qp2pce;
520
522 Teuchos::SerialDenseMatrix<ordinal_type,value_type> pce2qp;
523
525 Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> > qp2pce_k;
526
528 Teuchos::Array< Teuchos::SerialDenseMatrix<ordinal_type,value_type> > pce2qp_k;
529
531 Teuchos::BLAS<ordinal_type,value_type> blas;
532
533 };
534
535}
536
537#endif
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::const_iterator const_set_iterator
std::map< point_type, std::pair< value_type, ordinal_type >, point_compare_type > point_set_type
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.
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.
const_iterator end() const
Iterator to end of point set.
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.
const_set_iterator set_end() const
Iterator to end of point set.
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.
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.
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.
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.
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.
Teuchos::Array< ordinal_type > perm
Permutation array when reordering for PST.
Top-level namespace for Stokhos classes and functions.
LexographicLess< TensorProductElement< ordinal_type, value_type >, FloatingPointLess< value_type > > type