Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GramSchmidtBasisImp.hpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Teuchos_BLAS.hpp"
45
46template <typename ordinal_type, typename value_type>
49 const Teuchos::RCP<const OrthogPolyBasis<ordinal_type,value_type> >& basis_,
50 const Teuchos::Array< Teuchos::Array<value_type> >& points,
51 const Teuchos::Array<value_type>& weights_,
52 const value_type& sparse_tol_) :
53 name("Gram Schmidt Basis"),
54 basis(basis_),
55 weights(weights_),
56 sparse_tol(sparse_tol_),
57 p(basis->order()),
58 d(basis->dimension()),
59 sz(basis->size()),
60 norms(sz),
61 gs_mat(sz,sz),
62 basis_vals_tmp(sz)
63{
64 // Get quadrature data
65 ordinal_type nqp = weights.size();
66 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
67 for (ordinal_type k=0; k<nqp; k++) {
68 values[k].resize(sz);
69 basis->evaluateBases(points[k], values[k]);
70 }
71
72 // Compute all inner products
73 Teuchos::SerialDenseMatrix<ordinal_type, value_type> inner_product(sz,sz);
74 inner_product.putScalar(0.0);
75 for (ordinal_type i=0; i<sz; i++) {
76 for (ordinal_type j=0; j<=i; j++) {
77 value_type t = 0.0;
78 for (ordinal_type k=0; k<nqp; k++)
79 t += weights[k]*values[k][i]*values[k][j];
80 inner_product(i,j) = t;
81 }
82 }
83
84 // Classical Gram-Schmidt algorithm:
85 // u_i = v_i - \sum_{j<i} (v_i,u_j)/(u_j,u_j) u_j
86 // u_j = \sum_{k<=i} a_{jk} v_k
87 // => u_i = v_i - \sum_{j<i}\sum_{k<=j} (v_i,u_j)/(u_j,u_j)*a_{jk}*v_k
88 for (ordinal_type i=0; i<sz; i++) {
89
90 // a_{ii} = 1.0
91 gs_mat(i,i) = 1.0;
92
93 for (ordinal_type j=0; j<i; j++) {
94
95 // compute t = (v_i,u_j)/(u_j,u_j)
96 value_type t = 0.0;
97 for (ordinal_type k=0; k<=j; k++)
98 t += gs_mat(j,k)*inner_product(i,k);
99 t /= norms[j];
100
101 // substract contribution to a_{ik}: t*a_{jk}
102 for (ordinal_type k=0; k<=j; k++)
103 gs_mat(i,k) -= t*gs_mat(j,k);
104 }
105
106 // compute (u_i,u_i) = \sum_{j,k<=i} a_{ij}*a_{ik}*(v_j,v_k)
107 value_type nrm = 0.0;
108 for (ordinal_type j=0; j<=i; j++) {
109 for (ordinal_type k=0; k<=j; k++)
110 nrm += gs_mat(i,j)*gs_mat(i,k)*inner_product(j,k);
111 for (ordinal_type k=j+1; k<=i; k++)
112 nrm += gs_mat(i,j)*gs_mat(i,k)*inner_product(k,j);
113 }
114 norms[i] = nrm;
115
116 }
117
118 basis_values.resize(nqp);
119 for (ordinal_type k=0; k<nqp; k++) {
120 basis_values[k].resize(sz);
121 for (ordinal_type i=0; i<sz; i++) {
122 value_type t = 0.0;
123 for (ordinal_type j=0; j<=i; j++)
124 t += gs_mat(i,j)*values[k][j];
125 basis_values[k][i] = t;
126 }
127 }
128}
129
130template <typename ordinal_type, typename value_type>
133{
134}
135
136template <typename ordinal_type, typename value_type>
139order() const
140{
141 return p;
142}
143
144template <typename ordinal_type, typename value_type>
147dimension() const
148{
149 return d;
150}
151
152template <typename ordinal_type, typename value_type>
155size() const
156{
157 return sz;
158}
159
160template <typename ordinal_type, typename value_type>
161const Teuchos::Array<value_type>&
163norm_squared() const
164{
165 return norms;
166}
167
168template <typename ordinal_type, typename value_type>
169const value_type&
171norm_squared(ordinal_type i) const
172{
173 return norms[i];
174}
175
176template <typename ordinal_type, typename value_type>
177Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
180{
181 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
182 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>);
183 ordinal_type nqp = weights.size();
184 for (ordinal_type j=0; j<sz; j++) {
185 for (ordinal_type i=0; i<sz; i++) {
186 for (ordinal_type k=0; k<sz; k++) {
187 value_type t = 0.0;
188 for (ordinal_type l=0; l<nqp; l++)
189 t +=
190 weights[l]*basis_values[l][i]*basis_values[l][j]*basis_values[l][k];
191 if (std::abs(t) > sparse_tol)
192 Cijk->add_term(i,j,k,t);
193 }
194 }
195 }
196
197 Cijk->fillComplete();
198
199 return Cijk;
200}
201
202template <typename ordinal_type, typename value_type>
203Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
206{
207 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
208 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>);
209 ordinal_type nqp = weights.size();
210 for (ordinal_type j=0; j<sz; j++) {
211 for (ordinal_type i=0; i<sz; i++) {
212 for (ordinal_type k=0; k<d+1; k++) {
213 value_type t = 0.0;
214 for (ordinal_type l=0; l<nqp; l++)
215 t +=
216 weights[l]*basis_values[l][i]*basis_values[l][j]*basis_values[l][k];
217 if (std::abs(t) > sparse_tol)
218 Cijk->add_term(i,j,k,t);
219 }
220 }
221 }
222
223 Cijk->fillComplete();
224
225 return Cijk;
226}
227
228template <typename ordinal_type, typename value_type>
231evaluateZero(ordinal_type i) const
232{
233 value_type z = 0.0;
234 for (ordinal_type j=0; j<sz; j++)
235 z += gs_mat(i,j)*basis->evaluateZero(j);
236
237 return z;
238}
239
240template <typename ordinal_type, typename value_type>
241void
243evaluateBases(const Teuchos::ArrayView<const value_type>& point,
244 Teuchos::Array<value_type>& basis_vals) const
245{
246 basis->evaluateBases(point, basis_vals_tmp);
247 for (ordinal_type i=0; i<sz; i++) {
248 value_type t = 0.0;
249 for (ordinal_type j=0; j<sz; j++)
250 t += gs_mat(i,j)*basis_vals_tmp[j];
251 basis_vals[i] = t;
252 }
253}
254
255template <typename ordinal_type, typename value_type>
256void
258print(std::ostream& os) const
259{
260 os << "Gram-Schmidt basis of order " << p << ", dimension " << d
261 << ", and size " << sz << ". Matrix coefficients:\n";
262 os << printMat(gs_mat) << std::endl;
263 os << "Basis vector norms (squared):\n\t";
264 for (ordinal_type i=0; i<sz; i++)
265 os << norms[i] << " ";
266 os << "\n";
267 os << "Underlying basis:\n";
268 os << *basis;
269}
270
271template <typename ordinal_type, typename value_type>
272const std::string&
274getName() const
275{
276 return name;
277}
278
279template <typename ordinal_type, typename value_type>
280void
282transformCoeffs(const value_type *in, value_type *out) const
283{
284 Teuchos::BLAS<ordinal_type, value_type> blas;
285 for (ordinal_type i=0; i<sz; i++)
286 out[i] = in[i];
287 blas.TRSM(Teuchos::LEFT_SIDE, Teuchos::LOWER_TRI, Teuchos::TRANS,
288 Teuchos::UNIT_DIAG, sz, 1, 1.0, gs_mat.values(), sz, out, sz);
289}
Transforms a non-orthogonal multivariate basis to an orthogonal one using the Gram-Schmit procedure.
virtual value_type evaluateZero(ordinal_type i) const =0
Evaluate basis polynomial i at zero.
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const =0
Compute triple product tensor.
virtual ordinal_type order() const =0
Return order of basis.
virtual ordinal_type size() const =0
Return total size of basis.
virtual ordinal_type dimension() const =0
Return dimension of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const =0
Compute linear triple product tensor where k = 0,1.
virtual const std::string & getName() const =0
Return string name of basis.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const =0
Evaluate basis polynomials at given point point.
virtual void print(std::ostream &os) const =0
Print basis to stream os.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)