Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MonomialGramSchmidtPCEBasisImp.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 "Stokhos_SDMUtils.hpp"
44
45template <typename ordinal_type, typename value_type>
48 ordinal_type max_p,
49 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
50 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
51 const Teuchos::ParameterList& params) :
52 GSReducedPCEBasisBase<ordinal_type,value_type>(max_p, pce, quad, params),
53 name("Monomial Gram Schmidt PCE Basis")
54{
55 this->setup(max_p, pce, quad);
56}
57
58template <typename ordinal_type, typename value_type>
61{
62}
63
64template <typename ordinal_type, typename value_type>
65const std::string&
67getName() const
68{
69 return name;
70}
71
72template <typename ordinal_type, typename value_type>
73ordinal_type
76 ordinal_type max_p,
77 value_type threshold,
78 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
79 const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& F,
80 const Teuchos::Array<value_type>& weights,
81 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >& terms_,
82 Teuchos::Array<ordinal_type>& num_terms_,
83 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& Qp_,
84 Teuchos::SerialDenseMatrix<ordinal_type,value_type>& Q_)
85{
86 // Compute basis terms -- 2-D array giving powers for each linear index
87 ordinal_type max_sz;
88 CPBUtils::compute_terms(max_p, this->d, max_sz, terms_, num_terms_);
89
90 // Compute B matrix -- monomials in F
91 // for i=0,...,nqp-1
92 // for j=0,...,sz-1
93 // B(i,j) = F(i,1)^terms_[j][1] * ... * F(i,d)^terms_[j][d]
94 // where sz is the total size of a basis up to order p and terms_[j]
95 // is an array of powers for each term in the total-order basis
96 ordinal_type nqp = weights.size();
97 SDM B(nqp, max_sz);
98 for (ordinal_type i=0; i<nqp; i++) {
99 for (ordinal_type j=0; j<max_sz; j++) {
100 B(i,j) = 1.0;
101 for (ordinal_type k=0; k<this->d; k++)
102 B(i,j) *= std::pow(F(i,k), terms_[j][k]);
103 }
104 }
105
106 // Rescale columns of B to have unit norm
107 for (ordinal_type j=0; j<max_sz; j++) {
108 value_type nrm = 0.0;
109 for (ordinal_type i=0; i<nqp; i++)
110 nrm += B(i,j)*B(i,j)*weights[i];
111 nrm = std::sqrt(nrm);
112 for (ordinal_type i=0; i<nqp; i++)
113 B(i,j) /= nrm;
114 }
115
116 // Compute our new basis -- each column of Q is the new basis evaluated
117 // at the original quadrature points. Constraint pivoting so first d+1
118 // columns and included in Q.
119 SDM R;
120 Teuchos::Array<ordinal_type> piv(max_sz);
121 for (int i=0; i<this->d+1; i++)
122 piv[i] = 1;
124 ordinal_type sz_ = SOF::createOrthogonalBasis(
125 this->orthogonalization_method, threshold, this->verbose, B, weights,
126 Q_, R, piv);
127
128 // Compute Qp = A^T*W*Q
129 SDM tmp(nqp, sz_);
130 Qp_.reshape(this->pce_sz, sz_);
131 for (ordinal_type i=0; i<nqp; i++)
132 for (ordinal_type j=0; j<sz_; j++)
133 tmp(i,j) = Q_(i,j)*weights[i];
134 ordinal_type ret =
135 Qp_.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, A, tmp, 0.0);
136 TEUCHOS_ASSERT(ret == 0);
137
138 // It isn't clear that Qp is orthogonal, but if you derive the projection
139 // matrix from the original space to the reduced, you end up with
140 // Q^T*W*A = Qp^T
141
142 return sz_;
143}
UnitTestSetup< Kokkos::Cuda > setup
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual const std::string & getName() const
Return string name of basis.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
MonomialGramSchmidtPCEBasis(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::ParameterList &params=Teuchos::ParameterList())
Constructor.
virtual ordinal_type buildReducedBasis(ordinal_type max_p, value_type threshold, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights, Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > &terms_, Teuchos::Array< ordinal_type > &num_terms_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Qp_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q_)
Build the reduced basis, parameterized by total order max_p.
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Encapsulate various orthogonalization (ie QR) methods.
Abstract base class for quadrature methods.