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 // $Id$
2 // $Source$
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Stokhos Package
7 // Copyright (2009) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include "Teuchos_UnitTestHarness.hpp"
45 #include "Teuchos_TestingHelpers.hpp"
46 #include "Teuchos_UnitTestRepository.hpp"
47 #include "Teuchos_GlobalMPISession.hpp"
48 
49 #include "Stokhos.hpp"
50 #include "Stokhos_UnitTestHelpers.hpp"
51 #ifdef HAVE_STOKHOS_FORUQTK
52 #include "Stokhos_gaussq.h"
53 #endif
54 
55 namespace HermiteBasisUnitTest {
56 
57  // Common setup for unit tests
58  template <typename OrdinalType, typename ValueType>
59  struct UnitTestSetup {
60  ValueType rtol, atol;
61  OrdinalType p;
62  Stokhos::HermiteBasis<OrdinalType,ValueType> basis;
63 
64  UnitTestSetup() : rtol(1e-12), atol(1e-7), p(10), basis(p,true) {}
65 
66  };
67 
68  UnitTestSetup<int,double> setup;
69 
70  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Order ) {
71  int order = setup.basis.order();
72  TEUCHOS_TEST_EQUALITY(order, setup.p, out, success);
73  }
74 
75  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Size ) {
76  int size = setup.basis.size();
77  TEUCHOS_TEST_EQUALITY(size, setup.p+1, out, success);
78  }
79 
80  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm ) {
81  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
82  Teuchos::Array<double> n2(setup.p+1);
83  n2[0] = 1.0;
84  for (int i=1; i<=setup.p; i++)
85  n2[i] = 1.0;
86  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
87  out);
88  }
89 
90  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Norm2 ) {
91  Teuchos::Array<double> n1(setup.p+1);
92  Teuchos::Array<double> n2(setup.p+1);
93  n1[0] = setup.basis.norm_squared(0);
94  n2[0] = 1.0;
95  for (int i=1; i<=setup.p; i++) {
96  n1[i] = setup.basis.norm_squared(i);
97  n2[i] = 1.0;
98  }
99  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
100  out);
101  }
102 
103  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadNorm ) {
104  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
105  Teuchos::Array<double> n2(setup.p+1);
106  Teuchos::Array<double> x, w;
107  Teuchos::Array< Teuchos::Array<double> > v;
108  setup.basis.getQuadPoints(2*setup.p, x, w, v);
109  int nqp = w.size();
110  for (int i=0; i<=setup.p; i++) {
111  n2[i] = 0;
112  for (int j=0; j<nqp; j++)
113  n2[i] += w[j]*v[j][i]*v[j][i];
114  }
115  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
116  out);
117  }
118 
119  // Tests basis is orthogonal
120  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadOrthog ) {
121  const Teuchos::Array<double>& norms = setup.basis.norm_squared();
122  Teuchos::Array<double> x, w;
123  Teuchos::Array< Teuchos::Array<double> > v;
124  setup.basis.getQuadPoints(2*setup.p, x, w, v);
125  int nqp = w.size();
126  Teuchos::SerialDenseMatrix<int,double> mat(setup.p+1, setup.p+1);
127  for (int i=0; i<=setup.p; i++) {
128  for (int j=0; j<=setup.p; j++) {
129  for (int k=0; k<nqp; k++)
130  mat(i,j) += w[k]*v[k][i]*v[k][j];
131  mat(i,j) /= std::sqrt(norms[i]*norms[j]);
132  }
133  mat(i,i) -= 1.0;
134  }
135  success = mat.normInf() < setup.atol;
136  if (!success) {
137  out << "\n Error, mat.normInf() < atol = " << mat.normInf()
138  << " < " << setup.atol << ": failed!\n";
139  out << "mat = " << printMat(mat) << std::endl;
140  }
141  }
142 
143  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, TripleProduct ) {
144  Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
145  setup.basis.computeTripleProductTensor();
146 
147  Teuchos::Array<double> x, w;
148  Teuchos::Array< Teuchos::Array<double> > v;
149  setup.basis.getQuadPoints(3*setup.p, x, w, v);
150 
151  success = true;
152  for (int i=0; i<=setup.p; i++) {
153  for (int j=0; j<=setup.p; j++) {
154  for (int k=0; k<=setup.p; k++) {
155  double c = 0.0;
156  int nqp = w.size();
157  for (int qp=0; qp<nqp; qp++)
158  c += w[qp]*v[qp][i]*v[qp][j]*v[qp][k];
159  std::stringstream ss;
160  ss << "Cijk(" << i << "," << j << "," << k << ")";
161  success = success &&
162  Stokhos::compareValues((*Cijk)(i,j,k), ss.str(), c, "c",
163  setup.rtol, setup.atol, out);
164  }
165  }
166  }
167  }
168 
169  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct ) {
170  Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
171  setup.basis.computeDerivDoubleProductTensor();
172 
173  Teuchos::Array<double> x, w;
174  Teuchos::Array< Teuchos::Array<double> > v, val, deriv;
175  setup.basis.getQuadPoints(2*setup.p, x, w, v);
176  int nqp = w.size();
177  val.resize(nqp);
178  deriv.resize(nqp);
179  for (int i=0; i<nqp; i++) {
180  val[i].resize(setup.p+1);
181  deriv[i].resize(setup.p+1);
182  setup.basis.evaluateBasesAndDerivatives(x[i], val[i], deriv[i]);
183  }
184 
185  success = true;
186  for (int i=0; i<=setup.p; i++) {
187  for (int j=0; j<=setup.p; j++) {
188  double b = 0.0;
189  for (int qp=0; qp<nqp; qp++)
190  b += w[qp]*deriv[qp][i]*val[qp][j];
191  std::stringstream ss;
192  ss << "Bij(" << i << "," << j << ")";
193  success = success &&
194  Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
195  setup.rtol, setup.atol, out);
196  }
197  }
198  }
199 
200  // TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, DerivDoubleProduct2 ) {
201  // Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
202  // setup.basis.computeDerivDoubleProductTensor();
203  // const Teuchos::Array<double>& n = setup.basis.norm_squared();
204 
205  // success = true;
206  // for (int i=0; i<=setup.p; i++) {
207  // for (int j=0; j<=setup.p; j++) {
208  // double b = 0.0;
209  // if (j == i-1)
210  // b = i*n[j];
211  // std::stringstream ss;
212  // ss << "Bij(" << i << "," << j << ")";
213  // success = success &&
214  // Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
215  // setup.rtol, setup.atol, out);
216  // }
217  // }
218  // }
219 
220  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBases ) {
221  double x = 0.1234;
222  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
223  setup.basis.evaluateBases(x, v1);
224 
225  // evaluate bases using formula
226  // He_0(x) = 1
227  // He_1(x) = x
228  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
229  v2[0] = 1.0;
230  if (setup.p >= 1)
231  v2[1] = x;
232  for (int i=2; i<=setup.p; i++)
233  v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
234  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
235  out);
236  }
237 
238  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, EvaluateBasesAndDerivatives ) {
239  double x = 0.1234;
240  Teuchos::Array<double> val1(setup.p+1), deriv1(setup.p+1),
241  val2(setup.p+1), deriv2(setup.p+1);
242  setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
243 
244  // evaluate bases and derivatives using formula:
245  // He_0(x) = 1
246  // He_1(x) = x
247  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
248  val2[0] = 1.0;
249  if (setup.p >= 1)
250  val2[1] = x;
251  for (int i=2; i<=setup.p; i++)
252  val2[i] = (x*val2[i-1] - std::sqrt(i-1.0)*val2[i-2])/std::sqrt(1.0*i);
253 
254  deriv2[0] = 0.0;
255  if (setup.p >= 1)
256  deriv2[1] = 1.0;
257  for (int i=2; i<=setup.p; i++)
258  deriv2[i] = (val2[i-1] + x*deriv2[i-1] - std::sqrt(i-1.0)*deriv2[i-2]) /
259  std::sqrt(1.0*i);
260  success = Stokhos::compareArrays(val1, "val1", val2, "val2",
261  setup.rtol, setup.atol, out);
262  success = success &&
263  Stokhos::compareArrays(deriv1, "deriv1", deriv2, "deriv2",
264  setup.rtol, setup.atol, out);
265  }
266 
267  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Evaluate ) {
268  double x = 0.1234;
269  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
270  for (int i=0; i<=setup.p; i++)
271  v1[i] = setup.basis.evaluate(x, i);
272 
273  // evaluate bases using formula
274  // He_0(x) = 1
275  // He_1(x) = x
276  // He_i(x) = (x*He_{i-1}(x) - sqrt(i-1)*He_{i-2}(x))/sqrt(i), i=2,3,...
277  v2[0] = 1.0;
278  if (setup.p >= 1)
279  v2[1] = x;
280  for (int i=2; i<=setup.p; i++)
281  v2[i] = (x*v2[i-1] - std::sqrt(i-1.0)*v2[i-2])/std::sqrt(1.0*i);
282  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
283  out);
284  }
285 
286  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, Recurrence ) {
287  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
288  d1(setup.p+1);
289  Teuchos::Array<double> a2(setup.p+1), b2(setup.p+1), c2(setup.p+1),
290  d2(setup.p+1);
291  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
292 
293  // compute coefficients using formula
294  a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
295  for (int i=1; i<=setup.p; i++) {
296  a2[i] = 0.0;
297  b2[i] = std::sqrt(1.0*i);
298  c2[i] = 1.0;
299  d2[i] = std::sqrt(1.0*i);
300  }
301  success = true;
302  success = success &&
303  Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol, out);
304  success = success &&
305  Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol, out);
306  success = success &&
307  Stokhos::compareArrays(c1, "c1", c2, "c2", setup.rtol, setup.atol, out);
308  success = success &&
309  Stokhos::compareArrays(d1, "d1", d2, "d2", setup.rtol, setup.atol, out);
310  }
311 
312  // Tests alpha coefficients satisfy
313  // alpha_k = delta_k * (t*psi_k,psi_k)/(psi_k,psi_k)
314  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceAlpha ) {
315  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
316  d1(setup.p+1);
317  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
318 
319  Teuchos::Array<double> a2(setup.p+1);
320  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
321  Teuchos::Array<double> x, w;
322  Teuchos::Array< Teuchos::Array<double> > v;
323  setup.basis.getQuadPoints(2*setup.p, x, w, v);
324  int nqp = w.size();
325  for (int i=0; i<=setup.p; i++) {
326  a2[i] = 0;
327  for (int j=0; j<nqp; j++)
328  a2[i] += w[j]*x[j]*v[j][i]*v[j][i];
329  a2[i] = a2[i]*c1[i]/n1[i];
330  }
331  success = Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol,
332  out);
333  }
334 
335  // Tests beta coefficients satisfy
336  // beta_k =
337  // gamma_k * delta_k/delta_{k-1} * (psi_k,psi_k)/(psi_{k-1},psi_{k-1})
338  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, RecurrenceBeta ) {
339  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
340  d1(setup.p+1);
341  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
342 
343  Teuchos::Array<double> b2(setup.p+1);
344  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
345  b2[0] = b1[0];
346  for (int i=1; i<=setup.p; i++) {
347  b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
348  }
349  success = Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol,
350  out);
351  }
352 
353 #ifdef HAVE_STOKHOS_DAKOTA
354  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsDakota ) {
355  int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
356  Teuchos::Array<double> x1, w1;
357  Teuchos::Array< Teuchos::Array<double> > v1;
358  setup.basis.getQuadPoints(setup.p, x1, w1, v1);
359 
360  Teuchos::Array<double> x2(n), w2(n);
361  Teuchos::Array< Teuchos::Array<double> > v2(n);
362  webbur::hermite_compute(n, &x2[0], &w2[0]);
363 
364  for (int i=0; i<n; i++) {
365  w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
366  x2[i] *= std::sqrt(2.0);
367  v2[i].resize(setup.p+1);
368  setup.basis.evaluateBases(x2[i], v2[i]);
369  }
370  success = true;
371  success = success &&
372  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
373  success = success &&
374  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
375  for (int i=0; i<n; i++) {
376  std::stringstream ss1, ss2;
377  ss1 << "v1[" << i << "]";
378  ss2 << "v2[" << i << "]";
379  success = success &&
380  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
381  setup.rtol, setup.atol, out);
382  }
383  }
384 #endif
385 
386 #ifdef HAVE_STOKHOS_FORUQTK
387  TEUCHOS_UNIT_TEST( Stokhos_NormalizedHermiteBasis, QuadPointsForUQTK ) {
388  int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
389  Teuchos::Array<double> x1, w1;
390  Teuchos::Array< Teuchos::Array<double> > v1;
391  setup.basis.getQuadPoints(setup.p, x1, w1, v1);
392 
393  Teuchos::Array<double> x2(n), w2(n);
394  Teuchos::Array< Teuchos::Array<double> > v2(n);
395  int kind = 4;
396  int kpts = 0;
397  double endpts[2] = {0.0, 0.0};
398  Teuchos::Array<double> b(n);
399  double alpha = 0.0;
400  double beta = 0.0;
401  GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
402 
403  for (int i=0; i<n; i++) {
404  w2[i] *= 0.5/std::sqrt(std::atan(1.0)); // 1/sqrt(pi)
405  x2[i] *= std::sqrt(2.0);
406  v2[i].resize(setup.p+1);
407  setup.basis.evaluateBases(x2[i], v2[i]);
408  }
409  success = true;
410  success = success &&
411  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
412  success = success &&
413  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
414  for (int i=0; i<n; i++) {
415  std::stringstream ss1, ss2;
416  ss1 << "v1[" << i << "]";
417  ss2 << "v2[" << i << "]";
418  success = success &&
419  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
420  setup.rtol, setup.atol, out);
421  }
422  }
423 #endif
424 
425 }
426 
427 int main( int argc, char* argv[] ) {
428  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
429  return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
430 }
Sacado::UQ::sqrt
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:878
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
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:170
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:101
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
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:1191
Teuchos_UnitTestHarness.hpp
HermiteBasisUnitTest::setup
UnitTestSetup< int, double > setup
Definition: Stokhos_HermiteBasisUnitTest.cpp:68
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:1101
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:57
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:62
Teuchos_GlobalMPISession.hpp
Teuchos::Array::size
size_type size() const
true
true
Definition: Stokhos_LanczosUnitTest.cpp:267
HermiteBasisUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Stokhos_HermiteBasis, Order)
Definition: Stokhos_HermiteBasisUnitTest.cpp:70
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 by   doxygen 1.8.5