Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
46 #include "Stokhos_ProductBasis.hpp"
49 #include "Teuchos_BLAS.hpp"
50 
51 namespace Stokhos {
52 
57  template <typename ordinal_t,
58  typename value_t,
59  typename point_compare_type =
60  typename DefaultPointCompare<ordinal_t,value_t>::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;
68  typedef typename base_type::point_type point_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 
102 
103  // Make sure order of each 1-D basis is large enough for
104  // supplied set of coefficients
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
120  multiindex_type n(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();
152  point_type point(dim);
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 
251  const_set_iterator it = points.find(point);
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,
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,
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 
314  const value_type& alpha,
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());
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());
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 
366  void apply_pst(
368  const value_type& alpha,
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)
442  Teuchos::View, &tmp1[0], n, n, nk*p);
443 
444  // y = x'
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
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 
498  bool use_pst;
499 
502 
505 
508 
511 
513  bool reorder;
514 
517 
520 
523 
526 
529 
532 
533  };
534 
535 }
536 
537 #endif
ScalarType * values() const
std::map< point_type, std::pair< value_type, ordinal_type >, point_compare_type > point_set_type
const point_type & point(ordinal_type n) const
Get point for given index.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > qp2pce
Matrix mapping points to coefficients.
Container storing a term in a generalized tensor product.
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > qp2pce_k
Matrix mapping points to coefficients for each dimension for PST.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
An operator interface for building pseudo-spectral approximations.
set_iterator set_begin()
Iterator to begining of point set.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
const_iterator end() const
Iterator to end of point set.
Teuchos::BLAS< ordinal_type, value_type > blas
BLAS wrappers.
PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > base_type
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const =0
Get orders of each coordinate polynomial given an index i.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
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.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
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 resize(size_type new_size, const value_type &x=value_type())
const_iterator begin() const
Return iterator for first element in the set.
const_set_iterator set_begin() const
Iterator to begining of point set.
const_iterator begin() const
Iterator to begining of point set.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > pce2qp
Matrix mapping coefficients to points.
iterator end()
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.
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
OrdinalType numCols() const
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.
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.
int reshape(OrdinalType numRows, OrdinalType numCols)
size_type size() const
A comparison functor implementing a strict weak ordering based lexographic ordering.
Iterator class for iterating over elements of the index set.
An operator for building pseudo-spectral coefficients using tensor-product quadrature.
#define TEUCHOS_ASSERT(assertion_test)
const_set_iterator set_end() const
Iterator to end of point set.
Teuchos::Array< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > pce2qp_k
Matrix mapping coefficients to points for each dimension for PST.
const_iterator end() const
Return iterator for end of the index set.
iterator begin()
int n
ordinal_type index(const point_type &point) const
Get point index for given point.
Teuchos::Array< ordinal_type > perm
Permutation array when reordering for PST.
virtual ordinal_type size() const =0
Return total size of basis.
OrdinalType stride() const
virtual MultiIndex< ordinal_type > getMaxOrders() const =0
Return maximum order allowable for each coordinate basis.
OrdinalType numRows() const