Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
46 template <typename ordinal_type, typename value_type>
49  const Teuchos::RCP<const OrthogPolyBasis<ordinal_type,value_type> >& basis_,
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();
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
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 
130 template <typename ordinal_type, typename value_type>
133 {
134 }
135 
136 template <typename ordinal_type, typename value_type>
139 order() const
140 {
141  return p;
142 }
143 
144 template <typename ordinal_type, typename value_type>
147 dimension() const
148 {
149  return d;
150 }
151 
152 template <typename ordinal_type, typename value_type>
155 size() const
156 {
157  return sz;
158 }
159 
160 template <typename ordinal_type, typename value_type>
163 norm_squared() const
164 {
165  return norms;
166 }
167 
168 template <typename ordinal_type, typename value_type>
169 const value_type&
171 norm_squared(ordinal_type i) const
172 {
173  return norms[i];
174 }
175 
176 template <typename ordinal_type, typename value_type>
180 {
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 
202 template <typename ordinal_type, typename value_type>
206 {
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 
228 template <typename ordinal_type, typename value_type>
231 evaluateZero(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 
240 template <typename ordinal_type, typename value_type>
241 void
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 
255 template <typename ordinal_type, typename value_type>
256 void
258 print(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 
271 template <typename ordinal_type, typename value_type>
272 const std::string&
274 getName() const
275 {
276  return name;
277 }
278 
279 template <typename ordinal_type, typename value_type>
280 void
282 transformCoeffs(const value_type *in, value_type *out) const
283 {
285  for (ordinal_type i=0; i<sz; i++)
286  out[i] = in[i];
288  Teuchos::UNIT_DIAG, sz, 1, 1.0, gs_mat.values(), sz, out, sz);
289 }
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
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 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 ordinal_type dimension() const =0
Return dimension of basis.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual value_type evaluateZero(ordinal_type i) const =0
Evaluate basis polynomial i at zero.
Transforms a non-orthogonal multivariate basis to an orthogonal one using the Gram-Schmit procedure...
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual void print(std::ostream &os) const =0
Print basis to stream os.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const =0
Compute linear triple product tensor where k = 0,1.
virtual ordinal_type size() const =0
Return total size of basis.
virtual const std::string & getName() const =0
Return string name of basis.