42 #include "Teuchos_Assert.hpp" 43 #include "Teuchos_BLAS.hpp" 44 #include "Teuchos_TimeMonitor.hpp" 46 template <
typename ordinal_type,
typename value_type>
53 bool limit_integration_order_) :
57 limit_integration_order(limit_integration_order_),
60 const_cast<
value_type*>(quad->getQuadWeights().getRawPtr()),
64 lanczos_vecs(nqp, p+1),
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();
74 pce_vals[i] =
pce->evaluate(quad_points[i], basis_values[i]);
82 template <
typename ordinal_type,
typename value_type>
88 template <
typename ordinal_type,
typename value_type>
92 Teuchos::Array<value_type>& quad_points,
93 Teuchos::Array<value_type>& quad_weights,
94 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const 96 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 97 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::LanczosPCEBasis -- compute Gauss points");
106 if (limit_integration_order && quad_order > 2*this->p)
107 quad_order = 2*this->p;
114 if (quad_weights.size() < num_points) {
116 quad_weights.resize(num_points);
117 quad_points.resize(num_points);
118 quad_values.resize(num_points);
121 quad_points[i] = quad_points[0];
122 quad_values[i].resize(this->p+1);
123 this->evaluateBases(quad_points[i], quad_values[i]);
128 template <
typename ordinal_type,
typename value_type>
129 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
137 template <
typename ordinal_type,
typename value_type>
145 template <
typename ordinal_type,
typename value_type>
150 Teuchos::BLAS<ordinal_type, value_type> blas;
152 blas.GEMV(Teuchos::NO_TRANS, sz, this->p+1,
153 value_type(1.0), fromStieltjesMat.values(), sz,
157 template <
typename ordinal_type,
typename value_type>
161 Teuchos::Array<value_type>& alpha,
162 Teuchos::Array<value_type>& beta,
163 Teuchos::Array<value_type>& delta,
164 Teuchos::Array<value_type>& gamma)
const 166 Teuchos::Array<value_type> nrm(n);
172 Teuchos::RCP<matrix_type> lv;
174 lv = Teuchos::rcp(&lanczos_vecs,
false);
179 lanczos_type::computeNormalized(n, vs, A, u0, *lv, alpha, beta, nrm);
181 lanczos_type::compute(n, vs, A, u0, *lv, alpha, beta, nrm);
192 return this->normalize;
195 template <
typename ordinal_type,
typename value_type>
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();
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);
218 new_pce.resize(this->p+1);
221 u[i] = (*pce)[i]*pce->basis()->norm_squared(i);
222 new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, fromStieltjesMat, u,
225 new_pce[i] /= this->norms[i];
228 template <
typename ordinal_type,
typename value_type>
234 limit_integration_order(basis.limit_integration_order),
236 pce_weights(basis.pce_weights),
237 pce_vals(basis.pce_vals),
239 lanczos_vecs(nqp, p+1),
vector_type pce_vals
Values of PCE at quadrature points.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
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.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
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.
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...
Abstract base class for quadrature methods.
lanczos_type::matrix_type matrix_type
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.
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
lanczos_type::vector_type vector_type
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 transformCoeffsFromLanczos(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
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.
vector_type u0
Initial Lanczos vector.
~LanczosPCEBasis()
Destructor.
virtual void setup()
Setup basis after computing recurrence coefficients.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > quad
Quadrature object.