42 #include "Teuchos_Assert.hpp" 43 #include "Teuchos_BLAS.hpp" 44 #include "Teuchos_LAPACK.hpp" 45 #include "Teuchos_TimeMonitor.hpp" 46 #include "Teuchos_SerialDenseHelpers.hpp" 48 template <
typename ordinal_type,
typename value_type>
55 bool limit_integration_order_) :
57 limit_integration_order(limit_integration_order_),
58 pce_sz(pce.basis()->size()),
59 pce_norms(pce.basis()->norm_squared()),
62 basis_vecs(pce_sz, p+1),
76 Teuchos::Array<value_type> pce_vals(nqp);
78 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
80 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
83 pce_vals[i] = normalized_pce.
evaluate(quad_points[i], basis_values[i]);
92 for (
typename Cijk_type::k_iterator k_it = Cijk.
k_begin();
93 k_it != Cijk.
k_end(); ++k_it) {
95 for (
typename Cijk_type::kj_iterator j_it = Cijk.
j_begin(k_it);
96 j_it != Cijk.
j_end(k_it); ++j_it) {
99 for (
typename Cijk_type::kji_iterator i_it = Cijk.
i_begin(j_it);
100 i_it != Cijk.
i_end(j_it); ++i_it) {
111 vector_type u0 = Teuchos::getCol(Teuchos::View, K, 0);
116 vector_type u = Teuchos::getCol(Teuchos::View, K, k);
117 vector_type up = Teuchos::getCol(Teuchos::View, K, k-1);
118 u.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, up, 0.0);
131 std::cout << K << std::endl << std::endl;
136 Teuchos::Array<value_type> tau(
pce_sz);
138 lapack.GEQRF(
pce_sz,
pce_sz, K.values(), K.stride(), &tau[0],
139 &ws_size_query, -1, &info);
140 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
141 "GEQRF returned value " << info);
143 Teuchos::Array<value_type> work(ws_size);
144 lapack.GEQRF(
pce_sz,
pce_sz, K.values(), K.stride(), &tau[0],
145 &work[0], ws_size, &info);
146 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
147 "GEQRF returned value " << info);
151 &ws_size_query, -1, &info);
152 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
153 "ORGQR returned value " << info);
155 work.resize(ws_size);
157 &work[0], ws_size, &info);
158 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
159 "ORGQR returned value " << info);
169 prod.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, K, A, 0.0);
171 T.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, prod, K, 0.0);
189 u[i] = normalized_pce[i];
196 template <
typename ordinal_type,
typename value_type>
202 template <
typename ordinal_type,
typename value_type>
206 Teuchos::Array<value_type>& quad_points,
207 Teuchos::Array<value_type>& quad_weights,
208 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const 210 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 211 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::MonoProjPCEBasis -- compute Gauss points");
216 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
220 if (limit_integration_order && quad_order > 2*this->p)
221 quad_order = 2*this->p;
228 if (quad_weights.size() < num_points) {
230 quad_weights.resize(num_points);
231 quad_points.resize(num_points);
232 quad_values.resize(num_points);
235 quad_points[i] = quad_points[0];
236 quad_values[i].resize(this->p+1);
237 evaluateBases(quad_points[i], quad_values[i]);
242 template <
typename ordinal_type,
typename value_type>
243 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
251 template <
typename ordinal_type,
typename value_type>
259 template <
typename ordinal_type,
typename value_type>
265 Teuchos::BLAS<ordinal_type, value_type> blas;
266 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
272 out[i] /= pce_norms[i];
275 template <
typename ordinal_type,
typename value_type>
279 Teuchos::Array<value_type>& alpha,
280 Teuchos::Array<value_type>& beta,
281 Teuchos::Array<value_type>& delta,
282 Teuchos::Array<value_type>& gamma)
const 291 std::cout <<
"i = " << i <<
" alpha = " << alpha[i] <<
" beta = " << beta[i]
292 <<
" gamma = " << gamma[i] << std::endl;
298 template <
typename ordinal_type,
typename value_type>
302 limit_integration_order(basis.limit_integration_order),
303 pce_sz(basis.pce_sz),
304 pce_norms(basis.pce_norms),
307 basis_vecs(basis.basis_vecs),
308 new_pce(basis.new_pce)
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Teuchos::Array< value_type > norms
Norms.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
virtual ordinal_type size() const =0
Get number of quadrature points.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
matrix_type basis_vecs
Basis vectors.
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...
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 const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const =0
Get quadrature points.
void transformCoeffs(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
Abstract base class for quadrature methods.
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.
~MonoProjPCEBasis()
Destructor.
Teuchos::Array< value_type > pce_norms
Basis norms.
Stokhos::Sparse3Tensor< int, double > Cijk_type
virtual void setup()
Setup basis after computing recurrence coefficients.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > matrix_type
ordinal_type pce_sz
Size of PC expansion.
Teuchos::Array< value_type > b
Stores full set of beta coefficients.
Teuchos::Array< value_type > a
Stores full set of alpha coefficients.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
k_iterator k_end() const
Iterator pointing to last k entry.
ordinal_type p
Order of basis.
MonoProjPCEBasis(ordinal_type p, const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce, const Stokhos::Quadrature< ordinal_type, value_type > &quad, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, bool limit_integration_order=false)
Constructor.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
vector_type new_pce
Projection of pce in new basis.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
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.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
k_iterator k_begin() const
Iterator pointing to first k entry.