Intrepid
test_01.cpp
Go to the documentation of this file.
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
58using namespace std;
59using 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 */
92template<class Scalar>
93Scalar 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
105template<class Scalar>
106Scalar * 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------------------------------------------------------------------*/
193int 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}
Header file for a set of functions providing orthogonal polynomial polynomial calculus and interpolat...
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin,...