48template <
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();
60 ordinal_type pce_sz = pce_basis->size();
63 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
64 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
65 if (prod_basis != Teuchos::null)
66 coord_bases = prod_basis->getCoordinateBases();
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);
75 for (ordinal_type i=0; i<pce.size(); i++) {
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)));
113 ordinal_type d = coordinate_bases.size();
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();
136 ordinal_type nqp = weights.size();
137 SDM Psi(pce_sz, nqp);
138 for (ordinal_type i=0; i<pce_sz; i++)
139 for (ordinal_type k=0; k<nqp; k++)
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++) {
148 ordinal_type jdx = 0;
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;
184 value_type rank_threshold =
params.get(
"Rank Threshold", 1.0e-12);
185 ordinal_type rank =
svd_threshold(rank_threshold,
A, sigma, U, Vt);
186 Ainv.shape(sz, pce_sz);
187 TEUCHOS_ASSERT(rank == Vt.numRows());
188 for (ordinal_type i=0; i<Vt.numRows(); i++)
189 for (ordinal_type
j=0;
j<Vt.numCols();
j++)
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) = [";
199 for (ordinal_type i=0; i<rank; i++)
200 std::cout << sigma[i] <<
" ";
201 std::cout <<
"]" << std::endl;
204 SDM SVt(rank, Vt.numCols());
205 for (ordinal_type i=0; i<Vt.numRows(); i++)
206 for (ordinal_type
j=0;
j<Vt.numCols();
j++)
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;
220 for (ordinal_type i=0; i<sz; i++)
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);
231 for (ordinal_type i=0; i<pce_sz; i++)
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;
240template <
typename ordinal_type,
typename value_type>
246template <
typename ordinal_type,
typename value_type>
251 return tensor_lanczos_basis->order();
254template <
typename ordinal_type,
typename value_type>
259 return tensor_lanczos_basis->dimension();
262template <
typename ordinal_type,
typename value_type>
267 return tensor_lanczos_basis->size();
270template <
typename ordinal_type,
typename value_type>
271const Teuchos::Array<value_type>&
275 return tensor_lanczos_basis->norm_squared();
278template <
typename ordinal_type,
typename value_type>
283 return tensor_lanczos_basis->norm_squared(i);
286template <
typename ordinal_type,
typename value_type>
287Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
292 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
293 tensor_lanczos_basis->computeTripleProductTensor();
298template <
typename ordinal_type,
typename value_type>
299Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
304 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
305 tensor_lanczos_basis->computeLinearTripleProductTensor();
310template <
typename ordinal_type,
typename value_type>
315 return tensor_lanczos_basis->evaluateZero(i);
318template <
typename ordinal_type,
typename value_type>
321evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
322 Teuchos::Array<value_type>& basis_vals)
const
324 return tensor_lanczos_basis->evaluateBases(point, basis_vals);
327template <
typename ordinal_type,
typename value_type>
330print(std::ostream& os)
const
332 tensor_lanczos_basis->print(os);
335template <
typename ordinal_type,
typename value_type>
343template <
typename ordinal_type,
typename value_type>
346term(ordinal_type i)
const
348 return tensor_lanczos_basis->term(i);
351template <
typename ordinal_type,
typename value_type>
356 return tensor_lanczos_basis->index(term);
359template <
typename ordinal_type,
typename value_type>
360Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
364 return tensor_lanczos_basis->getCoordinateBases();
367template <
typename ordinal_type,
typename value_type>
372 return tensor_lanczos_basis->getMaxOrders();
375template <
typename ordinal_type,
typename value_type>
379 ordinal_type ncol,
bool transpose)
const
381 ordinal_type pce_sz = A.numRows();
382 ordinal_type sz = A.numCols();
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);
399template <
typename ordinal_type,
typename value_type>
403 ordinal_type ncol,
bool transpose)
const
405 ordinal_type pce_sz = A.numRows();
406 ordinal_type sz = A.numCols();
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);
423template <
typename ordinal_type,
typename value_type>
424Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
431template <
typename ordinal_type,
typename value_type>
436 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
438 ordinal_type dim = basis->dimension();
439 value_type tol = 1.0e-15;
442 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(basis);
443 if (prod_basis == Teuchos::null)
449 Teuchos::Array<ordinal_type> dependent_dims;
450 tmp_pce.reset(basis);
451 for (ordinal_type i=0; i<dim; i++) {
452 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
454 for (ordinal_type
j=1;
j<=p;
j++)
455 tmp_pce.term(i,
j) = pce.
term(i,
j);
456 value_type nrm = tmp_pce.two_norm();
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)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
Teuchos::ParameterList params
Algorithm parameters.
SDM Phi
Values of transformed basis at quadrature points.
virtual const std::string & getName() const
Return string name of basis.
SDM A
Transition matrix from reduced basis to original.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
ordinal_type isInvariant(const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce) const
Teuchos::RCP< Stokhos::CompletePolynomialBasis< ordinal_type, value_type > > tensor_lanczos_basis
Product Lanczos basis.
virtual ~ProductLanczosPCEBasis()
Destructor.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
ordinal_type order() const
Return order of basis.
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.
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 const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
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.
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.
ordinal_type dimension() const
Return dimension of basis.
SDM Ainv
Projection matrix from original matrix to reduced.
virtual ordinal_type size() const
Return total size of basis.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
static Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > create(Teuchos::ParameterList &sgParams)
Generate quadrature object.
Abstract base class for quadrature methods.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
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)