Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_QuadraturePseudoSpectralOperator.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_QUADRATURE_PSEUDO_SPECTRAL_OPERATOR_HPP
43#define STOKHOS_QUADRATURE_PSEUDO_SPECTRAL_OPERATOR_HPP
44
47#include "Teuchos_SerialDenseMatrix.hpp"
48#include "Teuchos_BLAS.hpp"
49
50namespace Stokhos {
51
56 template <typename ordinal_t,
57 typename value_t,
58 typename point_compare_type =
61 public PseudoSpectralOperator<ordinal_t,value_t,point_compare_type> {
62 public:
63
64 typedef ordinal_t ordinal_type;
65 typedef value_t value_type;
66 typedef PseudoSpectralOperator<ordinal_type,value_type,point_compare_type> base_type;
70 typedef typename base_type::iterator iterator;
74
75 typedef MultiIndex<ordinal_type> multiindex_type;
76
79 const OrthogPolyBasis<ordinal_type,value_type>& basis,
80 const Quadrature<ordinal_type,value_type>& quad,
81 const point_compare_type& point_compare = point_compare_type()) :
82 coeff_sz(basis.size()),
83 points(point_compare) {
84
85 const Teuchos::Array<value_type>& quad_weights = quad.getQuadWeights();
86 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
87 quad.getQuadPoints();
88 const Teuchos::Array< Teuchos::Array<value_type> > & quad_vals =
89 quad.getBasisAtQuadPoints();
90
91 ordinal_type nqp = quad.size();
92 ordinal_type npc = basis.size();
93 ordinal_type dim = basis.dimension();
94
95 // Generate point set
96 point_type thePoint(dim);
97 for (ordinal_type i=0; i<nqp; ++i) {
98 for (ordinal_type k=0; k<dim; k++) {
99 thePoint[k] = quad_points[i][k];
100 }
101 points[thePoint] = std::make_pair(quad_weights[i],i);
102 }
103
104 // Generate quadrature operator
105 // This has to happen before we change the global index since
106 // the point set likely will reorder the points
107 qp2pce.reshape(npc,nqp);
108 pce2qp.reshape(nqp,npc);
109 typename point_set_type::iterator di = points.begin();
110 typename point_set_type::iterator di_end = points.end();
111 ordinal_type jdx = 0;
112 for (; di != di_end; ++di) {
113 ordinal_type j = di->second.second;
114 for (ordinal_type i=0; i<npc; i++) {
115 qp2pce(i,jdx) = quad_weights[j]*quad_vals[j][i] /
116 basis.norm_squared(i);
117 pce2qp(jdx,i) = quad_vals[j][i];
118 }
119 ++jdx;
120 }
121
122 // Generate linear ordering of points
123 point_map.resize(nqp);
124 ordinal_type idx=0;
125 di = points.begin();
126 while (di != di_end) {
127 di->second.second = idx;
128 point_map[idx] = di->first;
129 ++idx;
130 ++di;
131 }
132
133
134 }
135
138
140 ordinal_type point_size() const { return points.size(); }
141
143 ordinal_type coeff_size() const { return coeff_sz; }
144
146 iterator begin() { return point_map.begin(); }
147
149 iterator end() { return point_map.end(); }
150
152 const_iterator begin() const { return point_map.begin(); }
153
155 const_iterator end() const { return point_map.end(); }
156
158 set_iterator set_begin() { return points.begin(); }
159
161 set_iterator set_end() { return points.end(); }
162
164 const_set_iterator set_begin() const { return points.begin(); }
165
167 const_set_iterator set_end() const { return points.end(); }
168
172 TEUCHOS_TEST_FOR_EXCEPTION(
173 it == points.end(), std::logic_error, "Invalid term " << point);
174 return it->second.second;
175 }
176
178 const point_type& point(ordinal_type n) const { return point_map[n]; }
179
181
188 virtual void transformQP2PCE(
189 const value_type& alpha,
190 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
191 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
192 const value_type& beta,
193 bool trans = false) const {
194 ordinal_type ret;
195 if (trans)
196 ret = result.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, alpha,
197 input, qp2pce, beta);
198 else
199 ret = result.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha,
200 qp2pce, input, beta);
201 TEUCHOS_ASSERT(ret == 0);
202 }
203
205
212 virtual void transformPCE2QP(
213 const value_type& alpha,
214 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
215 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result,
216 const value_type& beta,
217 bool trans = false) const {
218 if (trans) {
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(),
224 input.stride(), pce2qp.values(), pce2qp.stride(), beta,
225 result.values(), result.stride());
226 }
227 else {
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());
235 }
236 }
237
238 protected:
239
242
245
248
250 Teuchos::SerialDenseMatrix<ordinal_type,value_type> qp2pce;
251
253 Teuchos::SerialDenseMatrix<ordinal_type,value_type> pce2qp;
254
256 Teuchos::BLAS<ordinal_type,value_type> blas;
257
258 };
259
260}
261
262#endif
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
An operator for building pseudo-spectral coefficients using an arbitrary quadrature rule.
const_set_iterator set_begin() const
Iterator to begining of point set.
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.
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.
Teuchos::BLAS< ordinal_type, value_type > blas
BLAS wrappers.
const_set_iterator set_end() const
Iterator to end of point set.
ordinal_type index(const point_type &point) const
Get point index for given point.
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.
const_iterator end() const
Iterator to end of point set.
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 > pce2qp
Matrix mapping coefficients to points.
PseudoSpectralOperator< ordinal_type, value_type, point_compare_type > base_type
set_iterator set_begin()
Iterator to begining of point set.
const_iterator begin() const
Iterator to begining of point set.
Container storing a term in a generalized tensor product.
Top-level namespace for Stokhos classes and functions.
LexographicLess< TensorProductElement< ordinal_type, value_type >, FloatingPointLess< value_type > > type