Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_LanczosProjPCEBasisImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Teuchos_Assert.hpp"
43#include "Teuchos_BLAS.hpp"
44#include "Teuchos_TimeMonitor.hpp"
45
46template <typename ordinal_type, typename value_type>
49 ordinal_type p,
50 const Teuchos::RCP< const Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce_,
51 const Teuchos::RCP< const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk,
52 bool normalize,
53 bool limit_integration_order_) :
54 RecurrenceBasis<ordinal_type, value_type>("Lanczos-proj PCE", p, normalize),
55 pce(pce_),
56 limit_integration_order(limit_integration_order_),
57 pce_sz(pce->basis()->size()),
58 Cijk_matrix(pce_sz,pce_sz),
59 weights(Teuchos::Copy,
60 const_cast<value_type*>(pce->basis()->norm_squared().getRawPtr()),
61 pce_sz),
62 u0(pce_sz),
63 lanczos_vecs(pce_sz, p+1),
64 new_pce(p+1)
65{
66 u0[0] = value_type(1);
67
68 pce_norms = pce->basis()->norm_squared();
69 for (ordinal_type i=0; i<pce_sz; i++) {
70 pce_norms[i] = std::sqrt(pce_norms[i]);
71 weights[i] = value_type(1);
72 }
73
74 // Compute matrix -- For the matrix to be symmetric, the original basis
75 // must be normalized. However we don't want to require this, so we
76 // rescale the pce coefficients for a normalized basis
78 for (typename Cijk_type::k_iterator k_it = Cijk->k_begin();
79 k_it != Cijk->k_end(); ++k_it) {
80 ordinal_type k = index(k_it);
81 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
82 j_it != Cijk->j_end(k_it); ++j_it) {
83 ordinal_type j = index(j_it);
84 value_type val = 0;
85 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
86 i_it != Cijk->i_end(j_it); ++i_it) {
87 ordinal_type i = index(i_it);
88 value_type c = value(i_it);
89 val += (*pce)[i]*c / (pce_norms[j]*pce_norms[k]);
90 }
91 Cijk_matrix(k,j) = val;
92 }
93 }
94
95 // Setup of rest of recurrence basis
96 this->setup();
97
98
99}
100
101template <typename ordinal_type, typename value_type>
106
107template <typename ordinal_type, typename value_type>
108void
110getQuadPoints(ordinal_type quad_order,
111 Teuchos::Array<value_type>& quad_points,
112 Teuchos::Array<value_type>& quad_weights,
113 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
114{
115#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
116 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::LanczosPCEBasis -- compute Gauss points");
117#endif
118
119 // Call base class
120 ordinal_type num_points =
121 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
122
123 // We can't always reliably generate quadrature points of order > 2*p
124 // when using sparse grids for the underlying quadrature
125 if (limit_integration_order && quad_order > 2*this->p)
126 quad_order = 2*this->p;
128 quad_points,
129 quad_weights,
130 quad_values);
131
132 // Fill in the rest of the points with zero weight
133 if (quad_weights.size() < num_points) {
134 ordinal_type old_size = quad_weights.size();
135 quad_weights.resize(num_points);
136 quad_points.resize(num_points);
137 quad_values.resize(num_points);
138 for (ordinal_type i=old_size; i<num_points; i++) {
139 quad_weights[i] = value_type(0);
140 quad_points[i] = quad_points[0];
141 quad_values[i].resize(this->p+1);
142 this->evaluateBases(quad_points[i], quad_values[i]);
143 }
144 }
145}
146
147template <typename ordinal_type, typename value_type>
148Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
150cloneWithOrder(ordinal_type p) const
151{
152 return
154 p,*this));
155}
156
157template <typename ordinal_type, typename value_type>
158value_type
160getNewCoeffs(ordinal_type i) const
161{
162 return new_pce[i];
163}
164
165template <typename ordinal_type, typename value_type>
166void
168transformCoeffsFromLanczos(const value_type *in, value_type *out) const
169{
170 // Transform coefficients to normalized basis
171 Teuchos::BLAS<ordinal_type, value_type> blas;
172 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
173 value_type(1.0), lanczos_vecs.values(), pce_sz,
174 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
175
176 // Transform from normalized to original
177 for (ordinal_type i=0; i<pce_sz; i++)
178 out[i] /= pce_norms[i];
179}
180
181template <typename ordinal_type, typename value_type>
182bool
184computeRecurrenceCoefficients(ordinal_type n,
185 Teuchos::Array<value_type>& alpha,
186 Teuchos::Array<value_type>& beta,
187 Teuchos::Array<value_type>& delta,
188 Teuchos::Array<value_type>& gamma) const
189{
190 Teuchos::Array<value_type> nrm(n);
191 vectorspace_type vs(weights);
192 operator_type A(Cijk_matrix);
193
194 // Create space to store lanczos vectors -- use lanczos_vecs if
195 // we are requesting p+1 vectors
196 Teuchos::RCP<matrix_type> lv;
197 if (n == this->p+1)
198 lv = Teuchos::rcp(&lanczos_vecs, false);
199 else
200 lv = Teuchos::rcp(new matrix_type(pce_sz,n));
201
202 if (this->normalize)
203 lanczos_type::computeNormalized(n, vs, A, u0, *lv, alpha, beta, nrm);
204 else
205 lanczos_type::compute(n, vs, A, u0, *lv, alpha, beta, nrm);
206
207 for (ordinal_type i=0; i<n; i++) {
208 delta[i] = value_type(1.0);
209 }
210 if (this->normalize)
211 gamma = beta;
212 else
213 for (ordinal_type i=0; i<n; i++)
214 gamma[i] = value_type(1.0);
215
216 /*
217 matrix_type slv(pce_sz, n);
218 for (ordinal_type j=0; j<n; j++)
219 for (ordinal_type i=0; i<pce_sz; i++)
220 slv(i,j) = (*lv)(i,j) * weights[i];
221 matrix_type prod(n,n);
222 prod.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *lv, slv, 0.0);
223 for (ordinal_type j=0; j<n; j++) {
224 for (ordinal_type i=0; i<n; i++)
225 prod(i,j) /= std::sqrt(nrm[i]*nrm[j]);
226 prod(j,j) -= 1.0;
227 }
228 std::cout << "orthogonalization error = " << prod.normInf() << std::endl;
229 */
230
231 return this->normalize;
232}
233
234template <typename ordinal_type, typename value_type>
235void
237setup()
238{
240
241 // Project original PCE into the new basis
242 vector_type u(pce_sz);
243 for (ordinal_type i=0; i<pce_sz; i++)
244 u[i] = (*pce)[i]*pce_norms[i];
245 new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, lanczos_vecs, u,
246 0.0);
247 for (ordinal_type i=0; i<=this->p; i++)
248 new_pce[i] /= this->norms[i];
249}
250
251template <typename ordinal_type, typename value_type>
253LanczosProjPCEBasis(ordinal_type p, const LanczosProjPCEBasis& basis) :
254 RecurrenceBasis<ordinal_type, value_type>("Lanczos-proj PCE", p, false),
255 pce(basis.pce),
256 limit_integration_order(basis.limit_integration_order),
257 pce_sz(basis.pce_sz),
258 pce_norms(basis.pce_norms),
259 Cijk_matrix(basis.Cijk_matrix),
260 weights(basis.weights),
261 u0(basis.u0),
262 lanczos_vecs(pce_sz, p+1),
263 new_pce()
264{
265 this->setup();
266}
Copy
expr val()
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
virtual void setup()
Setup basis after computing recurrence coefficients.
WeightedVectorSpace< ordinal_type, value_type > vectorspace_type
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.
vector_type weights
Weighting vector used in inner-products.
matrix_type Cijk_matrix
Triple-product matrix used in generating lanczos vectors.
void transformCoeffsFromLanczos(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > pce
PCE Lanczos procedure is based on.
LanczosProjPCEBasis(ordinal_type p, const Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, bool normalize, bool limit_integration_order=false)
Constructor.
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.
DenseOperator< ordinal_type, value_type > operator_type
ordinal_type pce_sz
Size of PC expansion.
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.
vector_type u0
Initial Lanczos vector.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
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.