Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_NormalizedLegendreBasisUnitTest.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 LegendreBasisUnitTest {
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::LegendreBasis<OrdinalType,ValueType> basis;
29 
30  UnitTestSetup() : rtol(1e-12), atol(1e-10), p(10), basis(p,true) {}
31 
32  };
33 
34  UnitTestSetup<int,double> setup;
35 
36  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Order ) {
37  int order = setup.basis.order();
38  TEUCHOS_TEST_EQUALITY(order, setup.p, out, success);
39  }
40 
41  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Size ) {
42  int size = setup.basis.size();
43  TEUCHOS_TEST_EQUALITY(size, setup.p+1, out, success);
44  }
45 
46  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Norm ) {
47  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
48  Teuchos::Array<double> n2(setup.p+1);
49  for (int i=0; i<=setup.p; i++)
50  n2[i] = 1.0;
51  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
52  out);
53  }
54 
55  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Norm2 ) {
56  Teuchos::Array<double> n1(setup.p+1);
57  Teuchos::Array<double> n2(setup.p+1);
58  for (int i=0; i<=setup.p; i++) {
59  n1[i] = setup.basis.norm_squared(i);
60  n2[i] = 1.0;
61  }
62  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
63  out);
64  }
65 
66  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadNorm ) {
67  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
68  Teuchos::Array<double> n2(setup.p+1);
69  Teuchos::Array<double> x, w;
70  Teuchos::Array< Teuchos::Array<double> > v;
71  setup.basis.getQuadPoints(2*setup.p, x, w, v);
72  int nqp = w.size();
73  for (int i=0; i<=setup.p; i++) {
74  n2[i] = 0;
75  for (int j=0; j<nqp; j++)
76  n2[i] += w[j]*v[j][i]*v[j][i];
77  }
78  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
79  out);
80  }
81 
82  // Tests basis is orthogonal
83  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadOrthog ) {
84  const Teuchos::Array<double>& norms = setup.basis.norm_squared();
85  Teuchos::Array<double> x, w;
86  Teuchos::Array< Teuchos::Array<double> > v;
87  setup.basis.getQuadPoints(2*setup.p, x, w, v);
88  int nqp = w.size();
89  Teuchos::SerialDenseMatrix<int,double> mat(setup.p+1, setup.p+1);
90  for (int i=0; i<=setup.p; i++) {
91  for (int j=0; j<=setup.p; j++) {
92  for (int k=0; k<nqp; k++)
93  mat(i,j) += w[k]*v[k][i]*v[k][j];
94  mat(i,j) /= std::sqrt(norms[i]*norms[j]);
95  }
96  mat(i,i) -= 1.0;
97  }
98  success = mat.normInf() < setup.atol;
99  if (!success) {
100  out << "\n Error, mat.normInf() < atol = " << mat.normInf()
101  << " < " << setup.atol << ": failed!\n";
102  out << "mat = " << printMat(mat) << std::endl;
103  }
104  }
105 
106  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, TripleProduct ) {
107  Teuchos::RCP< Stokhos::Dense3Tensor<int, double> > Cijk =
108  setup.basis.computeTripleProductTensor();
109 
110  Teuchos::Array<double> x, w;
111  Teuchos::Array< Teuchos::Array<double> > v;
112  setup.basis.getQuadPoints(3*setup.p, x, w, v);
113 
114  success = true;
115  for (int i=0; i<=setup.p; i++) {
116  for (int j=0; j<=setup.p; j++) {
117  for (int k=0; k<=setup.p; k++) {
118  double c = 0.0;
119  int nqp = w.size();
120  for (int qp=0; qp<nqp; qp++)
121  c += w[qp]*v[qp][i]*v[qp][j]*v[qp][k];
122  std::stringstream ss;
123  ss << "Cijk(" << i << "," << j << "," << k << ")";
124  success = success &&
125  Stokhos::compareValues((*Cijk)(i,j,k), ss.str(), c, "c",
126  setup.rtol, setup.atol, out);
127  }
128  }
129  }
130  }
131 
132  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, DerivDoubleProduct ) {
133  Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
134  setup.basis.computeDerivDoubleProductTensor();
135 
136  Teuchos::Array<double> x, w;
137  Teuchos::Array< Teuchos::Array<double> > v, val, deriv;
138  setup.basis.getQuadPoints(3*setup.p, x, w, v);
139  int nqp = w.size();
140  val.resize(nqp);
141  deriv.resize(nqp);
142  for (int i=0; i<nqp; i++) {
143  val[i].resize(setup.p+1);
144  deriv[i].resize(setup.p+1);
145  setup.basis.evaluateBasesAndDerivatives(x[i], val[i], deriv[i]);
146  }
147 
148  success = true;
149  for (int i=0; i<=setup.p; i++) {
150  for (int j=0; j<=setup.p; j++) {
151  double b = 0.0;
152  for (int qp=0; qp<nqp; qp++)
153  b += w[qp]*deriv[qp][i]*val[qp][j];
154  std::stringstream ss;
155  ss << "Bij(" << i << "," << j << ")";
156  success = success &&
157  Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
158  setup.rtol, setup.atol, out);
159  }
160  }
161  }
162 
163  // TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, DerivDoubleProduct2 ) {
164  // Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
165  // setup.basis.computeDerivDoubleProductTensor();
166  // const Teuchos::Array<double>& n = setup.basis.norm_squared();
167 
168  // Teuchos::Array< Teuchos::Array<double> > deriv_coeffs(setup.p+1);
169  // deriv_coeffs[0].resize(setup.p+1,0.0);
170  // if (setup.p >= 1) {
171  // deriv_coeffs[1].resize(setup.p+1,0.0);
172  // deriv_coeffs[1][0] = 1.0;
173  // }
174  // if (setup.p >= 2) {
175  // deriv_coeffs[2].resize(setup.p+1,0.0);
176  // deriv_coeffs[2][1] = 3.0;
177  // }
178  // for (int k=3; k<=setup.p; k++) {
179  // deriv_coeffs[k].resize(setup.p+1,0.0);
180  // deriv_coeffs[k][0] = 1.0/3.0*deriv_coeffs[k-1][1];
181  // for (int i=1; i<=k-3; i++)
182  // deriv_coeffs[k][i] = i/(2.0*i - 1.0)*deriv_coeffs[k-1][i-1] +
183  // (i+1.0)/(2.0*i + 3.0)*deriv_coeffs[k-1][i+1];
184  // deriv_coeffs[k][k-2] = (k-2.0)/(2.0*k-5.0)*deriv_coeffs[k-1][k-3];
185  // deriv_coeffs[k][k-1] = (k-1.0)/(2.0*k-3.0)*deriv_coeffs[k-1][k-2] + k;
186  // }
187 
188  // success = true;
189  // for (int i=0; i<=setup.p; i++) {
190  // for (int j=0; j<=setup.p; j++) {
191  // double b = deriv_coeffs[i][j]*n[j];
192  // std::stringstream ss;
193  // ss << "Bij(" << i << "," << j << ")";
194  // success = success &&
195  // Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
196  // setup.rtol, setup.atol, out);
197  // }
198  // }
199  // }
200 
201  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, EvaluateBases ) {
202  double x = 0.1234;
203  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
204  setup.basis.evaluateBases(x, v1);
205 
206  double b1, b2;
207  b2 = std::sqrt(1.0/3.0);
208  b1 = b2;
209  v2[0] = 1.0;
210  if (setup.p >= 1)
211  v2[1] = x/b1;
212  for (int i=2; i<=setup.p; i++) {
213  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
214  v2[i] = (x*v2[i-1] - b2*v2[i-2]) / b1;
215  b2 = b1;
216  }
217  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
218  out);
219  }
220 
221  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, EvaluateBasesAndDerivatives ) {
222  double x = 0.1234;
223  Teuchos::Array<double> val1(setup.p+1), deriv1(setup.p+1),
224  val2(setup.p+1), deriv2(setup.p+1);
225  setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
226 
227  // evaluate bases and derivatives using formula:
228  // P_0(x) = 1
229  // P_1(x) = sqrt(3)*x
230  // P_i(x) = (x*P_{i-1}(x) - b[i-1]*P_{i-2}(x))/b[i], i=2,3,...
231  // where b[i] = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)))
232  // P_0'(x) = 0
233  // P_1'(x) = sqrt(3)
234  // P_i'(x) = (P_{i-1}(x) + x*P_{i-1}'(x) - b[i-1]*P_{i-2}')/b[i], i=2,3,...
235  double b1, b2;
236  b2 = std::sqrt(1.0/3.0);
237  b1 = b2;
238  val2[0] = 1.0;
239  if (setup.p >= 1)
240  val2[1] = x/b1;
241  for (int i=2; i<=setup.p; i++) {
242  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
243  val2[i] = (x*val2[i-1] - b2*val2[i-2])/b1;
244  b2 = b1;
245  }
246 
247  b2 = std::sqrt(1.0/3.0);
248  b1 = b2;
249  deriv2[0] = 0.0;
250  if (setup.p >= 1)
251  deriv2[1] = 1.0/b1;
252  for (int i=2; i<=setup.p; i++) {
253  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
254  deriv2[i] = (val2[i-1] + x*deriv2[i-1] - b2*deriv2[i-2])/b1;
255  b2 = b1;
256  }
257  success = Stokhos::compareArrays(val1, "val1", val2, "val2",
258  setup.rtol, setup.atol, out);
259  success = success &&
260  Stokhos::compareArrays(deriv1, "deriv1", deriv2, "deriv2",
261  setup.rtol, setup.atol, out);
262 
263 
264  }
265 
266  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Evaluate ) {
267  double x = 0.1234;
268  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
269  for (int i=0; i<=setup.p; i++)
270  v1[i] = setup.basis.evaluate(x, i);
271 
272  double b1, b2;
273  b2 = std::sqrt(1.0/3.0);
274  b1 = b2;
275  v2[0] = 1.0;
276  if (setup.p >= 1)
277  v2[1] = x/b1;
278  for (int i=2; i<=setup.p; i++) {
279  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
280  v2[i] = (x*v2[i-1] - b2*v2[i-2])/b1;
281  b2 = b1;
282  }
283  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
284  out);
285  }
286 
287  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Recurrence ) {
288  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
289  d1(setup.p+1);
290  Teuchos::Array<double> a2(setup.p+1), b2(setup.p+1), c2(setup.p+1),
291  d2(setup.p+1);
292  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
293 
294  // compute coefficients using formula
295  a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
296  for (int i=1; i<=setup.p; i++) {
297  a2[i] = 0.0;
298  b2[i] = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
299  c2[i] = 1.0;
300  d2[i] = b2[i];
301  }
302  success = true;
303  success = success &&
304  Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol, out);
305  success = success &&
306  Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol, out);
307  success = success &&
308  Stokhos::compareArrays(c1, "c1", c2, "c2", setup.rtol, setup.atol, out);
309  success = success &&
310  Stokhos::compareArrays(d1, "d1", d2, "d2", setup.rtol, setup.atol, out);
311  }
312 
313  // Tests alpha coefficients satisfy
314  // alpha_k = delta_k * (t*psi_k,psi_k)/(psi_k,psi_k)
315  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, RecurrenceAlpha ) {
316  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
317  d1(setup.p+1);
318  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
319 
320  Teuchos::Array<double> a2(setup.p+1);
321  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
322  Teuchos::Array<double> x, w;
323  Teuchos::Array< Teuchos::Array<double> > v;
324  setup.basis.getQuadPoints(2*setup.p, x, w, v);
325  int nqp = w.size();
326  for (int i=0; i<=setup.p; i++) {
327  a2[i] = 0;
328  for (int j=0; j<nqp; j++)
329  a2[i] += w[j]*x[j]*v[j][i]*v[j][i];
330  a2[i] = a2[i]*c1[i]/n1[i];
331  }
332  success = Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol,
333  out);
334  }
335 
336  // Tests beta coefficients satisfy
337  // beta_k =
338  // gamma_k * delta_k/delta_{k-1} * (psi_k,psi_k)/(psi_{k-1},psi_{k-1})
339  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, RecurrenceBeta ) {
340  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
341  d1(setup.p+1);
342  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
343 
344  Teuchos::Array<double> b2(setup.p+1);
345  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
346  b2[0] = b1[0];
347  for (int i=1; i<=setup.p; i++) {
348  b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
349  }
350  success = Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol,
351  out);
352  }
353 
354 #ifdef HAVE_STOKHOS_DAKOTA
355  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadPointsDakota ) {
356  int n = static_cast<int>(std::ceil((2*setup.p+1)/2.0));
357  Teuchos::Array<double> x1, w1;
358  Teuchos::Array< Teuchos::Array<double> > v1;
359  setup.basis.getQuadPoints(2*setup.p, x1, w1, v1);
360 
361  Teuchos::Array<double> x2(n), w2(n);
362  Teuchos::Array< Teuchos::Array<double> > v2(n);
363  webbur::legendre_compute(n, &x2[0], &w2[0]);
364 
365  for (int i=0; i<n; i++) {
366  w2[i] *= 0.5; // measure = 1/2
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_NormalizedLegendreBasis, 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 = 1;
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; // measure = 1/2
405  v2[i].resize(setup.p+1);
406  setup.basis.evaluateBases(x2[i], v2[i]);
407  }
408  success = true;
409  success = success &&
410  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
411  success = success &&
412  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
413  for (int i=0; i<n; i++) {
414  std::stringstream ss1, ss2;
415  ss1 << "v1[" << i << "]";
416  ss2 << "v2[" << i << "]";
417  success = success &&
418  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
419  setup.rtol, setup.atol, out);
420  }
421  }
422 #endif
423 
424 }
425 
426 int main( int argc, char* argv[] ) {
427  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
428  return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
429 }
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
LegendreBasisUnitTest::UnitTestSetup::basis
Stokhos::LegendreBasis< OrdinalType, ValueType > basis
Definition: Stokhos_LegendreBasisUnitTest.cpp:28
Teuchos_UnitTestHarness.hpp
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
Stokhos_UnitTestHelpers.hpp
LegendreBasisUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Stokhos_LegendreBasis, Order)
Definition: Stokhos_LegendreBasisUnitTest.cpp:36
LegendreBasisUnitTest::setup
UnitTestSetup< int, double > setup
Definition: Stokhos_LegendreBasisUnitTest.cpp:34
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
Teuchos_GlobalMPISession.hpp
Teuchos::Array::size
size_type size() const
true
true
Definition: Stokhos_LanczosUnitTest.cpp:233
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