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 // $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 
48 
49 #include "Stokhos.hpp"
51 #ifdef HAVE_STOKHOS_FORUQTK
52 #include "Stokhos_gaussq.h"
53 #endif
54 
55 namespace LegendreBasisUnitTest {
56 
57  // Common setup for unit tests
58  template <typename OrdinalType, typename ValueType>
59  struct UnitTestSetup {
60  ValueType rtol, atol;
61  OrdinalType p;
63 
64  UnitTestSetup() : rtol(1e-12), atol(1e-10), p(10), basis(p,true) {}
65 
66  };
67 
69 
70  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Order ) {
71  int order = setup.basis.order();
72  TEUCHOS_TEST_EQUALITY(order, setup.p, out, success);
73  }
74 
75  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Size ) {
76  int size = setup.basis.size();
77  TEUCHOS_TEST_EQUALITY(size, setup.p+1, out, success);
78  }
79 
80  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Norm ) {
81  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
83  for (int i=0; i<=setup.p; i++)
84  n2[i] = 1.0;
85  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
86  out);
87  }
88 
89  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Norm2 ) {
92  for (int i=0; i<=setup.p; i++) {
93  n1[i] = setup.basis.norm_squared(i);
94  n2[i] = 1.0;
95  }
96  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
97  out);
98  }
99 
100  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadNorm ) {
101  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
105  setup.basis.getQuadPoints(2*setup.p, x, w, v);
106  int nqp = w.size();
107  for (int i=0; i<=setup.p; i++) {
108  n2[i] = 0;
109  for (int j=0; j<nqp; j++)
110  n2[i] += w[j]*v[j][i]*v[j][i];
111  }
112  success = Stokhos::compareArrays(n1, "n1", n2, "n2", setup.rtol, setup.atol,
113  out);
114  }
115 
116  // Tests basis is orthogonal
117  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadOrthog ) {
118  const Teuchos::Array<double>& norms = setup.basis.norm_squared();
121  setup.basis.getQuadPoints(2*setup.p, x, w, v);
122  int nqp = w.size();
124  for (int i=0; i<=setup.p; i++) {
125  for (int j=0; j<=setup.p; j++) {
126  for (int k=0; k<nqp; k++)
127  mat(i,j) += w[k]*v[k][i]*v[k][j];
128  mat(i,j) /= std::sqrt(norms[i]*norms[j]);
129  }
130  mat(i,i) -= 1.0;
131  }
132  success = mat.normInf() < setup.atol;
133  if (!success) {
134  out << "\n Error, mat.normInf() < atol = " << mat.normInf()
135  << " < " << setup.atol << ": failed!\n";
136  out << "mat = " << printMat(mat) << std::endl;
137  }
138  }
139 
140  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, TripleProduct ) {
142  setup.basis.computeTripleProductTensor();
143 
146  setup.basis.getQuadPoints(3*setup.p, x, w, v);
147 
148  success = true;
149  for (int i=0; i<=setup.p; i++) {
150  for (int j=0; j<=setup.p; j++) {
151  for (int k=0; k<=setup.p; k++) {
152  double c = 0.0;
153  int nqp = w.size();
154  for (int qp=0; qp<nqp; qp++)
155  c += w[qp]*v[qp][i]*v[qp][j]*v[qp][k];
156  std::stringstream ss;
157  ss << "Cijk(" << i << "," << j << "," << k << ")";
158  success = success &&
159  Stokhos::compareValues((*Cijk)(i,j,k), ss.str(), c, "c",
160  setup.rtol, setup.atol, out);
161  }
162  }
163  }
164  }
165 
166  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, DerivDoubleProduct ) {
168  setup.basis.computeDerivDoubleProductTensor();
169 
172  setup.basis.getQuadPoints(3*setup.p, x, w, v);
173  int nqp = w.size();
174  val.resize(nqp);
175  deriv.resize(nqp);
176  for (int i=0; i<nqp; i++) {
177  val[i].resize(setup.p+1);
178  deriv[i].resize(setup.p+1);
179  setup.basis.evaluateBasesAndDerivatives(x[i], val[i], deriv[i]);
180  }
181 
182  success = true;
183  for (int i=0; i<=setup.p; i++) {
184  for (int j=0; j<=setup.p; j++) {
185  double b = 0.0;
186  for (int qp=0; qp<nqp; qp++)
187  b += w[qp]*deriv[qp][i]*val[qp][j];
188  std::stringstream ss;
189  ss << "Bij(" << i << "," << j << ")";
190  success = success &&
191  Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
192  setup.rtol, setup.atol, out);
193  }
194  }
195  }
196 
197  // TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, DerivDoubleProduct2 ) {
198  // Teuchos::RCP< Teuchos::SerialDenseMatrix<int, double> > Bij =
199  // setup.basis.computeDerivDoubleProductTensor();
200  // const Teuchos::Array<double>& n = setup.basis.norm_squared();
201 
202  // Teuchos::Array< Teuchos::Array<double> > deriv_coeffs(setup.p+1);
203  // deriv_coeffs[0].resize(setup.p+1,0.0);
204  // if (setup.p >= 1) {
205  // deriv_coeffs[1].resize(setup.p+1,0.0);
206  // deriv_coeffs[1][0] = 1.0;
207  // }
208  // if (setup.p >= 2) {
209  // deriv_coeffs[2].resize(setup.p+1,0.0);
210  // deriv_coeffs[2][1] = 3.0;
211  // }
212  // for (int k=3; k<=setup.p; k++) {
213  // deriv_coeffs[k].resize(setup.p+1,0.0);
214  // deriv_coeffs[k][0] = 1.0/3.0*deriv_coeffs[k-1][1];
215  // for (int i=1; i<=k-3; i++)
216  // deriv_coeffs[k][i] = i/(2.0*i - 1.0)*deriv_coeffs[k-1][i-1] +
217  // (i+1.0)/(2.0*i + 3.0)*deriv_coeffs[k-1][i+1];
218  // deriv_coeffs[k][k-2] = (k-2.0)/(2.0*k-5.0)*deriv_coeffs[k-1][k-3];
219  // deriv_coeffs[k][k-1] = (k-1.0)/(2.0*k-3.0)*deriv_coeffs[k-1][k-2] + k;
220  // }
221 
222  // success = true;
223  // for (int i=0; i<=setup.p; i++) {
224  // for (int j=0; j<=setup.p; j++) {
225  // double b = deriv_coeffs[i][j]*n[j];
226  // std::stringstream ss;
227  // ss << "Bij(" << i << "," << j << ")";
228  // success = success &&
229  // Stokhos::compareValues((*Bij)(i,j), ss.str(), b, "b",
230  // setup.rtol, setup.atol, out);
231  // }
232  // }
233  // }
234 
235  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, EvaluateBases ) {
236  double x = 0.1234;
237  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
238  setup.basis.evaluateBases(x, v1);
239 
240  double b1, b2;
241  b2 = std::sqrt(1.0/3.0);
242  b1 = b2;
243  v2[0] = 1.0;
244  if (setup.p >= 1)
245  v2[1] = x/b1;
246  for (int i=2; i<=setup.p; i++) {
247  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
248  v2[i] = (x*v2[i-1] - b2*v2[i-2]) / b1;
249  b2 = b1;
250  }
251  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
252  out);
253  }
254 
255  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, EvaluateBasesAndDerivatives ) {
256  double x = 0.1234;
257  Teuchos::Array<double> val1(setup.p+1), deriv1(setup.p+1),
258  val2(setup.p+1), deriv2(setup.p+1);
259  setup.basis.evaluateBasesAndDerivatives(x, val1, deriv1);
260 
261  // evaluate bases and derivatives using formula:
262  // P_0(x) = 1
263  // P_1(x) = sqrt(3)*x
264  // P_i(x) = (x*P_{i-1}(x) - b[i-1]*P_{i-2}(x))/b[i], i=2,3,...
265  // where b[i] = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)))
266  // P_0'(x) = 0
267  // P_1'(x) = sqrt(3)
268  // P_i'(x) = (P_{i-1}(x) + x*P_{i-1}'(x) - b[i-1]*P_{i-2}')/b[i], i=2,3,...
269  double b1, b2;
270  b2 = std::sqrt(1.0/3.0);
271  b1 = b2;
272  val2[0] = 1.0;
273  if (setup.p >= 1)
274  val2[1] = x/b1;
275  for (int i=2; i<=setup.p; i++) {
276  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
277  val2[i] = (x*val2[i-1] - b2*val2[i-2])/b1;
278  b2 = b1;
279  }
280 
281  b2 = std::sqrt(1.0/3.0);
282  b1 = b2;
283  deriv2[0] = 0.0;
284  if (setup.p >= 1)
285  deriv2[1] = 1.0/b1;
286  for (int i=2; i<=setup.p; i++) {
287  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
288  deriv2[i] = (val2[i-1] + x*deriv2[i-1] - b2*deriv2[i-2])/b1;
289  b2 = b1;
290  }
291  success = Stokhos::compareArrays(val1, "val1", val2, "val2",
292  setup.rtol, setup.atol, out);
293  success = success &&
294  Stokhos::compareArrays(deriv1, "deriv1", deriv2, "deriv2",
295  setup.rtol, setup.atol, out);
296 
297 
298  }
299 
300  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Evaluate ) {
301  double x = 0.1234;
302  Teuchos::Array<double> v1(setup.p+1), v2(setup.p+1);
303  for (int i=0; i<=setup.p; i++)
304  v1[i] = setup.basis.evaluate(x, i);
305 
306  double b1, b2;
307  b2 = std::sqrt(1.0/3.0);
308  b1 = b2;
309  v2[0] = 1.0;
310  if (setup.p >= 1)
311  v2[1] = x/b1;
312  for (int i=2; i<=setup.p; i++) {
313  b1 = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
314  v2[i] = (x*v2[i-1] - b2*v2[i-2])/b1;
315  b2 = b1;
316  }
317  success = Stokhos::compareArrays(v1, "v1", v2, "v2", setup.rtol, setup.atol,
318  out);
319  }
320 
321  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, Recurrence ) {
322  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
323  d1(setup.p+1);
324  Teuchos::Array<double> a2(setup.p+1), b2(setup.p+1), c2(setup.p+1),
325  d2(setup.p+1);
326  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
327 
328  // compute coefficients using formula
329  a2[0] = 0.0; b2[0] = 1.0; c2[0] = 1.0; d2[0] = 1.0;
330  for (int i=1; i<=setup.p; i++) {
331  a2[i] = 0.0;
332  b2[i] = std::sqrt((i*i)/((2.0*i+1.0)*(2.0*i-1.0)));
333  c2[i] = 1.0;
334  d2[i] = b2[i];
335  }
336  success = true;
337  success = success &&
338  Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol, out);
339  success = success &&
340  Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol, out);
341  success = success &&
342  Stokhos::compareArrays(c1, "c1", c2, "c2", setup.rtol, setup.atol, out);
343  success = success &&
344  Stokhos::compareArrays(d1, "d1", d2, "d2", setup.rtol, setup.atol, out);
345  }
346 
347  // Tests alpha coefficients satisfy
348  // alpha_k = delta_k * (t*psi_k,psi_k)/(psi_k,psi_k)
349  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, RecurrenceAlpha ) {
350  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
351  d1(setup.p+1);
352  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
353 
355  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
358  setup.basis.getQuadPoints(2*setup.p, x, w, v);
359  int nqp = w.size();
360  for (int i=0; i<=setup.p; i++) {
361  a2[i] = 0;
362  for (int j=0; j<nqp; j++)
363  a2[i] += w[j]*x[j]*v[j][i]*v[j][i];
364  a2[i] = a2[i]*c1[i]/n1[i];
365  }
366  success = Stokhos::compareArrays(a1, "a1", a2, "a2", setup.rtol, setup.atol,
367  out);
368  }
369 
370  // Tests beta coefficients satisfy
371  // beta_k =
372  // gamma_k * delta_k/delta_{k-1} * (psi_k,psi_k)/(psi_{k-1},psi_{k-1})
373  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, RecurrenceBeta ) {
374  Teuchos::Array<double> a1(setup.p+1), b1(setup.p+1), c1(setup.p+1),
375  d1(setup.p+1);
376  setup.basis.getRecurrenceCoefficients(a1, b1, c1, d1);
377 
379  const Teuchos::Array<double>& n1 = setup.basis.norm_squared();
380  b2[0] = b1[0];
381  for (int i=1; i<=setup.p; i++) {
382  b2[i] = d1[i]*(c1[i]/c1[i-1])*(n1[i]/n1[i-1]);
383  }
384  success = Stokhos::compareArrays(b1, "b1", b2, "b2", setup.rtol, setup.atol,
385  out);
386  }
387 
388 #ifdef HAVE_STOKHOS_DAKOTA
389  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadPointsDakota ) {
390  int n = static_cast<int>(std::ceil((2*setup.p+1)/2.0));
391  Teuchos::Array<double> x1, w1;
393  setup.basis.getQuadPoints(2*setup.p, x1, w1, v1);
394 
395  Teuchos::Array<double> x2(n), w2(n);
397  webbur::legendre_compute(n, &x2[0], &w2[0]);
398 
399  for (int i=0; i<n; i++) {
400  w2[i] *= 0.5; // measure = 1/2
401  v2[i].resize(setup.p+1);
402  setup.basis.evaluateBases(x2[i], v2[i]);
403  }
404  success = true;
405  success = success &&
406  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
407  success = success &&
408  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
409  for (int i=0; i<n; i++) {
410  std::stringstream ss1, ss2;
411  ss1 << "v1[" << i << "]";
412  ss2 << "v2[" << i << "]";
413  success = success &&
414  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
415  setup.rtol, setup.atol, out);
416  }
417  }
418 #endif
419 
420 #ifdef HAVE_STOKHOS_FORUQTK
421  TEUCHOS_UNIT_TEST( Stokhos_NormalizedLegendreBasis, QuadPointsForUQTK ) {
422  int n = static_cast<int>(std::ceil((setup.p+1)/2.0));
423  Teuchos::Array<double> x1, w1;
425  setup.basis.getQuadPoints(setup.p, x1, w1, v1);
426 
427  Teuchos::Array<double> x2(n), w2(n);
429  int kind = 1;
430  int kpts = 0;
431  double endpts[2] = {0.0, 0.0};
433  double alpha = 0.0;
434  double beta = 0.0;
435  GAUSSQ_F77(&kind, &n, &alpha, &beta, &kpts, endpts, &b[0], &x2[0], &w2[0]);
436 
437  for (int i=0; i<n; i++) {
438  w2[i] *= 0.5; // measure = 1/2
439  v2[i].resize(setup.p+1);
440  setup.basis.evaluateBases(x2[i], v2[i]);
441  }
442  success = true;
443  success = success &&
444  Stokhos::compareArrays(x1, "x1", x2, "x2", setup.rtol, setup.atol, out);
445  success = success &&
446  Stokhos::compareArrays(w1, "w1", w2, "w2", setup.rtol, setup.atol, out);
447  for (int i=0; i<n; i++) {
448  std::stringstream ss1, ss2;
449  ss1 << "v1[" << i << "]";
450  ss2 << "v2[" << i << "]";
451  success = success &&
452  Stokhos::compareArrays(v1[i], ss1.str(), v2[i], ss2.str(),
453  setup.rtol, setup.atol, out);
454  }
455  }
456 #endif
457 
458 }
459 
460 int main( int argc, char* argv[] ) {
461  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
463 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
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 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)
static int runUnitTestsFromMain(int argc, char *argv[])
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
Stokhos::LegendreBasis< OrdinalType, ValueType > basis
void resize(size_type new_size, const value_type &x=value_type())
TEUCHOS_UNIT_TEST(Stokhos_LegendreBasis, Order)
UnitTestSetup< int, double > setup
int main(int argc, char **argv)
ScalarTraits< ScalarType >::magnitudeType normInf() const
expr val()
#define TEUCHOS_TEST_EQUALITY(v1, v2, out, success)
size_type size() const
int n
void GAUSSQ_F77(int *kind, int *n, double *alpha, double *beta, int *kpts, double *endpts, double *b, double *x, double *w)