48 template <
typename ordinal_type,
typename value_type>
55 const Teuchos::ParameterList& params_) :
56 name(
"Product Lanczos PCE Basis"),
59 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis = pce[0].basis();
64 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
65 if (prod_basis != Teuchos::null)
69 bool project =
params.get(
"Project",
true);
70 bool normalize =
params.get(
"Normalize",
true);
71 bool limit_integration_order =
params.get(
"Limit Integration Order",
false);
72 bool use_stieltjes =
params.get(
"Use Old Stieltjes Method",
false);
73 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double > > > coordinate_bases;
74 Teuchos::Array<int> is_invariant(pce.size(),-2);
81 if (prod_basis != Teuchos::null)
83 if (is_invariant[i] >= 0) {
84 coordinate_bases.push_back(coord_bases[is_invariant[i]]);
89 else if (is_invariant[i] != -1) {
91 coordinate_bases.push_back(
94 p, Teuchos::rcp(&(pce[i]),
false), quad,
false,
95 normalize, project, Cijk)));
99 coordinate_bases.push_back(
102 p, Teuchos::rcp(&(pce[i]),
false), Cijk,
103 normalize, limit_integration_order)));
105 coordinate_bases.push_back(
108 p, Teuchos::rcp(&(pce[i]),
false), quad,
109 normalize, limit_integration_order)));
120 params.get(
"Cijk Drop Tolerance", 1.0e-15),
121 params.get(
"Use Old Cijk Algorithm",
false)));
124 Teuchos::ParameterList sg_params;
125 sg_params.sublist(
"Basis").set< Teuchos::RCP< const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > >(
"Stochastic Galerkin Basis",
tensor_lanczos_basis);
126 sg_params.sublist(
"Quadrature") =
params.sublist(
"Reduced Quadrature");
131 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
132 const Teuchos::Array< Teuchos::Array<value_type> >& points =
133 quad->getQuadPoints();
134 const Teuchos::Array< Teuchos::Array<value_type> >& basis_vals =
135 quad->getBasisAtQuadPoints();
137 SDM Psi(pce_sz, nqp);
140 Psi(i,k) = basis_vals[k][i]*weights[k]/pce_basis->norm_squared(i);
144 Teuchos::Array<value_type> red_basis_vals(sz);
145 Teuchos::Array<value_type> pce_vals(d);
147 for (
int k=0; k<nqp; k++) {
149 for (
int j=0;
j<pce.size();
j++) {
152 if (is_invariant[
j] != -1) {
155 if (is_invariant[
j] >= 0)
156 pce_vals[jdx] = points[k][is_invariant[
j]];
158 pce_vals[jdx] = pce[
j].evaluate(points[k], basis_vals[k]);
165 for (
int i=0; i<sz; i++)
166 Phi(i,k) = red_basis_vals[i];
169 bool verbose =
params.get(
"Verbose",
false);
174 A.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, Psi,
Phi, 0.0);
175 TEUCHOS_ASSERT(ret == 0);
182 Teuchos::Array<value_type> sigma;
186 Ainv.shape(sz, pce_sz);
187 TEUCHOS_ASSERT(rank == Vt.numRows());
190 Vt(i,
j) = Vt(i,
j) / sigma[i];
191 ret =
Ainv.multiply(Teuchos::TRANS, Teuchos::TRANS, 1.0, Vt, U, 0.0);
192 TEUCHOS_ASSERT(ret == 0);
196 std::cout <<
"rank = " << rank << std::endl;
198 std::cout <<
"diag(S) = [";
200 std::cout << sigma[i] <<
" ";
201 std::cout <<
"]" << std::endl;
204 SDM SVt(rank, Vt.numCols());
207 SVt(i,
j) = Vt(i,
j) * sigma[i] * sigma[i];
209 SDM err_A(pce_sz,sz);
211 ret = err_A.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -1.0, U, SVt,
213 TEUCHOS_ASSERT(ret == 0);
214 std::cout <<
"||A - U*S*V^T||_infty = " << err_A.normInf() << std::endl;
222 ret = err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0,
Ainv,
A,
224 TEUCHOS_ASSERT(ret == 0);
225 std::cout <<
"||Ainv*A - I||_infty = " << err.normInf() << std::endl;
229 SDM err2(pce_sz,pce_sz);
233 ret = err2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0,
A,
Ainv, -1.0);
234 TEUCHOS_ASSERT(ret == 0);
235 std::cout <<
"||A*Ainv - I||_infty = " << err2.normInf() << std::endl;
240 template <
typename ordinal_type,
typename value_type>
246 template <
typename ordinal_type,
typename value_type>
251 return tensor_lanczos_basis->order();
254 template <
typename ordinal_type,
typename value_type>
259 return tensor_lanczos_basis->dimension();
262 template <
typename ordinal_type,
typename value_type>
267 return tensor_lanczos_basis->size();
270 template <
typename ordinal_type,
typename value_type>
271 const Teuchos::Array<value_type>&
275 return tensor_lanczos_basis->norm_squared();
278 template <
typename ordinal_type,
typename value_type>
283 return tensor_lanczos_basis->norm_squared(i);
286 template <
typename ordinal_type,
typename value_type>
287 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
292 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
293 tensor_lanczos_basis->computeTripleProductTensor();
298 template <
typename ordinal_type,
typename value_type>
299 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
304 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
305 tensor_lanczos_basis->computeLinearTripleProductTensor();
310 template <
typename ordinal_type,
typename value_type>
315 return tensor_lanczos_basis->evaluateZero(i);
318 template <
typename ordinal_type,
typename value_type>
322 Teuchos::Array<value_type>& basis_vals)
const 324 return tensor_lanczos_basis->evaluateBases(point, basis_vals);
327 template <
typename ordinal_type,
typename value_type>
332 tensor_lanczos_basis->print(os);
335 template <
typename ordinal_type,
typename value_type>
343 template <
typename ordinal_type,
typename value_type>
348 return tensor_lanczos_basis->term(i);
351 template <
typename ordinal_type,
typename value_type>
356 return tensor_lanczos_basis->index(term);
359 template <
typename ordinal_type,
typename value_type>
360 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
364 return tensor_lanczos_basis->getCoordinateBases();
367 template <
typename ordinal_type,
typename value_type>
372 return tensor_lanczos_basis->getMaxOrders();
375 template <
typename ordinal_type,
typename value_type>
384 SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
385 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
387 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, A, 0.0);
388 TEUCHOS_ASSERT(ret == 0);
391 SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
392 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
394 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, zbar, 0.0);
395 TEUCHOS_ASSERT(ret == 0);
399 template <
typename ordinal_type,
typename value_type>
408 SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
409 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
411 zbar.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, z, Ainv, 0.0);
412 TEUCHOS_ASSERT(ret == 0);
415 SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
416 SDM zbar(Teuchos::View, out, sz, sz, ncol);
418 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Ainv, z, 0.0);
419 TEUCHOS_ASSERT(ret == 0);
423 template <
typename ordinal_type,
typename value_type>
424 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
431 template <
typename ordinal_type,
typename value_type>
436 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
443 if (prod_basis == Teuchos::null)
449 Teuchos::Array<ordinal_type> dependent_dims;
450 tmp_pce.reset(basis);
452 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
455 tmp_pce.term(i,
j) = pce.
term(i,
j);
457 if (nrm > tol) dependent_dims.push_back(i);
462 if (dependent_dims.size() == 1)
463 return dependent_dims[0];
466 else if (dependent_dims.size() == 0)
virtual void transformToOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients to original basis from this basis.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
virtual void print(std::ostream &os) const
Print basis to stream os.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
SDM Phi
Values of transformed basis at quadrature points.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
virtual ordinal_type size() const
Return total size of basis.
SDM A
Transition matrix from reduced basis to original.
SDM Ainv
Projection matrix from original matrix to reduced.
virtual const std::string & getName() const
Return string name of basis.
ordinal_type order() const
Return order of basis.
virtual ~ProductLanczosPCEBasis()
Destructor.
ordinal_type svd_threshold(const scalar_type &rank_threshold, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, Teuchos::Array< scalar_type > &s, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &U, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Vt)
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
ProductLanczosPCEBasis(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::ParameterList ¶ms=Teuchos::ParameterList())
Constructor.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Abstract base class for quadrature methods.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
ordinal_type dimension() const
Return dimension of basis.
virtual void transformFromOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients from original basis to this basis.
Teuchos::ParameterList params
Algorithm parameters.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Teuchos::RCP< Stokhos::CompletePolynomialBasis< ordinal_type, value_type > > tensor_lanczos_basis
Product Lanczos basis.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
ordinal_type isInvariant(const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce) const
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
static Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > create(Teuchos::ParameterList &sgParams)
Generate quadrature object.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.