42#include "Teuchos_Assert.hpp"
43#include "Teuchos_BLAS.hpp"
44#include "Teuchos_TimeMonitor.hpp"
46template <
typename ordinal_type,
typename value_type>
52 bool limit_integration_order_) :
54 limit_integration_order(limit_integration_order_),
55 pce_sz(pce.basis()->size()),
58 basis_vecs(pce_sz, p+1),
66 for (ordinal_type i=0; i<
pce_sz; i++) {
71 for (
typename Cijk_type::k_iterator k_it = Cijk.
k_begin();
72 k_it != Cijk.
k_end(); ++k_it) {
73 ordinal_type k = index(k_it);
74 for (
typename Cijk_type::kj_iterator j_it = Cijk.
j_begin(k_it);
75 j_it != Cijk.
j_end(k_it); ++j_it) {
76 ordinal_type
j = index(j_it);
78 for (
typename Cijk_type::kji_iterator i_it = Cijk.
i_begin(j_it);
79 i_it != Cijk.
i_end(j_it); ++i_it) {
80 ordinal_type i = index(i_it);
90 ordinal_type ws_size, info;
91 value_type ws_size_query;
92 Teuchos::Array<value_type> tau(
pce_sz-1);
93 lapack.SYTRD(
'L',
pce_sz+1, A.values(), A.stride(), &
a[0], &
b[0], &tau[0],
94 &ws_size_query, -1, &info);
95 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
96 "SYTRD returned value " << info);
97 ws_size =
static_cast<ordinal_type
>(ws_size_query);
98 Teuchos::Array<value_type> work(ws_size);
99 lapack.SYTRD(
'L',
pce_sz+1, A.values(), A.stride(), &
a[0], &
b[0], &tau[0],
100 &work[0], ws_size, &info);
101 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
102 "SYTRD returned value " << info);
110 &ws_size_query, -1, &info);
111 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
112 "ORGQR returned value " << info);
113 ws_size =
static_cast<ordinal_type
>(ws_size_query);
114 work.resize(ws_size);
116 &work[0], ws_size, &info);
117 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
118 "ORGQR returned value " << info);
121 for (ordinal_type
j=0;
j<
p+1;
j++)
122 for (ordinal_type i=0; i<
pce_sz; i++)
130 for (ordinal_type i=0; i<
pce_sz; i++)
134 for (ordinal_type i=0; i<=
p; i++)
137 std::cout <<
new_pce << std::endl;
143 for (ordinal_type
j=0;
j<
p+1;
j++)
145 std::cout <<
"orthogonalization error = " << prod.normInf() << std::endl;
149template <
typename ordinal_type,
typename value_type>
155template <
typename ordinal_type,
typename value_type>
159 Teuchos::Array<value_type>& quad_points,
160 Teuchos::Array<value_type>& quad_weights,
161 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const
163#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
164 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::LanczosPCEBasis -- compute Gauss points");
168 ordinal_type num_points =
169 static_cast<ordinal_type
>(std::ceil((quad_order+1)/2.0));
173 if (limit_integration_order && quad_order > 2*this->p)
174 quad_order = 2*this->p;
181 if (quad_weights.size() < num_points) {
182 ordinal_type old_size = quad_weights.size();
183 quad_weights.resize(num_points);
184 quad_points.resize(num_points);
185 quad_values.resize(num_points);
186 for (ordinal_type i=old_size; i<num_points; i++) {
187 quad_weights[i] = value_type(0);
188 quad_points[i] = quad_points[0];
189 quad_values[i].resize(this->p+1);
190 evaluateBases(quad_points[i], quad_values[i]);
195template <
typename ordinal_type,
typename value_type>
196Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
205template <
typename ordinal_type,
typename value_type>
213template <
typename ordinal_type,
typename value_type>
218 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
219 value_type(1.0), basis_vecs.values(), pce_sz,
220 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
221 for (ordinal_type i=0; i<pce_sz; i++)
222 out[i] /= pce_norms[i];
225template <
typename ordinal_type,
typename value_type>
229 Teuchos::Array<value_type>& alpha,
230 Teuchos::Array<value_type>& beta,
231 Teuchos::Array<value_type>& delta,
232 Teuchos::Array<value_type>& gamma)
const
235 for (ordinal_type i=0; i<n; i++) {
238 delta[i] = value_type(1.0);
241 std::cout <<
"i = " << i <<
" alpha = " << alpha[i] <<
" beta = " << beta[i]
242 <<
" gamma = " << gamma[i] << std::endl;
248template <
typename ordinal_type,
typename value_type>
251 RecurrenceBasis<ordinal_type, value_type>(
"Householder Tridiagonalization PCE", p, false),
252 limit_integration_order(basis.limit_integration_order),
253 pce_sz(basis.pce_sz),
256 basis_vecs(basis.basis_vecs),
257 new_pce(basis.new_pce)
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
Teuchos::Array< value_type > b
Stores full set of beta coefficients.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
Teuchos::Array< value_type > a
Stores full set of alpha coefficients.
~HouseTriDiagPCEBasis()
Destructor.
vector_type new_pce
Projection of pce in new basis.
ordinal_type pce_sz
Size of PC expansion.
Teuchos::Array< value_type > pce_norms
Basis norms.
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::LAPACK< ordinal_type, value_type > lapack
LAPACK routines.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > matrix_type
HouseTriDiagPCEBasis(ordinal_type p, const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, bool limit_integration_order=false)
Constructor.
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 Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
void transformCoeffsFromHouse(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
matrix_type basis_vecs
Basis vectors.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
Teuchos::Array< value_type > norms
Norms.
ordinal_type p
Order of basis.
virtual void setup()
Setup basis after computing 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.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
k_iterator k_end() const
Iterator pointing to last k entry.