44#include "Teuchos_BLAS.hpp"
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"),
56 sparse_tol(sparse_tol_),
58 d(basis->dimension()),
66 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
67 for (ordinal_type k=0; k<nqp; k++) {
69 basis->evaluateBases(points[k], values[k]);
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++) {
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;
88 for (ordinal_type i=0; i<sz; i++) {
93 for (ordinal_type
j=0;
j<i;
j++) {
97 for (ordinal_type k=0; k<=
j; k++)
98 t += gs_mat(
j,k)*inner_product(i,k);
102 for (ordinal_type k=0; k<=
j; k++)
103 gs_mat(i,k) -= t*gs_mat(
j,k);
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);
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++) {
123 for (ordinal_type
j=0;
j<=i;
j++)
124 t += gs_mat(i,
j)*values[k][
j];
125 basis_values[k][i] = t;
130template <
typename ordinal_type,
typename value_type>
136template <
typename ordinal_type,
typename value_type>
144template <
typename ordinal_type,
typename value_type>
152template <
typename ordinal_type,
typename value_type>
160template <
typename ordinal_type,
typename value_type>
161const Teuchos::Array<value_type>&
168template <
typename ordinal_type,
typename value_type>
176template <
typename ordinal_type,
typename value_type>
177Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
181 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
182 Teuchos::rcp(
new Sparse3Tensor<ordinal_type, value_type>);
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++) {
188 for (ordinal_type l=0; l<nqp; l++)
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);
197 Cijk->fillComplete();
202template <
typename ordinal_type,
typename value_type>
203Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
207 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
208 Teuchos::rcp(
new Sparse3Tensor<ordinal_type, value_type>);
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++) {
214 for (ordinal_type l=0; l<nqp; l++)
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);
223 Cijk->fillComplete();
228template <
typename ordinal_type,
typename value_type>
234 for (ordinal_type
j=0;
j<sz;
j++)
235 z += gs_mat(i,
j)*basis->evaluateZero(
j);
240template <
typename ordinal_type,
typename value_type>
243evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
244 Teuchos::Array<value_type>& basis_vals)
const
246 basis->evaluateBases(point, basis_vals_tmp);
247 for (ordinal_type i=0; i<sz; i++) {
249 for (ordinal_type
j=0;
j<sz;
j++)
250 t += gs_mat(i,
j)*basis_vals_tmp[
j];
255template <
typename ordinal_type,
typename value_type>
258print(std::ostream& os)
const
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] <<
" ";
267 os <<
"Underlying basis:\n";
271template <
typename ordinal_type,
typename value_type>
279template <
typename ordinal_type,
typename value_type>
284 Teuchos::BLAS<ordinal_type, value_type> blas;
285 for (ordinal_type i=0; i<sz; 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);
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)