42#include "Teuchos_UnitTestHarness.hpp"
43#include "Teuchos_TestingHelpers.hpp"
44#include "Teuchos_UnitTestRepository.hpp"
45#include "Teuchos_GlobalMPISession.hpp"
55 template <
typename OrdinalType,
typename ValueType>
59 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > >
bases;
60 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> >
basis;
61 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> >
quad;
62 Teuchos::RCP<Cijk_type>
Cijk;
70 for (OrdinalType i=0; i<
d; i++)
83 Cijk =
basis->computeTripleProductTensor();
91 const Teuchos::Array<double>& weights =
setup.quad->getQuadWeights();
92 const Teuchos::Array< Teuchos::Array<double> > & values =
93 setup.quad->getBasisAtQuadPoints();
97 k_it!=
setup.Cijk->k_end(); ++k_it) {
98 int k = Stokhos::index(k_it);
100 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
101 int j = Stokhos::index(j_it);
103 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
104 int i = Stokhos::index(i_it);
105 double c = Stokhos::value(i_it);
108 int nqp = weights.size();
109 for (
int qp=0; qp<nqp; qp++)
110 c2 += weights[qp]*values[qp][i]*values[qp][
j]*values[qp][k];
113 double err = std::abs(c-c2);
117 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )"
118 <<
" = " <<
"rel_err( " << c <<
", " << c2 <<
" ) = " << err
119 <<
" <= " << tol <<
" : ";
120 if (s) out <<
"Passed.";
125 success = success && s;
132 Teuchos::RCP< Stokhos::Sparse3Tensor<int, double> > Cijk_quad =
136 Teuchos::Array< Teuchos::RCP<Stokhos::Dense3Tensor<int,double> > > Cijk_1d(
setup.d);
137 for (
int i=0; i<
setup.d; i++)
138 Cijk_1d[i] =
setup.bases[i]->computeTripleProductTensor();
142 for (
int i=0; i<
setup.sz; i++) {
144 for (
int k=0; k<
setup.sz; k++) {
147 for (
int l=0; l<
setup.d; l++)
148 c *= (*Cijk_1d[l])(terms_i[l],terms_j[l],terms_k[l]);
149 if (std::abs(c/
setup.basis->norm_squared(i)) >
setup.sparse_tol)
150 Cijk_quad->add_term(i,
j,k,c);
154 Cijk_quad->fillComplete();
157 int nnz =
setup.Cijk->num_entries();
158 int nnz_quad = Cijk_quad->num_entries();
159 success = (nnz == nnz_quad);
162 <<
"Check: nnz(C) = " << nnz <<
" == nnz(C_quad) = " << nnz_quad
165 k_it!=Cijk_quad->k_end(); ++k_it) {
166 int k = Stokhos::index(k_it);
168 j_it != Cijk_quad->j_end(k_it); ++j_it) {
169 int j = Stokhos::index(j_it);
171 i_it != Cijk_quad->i_end(j_it); ++i_it) {
172 int i = Stokhos::index(i_it);
173 double c =
setup.Cijk->getValue(i,
j,k);
174 double c2 = Stokhos::value(i_it);
177 double err = std::abs(c-c2);
181 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )"
182 <<
" = " <<
"rel_err( " << c <<
", " << c2 <<
" ) = " << err
183 <<
" <= " << tol <<
" : ";
184 if (s) out <<
"Passed.";
189 success = success && s;
202 c =
setup.Cijk->getValue(0, 0, 0);
206 success = success && s;
208 c =
setup.Cijk->getValue(9, 25, 4);
212 success = success && s;
214 c =
setup.Cijk->getValue(8, 25, 4);
218 success = success && s;
224 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
225 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
A multidimensional index.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
TEUCHOS_UNIT_TEST(Stokhos_Sparse3Tensor, Values)
Stokhos::Sparse3Tensor< int, double > Cijk_type
UnitTestSetup< int, double > setup
bool compareValues(const ValueType &a1, const std::string &a1_name, const ValueType &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Teuchos::Array< Teuchos::RCP< const Stokhos::OneDOrthogPolyBasis< OrdinalType, ValueType > > > bases
Teuchos::RCP< Cijk_type > Cijk
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
Bi-directional iterator for traversing a sparse array.