Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_HouseTriDiagPCEBasisImp.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,
52 bool limit_integration_order_) :
53 RecurrenceBasis<ordinal_type, value_type>("Householder Tridiagonalization PCE", p, true),
54 limit_integration_order(limit_integration_order_),
55 pce_sz(pce.basis()->size()),
56 a(pce_sz+1),
57 b(pce_sz),
58 basis_vecs(pce_sz, p+1),
59 new_pce(p+1)
60{
61 pce_norms = pce.basis()->norm_squared();
62
63 // Compute matrix -- rescale basis to unit norm
65 A(0,0) = 1.0;
66 for (ordinal_type i=0; i<pce_sz; i++) {
67 A(0,i+1) = 1.0/pce_sz;
68 A(i+1,0) = 1.0/pce_sz;
69 }
71 for (typename Cijk_type::k_iterator k_it = Cijk.k_begin();
72 k_it != Cijk.k_end(); ++k_it) {
73 ordinal_type k = index(k_it);
74 for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it);
75 j_it != Cijk.j_end(k_it); ++j_it) {
76 ordinal_type j = index(j_it);
77 value_type val = 0;
78 for (typename Cijk_type::kji_iterator i_it = Cijk.i_begin(j_it);
79 i_it != Cijk.i_end(j_it); ++i_it) {
80 ordinal_type i = index(i_it);
81 value_type c = value(i_it) / std::sqrt(pce_norms[i]*pce_norms[j]*pce_norms[k]);
82 val += std::sqrt(pce_norms[i])*pce[i]*c;
83 }
84 A(k+1,j+1) = val;
85 }
86 }
87
88 // Call LAPACK Householder tridiagonalization algorithm
89 // Householder vectors are stored in lower-triangular part
90 ordinal_type ws_size, info;
91 value_type ws_size_query;
92 Teuchos::Array<value_type> tau(pce_sz-1);
93 lapack.SYTRD('L', pce_sz+1, A.values(), A.stride(), &a[0], &b[0], &tau[0],
94 &ws_size_query, -1, &info);
95 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
96 "SYTRD returned value " << info);
97 ws_size = static_cast<ordinal_type>(ws_size_query);
98 Teuchos::Array<value_type> work(ws_size);
99 lapack.SYTRD('L', pce_sz+1, A.values(), A.stride(), &a[0], &b[0], &tau[0],
100 &work[0], ws_size, &info);
101 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
102 "SYTRD returned value " << info);
103
104 // Set sub-diagonal terms to zero
105 for (ordinal_type j=0; j<pce_sz; j++)
106 A(j+1,j) = 0.0;
107
108 // Now compute orthogonal matrix
109 lapack.ORGQR(pce_sz+1, pce_sz+1, pce_sz-1, A.values(), A.stride(), &tau[0],
110 &ws_size_query, -1, &info);
111 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
112 "ORGQR returned value " << info);
113 ws_size = static_cast<ordinal_type>(ws_size_query);
114 work.resize(ws_size);
115 lapack.ORGQR(pce_sz+1, pce_sz+1, pce_sz-1, A.values(), A.stride(), &tau[0],
116 &work[0], ws_size, &info);
117 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
118 "ORGQR returned value " << info);
119
120 // Get basis vectors
121 for (ordinal_type j=0; j<p+1; j++)
122 for (ordinal_type i=0; i<pce_sz; i++)
123 basis_vecs(i,j) = A(i+1,j+1);
124
125 // Setup of rest of recurrence basis
126 this->setup();
127
128 // Project original PCE into the new basis
130 for (ordinal_type i=0; i<pce_sz; i++)
131 u[i] = pce[i]*std::sqrt(pce_norms[i]);
132 new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, basis_vecs, u,
133 0.0);
134 for (ordinal_type i=0; i<=p; i++)
135 new_pce[i] /= this->norms[i];
136
137 std::cout << new_pce << std::endl;
138
139 // Test orthogonality of basis vectors
140 matrix_type prod(p+1,p+1);
141 prod.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, basis_vecs, basis_vecs,
142 0.0);
143 for (ordinal_type j=0; j<p+1; j++)
144 prod(j,j) -= 1.0;
145 std::cout << "orthogonalization error = " << prod.normInf() << std::endl;
146 //std::cout << prod << std::endl;
147}
148
149template <typename ordinal_type, typename value_type>
152{
153}
154
155template <typename ordinal_type, typename value_type>
156void
158getQuadPoints(ordinal_type quad_order,
159 Teuchos::Array<value_type>& quad_points,
160 Teuchos::Array<value_type>& quad_weights,
161 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
162{
163#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
164 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::LanczosPCEBasis -- compute Gauss points");
165#endif
166
167 // Call base class
168 ordinal_type num_points =
169 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
170
171 // We can't always reliably generate quadrature points of order > 2*p
172 // when using sparse grids for the underlying quadrature
173 if (limit_integration_order && quad_order > 2*this->p)
174 quad_order = 2*this->p;
176 quad_points,
177 quad_weights,
178 quad_values);
179
180 // Fill in the rest of the points with zero weight
181 if (quad_weights.size() < num_points) {
182 ordinal_type old_size = quad_weights.size();
183 quad_weights.resize(num_points);
184 quad_points.resize(num_points);
185 quad_values.resize(num_points);
186 for (ordinal_type i=old_size; i<num_points; i++) {
187 quad_weights[i] = value_type(0);
188 quad_points[i] = quad_points[0];
189 quad_values[i].resize(this->p+1);
190 evaluateBases(quad_points[i], quad_values[i]);
191 }
192 }
193}
194
195template <typename ordinal_type, typename value_type>
196Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
198cloneWithOrder(ordinal_type p) const
199{
200 return
202 p,*this));
203}
204
205template <typename ordinal_type, typename value_type>
206value_type
208getNewCoeffs(ordinal_type i) const
209{
210 return new_pce[i];
211}
212
213template <typename ordinal_type, typename value_type>
214void
216transformCoeffsFromHouse(const value_type *in, value_type *out) const
217{
218 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
219 value_type(1.0), basis_vecs.values(), pce_sz,
220 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
221 for (ordinal_type i=0; i<pce_sz; i++)
222 out[i] /= pce_norms[i];
223}
224
225template <typename ordinal_type, typename value_type>
226bool
228computeRecurrenceCoefficients(ordinal_type n,
229 Teuchos::Array<value_type>& alpha,
230 Teuchos::Array<value_type>& beta,
231 Teuchos::Array<value_type>& delta,
232 Teuchos::Array<value_type>& gamma) const
233{
234 // Get recurrence coefficients from the full set we already computed
235 for (ordinal_type i=0; i<n; i++) {
236 alpha[i] = a[i];
237 beta[i] = b[i];
238 delta[i] = value_type(1.0);
239 gamma[i] = b[i];
240
241 std::cout << "i = " << i << " alpha = " << alpha[i] << " beta = " << beta[i]
242 << " gamma = " << gamma[i] << std::endl;
243 }
244
245 return true;
246}
247
248template <typename ordinal_type, typename value_type>
250HouseTriDiagPCEBasis(ordinal_type p, const HouseTriDiagPCEBasis& basis) :
251 RecurrenceBasis<ordinal_type, value_type>("Householder Tridiagonalization PCE", p, false),
252 limit_integration_order(basis.limit_integration_order),
253 pce_sz(basis.pce_sz),
254 a(basis.a),
255 b(basis.b),
256 basis_vecs(basis.basis_vecs),
257 new_pce(basis.new_pce)
258{
259 this->setup();
260}
expr val()
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
Teuchos::Array< value_type > b
Stores full set of beta coefficients.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new basis.
Teuchos::Array< value_type > a
Stores full set of alpha coefficients.
vector_type new_pce
Projection of pce in new basis.
ordinal_type pce_sz
Size of PC expansion.
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.
Teuchos::LAPACK< ordinal_type, value_type > lapack
LAPACK routines.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > matrix_type
HouseTriDiagPCEBasis(ordinal_type p, const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, bool limit_integration_order=false)
Constructor.
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.
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.
void transformCoeffsFromHouse(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
Teuchos::Array< value_type > norms
Norms.
ordinal_type p
Order of basis.
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.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
k_iterator k_end() const
Iterator pointing to last k entry.