Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
44 template <typename ordinal_type, typename value_type>
47  ordinal_type max_p,
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 
63 template <typename ordinal_type, typename value_type>
64 void
67  ordinal_type max_p,
70 {
71  // Check for pce's that are constant and don't represent true random
72  // dimensions
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();
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);
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") {
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 
141 template <typename ordinal_type, typename value_type>
144 {
145 }
146 
147 template <typename ordinal_type, typename value_type>
150 order() const
151 {
152  return p;
153 }
154 
155 template <typename ordinal_type, typename value_type>
158 dimension() const
159 {
160  return d;
161 }
162 
163 template <typename ordinal_type, typename value_type>
166 size() const
167 {
168  return sz;
169 }
170 
171 template <typename ordinal_type, typename value_type>
175 {
176  return norms;
177 }
178 
179 template <typename ordinal_type, typename value_type>
180 const value_type&
183 {
184  return norms[i];
185 }
186 
187 template <typename ordinal_type, typename value_type>
191 
192 {
193  return Teuchos::null;
194 }
195 
196 template <typename ordinal_type, typename value_type>
200 
201 {
202  return Teuchos::null;
203 }
204 
205 template <typename ordinal_type, typename value_type>
209 {
210  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
211 }
212 
213 template <typename ordinal_type, typename value_type>
214 void
217  Teuchos::Array<value_type>& basis_vals) const
218 {
219  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
220 }
221 
222 template <typename ordinal_type, typename value_type>
223 void
225 print(std::ostream& os) const
226 {
227  os << "Gram-Schmidt basis of order " << p << ", dimension " << d
228  << ", and size " << sz << ". Matrix coefficients:\n";
229  os << 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 
236 template <typename ordinal_type, typename value_type>
237 void
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 
258 template <typename ordinal_type, typename value_type>
259 void
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 
280 template <typename ordinal_type, typename value_type>
284 {
285  return reduced_quad;
286 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
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.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
T & get(ParameterList &l, const std::string &name)
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)
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
bool isParameter(const std::string &name) const
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 void print(std::ostream &os) const
Print basis to stream os.
Abstract base class for quadrature methods.
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.
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 ordinal_type size() const
Return total size of basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
ordinal_type order() const
Return order of basis.
void push_back(const value_type &x)
size_type size() const
ordinal_type dimension() const
Return dimension of basis.
#define TEUCHOS_ASSERT(assertion_test)
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
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.