Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GSReducedPCEBasisBaseImp.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
43
44template <typename ordinal_type, typename value_type>
47 ordinal_type max_p,
48 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
49 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
50 const Teuchos::ParameterList& params_) :
51 params(params_),
52 pce_basis(pce[0].basis()),
53 pce_sz(pce_basis->size()),
54 p(max_p),
55 d(pce.size()),
56 verbose(params.get("Verbose", false)),
57 rank_threshold(params.get("Rank Threshold", 1.0e-12)),
58 orthogonalization_method(params.get("Orthogonalization Method",
59 "Householder"))
60{
61}
62
63template <typename ordinal_type, typename value_type>
64void
66setup(
67 ordinal_type max_p,
68 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
69 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad)
70{
71 // Check for pce's that are constant and don't represent true random
72 // dimensions
73 Teuchos::Array< const Stokhos::OrthogPolyApprox<ordinal_type, value_type>* > pce2;
74 for (ordinal_type i=0; i<pce.size(); i++) {
75 if (pce[i].standard_deviation() > 1.0e-15)
76 pce2.push_back(&pce[i]);
77 }
78 d = pce2.size();
79
80 // Get quadrature data
81 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
82 const Teuchos::Array< Teuchos::Array<value_type> >& points =
83 quad->getQuadPoints();
84 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
85 quad->getBasisAtQuadPoints();
86 ordinal_type nqp = weights.size();
87
88 // Original basis at quadrature points -- needed to transform expansions
89 // in this basis back to original
90 SDM A(nqp, pce_sz);
91 for (ordinal_type i=0; i<nqp; i++)
92 for (ordinal_type j=0; j<pce_sz; j++)
93 A(i,j) = basis_values[i][j];
94
95 // Compute norms of each pce for rescaling
96 Teuchos::Array<value_type> pce_norms(d, 0.0);
97 for (ordinal_type j=0; j<d; j++) {
98 for (ordinal_type i=0; i<pce_sz; i++)
99 pce_norms[j] += (*pce2[j])[i]*(*pce2[j])[i]*pce_basis->norm_squared(i);
100 pce_norms[j] = std::sqrt(pce_norms[j]);
101 }
102
103 // Compute F matrix -- PCEs evaluated at all quadrature points
104 // Since F is used in the reduced quadrature below as the quadrature points
105 // for this reduced basis, does scaling by the pce_norms mess up the points?
106 // No -- F essentially defines the random variables this basis is a function
107 // of, and thus they can be scaled in any way we want. Because we don't
108 // explicitly write the basis in terms of F, the scaling is implicit.
109 SDM F(nqp, d);
110 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
111 for (ordinal_type i=0; i<nqp; i++)
112 for (ordinal_type j=0; j<d; j++)
113 F(i,j) = pce2[j]->evaluate(points[i], basis_values[i]);
114
115 // Build the reduced basis
116 sz = buildReducedBasis(max_p, rank_threshold, A, F, weights, terms, num_terms,
117 Qp, Q);
118
119 // Compute reduced quadrature rule
120 Teuchos::ParameterList quad_params = params.sublist("Reduced Quadrature");
122 quad_params);
123 SDM Q2;
124 if (quad_params.isParameter("Reduced Quadrature Method") &&
125 quad_params.get<std::string>("Reduced Quadrature Method") == "Q2") {
126 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms2;
127 Teuchos::Array<ordinal_type> num_terms2;
128 value_type rank_threshold2 = quad_params.get("Q2 Rank Threshold",
129 rank_threshold);
130 SDM Qp2;
131 //ordinal_type sz2 =
132 buildReducedBasis(2*max_p, rank_threshold2, A, F, weights, terms2,
133 num_terms2, Qp2, Q2);
134 }
135 reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F, weights);
136
137 // Basis is orthonormal by construction
138 norms.resize(sz, 1.0);
139}
140
141template <typename ordinal_type, typename value_type>
144{
145}
146
147template <typename ordinal_type, typename value_type>
148ordinal_type
150order() const
151{
152 return p;
153}
154
155template <typename ordinal_type, typename value_type>
156ordinal_type
158dimension() const
159{
160 return d;
161}
162
163template <typename ordinal_type, typename value_type>
164ordinal_type
166size() const
167{
168 return sz;
169}
170
171template <typename ordinal_type, typename value_type>
172const Teuchos::Array<value_type>&
174norm_squared() const
175{
176 return norms;
177}
178
179template <typename ordinal_type, typename value_type>
180const value_type&
182norm_squared(ordinal_type i) const
183{
184 return norms[i];
185}
186
187template <typename ordinal_type, typename value_type>
188Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
191
192{
193 return Teuchos::null;
194}
195
196template <typename ordinal_type, typename value_type>
197Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
200
201{
202 return Teuchos::null;
203}
204
205template <typename ordinal_type, typename value_type>
206value_type
208evaluateZero(ordinal_type i) const
209{
210 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
211}
212
213template <typename ordinal_type, typename value_type>
214void
216evaluateBases(const Teuchos::ArrayView<const value_type>& point,
217 Teuchos::Array<value_type>& basis_vals) const
218{
219 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
220}
221
222template <typename ordinal_type, typename value_type>
223void
225print(std::ostream& os) const
226{
227 os << "Gram-Schmidt basis of order " << p << ", dimension " << d
228 << ", and size " << sz << ". Matrix coefficients:\n";
229 os << printMat(Qp) << std::endl;
230 os << "Basis vector norms (squared):\n\t";
231 for (ordinal_type i=0; i<sz; i++)
232 os << norms[i] << " ";
233 os << "\n";
234}
235
236template <typename ordinal_type, typename value_type>
237void
239transformToOriginalBasis(const value_type *in, value_type *out,
240 ordinal_type ncol, bool transpose) const
241{
242 if (transpose) {
243 SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
244 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
245 ordinal_type ret =
246 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
247 TEUCHOS_ASSERT(ret == 0);
248 }
249 else {
250 SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
251 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
252 ordinal_type ret =
253 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
254 TEUCHOS_ASSERT(ret == 0);
255 }
256}
257
258template <typename ordinal_type, typename value_type>
259void
261transformFromOriginalBasis(const value_type *in, value_type *out,
262 ordinal_type ncol, bool transpose) const
263{
264 if (transpose) {
265 SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
266 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
267 ordinal_type ret =
268 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
269 TEUCHOS_ASSERT(ret == 0);
270 }
271 else {
272 SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
273 SDM zbar(Teuchos::View, out, sz, sz, ncol);
274 ordinal_type ret =
275 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
276 TEUCHOS_ASSERT(ret == 0);
277 }
278}
279
280template <typename ordinal_type, typename value_type>
281Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
284{
285 return reduced_quad;
286}
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
ordinal_type dimension() const
Return dimension of basis.
GSReducedPCEBasisBase(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::ParameterList &params=Teuchos::ParameterList())
Constructor.
virtual void transformFromOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients from original basis to this basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
void setup(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad)
virtual ordinal_type size() const
Return total size of basis.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual void print(std::ostream &os) const
Print basis to stream os.
ordinal_type order() const
Return order of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual void transformToOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients to original basis from this basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Abstract base class for quadrature methods.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< const Stokhos::UserDefinedQuadrature< ordinal_type, value_type > > createReducedQuadrature(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights) const
Get reduced quadrature object.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)