Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_LanczosPCEBasisImp.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 "Teuchos_Assert.hpp"
43 #include "Teuchos_BLAS.hpp"
44 #include "Teuchos_TimeMonitor.hpp"
45 
46 template <typename ordinal_type, typename value_type>
49  ordinal_type p,
52  bool normalize,
53  bool limit_integration_order_) :
54  RecurrenceBasis<ordinal_type, value_type>("Lanczos PCE", p, normalize),
55  pce(pce_),
56  quad(quad_),
57  limit_integration_order(limit_integration_order_),
58  nqp(quad->size()),
59  pce_weights(Teuchos::Copy,
60  const_cast<value_type*>(quad->getQuadWeights().getRawPtr()),
61  nqp),
62  pce_vals(nqp),
63  u0(nqp),
64  lanczos_vecs(nqp, p+1),
65  fromStieltjesMat(),
66  new_pce()
67 {
68  // Evaluate PCE at quad points
69  const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
70  quad->getQuadPoints();
71  const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
72  quad->getBasisAtQuadPoints();
73  for (ordinal_type i=0; i<nqp; i++) {
74  pce_vals[i] = pce->evaluate(quad_points[i], basis_values[i]);
75  u0[i] = value_type(1);
76  }
77 
78  // Setup rest of basis
79  this->setup();
80 }
81 
82 template <typename ordinal_type, typename value_type>
85 {
86 }
87 
88 template <typename ordinal_type, typename value_type>
89 void
92  Teuchos::Array<value_type>& quad_points,
93  Teuchos::Array<value_type>& quad_weights,
94  Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
95 {
96 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
97  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::LanczosPCEBasis -- compute Gauss points");
98 #endif
99 
100  // Call base class
101  ordinal_type num_points =
102  static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
103 
104  // We can't always reliably generate quadrature points of order > 2*p
105  // when using sparse grids for the underlying quadrature
106  if (limit_integration_order && quad_order > 2*this->p)
107  quad_order = 2*this->p;
109  quad_points,
110  quad_weights,
111  quad_values);
112 
113  // Fill in the rest of the points with zero weight
114  if (quad_weights.size() < num_points) {
115  ordinal_type old_size = quad_weights.size();
116  quad_weights.resize(num_points);
117  quad_points.resize(num_points);
118  quad_values.resize(num_points);
119  for (ordinal_type i=old_size; i<num_points; i++) {
120  quad_weights[i] = value_type(0);
121  quad_points[i] = quad_points[0];
122  quad_values[i].resize(this->p+1);
123  this->evaluateBases(quad_points[i], quad_values[i]);
124  }
125  }
126 }
127 
128 template <typename ordinal_type, typename value_type>
132 {
134  p,*this));
135 }
136 
137 template <typename ordinal_type, typename value_type>
141 {
142  return new_pce[i];
143 }
144 
145 template <typename ordinal_type, typename value_type>
146 void
149 {
151  ordinal_type sz = fromStieltjesMat.numRows();
152  blas.GEMV(Teuchos::NO_TRANS, sz, this->p+1,
153  value_type(1.0), fromStieltjesMat.values(), sz,
154  in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
155 }
156 
157 template <typename ordinal_type, typename value_type>
158 bool
164  Teuchos::Array<value_type>& gamma) const
165 {
167  vectorspace_type vs(pce_weights);
168  operator_type A(pce_vals);
169 
170  // Create space to store lanczos vectors -- use lanczos_vecs if
171  // we are requesting p+1 vectors
173  if (n == this->p+1)
174  lv = Teuchos::rcp(&lanczos_vecs, false);
175  else
176  lv = Teuchos::rcp(new matrix_type(nqp,n));
177 
178  if (this->normalize)
179  lanczos_type::computeNormalized(n, vs, A, u0, *lv, alpha, beta, nrm);
180  else
181  lanczos_type::compute(n, vs, A, u0, *lv, alpha, beta, nrm);
182 
183  for (ordinal_type i=0; i<n; i++) {
184  delta[i] = value_type(1.0);
185  }
186  if (this->normalize)
187  gamma = beta;
188  else
189  for (ordinal_type i=0; i<n; i++)
190  gamma[i] = value_type(1.0);
191 
192  return this->normalize;
193 }
194 
195 template <typename ordinal_type, typename value_type>
196 void
199 {
201 
202  // Compute transformation matrix back to original basis
203  ordinal_type sz = pce->size();
204  fromStieltjesMat.shape(sz, this->p+1);
205  fromStieltjesMat.putScalar(0.0);
206  const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
207  quad->getBasisAtQuadPoints();
208  for (ordinal_type i=0; i<sz; i++) {
209  for (ordinal_type j=0; j<=this->p; j++) {
210  for (ordinal_type k=0; k<nqp; k++)
211  fromStieltjesMat(i,j) +=
212  pce_weights[k]*lanczos_vecs(k,j)*basis_values[k][i];
213  fromStieltjesMat(i,j) /= pce->basis()->norm_squared(i);
214  }
215  }
216 
217  // Project original PCE into the new basis
218  new_pce.resize(this->p+1);
219  vector_type u(sz);
220  for (ordinal_type i=0; i<sz; i++)
221  u[i] = (*pce)[i]*pce->basis()->norm_squared(i);
222  new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, fromStieltjesMat, u,
223  0.0);
224  for (ordinal_type i=0; i<=this->p; i++)
225  new_pce[i] /= this->norms[i];
226 }
227 
228 template <typename ordinal_type, typename value_type>
232  pce(basis.pce),
233  quad(basis.quad),
234  limit_integration_order(basis.limit_integration_order),
235  nqp(basis.nqp),
236  pce_weights(basis.pce_weights),
237  pce_vals(basis.pce_vals),
238  u0(basis.u0),
239  lanczos_vecs(nqp, p+1),
240  fromStieltjesMat(),
241  new_pce()
242 {
243  this->setup();
244 }
vector_type pce_vals
Values of PCE at quadrature points.
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points...
RawPointerConversionTraits< Container >::Ptr_t getRawPtr(const Container &c)
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Get Gauss quadrature points, weights, and values of basis at points.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract base class for quadrature methods.
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
LanczosPCEBasis(ordinal_type p, const Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool normalize, bool limit_integration_order)
Constructor.
ordinal_type nqp
Number of quadrature points.
void resize(size_type new_size, const value_type &x=value_type())
virtual void setup()
Setup basis after computing recurrence coefficients.
Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > pce
PCE Lanczos procedure is based on.
size_type size() const
virtual bool computeRecurrenceCoefficients(ordinal_type n, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Compute recurrence coefficients.
vector_type u0
Initial Lanczos vector.
void transformCoeffsFromLanczos(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
virtual void setup()
Setup basis after computing recurrence coefficients.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > quad
Quadrature object.