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),
69 for (ordinal_type i=0; i<
pce_sz; i++) {
75 ordinal_type nqp = quad.
size();
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 =
82 for (ordinal_type i=0; i<nqp; i++) {
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) {
94 ordinal_type k = index(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) {
97 ordinal_type
j = index(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) {
101 ordinal_type i = index(i_it);
111 vector_type u0 = Teuchos::getCol(Teuchos::View, K, 0);
113 for (ordinal_type i=1; i<
pce_sz; i++)
115 for (ordinal_type k=1; k<
pce_sz; k++) {
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;
134 ordinal_type ws_size, info;
135 value_type ws_size_query;
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);
142 ws_size =
static_cast<ordinal_type
>(ws_size_query);
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);
154 ws_size =
static_cast<ordinal_type
>(ws_size_query);
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);
162 for (ordinal_type
j=0;
j<
p+1;
j++)
163 for (ordinal_type i=0; i<
pce_sz; i++)
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);
177 for (ordinal_type i=0; i<
pce_sz-1; i++) {
188 for (ordinal_type i=0; i<
pce_sz; i++)
189 u[i] = normalized_pce[i];
192 for (ordinal_type i=0; i<=
p; i++)
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");
215 ordinal_type num_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) {
229 ordinal_type old_size = quad_weights.size();
230 quad_weights.resize(num_points);
231 quad_points.resize(num_points);
232 quad_values.resize(num_points);
233 for (ordinal_type i=old_size; i<num_points; i++) {
234 quad_weights[i] = value_type(0);
235 quad_points[i] = quad_points[0];
236 quad_values[i].resize(this->p+1);
237 evaluateBases(quad_points[i], quad_values[i]);
265 Teuchos::BLAS<ordinal_type, value_type> blas;
266 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
267 value_type(1.0), basis_vecs.values(), pce_sz,
268 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
271 for (ordinal_type i=0; i<pce_sz; i++)
272 out[i] /= pce_norms[i];