Intrepid
test_01.cpp
1 // @HEADER
3 //
4 // File: test_01.cpp
5 //
6 // For more information, please see: http://www.nektar.info
7 //
8 // The MIT License
9 //
10 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11 // Department of Aeronautics, Imperial College London (UK), and Scientific
12 // Computing and Imaging Institute, University of Utah (USA).
13 //
14 // License for the specific language governing rights and limitations under
15 // Permission is hereby granted, free of charge, to any person obtaining a
16 // copy of this software and associated documentation files (the "Software"),
17 // to deal in the Software without restriction, including without limitation
18 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
19 // and/or sell copies of the Software, and to permit persons to whom the
20 // Software is furnished to do so, subject to the following conditions:
21 //
22 // The above copyright notice and this permission notice shall be included
23 // in all copies or substantial portions of the Software.
24 //
25 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
26 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
27 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
28 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
29 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
30 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
31 // DEALINGS IN THE SOFTWARE.
32 //
33 // Description:
34 // This file is redistributed with the Intrepid package. It should be used
35 // in accordance with the above MIT license, at the request of the authors.
36 // This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
37 //
38 // Origin: Nektar++ library, http://www.nektar.info, downloaded on
39 // March 10, 2009.
40 //
42 
43 
51 #include "Intrepid_Polylib.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 
57 
58 using namespace std;
59 using namespace Intrepid;
60 
61 #define NPLOWER 5
62 #define NPUPPER 15
63 #define TEST_EPS 1000*INTREPID_TOL
64 
65 #define GAUSS_INT 1
66 #define GAUSS_RADAUM_INT 1
67 #define GAUSS_RADAUP_INT 1
68 #define GAUSS_LOBATTO_INT 1
69 #define GAUSS_DIFF 1
70 #define GAUSS_RADAUM_DIFF 1
71 #define GAUSS_RADAUP_DIFF 1
72 #define GAUSS_LOBATTO_DIFF 1
73 #define GAUSS_INTERP 1
74 #define GAUSS_RADAUM_INTERP 1
75 #define GAUSS_RADAUP_INTERP 1
76 #define GAUSS_LOBATTO_INTERP 1
77 
78 
79 #define INTREPID_TEST_COMMAND( S ) \
80 { \
81  try { \
82  S ; \
83  } \
84  catch (const std::logic_error & err) { \
85  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
86  *outStream << err.what() << '\n'; \
87  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
88  }; \
89 }
90 
91 /* local routines */
92 template<class Scalar>
93 Scalar ddot(int n, Scalar *x, int incx, Scalar *y, int incy)
94 {
95  Scalar sum = 0.;
96 
97  while (n--) {
98  sum += (*x) * (*y);
99  x += incx;
100  y += incy;
101  }
102  return sum;
103 }
104 
105 template<class Scalar>
106 Scalar * dvector(int nl,int nh)
107 {
108  Scalar * v;
109 
110  v = (Scalar *)malloc((unsigned) (nh-nl+1)*sizeof(Scalar));
111  return v-nl;
112 }
113 
114 
115 /* -------------------------------------------------------------------
116  This is a routine to test the integration, differentiation and
117  interpolation routines in the polylib.c.
118 
119  First, it performs the integral
120 
121  /1 alpha beta alpha,beta
122  | (1-x) (1+x) P (x) dx = 0
123  /-1 n
124 
125  for all -0.5 <= alpha <= 5 (increments of 0.5)
126  -0.5 <= beta <= 5 (increments of 0.5)
127 
128  using np points where
129  NPLOWER <= np <= NPUPPER
130  2 <= n <= 2*np - delta
131 
132  delta = 1 (gauss), 2(radau), 3(lobatto).
133  The integral is evaluated and if it is larger that EPS then the
134  value of alpha,beta,np,n and the integral is printed to the screen.
135 
136  After every alpha value the statement
137  "finished checking all beta values for alpha = #"
138  is printed
139 
140  The routine then evaluates the derivate of
141 
142  d n n-1
143  -- x = n x
144  dx
145 
146  for all -0.5 <= alpha <= 5 (increments of 0.5)
147  -0.5 <= beta <= 5 (increments of 0.5)
148 
149  using np points where
150  NPLOWER <= np <= NPUPPER
151  2 <= n <= np - 1
152 
153  The error is check in a pointwise sense and if it is larger than
154  EPS then the value of alpha,beta,np,n and the error is printed to
155  the screen. After every alpha value the statement
156  "finished checking all beta values for alpha = #"
157  is printed
158 
159  Finally the routine evaluates the interpolation of
160 
161  n n
162  z to x
163 
164  where z are the quadrature zeros and x are the equispaced points
165 
166  2*i
167  x = ----- - 1.0 (0 <= i <= np-1)
168  i (np-1)
169 
170 
171  for all -0.5 <= alpha <= 5 (increments of 0.5)
172  -0.5 <= beta <= 5 (increments of 0.5)
173 
174  using np points where
175  NPLOWER <= np <= NPUPPER
176  2 <= n <= np - 1
177 
178  The error is check in a pointwise sense and if it is larger than
179  EPS then the value of alpha,beta,np,n and the error is printed to
180  the screen. After every alpha value the statement
181  "finished checking all beta values for alpha = #"
182  is printed
183 
184  The above checks are performed for all the Gauss, Gauss-Radau and
185  Gauss-Lobatto points. If you want to disable any routine then set
186  GAUSS_INT, GAUSS_RADAU_INT, GAUSS_LOBATTO_INT = 0
187  for the integration rouintes
188  GAUSS_DIFF,GAUSS_RADAU_DIFF, GAUSS_LOBATTO_DIFF = 0
189  for the differentiation routines
190  GAUSS_INTERP,GAUSS_RADAU_INTERP, GAUSS_LOBATTO_INTERP = 0
191  for the interpolation routines.
192 ------------------------------------------------------------------*/
193 int main(int argc, char *argv[]) {
194 
195  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
196 
197  // This little trick lets us print to std::cout only if
198  // a (dummy) command-line argument is provided.
199  int iprint = argc - 1;
200  Teuchos::RCP<std::ostream> outStream;
201  Teuchos::oblackholestream bhs; // outputs nothing
202  if (iprint > 0)
203  outStream = Teuchos::rcp(&std::cout, false);
204  else
205  outStream = Teuchos::rcp(&bhs, false);
206 
207  // Save the format state of the original std::cout.
208  Teuchos::oblackholestream oldFormatState;
209  oldFormatState.copyfmt(std::cout);
210 
211  *outStream \
212  << "===============================================================================\n" \
213  << "| |\n" \
214  << "| Unit Test (IntrepidPolylib) |\n" \
215  << "| |\n" \
216  << "| 1) the original Polylib tests |\n" \
217  << "| |\n" \
218  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
219  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
220  << "| |\n" \
221  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
222  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
223  << "| |\n" \
224  << "===============================================================================\n";
225 
226 
227  int errorFlag = 0;
228  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
229  int endThrowNumber = beginThrowNumber + 1;
230 
231  typedef IntrepidPolylib ipl;
232  IntrepidPolylib iplib;
233 
234  *outStream \
235  << "\n"
236  << "===============================================================================\n"\
237  << "| TEST 1: exceptions |\n"\
238  << "===============================================================================\n";
239 
240  try{
241  INTREPID_TEST_COMMAND( iplib.gammaF((double)0.33) );
242  }
243  catch (const std::logic_error & err) {
244  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
245  *outStream << err.what() << '\n';
246  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
247  errorFlag = -1000;
248  };
249 
250  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
251  errorFlag = -1000;
252 
253  *outStream \
254  << "\n"
255  << "===============================================================================\n"\
256  << "| TEST 2: correctness of math operations |\n"\
257  << "===============================================================================\n";
258 
259  outStream->precision(20);
260 
261  try {
262  { // start scope
263  int np,n,i;
264  double *z,*w,*p,sum=0,alpha,beta,*d;
265 
266  z = dvector<double>(0,NPUPPER-1);
267  w = dvector<double>(0,NPUPPER-1);
268  p = dvector<double>(0,NPUPPER-1);
269 
270  d = dvector<double>(0,NPUPPER*NPUPPER-1);
271 
272 #if GAUSS_INT
273  /* Gauss Integration */
274  *outStream << "Begin checking Gauss integration\n";
275  alpha = -0.5;
276  while(alpha <= 5.0){
277  beta = -0.5;
278  while(beta <= 5.0){
279 
280  for(np = NPLOWER; np <= NPUPPER; ++np){
281  ipl::zwgj(z,w,np,alpha,beta);
282  for(n = 2; n < 2*np-1; ++n){
283  ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta);
284  sum = ddot(np,w,1,p,1);
285  if(std::abs(sum)>TEST_EPS) {
286  errorFlag = -1000;
287  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
288  ", np = " << np << ", n = " << n << " integral was " << sum << "\n";
289  }
290  }
291  }
292 
293  beta += 0.5;
294  }
295  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
296  alpha += 0.5;
297  }
298  *outStream << "Finished checking Gauss Integration\n";
299 #endif
300 
301 #if GAUSS_RADAUM_INT
302  /* Gauss Radau (z = -1) Integration */
303  *outStream << "Begin checking Gauss-Radau (z = -1) integration\n";
304  alpha = -0.5;
305  while(alpha <= 5.0){
306  beta = -0.5;
307  while(beta <= 5.0){
308 
309  for(np = NPLOWER; np <= NPUPPER; ++np){
310  ipl::zwgrjm(z,w,np,alpha,beta);
311  for(n = 2; n < 2*np-2; ++n){
312  ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta);
313  sum = ddot(np,w,1,p,1);
314  if(std::abs(sum)>TEST_EPS) {
315  errorFlag = -1000;
316  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
317  ", np = " << np << ", n = " << n << " integral was " << sum << "\n";
318  }
319  }
320  }
321 
322  beta += 0.5;
323  }
324  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
325  alpha += 0.5;
326  }
327  *outStream << "Finished checking Gauss-Radau (z = -1) Integration\n";
328 #endif
329 
330 #if GAUSS_RADAUP_INT
331  /* Gauss Radau (z = +1) Integration */
332  *outStream << "Begin checking Gauss-Radau (z = +1) integration\n";
333  alpha = -0.5;
334  while(alpha <= 5.0){
335  beta = -0.5;
336  while(beta <= 5.0){
337 
338  for(np = NPLOWER; np <= NPUPPER; ++np){
339  ipl::zwgrjp(z,w,np,alpha,beta);
340  for(n = 2; n < 2*np-2; ++n){
341  ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta);
342  sum = ddot(np,w,1,p,1);
343  if(std::abs(sum)>TEST_EPS) {
344  errorFlag = -1000;
345  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
346  ", np = " << np << ", n = " << n << " integral was " << sum << "\n";
347  }
348  }
349  }
350 
351  beta += 0.5;
352  }
353  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
354  alpha += 0.5;
355  }
356  *outStream << "Finished checking Gauss-Radau (z = +1) Integration\n";
357 #endif
358 
359 #if GAUSS_LOBATTO_INT
360  /* Gauss Lobatto Integration */
361  *outStream << "Begin checking Gauss-Lobatto integration\n";
362  alpha = -0.5;
363  while(alpha <= 5.0){
364  beta = -0.5;
365  while(beta <= 5.0){
366 
367  for(np = NPLOWER; np <= NPUPPER; ++np){
368  ipl::zwglj(z,w,np,alpha,beta);
369  for(n = 2; n < 2*np-3; ++n){
370  ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta);
371  sum = ddot(np,w,1,p,1);
372  if(std::abs(sum)>TEST_EPS) {
373  errorFlag = -1000;
374  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
375  ", np = " << np << ", n = " << n << " integral was " << sum << "\n";
376  }
377  }
378  }
379 
380  beta += 0.5;
381  }
382  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
383  alpha += 0.5;
384  }
385  *outStream << "Finished checking Gauss-Lobatto Integration\n";
386 #endif
387 
388 #if GAUSS_DIFF
389  *outStream << "Begin checking differentiation through Gauss points\n";
390  alpha = -0.5;
391  while(alpha <= 5.0){
392  beta = -0.5;
393  while(beta <= 5.0){
394 
395  for(np = NPLOWER; np <= NPUPPER; ++np){
396  ipl::zwgj(z,w,np,alpha,beta);
397  for(n = 2; n < np-1; ++n){
398  ipl::Dgj(d,z,np,alpha,beta);
399 
400  for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
401  sum = 0;
402  for(i = 0; i < np; ++i)
403  sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
404  sum /= np;
405  if(std::abs(sum)>TEST_EPS) {
406  errorFlag = -1000;
407  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
408  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
409  }
410  }
411  }
412  beta += 0.5;
413  }
414  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
415  alpha += 0.5;
416  }
417  *outStream << "Finished checking Gauss Jacobi differentiation\n";
418 #endif
419 
420 #if GAUSS_RADAUM_DIFF
421  *outStream << "Begin checking differentiation through Gauss-Radau (z=-1) points\n";
422  alpha = -0.5;
423  while(alpha <= 5.0){
424  beta = -0.5;
425  while(beta <= 5.0){
426 
427  for(np = NPLOWER; np <= NPUPPER; ++np){
428  ipl::zwgrjm(z,w,np,alpha,beta);
429  for(n = 2; n < np-1; ++n){
430  ipl::Dgrjm(d,z,np,alpha,beta);
431 
432  for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
433  sum = 0;
434  for(i = 0; i < np; ++i)
435  sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
436  sum /= np;
437  if(std::abs(sum)>TEST_EPS) {
438  errorFlag = -1000;
439  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
440  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
441  }
442  }
443  }
444  beta += 0.5;
445  }
446  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
447  alpha += 0.5;
448  }
449  *outStream << "Finished checking Gauss-Radau (z=-1) differentiation\n";
450 #endif
451 
452 #if GAUSS_RADAUP_DIFF
453  *outStream << "Begin checking differentiation through Gauss-Radau (z=+1) points\n";
454  alpha = -0.5;
455  while(alpha <= 5.0){
456  beta = -0.5;
457  while(beta <= 5.0){
458 
459  for(np = NPLOWER; np <= NPUPPER; ++np){
460  ipl::zwgrjp(z,w,np,alpha,beta);
461  for(n = 2; n < np-1; ++n){
462  ipl::Dgrjp(d,z,np,alpha,beta);
463 
464  for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
465  sum = 0;
466  for(i = 0; i < np; ++i)
467  sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
468  sum /= np;
469  if(std::abs(sum)>TEST_EPS) {
470  errorFlag = -1000;
471  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
472  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
473  }
474  }
475  }
476  beta += 0.5;
477  }
478  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
479  alpha += 0.5;
480  }
481  *outStream << "Finished checking Gauss-Radau (z=+1) differentiation\n";
482 #endif
483 
484 #if GAUSS_RADAUP_DIFF
485  *outStream << "Begin checking differentiation through Gauss-Lobatto (z=+1) points\n";
486  alpha = -0.5;
487  while(alpha <= 5.0){
488  beta = -0.5;
489  while(beta <= 5.0){
490 
491  for(np = NPLOWER; np <= NPUPPER; ++np){
492  ipl::zwglj(z,w,np,alpha,beta);
493  for(n = 2; n < np-1; ++n){
494  ipl::Dglj(d,z,np,alpha,beta);
495 
496  for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
497  sum = 0;
498  for(i = 0; i < np; ++i)
499  sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
500  sum /= np;
501  if(std::abs(sum)>TEST_EPS) {
502  errorFlag = -1000;
503  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
504  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
505  }
506  }
507  }
508  beta += 0.5;
509  }
510  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
511  alpha += 0.5;
512  }
513  *outStream << "Finished checking Gauss-Lobatto differentiation\n";
514 #endif
515 
516 #if GAUSS_INTERP
517  *outStream << "Begin checking interpolation through Gauss points\n";
518  alpha = -0.5;
519  while(alpha <= 5.0){
520  beta = -0.5;
521  while(beta <= 5.0){
522 
523  for(np = NPLOWER; np <= NPUPPER; ++np){
524  ipl::zwgj(z,w,np,alpha,beta);
525  for(n = 2; n < np-1; ++n){
526  for(i = 0; i < np; ++i) {
527  w[i] = 2.0*i/(double)(np-1)-1.0;
528  p[i] = std::pow(z[i],n);
529  }
530  ipl::Imgj(d,z,w,np,np,alpha,beta);
531  sum = 0;
532  for(i = 0; i < np; ++i)
533  sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
534  sum /= np;
535  if(std::abs(sum)>TEST_EPS) {
536  errorFlag = -1000;
537  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
538  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
539  }
540  }
541  }
542  beta += 0.5;
543  }
544  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
545  alpha += 0.5;
546  }
547  *outStream << "Finished checking Gauss Jacobi interpolation\n";
548 #endif
549 
550 #if GAUSS_RADAUM_INTERP
551  *outStream << "Begin checking interpolation through Gauss-Radau (z=-1) points\n";
552  alpha = -0.5;
553  while(alpha <= 5.0){
554  beta = -0.5;
555  while(beta <= 5.0){
556 
557  for(np = NPLOWER; np <= NPUPPER; ++np){
558  ipl::zwgrjm(z,w,np,alpha,beta);
559  for(n = 2; n < np-1; ++n){
560  for(i = 0; i < np; ++i) {
561  w[i] = 2.0*i/(double)(np-1)-1.0;
562  p[i] = std::pow(z[i],n);
563  }
564  ipl::Imgrjm(d,z,w,np,np,alpha,beta);
565  sum = 0;
566  for(i = 0; i < np; ++i)
567  sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
568  sum /= np;
569  if(std::abs(sum)>TEST_EPS) {
570  errorFlag = -1000;
571  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
572  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
573  }
574  }
575  }
576  beta += 0.5;
577  }
578  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
579  alpha += 0.5;
580  }
581  *outStream << "Finished checking Gauss-Radau (z=-1) interpolation\n";
582 #endif
583 
584 #if GAUSS_RADAUP_INTERP
585  *outStream << "Begin checking interpolation through Gauss-Radau (z=+1) points\n";
586  alpha = -0.5;
587  while(alpha <= 5.0){
588  beta = -0.5;
589  while(beta <= 5.0){
590 
591  for(np = NPLOWER; np <= NPUPPER; ++np){
592  ipl::zwgrjp(z,w,np,alpha,beta);
593  for(n = 2; n < np-1; ++n){
594  for(i = 0; i < np; ++i) {
595  w[i] = 2.0*i/(double)(np-1)-1.0;
596  p[i] = std::pow(z[i],n);
597  }
598  ipl::Imgrjp(d,z,w,np,np,alpha,beta);
599  sum = 0;
600  for(i = 0; i < np; ++i)
601  sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
602  sum /= np;
603  if(std::abs(sum)>TEST_EPS) {
604  errorFlag = -1000;
605  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
606  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
607  }
608  }
609  }
610  beta += 0.5;
611  }
612  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
613  alpha += 0.5;
614  }
615  *outStream << "Finished checking Gauss-Radau (z=+1) interpolation\n";
616 #endif
617 
618 #if GAUSS_LOBATTO_INTERP
619  *outStream << "Begin checking interpolation through Gauss-Lobatto points\n";
620  alpha = -0.5;
621  while(alpha <= 5.0){
622  beta = -0.5;
623  while(beta <= 5.0){
624 
625  for(np = NPLOWER; np <= NPUPPER; ++np){
626  ipl::zwglj(z,w,np,alpha,beta);
627  for(n = 2; n < np-1; ++n){
628  for(i = 0; i < np; ++i) {
629  w[i] = 2.0*i/(double)(np-1)-1.0;
630  p[i] = std::pow(z[i],n);
631  }
632  ipl::Imglj(d,z,w,np,np,alpha,beta);
633  sum = 0;
634  for(i = 0; i < np; ++i)
635  sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
636  sum /= np;
637  if(std::abs(sum)>TEST_EPS) {
638  errorFlag = -1000;
639  *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta <<
640  ", np = " << np << ", n = " << n << " difference " << sum << "\n";
641  }
642  }
643  }
644  beta += 0.5;
645  }
646  *outStream << "finished checking all beta values for alpha = " << alpha << "\n";
647  alpha += 0.5;
648  }
649  *outStream << "Finished checking Gauss-Lobatto interpolation\n";
650 #endif
651 
652  free(z); free(w); free(p); free(d);
653 
654  } // end scope
655 
656  /******************************************/
657  *outStream << "\n";
658  }
659  catch (const std::logic_error & err) {
660  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
661  *outStream << err.what() << '\n';
662  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
663  errorFlag = -1000;
664  };
665 
666 
667  if (errorFlag != 0)
668  std::cout << "End Result: TEST FAILED\n";
669  else
670  std::cout << "End Result: TEST PASSED\n";
671 
672  // reset format state of std::cout
673  std::cout.copyfmt(oldFormatState);
674 
675  return errorFlag;
676 }
Intrepid::IntrepidPolylib
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin, Aeronautics, Imperial College London, modified and redistributed by D. Ridzal.
Definition: Intrepid_Polylib.hpp:223
Intrepid::IntrepidPolylib::gammaF
static Scalar gammaF(const Scalar x)
Calculate the Gamma function , , for integer values x and halves.
Definition: Intrepid_PolylibDef.hpp:814
Intrepid_Polylib.hpp
Header file for a set of functions providing orthogonal polynomial polynomial calculus and interpolat...
Generated by   doxygen 1.8.5