Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
46 template <typename ordinal_type, typename value_type>
49  ordinal_type p,
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 
101 template <typename ordinal_type, typename value_type>
104 {
105 }
106 
107 template <typename ordinal_type, typename value_type>
108 void
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 
147 template <typename ordinal_type, typename value_type>
151 {
152  return
154  p,*this));
155 }
156 
157 template <typename ordinal_type, typename value_type>
161 {
162  return new_pce[i];
163 }
164 
165 template <typename ordinal_type, typename value_type>
166 void
169 {
170  // Transform coefficients to normalized basis
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 
181 template <typename ordinal_type, typename value_type>
182 bool
188  Teuchos::Array<value_type>& gamma) const
189 {
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
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 
234 template <typename ordinal_type, typename value_type>
235 void
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 
251 template <typename ordinal_type, typename value_type>
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 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
matrix_type Cijk_matrix
Triple-product matrix used in generating lanczos vectors.
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.
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
virtual void setup()
Setup basis after computing recurrence coefficients.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
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...
RawPointerConversionTraits< Container >::Ptr_t getRawPtr(const Container &c)
Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > pce
PCE Lanczos procedure is based on.
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.
Teuchos::Array< value_type > pce_norms
Basis norms.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void transformCoeffsFromLanczos(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
vector_type u0
Initial Lanczos vector.
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
vector_type weights
Weighting vector used in inner-products.
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.
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.
void resize(size_type new_size, const value_type &x=value_type())
Stokhos::Sparse3Tensor< int, double > Cijk_type
virtual void setup()
Setup basis after computing recurrence coefficients.
expr val()
size_type size() const
ordinal_type pce_sz
Size of PC expansion.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...