42#ifndef STOKHOS_UNIT_TEST_HELPERS_HPP
43#define STOKHOS_UNIT_TEST_HELPERS_HPP
48#include "Teuchos_SerialDenseMatrix.hpp"
49#include "Teuchos_FancyOStream.hpp"
53 template<
class OrdinalType,
class ValueType>
55 const std::string& a1_name,
57 const std::string& a2_name,
58 const ValueType& rel_tol,
const ValueType& abs_tol,
59 Teuchos::FancyOStream& out)
63 out <<
"Comparing " << a1_name <<
" == " << a2_name <<
" ... ";
65 const OrdinalType n = a1.
size();
69 out <<
"\nError, "<<a1_name<<
".size() = "<<a1.
size()<<
" == "
70 << a2_name<<
".size() = "<<a2.
size()<<
" : failed!\n";
75 for( OrdinalType i = 0; i < n; ++i ) {
76 ValueType nrm = std::sqrt(a1.
basis()->norm_squared(i));
77 ValueType err = std::abs(a1[i] - a2[i]) / nrm;
79 abs_tol + rel_tol*std::max(std::abs(a1[i]),std::abs(a2[i]))/nrm;
82 <<
"\nError, relErr("<<a1_name<<
"["<<i<<
"],"
83 <<a2_name<<
"["<<i<<
"]) = relErr("<<a1[i]<<
","<<a2[i]<<
") = "
84 <<err<<
" <= tol = "<<tol<<
": failed!\n";
93 << a1_name <<
" = " << a1 << std::endl
94 << a2_name <<
" = " << a2 << std::endl;
100 template<
class ValueType>
102 const std::string& a1_name,
104 const std::string& a2_name,
105 const ValueType& rel_tol,
const ValueType& abs_tol,
106 Teuchos::FancyOStream& out)
110 ValueType err = std::abs(a1 - a2);
111 ValueType tol = abs_tol + rel_tol*std::max(std::abs(a1),std::abs(a2));
113 out <<
"\nError, relErr(" << a1_name <<
","
114 << a2_name <<
") = relErr(" << a1 <<
"," << a2 <<
") = "
115 << err <<
" <= tol = " << tol <<
": failed!\n";
122 template<
class ValueType,
class VectorType1,
class VectorType2>
124 const std::string& a1_name,
125 const VectorType2&a2,
126 const std::string& a2_name,
127 const ValueType& rel_tol,
const ValueType& abs_tol,
128 Teuchos::FancyOStream& out)
130 using value_t_1 =
typename VectorType1::value_type;
131 using value_t_2 =
typename VectorType2::value_type;
132 static_assert(std::is_same<value_t_1,value_t_2>::value,
"Inconsistent types.");
136 out <<
"Comparing " << a1_name <<
" == " << a2_name <<
" ... ";
138 const int n = a1.size();
141 if (a2.size() != n) {
142 out <<
"\nError, "<<a1_name<<
".size() = "<<a1.size()<<
" == "
143 << a2_name<<
".size() = "<<a2.size()<<
" : failed!\n";
148 for(
int i = 0; i < n; ++i ) {
149 ValueType err = abs(a1.coeff(i) - a2.coeff(i));
151 abs_tol + rel_tol*std::max(abs(a1.fastAccessCoeff(i)),
152 abs(a2.fastAccessCoeff(i)));
155 <<
"\nError, relErr("<<a1_name<<
"["<<i<<
"],"
156 <<a2_name<<
"["<<i<<
"]) = relErr("<<a1.coeff(i)<<
","<<a2.coeff(i)
157 <<
") = "<<err<<
" <= tol = "<<tol<<
": failed!\n";
166 << a1_name <<
" = " << a1 << std::endl
167 << a2_name <<
" = " << a2 << std::endl;
173 template<
class Array1,
class Array2,
class ValueType>
175 const Array2& a2,
const std::string& a2_name,
176 const ValueType& rel_tol,
177 const ValueType& abs_tol,
178 Teuchos::FancyOStream& out)
183 out <<
"Comparing " << a1_name <<
" == " << a2_name <<
" ... ";
185 const int n = a1.size();
188 if (as<int>(a2.size()) != n) {
189 out <<
"\nError, "<<a1_name<<
".size() = "<<a1.size()<<
" == "
190 << a2_name<<
".size() = "<<a2.size()<<
" : failed!\n";
195 for(
int i = 0; i < n; ++i ) {
196 ValueType err = std::abs(a1[i] - a2[i]);
198 abs_tol + rel_tol*std::max(std::abs(a1[i]),std::abs(a2[i]));
200 out <<
"\nError, relErr(" << a1_name <<
"[" << i <<
"]," << a2_name
201 <<
"[" << i <<
"]) = relErr(" << a1[i] <<
"," <<a2[i] <<
") = "
202 << err <<
" <= tol = " << tol <<
": failed!\n";
213 template<
class ordinal_type,
class scalar_type>
215 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& a1,
216 const std::string& a1_name,
217 const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& a2,
218 const std::string& a2_name,
219 const scalar_type& rel_tol,
220 const scalar_type& abs_tol,
221 Teuchos::FancyOStream& out)
226 out <<
"Comparing " << a1_name <<
" == " << a2_name <<
" ... ";
228 const ordinal_type m = a1.numRows();
229 const ordinal_type n = a1.numCols();
232 if (a2.numRows() != m) {
233 out <<
"\nError, "<<a1_name<<
".numRows() = "<<a1.numRows()<<
" == "
234 << a2_name<<
".numRows() = "<<a2.numRows()<<
" : failed!\n";
239 if (a2.numCols() != n) {
240 out <<
"\nError, "<<a1_name<<
".numCols() = "<<a1.numCols()<<
" == "
241 << a2_name<<
".numCols() = "<<a2.numCols()<<
" : failed!\n";
246 for (ordinal_type i=0; i<m; i++) {
247 for (ordinal_type
j=0;
j<n;
j++) {
248 scalar_type err = std::abs(a1(i,
j) - a2(i,
j));
250 abs_tol + rel_tol*std::max(std::abs(a1(i,
j)),std::abs(a2(i,
j)));
252 out <<
"\nError, relErr(" << a1_name <<
"(" << i <<
"," <<
j <<
"), "
253 << a2_name <<
"(" << i <<
"," <<
j <<
")) = relErr("
254 << a1(i,
j) <<
", " <<a2(i,
j) <<
") = "
255 << err <<
" <= tol = " << tol <<
": failed!\n";
267 template<
class ordinal_type,
class scalar_type>
270 const std::string& cijk1_name,
272 const std::string& cijk2_name,
273 const scalar_type& rel_tol,
274 const scalar_type& abs_tol,
275 Teuchos::FancyOStream& out)
280 out <<
"Comparing " << cijk1_name <<
" == " << cijk2_name <<
" ... ";
287 for (
typename Cijk_type::k_iterator k_it=Cijk2.
k_begin();
288 k_it!=Cijk2.
k_end(); ++k_it) {
289 ordinal_type k = Stokhos::index(k_it);
290 for (
typename Cijk_type::kj_iterator j_it = Cijk2.
j_begin(k_it);
291 j_it != Cijk2.
j_end(k_it); ++j_it) {
292 ordinal_type
j = Stokhos::index(j_it);
293 for (
typename Cijk_type::kji_iterator i_it = Cijk2.
i_begin(j_it);
294 i_it != Cijk2.
i_end(j_it); ++i_it) {
295 ordinal_type i = Stokhos::index(i_it);
297 scalar_type c2 = Stokhos::value(i_it);
299 double tol = abs_tol + c2*rel_tol;
300 double err = std::abs(c1-c2);
304 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )"
305 <<
" = " <<
"rel_err( " << c1 <<
", " << c2 <<
" ) = " << err
306 <<
" <= " << tol <<
" : ";
307 if (s) out <<
"Passed.";
312 success = success && s;
320 template<
class ordinal_type,
class scalar_type>
324 const scalar_type& sparse_tol,
325 const scalar_type& rel_tol,
326 const scalar_type& abs_tol,
327 Teuchos::FancyOStream& out,
331 ordinal_type sz = basis.
size();
332 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,scalar_type> > > bases = basis.
getCoordinateBases();
334 Teuchos::Array< Teuchos::RCP<Stokhos::Dense3Tensor<ordinal_type,scalar_type> > > Cijk_1d(d);
335 for (ordinal_type i=0; i<d; i++)
336 Cijk_1d[i] = bases[i]->computeTripleProductTensor();
337 for (ordinal_type
j=0;
j<sz;
j++) {
339 for (ordinal_type i=0; i<sz; i++) {
341 for (ordinal_type k=0; k<sz; k++) {
345 for (ordinal_type l=0; l<d; l++)
346 c *= (*Cijk_1d[l])(terms_i[l],terms_j[l],terms_k[l]);
359 rel_tol, abs_tol, out);
364 template <
typename operator_type1,
typename operator_type2,
365 typename scalar_type>
367 const operator_type2& op2,
368 const scalar_type& rel_tol,
369 const scalar_type& abs_tol,
370 Teuchos::FancyOStream& out) {
373 typedef typename operator_type1::const_set_iterator point_iterator_type;
376 TEUCHOS_TEST_EQUALITY(op1.point_size(), op2.point_size(), out, success);
377 if (!success)
return false;
380 point_iterator_type pi1 = op1.set_begin();
381 point_iterator_type pi2 = op2.set_begin();
382 while (pi1 != op1.set_end() || pi2 != op2.set_end()) {
383 std::stringstream ss1, ss2;
388 pi2->first, ss2.str(),
389 rel_tol, abs_tol, out);
391 std::stringstream ss3, ss4;
392 ss3 << pi1->second.first;
393 ss4 << pi2->second.first;
396 pi2->second.first, ss4.str(),
397 rel_tol, abs_tol, out);
405 template <
typename operator_type1,
typename operator_type2,
406 typename func_type1,
typename func_type2,
407 typename scalar_type>
409 const operator_type2& op2,
410 const func_type1& func1,
411 const func_type2& func2,
412 const scalar_type& rel_tol,
413 const scalar_type& abs_tol,
414 Teuchos::FancyOStream& out) {
417 typedef typename operator_type1::ordinal_type ordinal_type;
418 typedef typename operator_type1::value_type value_type;
419 typedef typename operator_type1::const_iterator point_iterator_type;
422 TEUCHOS_TEST_EQUALITY(op1.coeff_size(), op2.coeff_size(), out, success);
426 ordinal_type point_sz1 = op1.point_size();
427 ordinal_type point_sz2 = op2.point_size();
428 ordinal_type coeff_sz = op1.coeff_size();
431 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
f(point_sz1,2);
432 ordinal_type idx = 0;
433 for (point_iterator_type pi = op1.begin(); pi != op1.end(); ++pi) {
434 f(idx,0) = func1(*pi);
435 f(idx,1) = func2(*pi);
439 Teuchos::SerialDenseMatrix<ordinal_type,value_type> f2(point_sz2,2);
441 for (point_iterator_type pi = op2.begin(); pi != op2.end(); ++pi) {
442 f2(idx,0) = func1(*pi);
443 f2(idx,1) = func2(*pi);
448 if (point_sz1 == point_sz2)
453 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x(coeff_sz,2);
454 op1.transformQP2PCE(1.0,
f, x, 0.0);
456 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x2(coeff_sz,2);
457 op2.transformQP2PCE(1.0, f2, x2, 0.0);
466 template <
typename operator_type1,
typename operator_type2,
467 typename func_type1,
typename func_type2,
468 typename scalar_type>
470 const operator_type2& op2,
471 const func_type1& func1,
472 const func_type2& func2,
473 const scalar_type& rel_tol,
474 const scalar_type& abs_tol,
475 Teuchos::FancyOStream& out) {
478 typedef typename operator_type1::ordinal_type ordinal_type;
479 typedef typename operator_type1::value_type value_type;
480 typedef typename operator_type1::const_iterator point_iterator_type;
483 TEUCHOS_TEST_EQUALITY(op1.coeff_size(), op2.coeff_size(), out, success);
487 ordinal_type point_sz1 = op1.point_size();
488 ordinal_type point_sz2 = op2.point_size();
489 ordinal_type coeff_sz = op1.coeff_size();
492 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
f(2,point_sz1);
493 ordinal_type idx = 0;
494 for (point_iterator_type pi = op1.begin(); pi != op1.end(); ++pi) {
495 f(0,idx) = func1(*pi);
496 f(1,idx) = func2(*pi);
500 Teuchos::SerialDenseMatrix<ordinal_type,value_type> f2(2,point_sz2);
502 for (point_iterator_type pi = op2.begin(); pi != op2.end(); ++pi) {
503 f2(0,idx) = func1(*pi);
504 f2(1,idx) = func2(*pi);
509 if (point_sz1 == point_sz2)
514 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x(2,coeff_sz);
515 op1.transformQP2PCE(1.0,
f, x, 0.0,
true);
517 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x2(2,coeff_sz);
518 op2.transformQP2PCE(1.0, f2, x2, 0.0,
true);
527 template <
typename basis_type,
typename operator_type,
typename scalar_type>
529 const operator_type& op,
530 const scalar_type& rel_tol,
531 const scalar_type& abs_tol,
532 Teuchos::FancyOStream& out) {
533 typedef typename operator_type::ordinal_type ordinal_type;
534 typedef typename operator_type::value_type value_type;
536 ordinal_type coeff_sz = op.coeff_size();
537 ordinal_type point_sz = op.point_size();
538 Teuchos::SerialDenseMatrix<ordinal_type,value_type> eye(coeff_sz,coeff_sz);
540 for (ordinal_type i=0; i<coeff_sz; ++i)
544 Teuchos::SerialDenseMatrix<ordinal_type,value_type>
f(point_sz,coeff_sz);
545 op.transformPCE2QP(1.0, eye,
f, 0.0);
548 Teuchos::SerialDenseMatrix<ordinal_type,value_type> x(coeff_sz,coeff_sz);
549 op.transformQP2PCE(1.0,
f, x, 0.0);
552 for (ordinal_type i=0; i<coeff_sz; ++i)
556 Teuchos::SerialDenseMatrix<ordinal_type,value_type> z(coeff_sz,coeff_sz);
559 out <<
"Discrete orthogonality error = " << x.normInf() << std::endl;
Legendre polynomial basis.
A multidimensional index.
ordinal_type order() const
Compute total order of index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
ordinal_type size() const
Return size.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual ordinal_type size() const =0
Return total size of basis.
virtual ordinal_type dimension() const =0
Return dimension of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const =0
Get orders of each coordinate polynomial given an index i.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
Get Cijk value for a given i, j, k indices.
void fillComplete()
Signal all terms have been added.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
ordinal_type num_entries() const
Return number of non-zero entries.
k_iterator k_end() const
Iterator pointing to last k entry.
void add_term(ordinal_type i, ordinal_type j, ordinal_type k, const value_type &c)
Add new term for given (i,j,k)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Top-level namespace for Stokhos classes and functions.
bool testSparse3Tensor(const Stokhos::Sparse3Tensor< ordinal_type, scalar_type > &Cijk, const Stokhos::ProductBasis< ordinal_type, scalar_type > &basis, const scalar_type &sparse_tol, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out, bool linear=false)
bool compareSDM(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &a1, const std::string &a1_name, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &a2, const std::string &a2_name, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
bool testPseudoSpectralDiscreteOrthogonality(const basis_type &basis, const operator_type &op, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
bool testPseudoSpectralApplyTrans(const operator_type1 &op1, const operator_type2 &op2, const func_type1 &func1, const func_type2 &func2, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
bool testPseudoSpectralPoints(const operator_type1 &op1, const operator_type2 &op2, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
bool compareSparse3Tensor(const Stokhos::Sparse3Tensor< ordinal_type, scalar_type > &Cijk1, const std::string &cijk1_name, const Stokhos::Sparse3Tensor< ordinal_type, scalar_type > &Cijk2, const std::string &cijk2_name, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
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)
bool compareArrays(const Array1 &a1, const std::string &a1_name, const Array2 &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
bool testPseudoSpectralApply(const operator_type1 &op1, const operator_type2 &op2, const func_type1 &func1, const func_type2 &func2, const scalar_type &rel_tol, const scalar_type &abs_tol, Teuchos::FancyOStream &out)
bool compareVecs(const VectorType1 &a1, const std::string &a1_name, const VectorType2 &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
bool comparePCEs(const PCEType &a1, const std::string &a1_name, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)