Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_NormalizedHermiteBasisUnitTest.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_UnitTestHarness.hpp"
11 #include "Teuchos_TestingHelpers.hpp"
12 #include "Teuchos_UnitTestRepository.hpp"
13 #include "Teuchos_GlobalMPISession.hpp"
14 
15 #include "Stokhos.hpp"
16 #include "Stokhos_UnitTestHelpers.hpp"
17 #ifdef HAVE_STOKHOS_FORUQTK
18 #include "Stokhos_gaussq.h"
19 #endif
20 
21 namespace HermiteBasisUnitTest {
22 
23  // Common setup for unit tests
24  template <typename OrdinalType, typename ValueType>
25  struct UnitTestSetup {
26  ValueType rtol, atol;
27  OrdinalType p;
28  Stokhos::HermiteBasis<OrdinalType,ValueType> basis;
29 
30  UnitTestSetup() : rtol(1e-12), atol(1e-7), p(10), basis(p,true) {}
31 
32  };
33 
34  UnitTestSetup<int,double> setup;
35 
36  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Order ) {
37  int order = setup.basis.order();
38  TEUCHOS_TEST_EQUALITY(order, setup.p, out, success);
39  }
40 
41  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Size ) {
42  int size = setup.basis.size();
43  TEUCHOS_TEST_EQUALITY(size, setup.p+1, out, success);
44  }
45 
46  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm ) {
47  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
48  Teuchos::Array<double> n2(setup.p+1);
49  n2[0] = 1.0;
50  for (int i=1; i<=setup.p; i++)
51  n2[i] = 1.0;
52  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
53  out);
54  }
55 
56  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm2 ) {
57  Teuchos::Array<double> n1(setup.p+1);
58  Teuchos::Array<double> n2(setup.p+1);
59  n1[0] = setup.basis.norm_squared(0);
60  n2[0] = 1.0;
61  for (int i=1; i<=setup.p; i++) {
62  n1[i] = setup.basis.norm_squared(i);
63  n2[i] = 1.0;
64  }
65  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
66  out);
67  }
68 
69  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadNorm ) {
70  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
71  Teuchos::Array<double> n2(setup.p+1);
72  Teuchos::Array<double> x, w;
73  Teuchos::Array< Teuchos::Array<double> > v;
74  setup.basis.getQuadPoints(2*setup.p, x, w, v);
75  int nqp = w.size();
76  for (int i=0; i<=setup.p; i++) {
77  n2[i] = 0;
78  for (int j=0; j<nqp; j++)
79  n2[i] += w[j]*v[j][i]*v[j][i];
80  }
81  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
82  out);
83  }
84 
85  // Tests basis is orthogonal
86  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadOrthog ) {
87  const Teuchos::Array<double>& norms = setup.basis.norm_squared();
88  Teuchos::Array<double> x, w;
89  Teuchos::Array< Teuchos::Array<double> > v;
90  setup.basis.getQuadPoints(2*setup.p, x, w, v);
91  int nqp = w.size();
92  Teuchos::SerialDenseMatrix<int,double> mat(setup.p+1, setup.p+1);
93  for (int i=0; i<=setup.p; i++) {
94  for (int j=0; j<=setup.p; j++) {
95  for (int k=0; k<nqp; k++)
96  mat(i,j) += w[k]*v[k][i]*v[k][j];
97  mat(i,j) /= std::sqrt(norms[i]*norms[j]);
98  }
99  mat(i,i) -= 1.0;
100  }
101  success = mat.normInf() < setup.atol;
102  if (!success) {
103  out << "\n Error, mat.normInf() < atol = " << mat.normInf()
104  << " < " << setup.atol << ": failed!\n";
105  out << "mat = " << printMat(mat) << std::endl;
106  }
107  }
108 
109  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, TripleProduct ) {
110  Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
111  setup.basis.computeTripleProductTensor();
112 
113  Teuchos::Array<double> x, w;
114  Teuchos::Array< Teuchos::Array<double> > v;
115  setup.basis.getQuadPoints(3*setup.p, x, w, v);
116 
117  success = true;
118  for (int i=0; i<=setup.p; i++) {
119  for (int j=0; j<=setup.p; j++) {
120  for (int k=0; k<=setup.p; k++) {
121  double c = 0.0;
122  int nqp = w.size();
123  for (int qp=0; qp<nqp; qp++)
124  c += w[qp]*v[qp][i]*v[qp][j]*v[qp][k];
125  std::stringstream ss;
126  ss << "Cijk(" << i << "," << j << "," << k << ")";
127  success = success &&
128  Stokhos::compareValues((*Cijk)(i,j,k), ss.str(), c, "c",
129  setup.rtol, setup.atol, out);
130  }
131  }
132  }
133  }
134 
135  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct ) {
136  Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
137  setup.basis.computeDerivDoubleProductTensor();
138 
139  Teuchos::Array<double> x, w;
140  Teuchos::Array< Teuchos::Array<double> > v, val, deriv;
141  setup.basis.getQuadPoints(2*setup.p, x, w, v);
142  int nqp = w.size();
143  val.resize(nqp);
144  deriv.resize(nqp);
145  for (int i=0; i<nqp; i++) {
146  val[i].resize(setup.p+1);
147  deriv[i].resize(setup.p+1);
148  setup.basis.evaluateBasesAndDerivatives(x[i], val[i], deriv[i]);
149  }
150 
151  success = true;
152  for (int i=0; i<=setup.p; i++) {
153  for (int j=0; j<=setup.p; j++) {
154  double b = 0.0;
155  for (int qp=0; qp<nqp; qp++)
156  b += w[qp]*deriv[qp][i]*val[qp][j];
157  std::stringstream ss;
158  ss << "Bij(" << i << "," << j << ")";
159  success = success &&
160  Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
161  setup.rtol, setup.atol, out);
162  }
163  }
164  }
165 
166  // TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct2 ) {
167  // Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
168  // setup.basis.computeDerivDoubleProductTensor();
169  // const Teuchos::Array<double>& n = setup.basis.norm_squared();
170 
171  // success = true;
172  // for (int i=0; i<=setup.p; i++) {
173  // for (int j=0; j<=setup.p; j++) {
174  // double b = 0.0;
175  // if (j == i-1)
176  // b = i*n[j];
177  // std::stringstream ss;
178  // ss << "Bij(" << i << "," << j << ")";
179  // success = success &&
180  // Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
181  // setup.rtol, setup.atol, out);
182  // }
183  // }
184  // }
185 
186  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBases ) {
187  double x = 0.1234;
188  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
189  setup.basis.evaluateBases(x, v1);
190 
191  // evaluate bases using formula
192  // He_0(x) = 1
193  // He_1(x) = x
194  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
195  v2[0] = 1.0;
196  if (setup.p >= 1)
197  v2[1] = x;
198  for (int i=2; i<=setup.p; i++)
199  v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
200  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
201  out);
202  }
203 
204  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBasesAndDerivatives ) {
205  double x = 0.1234;
206  Teuchos::Array<double> val1(setup.p+1), deriv1(setup.p+1),
207  val2(setup.p+1), deriv2(setup.p+1);
208  setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
209 
210  // evaluate bases and derivatives using formula:
211  // He_0(x) = 1
212  // He_1(x) = x
213  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
214  val2[0] = 1.0;
215  if (setup.p >= 1)
216  val2[1] = x;
217  for (int i=2; i<=setup.p; i++)
218  val2[i] = (x*val2[i-1] - std::sqrt(i-1.0)*val2[i-2])/std::sqrt(1.0*i);
219 
220  deriv2[0] = 0.0;
221  if (setup.p >= 1)
222  deriv2[1] = 1.0;
223  for (int i=2; i<=setup.p; i++)
224  deriv2[i] = (val2[i-1] + x*deriv2[i-1] - std::sqrt(i-1.0)*deriv2[i-2]) /
225  std::sqrt(1.0*i);
226  success = Stokhos::compareArrays(val1, "val1", val2, "val2",
227  setup.rtol, setup.atol, out);
228  success = success &&
229  Stokhos::compareArrays(deriv1, "deriv1", deriv2, "deriv2",
230  setup.rtol, setup.atol, out);
231  }
232 
233  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Evaluate ) {
234  double x = 0.1234;
235  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
236  for (int i=0; i<=setup.p; i++)
237  v1[i] = setup.basis.evaluate(x, i);
238 
239  // evaluate bases using formula
240  // He_0(x) = 1
241  // He_1(x) = x
242  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
243  v2[0] = 1.0;
244  if (setup.p >= 1)
245  v2[1] = x;
246  for (int i=2; i<=setup.p; i++)
247  v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
248  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
249  out);
250  }
251 
252  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Recurrence ) {
253  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
254  d1(setup.p+1);
255  Teuchos::Array<double> a2(setup.p+1), b2(setup.p+1), c2(setup.p+1),
256  d2(setup.p+1);
257  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
258 
259  // compute coefficients using formula
260  a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
261  for (int i=1; i<=setup.p; i++) {
262  a2[i] = 0.0;
263  b2[i] = std::sqrt(1.0*i);
264  c2[i] = 1.0;
265  d2[i] = std::sqrt(1.0*i);
266  }
267  success = true;
268  success = success &&
269  Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol, out);
270  success = success &&
271  Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol, out);
272  success = success &&
273  Stokhos::compareArrays(c1, "c1", c2, "c2", setup.rtol, setup.atol, out);
274  success = success &&
275  Stokhos::compareArrays(d1, "d1", d2, "d2", setup.rtol, setup.atol, out);
276  }
277 
278  // Tests alpha coefficients satisfy
279  // alpha_k = delta_k * (t*psi_k,psi_k)/(psi_k,psi_k)
280  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceAlpha ) {
281  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
282  d1(setup.p+1);
283  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
284 
285  Teuchos::Array<double> a2(setup.p+1);
286  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
287  Teuchos::Array<double> x, w;
288  Teuchos::Array< Teuchos::Array<double> > v;
289  setup.basis.getQuadPoints(2*setup.p, x, w, v);
290  int nqp = w.size();
291  for (int i=0; i<=setup.p; i++) {
292  a2[i] = 0;
293  for (int j=0; j<nqp; j++)
294  a2[i] += w[j]*x[j]*v[j][i]*v[j][i];
295  a2[i] = a2[i]*c1[i]/n1[i];
296  }
297  success = Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol,
298  out);
299  }
300 
301  // Tests beta coefficients satisfy
302  // beta_k =
303  // gamma_k * delta_k/delta_{k-1} * (psi_k,psi_k)/(psi_{k-1},psi_{k-1})
304  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceBeta ) {
305  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
306  d1(setup.p+1);
307  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
308 
309  Teuchos::Array<double> b2(setup.p+1);
310  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
311  b2[0] = b1[0];
312  for (int i=1; i<=setup.p; i++) {
313  b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
314  }
315  success = Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol,
316  out);
317  }
318 
319 #ifdef HAVE_STOKHOS_DAKOTA
320  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsDakota ) {
321  int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
322  Teuchos::Array<double> x1, w1;
323  Teuchos::Array< Teuchos::Array<double> > v1;
324  setup.basis.getQuadPoints(setup.p, x1, w1, v1);
325 
326  Teuchos::Array<double> x2(n), w2(n);
327  Teuchos::Array< Teuchos::Array<double> > v2(n);
328  webbur::hermite_compute(n, &x2[0], &w2[0]);
329 
330  for (int i=0; i<n; i++) {
331  w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
332  x2[i] *= std::sqrt(2.0);
333  v2[i].resize(setup.p+1);
334  setup.basis.evaluateBases(x2[i], v2[i]);
335  }
336  success = true;
337  success = success &&
338  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
339  success = success &&
340  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
341  for (int i=0; i<n; i++) {
342  std::stringstream ss1, ss2;
343  ss1 << "v1[" << i << "]";
344  ss2 << "v2[" << i << "]";
345  success = success &&
346  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
347  setup.rtol, setup.atol, out);
348  }
349  }
350 #endif
351 
352 #ifdef HAVE_STOKHOS_FORUQTK
353  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsForUQTK ) {
354  int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
355  Teuchos::Array<double> x1, w1;
356  Teuchos::Array< Teuchos::Array<double> > v1;
357  setup.basis.getQuadPoints(setup.p, x1, w1, v1);
358 
359  Teuchos::Array<double> x2(n), w2(n);
360  Teuchos::Array< Teuchos::Array<double> > v2(n);
361  int kind = 4;
362  int kpts = 0;
363  double endpts[2] = {0.0, 0.0};
364  Teuchos::Array<double> b(n);
365  double alpha = 0.0;
366  double beta = 0.0;
367  GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
368 
369  for (int i=0; i<n; i++) {
370  w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
371  x2[i] *= std::sqrt(2.0);
372  v2[i].resize(setup.p+1);
373  setup.basis.evaluateBases(x2[i], v2[i]);
374  }
375  success = true;
376  success = success &&
377  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
378  success = success &&
379  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
380  for (int i=0; i<n; i++) {
381  std::stringstream ss1, ss2;
382  ss1 << "v1[" << i << "]";
383  ss2 << "v2[" << i << "]";
384  success = success &&
385  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
386  setup.rtol, setup.atol, out);
387  }
388  }
389 #endif
390 
391 }
392 
393 int main( int argc, char* argv[] ) {
394  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
395  return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
396 }
Sacado::UQ::sqrt
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:847
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:252
printMat
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Stokhos::compareArrays
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)
Definition: Stokhos_UnitTestHelpers.hpp:142
Stokhos::compareValues
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)
Definition: Stokhos_UnitTestHelpers.hpp:69
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Teuchos_UnitTestRepository.hpp
Teuchos::UnitTestRepository::runUnitTestsFromMain
static int runUnitTestsFromMain(int argc, char *argv[])
Sacado::UQ::ceil
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1160
Teuchos_UnitTestHarness.hpp
HermiteBasisUnitTest::setup
UnitTestSetup< int, double > setup
Definition: Stokhos_HermiteBasisUnitTest.cpp:34
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Stokhos_UnitTestHelpers.hpp
Sacado::UQ::atan
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1070
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:25
Teuchos::SerialDenseMatrix::normInf
ScalarTraits< ScalarType >::magnitudeType normInf() const
Teuchos_TestingHelpers.hpp
val
expr val()
TEUCHOS_TEST_EQUALITY
#define TEUCHOS_TEST_EQUALITY(v1, v2, out, success)
Stokhos_gaussq.h
HermiteBasisUnitTest::UnitTestSetup::basis
Stokhos::HermiteBasis< OrdinalType, ValueType > basis
Definition: Stokhos_HermiteBasisUnitTest.cpp:28
Teuchos_GlobalMPISession.hpp
Teuchos::Array::size
size_type size() const
true
true
Definition: Stokhos_LanczosUnitTest.cpp:233
HermiteBasisUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Stokhos_HermiteBasis, Order)
Definition: Stokhos_HermiteBasisUnitTest.cpp:36
n
int n
GAUSSQ_F77
void GAUSSQ_F77(int *kind, int *n, double *alpha, double *beta, int *kpts, double *endpts, double *b, double *x, double *w)
Generated on Fri Nov 22 2024 09:18:13 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5