Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_ProductLanczosGramSchmidtPCEBasisImp.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 #include "Stokhos_SDMUtils.hpp"
47 #include "Stokhos_ProductBasis.hpp"
49 
50 template <typename ordinal_type, typename value_type>
53  ordinal_type max_p,
54  const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
55  const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
56  const Teuchos::RCP< const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk,
57  const Teuchos::ParameterList& params_) :
58  name("Product Lanczos Gram-Schmidt PCE Basis"),
59  params(params_),
60  p(max_p)
61 {
62  Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis = pce[0].basis();
63  pce_sz = pce_basis->size();
64 
65  // Check if basis is a product basis
66  Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
67  Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
68  if (prod_basis != Teuchos::null)
69  coord_bases = prod_basis->getCoordinateBases();
70 
71  // Build Lanczos basis for each pce
72  bool project = params.get("Project", true);
73  bool normalize = params.get("Normalize", true);
74  bool limit_integration_order = params.get("Limit Integration Order", false);
75  bool use_stieltjes = params.get("Use Old Stieltjes Method", false);
76  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double > > > coordinate_bases;
77  Teuchos::Array<int> is_invariant(pce.size(),-2);
78  for (ordinal_type i=0; i<pce.size(); i++) {
79 
80  // Check for pce's lying in invariant subspaces, which are pce's that
81  // depend on only a single dimension. In this case use the corresponding
82  // original coordinate basis. Convention is: -2 -- not invariant, -1 --
83  // constant, i >= 0 pce depends only on dimension i
84  if (prod_basis != Teuchos::null)
85  is_invariant[i] = isInvariant(pce[i]);
86  if (is_invariant[i] >= 0) {
87  coordinate_bases.push_back(coord_bases[is_invariant[i]]);
88  }
89 
90  // Exclude constant pce's from the basis since they don't represent
91  // stochastic dimensions
92  else if (is_invariant[i] != -1) {
93  if (use_stieltjes) {
94  coordinate_bases.push_back(
95  Teuchos::rcp(
97  p, Teuchos::rcp(&(pce[i]),false), quad, false,
98  normalize, project, Cijk)));
99  }
100  else {
101  if (project)
102  coordinate_bases.push_back(
103  Teuchos::rcp(
105  p, Teuchos::rcp(&(pce[i]),false), Cijk,
106  normalize, limit_integration_order)));
107  else
108  coordinate_bases.push_back(
109  Teuchos::rcp(
111  p, Teuchos::rcp(&(pce[i]),false), quad,
112  normalize, limit_integration_order)));
113  }
114  }
115  }
116  d = coordinate_bases.size();
117 
118  // Build tensor product basis
120  Teuchos::rcp(
122  coordinate_bases,
123  params.get("Cijk Drop Tolerance", 1.0e-15),
124  params.get("Use Old Cijk Algorithm", false)));
125 
126  // Build Psi matrix -- Psi_ij = Psi_i(x^j)*w_j/<Psi_i^2>
127  const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
128  const Teuchos::Array< Teuchos::Array<value_type> >& points =
129  quad->getQuadPoints();
130  const Teuchos::Array< Teuchos::Array<value_type> >& basis_vals =
131  quad->getBasisAtQuadPoints();
132  ordinal_type nqp = weights.size();
133  SDM Psi(pce_sz, nqp);
134  for (ordinal_type i=0; i<pce_sz; i++)
135  for (ordinal_type k=0; k<nqp; k++)
136  Psi(i,k) = basis_vals[k][i]*weights[k]/pce_basis->norm_squared(i);
137 
138  // Build Phi matrix -- Phi_ij = Phi_i(y(x^j))
139  sz = tensor_lanczos_basis->size();
140  Teuchos::Array<value_type> red_basis_vals(sz);
141  Teuchos::Array<value_type> pce_vals(d);
142  SDM Phi(sz, nqp);
143  SDM F(nqp, d);
144  for (int k=0; k<nqp; k++) {
145  ordinal_type jdx = 0;
146  for (int j=0; j<pce.size(); j++) {
147 
148  // Exclude constant pce's
149  if (is_invariant[j] != -1) {
150 
151  // Use the identity mapping for invariant subspaces
152  if (is_invariant[j] >= 0)
153  pce_vals[jdx] = points[k][is_invariant[j]];
154  else
155  pce_vals[jdx] = pce[j].evaluate(points[k], basis_vals[k]);
156  F(k,jdx) = pce_vals[jdx];
157  jdx++;
158 
159  }
160 
161  }
162  tensor_lanczos_basis->evaluateBases(pce_vals, red_basis_vals);
163  for (int i=0; i<sz; i++)
164  Phi(i,k) = red_basis_vals[i];
165  }
166 
167  bool verbose = params.get("Verbose", false);
168 
169  // Compute matrix A mapping reduced space to original
170  SDM A(pce_sz, sz);
171  ordinal_type ret =
172  A.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, Psi, Phi, 0.0);
173  TEUCHOS_ASSERT(ret == 0);
174 
175  // Rescale columns of A to have unit norm
176  const Teuchos::Array<value_type>& basis_norms = pce_basis->norm_squared();
177  for (ordinal_type j=0; j<sz; j++) {
178  value_type nrm = 0.0;
179  for (ordinal_type i=0; i<pce_sz; i++)
180  nrm += A(i,j)*A(i,j)*basis_norms[i];
181  nrm = std::sqrt(nrm);
182  for (ordinal_type i=0; i<pce_sz; i++)
183  A(i,j) /= nrm;
184  }
185 
186  // Compute our new basis -- each column of Qp is the coefficients of the
187  // new basis in the original basis. Constraint pivoting so first d+1
188  // columns and included in Qp.
189  value_type rank_threshold = params.get("Rank Threshold", 1.0e-12);
190  std::string orthogonalization_method =
191  params.get("Orthogonalization Method", "Householder");
192  Teuchos::Array<value_type> w(pce_sz, 1.0);
193  SDM R;
194  Teuchos::Array<ordinal_type> piv(sz);
195  for (int i=0; i<d+1; i++)
196  piv[i] = 1;
198  sz = SOF::createOrthogonalBasis(
199  orthogonalization_method, rank_threshold, verbose, A, w, Qp, R, piv);
200 
201  // Original basis at quadrature points -- needed to transform expansions
202  // in this basis back to original
203  SDM B(nqp, pce_sz);
204  for (ordinal_type i=0; i<nqp; i++)
205  for (ordinal_type j=0; j<pce_sz; j++)
206  B(i,j) = basis_vals[i][j];
207 
208  // Evaluate new basis at original quadrature points
209  Q.reshape(nqp, sz);
210  ret = Q.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, B, Qp, 0.0);
211  TEUCHOS_ASSERT(ret == 0);
212 
213  // Compute reduced quadrature rule
215  params.sublist("Reduced Quadrature"));
216  Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2;
217  reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F, weights);
218 
219  // Basis is orthonormal by construction
220  norms.resize(sz, 1.0);
221 }
222 
223 template <typename ordinal_type, typename value_type>
226 {
227 }
228 
229 template <typename ordinal_type, typename value_type>
232 order() const
233 {
234  return p;
235 }
236 
237 template <typename ordinal_type, typename value_type>
240 dimension() const
241 {
242  return d;
243 }
244 
245 template <typename ordinal_type, typename value_type>
248 size() const
249 {
250  return sz;
251 }
252 
253 template <typename ordinal_type, typename value_type>
254 const Teuchos::Array<value_type>&
257 {
258  return norms;
259 }
260 
261 template <typename ordinal_type, typename value_type>
262 const value_type&
265 {
266  return norms[i];
267 }
268 
269 template <typename ordinal_type, typename value_type>
270 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
273 
274 {
275  return Teuchos::null;
276 }
277 
278 template <typename ordinal_type, typename value_type>
279 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
282 
283 {
284  return Teuchos::null;
285 }
286 
287 template <typename ordinal_type, typename value_type>
291 {
292  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
293 }
294 
295 template <typename ordinal_type, typename value_type>
296 void
298 evaluateBases(const Teuchos::ArrayView<const value_type>& point,
299  Teuchos::Array<value_type>& basis_vals) const
300 {
301  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
302 }
303 
304 template <typename ordinal_type, typename value_type>
305 void
307 print(std::ostream& os) const
308 {
309  os << "Product Lanczos Gram-Schmidt basis of order " << p << ", dimension "
310  << d
311  << ", and size " << sz << ". Matrix coefficients:\n";
312  os << printMat(Qp) << std::endl;
313  os << "Basis vector norms (squared):\n\t";
314  for (ordinal_type i=0; i<sz; i++)
315  os << norms[i] << " ";
316  os << "\n";
317 }
318 
319 template <typename ordinal_type, typename value_type>
320 const std::string&
322 getName() const
323 {
324  return name;
325 }
326 
327 template <typename ordinal_type, typename value_type>
328 void
331  ordinal_type ncol, bool transpose) const
332 {
333  if (transpose) {
334  SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
335  SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
336  ordinal_type ret =
337  z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
338  TEUCHOS_ASSERT(ret == 0);
339  }
340  else {
341  SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
342  SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
343  ordinal_type ret =
344  z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
345  TEUCHOS_ASSERT(ret == 0);
346  }
347 }
348 
349 template <typename ordinal_type, typename value_type>
350 void
353  ordinal_type ncol, bool transpose) const
354 {
355  if (transpose) {
356  SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
357  SDM zbar(Teuchos::View, out, ncol, ncol, sz);
358  ordinal_type ret =
359  zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
360  TEUCHOS_ASSERT(ret == 0);
361  }
362  else {
363  SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
364  SDM zbar(Teuchos::View, out, sz, sz, ncol);
365  ordinal_type ret =
366  zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
367  TEUCHOS_ASSERT(ret == 0);
368  }
369 }
370 
371 template <typename ordinal_type, typename value_type>
372 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
375 {
376  return reduced_quad;
377 }
378 
379 template <typename ordinal_type, typename value_type>
383 {
384  Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
385  pce.basis();
386  ordinal_type dim = basis->dimension();
387  value_type tol = 1.0e-15;
388 
389  // Check if basis is a product basis
390  Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(basis);
391  if (prod_basis == Teuchos::null)
392  return -2;
393 
394  // Build list of dimensions pce depends on by looping over each dimension,
395  // computing norm of pce with just that dimension -- note we don't include
396  // the constant term
397  Teuchos::Array<ordinal_type> dependent_dims;
398  tmp_pce.reset(basis);
399  for (ordinal_type i=0; i<dim; i++) {
400  ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
401  tmp_pce.init(0.0);
402  for (ordinal_type j=1; j<=p; j++)
403  tmp_pce.term(i,j) = pce.term(i,j);
404  value_type nrm = tmp_pce.two_norm();
405  if (nrm > tol) dependent_dims.push_back(i);
406  }
407 
408  // If dependent_dims has length 1, pce a function of a single variable,
409  // which is an invariant subspace
410  if (dependent_dims.size() == 1)
411  return dependent_dims[0];
412 
413  // If dependent_dims has length 0, pce is constant
414  else if (dependent_dims.size() == 0)
415  return -1;
416 
417  // Otherwise pce depends on more than one variable
418  return -2;
419 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
virtual const std::string & getName() const
Return string name of basis.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
SDM Qp
Coefficients of transformed basis in original basis.
SDM Q
Values of transformed basis at quadrature points.
virtual ordinal_type size() const
Return total size of basis.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
ordinal_type isInvariant(const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce) const
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
virtual void print(std::ostream &os) const
Print basis to stream os.
Abstract base class for quadrature methods.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
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.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
Teuchos::RCP< Stokhos::CompletePolynomialBasis< ordinal_type, value_type > > tensor_lanczos_basis
Product Lanczos basis.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
ProductLanczosGramSchmidtPCEBasis(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::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
Constructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
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.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Encapsulate various orthogonalization (ie QR) methods.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.