|
Intrepid
|
00001 // @HEADER 00003 // 00004 // File: test_01.cpp 00005 // 00006 // For more information, please see: http://www.nektar.info 00007 // 00008 // The MIT License 00009 // 00010 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), 00011 // Department of Aeronautics, Imperial College London (UK), and Scientific 00012 // Computing and Imaging Institute, University of Utah (USA). 00013 // 00014 // License for the specific language governing rights and limitations under 00015 // Permission is hereby granted, free of charge, to any person obtaining a 00016 // copy of this software and associated documentation files (the "Software"), 00017 // to deal in the Software without restriction, including without limitation 00018 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 00019 // and/or sell copies of the Software, and to permit persons to whom the 00020 // Software is furnished to do so, subject to the following conditions: 00021 // 00022 // The above copyright notice and this permission notice shall be included 00023 // in all copies or substantial portions of the Software. 00024 // 00025 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 00026 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 00027 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 00028 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 00029 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 00030 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 00031 // DEALINGS IN THE SOFTWARE. 00032 // 00033 // Description: 00034 // This file is redistributed with the Intrepid package. It should be used 00035 // in accordance with the above MIT license, at the request of the authors. 00036 // This file is NOT covered by the usual Intrepid/Trilinos LGPL license. 00037 // 00038 // Origin: Nektar++ library, http://www.nektar.info, downloaded on 00039 // March 10, 2009. 00040 // 00042 00043 00051 #include "Intrepid_Polylib.hpp" 00052 #include "Teuchos_oblackholestream.hpp" 00053 #include "Teuchos_RCP.hpp" 00054 #include "Teuchos_ScalarTraits.hpp" 00055 #include "Teuchos_GlobalMPISession.hpp" 00056 00057 00058 using namespace std; 00059 using namespace Intrepid; 00060 00061 #define NPLOWER 5 00062 #define NPUPPER 15 00063 #define TEST_EPS 1000*INTREPID_TOL 00064 00065 #define GAUSS_INT 1 00066 #define GAUSS_RADAUM_INT 1 00067 #define GAUSS_RADAUP_INT 1 00068 #define GAUSS_LOBATTO_INT 1 00069 #define GAUSS_DIFF 1 00070 #define GAUSS_RADAUM_DIFF 1 00071 #define GAUSS_RADAUP_DIFF 1 00072 #define GAUSS_LOBATTO_DIFF 1 00073 #define GAUSS_INTERP 1 00074 #define GAUSS_RADAUM_INTERP 1 00075 #define GAUSS_RADAUP_INTERP 1 00076 #define GAUSS_LOBATTO_INTERP 1 00077 00078 00079 #define INTREPID_TEST_COMMAND( S ) \ 00080 { \ 00081 try { \ 00082 S ; \ 00083 } \ 00084 catch (std::logic_error err) { \ 00085 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00086 *outStream << err.what() << '\n'; \ 00087 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00088 }; \ 00089 } 00090 00091 /* local routines */ 00092 template<class Scalar> 00093 Scalar ddot(int n, Scalar *x, int incx, Scalar *y, int incy) 00094 { 00095 register Scalar sum = 0.; 00096 00097 while (n--) { 00098 sum += (*x) * (*y); 00099 x += incx; 00100 y += incy; 00101 } 00102 return sum; 00103 } 00104 00105 template<class Scalar> 00106 Scalar * dvector(int nl,int nh) 00107 { 00108 Scalar * v; 00109 00110 v = (Scalar *)malloc((unsigned) (nh-nl+1)*sizeof(Scalar)); 00111 return v-nl; 00112 } 00113 00114 00115 /* ------------------------------------------------------------------- 00116 This is a routine to test the integration, differentiation and 00117 interpolation routines in the polylib.c. 00118 00119 First, it performs the integral 00120 00121 /1 alpha beta alpha,beta 00122 | (1-x) (1+x) P (x) dx = 0 00123 /-1 n 00124 00125 for all -0.5 <= alpha <= 5 (increments of 0.5) 00126 -0.5 <= beta <= 5 (increments of 0.5) 00127 00128 using np points where 00129 NPLOWER <= np <= NPUPPER 00130 2 <= n <= 2*np - delta 00131 00132 delta = 1 (gauss), 2(radau), 3(lobatto). 00133 The integral is evaluated and if it is larger that EPS then the 00134 value of alpha,beta,np,n and the integral is printed to the screen. 00135 00136 After every alpha value the statement 00137 "finished checking all beta values for alpha = #" 00138 is printed 00139 00140 The routine then evaluates the derivate of 00141 00142 d n n-1 00143 -- x = n x 00144 dx 00145 00146 for all -0.5 <= alpha <= 5 (increments of 0.5) 00147 -0.5 <= beta <= 5 (increments of 0.5) 00148 00149 using np points where 00150 NPLOWER <= np <= NPUPPER 00151 2 <= n <= np - 1 00152 00153 The error is check in a pointwise sense and if it is larger than 00154 EPS then the value of alpha,beta,np,n and the error is printed to 00155 the screen. After every alpha value the statement 00156 "finished checking all beta values for alpha = #" 00157 is printed 00158 00159 Finally the routine evaluates the interpolation of 00160 00161 n n 00162 z to x 00163 00164 where z are the quadrature zeros and x are the equispaced points 00165 00166 2*i 00167 x = ----- - 1.0 (0 <= i <= np-1) 00168 i (np-1) 00169 00170 00171 for all -0.5 <= alpha <= 5 (increments of 0.5) 00172 -0.5 <= beta <= 5 (increments of 0.5) 00173 00174 using np points where 00175 NPLOWER <= np <= NPUPPER 00176 2 <= n <= np - 1 00177 00178 The error is check in a pointwise sense and if it is larger than 00179 EPS then the value of alpha,beta,np,n and the error is printed to 00180 the screen. After every alpha value the statement 00181 "finished checking all beta values for alpha = #" 00182 is printed 00183 00184 The above checks are performed for all the Gauss, Gauss-Radau and 00185 Gauss-Lobatto points. If you want to disable any routine then set 00186 GAUSS_INT, GAUSS_RADAU_INT, GAUSS_LOBATTO_INT = 0 00187 for the integration rouintes 00188 GAUSS_DIFF,GAUSS_RADAU_DIFF, GAUSS_LOBATTO_DIFF = 0 00189 for the differentiation routines 00190 GAUSS_INTERP,GAUSS_RADAU_INTERP, GAUSS_LOBATTO_INTERP = 0 00191 for the interpolation routines. 00192 ------------------------------------------------------------------*/ 00193 int main(int argc, char *argv[]) { 00194 00195 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00196 00197 // This little trick lets us print to std::cout only if 00198 // a (dummy) command-line argument is provided. 00199 int iprint = argc - 1; 00200 Teuchos::RCP<std::ostream> outStream; 00201 Teuchos::oblackholestream bhs; // outputs nothing 00202 if (iprint > 0) 00203 outStream = Teuchos::rcp(&std::cout, false); 00204 else 00205 outStream = Teuchos::rcp(&bhs, false); 00206 00207 // Save the format state of the original std::cout. 00208 Teuchos::oblackholestream oldFormatState; 00209 oldFormatState.copyfmt(std::cout); 00210 00211 *outStream \ 00212 << "===============================================================================\n" \ 00213 << "| |\n" \ 00214 << "| Unit Test (IntrepidPolylib) |\n" \ 00215 << "| |\n" \ 00216 << "| 1) the original Polylib tests |\n" \ 00217 << "| |\n" \ 00218 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00219 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00220 << "| |\n" \ 00221 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00222 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00223 << "| |\n" \ 00224 << "===============================================================================\n"; 00225 00226 00227 int errorFlag = 0; 00228 int beginThrowNumber = Teuchos::TestForException_getThrowNumber(); 00229 int endThrowNumber = beginThrowNumber + 1; 00230 00231 typedef IntrepidPolylib ipl; 00232 IntrepidPolylib iplib; 00233 00234 *outStream \ 00235 << "\n" 00236 << "===============================================================================\n"\ 00237 << "| TEST 1: exceptions |\n"\ 00238 << "===============================================================================\n"; 00239 00240 try{ 00241 INTREPID_TEST_COMMAND( iplib.gammaF((double)0.33) ); 00242 } 00243 catch (std::logic_error err) { 00244 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00245 *outStream << err.what() << '\n'; 00246 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00247 errorFlag = -1000; 00248 }; 00249 00250 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) 00251 errorFlag = -1000; 00252 00253 *outStream \ 00254 << "\n" 00255 << "===============================================================================\n"\ 00256 << "| TEST 2: correctness of math operations |\n"\ 00257 << "===============================================================================\n"; 00258 00259 outStream->precision(20); 00260 00261 try { 00262 { // start scope 00263 int np,n,i; 00264 double *z,*w,*p,sum=0,alpha,beta,*d; 00265 00266 z = dvector<double>(0,NPUPPER-1); 00267 w = dvector<double>(0,NPUPPER-1); 00268 p = dvector<double>(0,NPUPPER-1); 00269 00270 d = dvector<double>(0,NPUPPER*NPUPPER-1); 00271 00272 #if GAUSS_INT 00273 /* Gauss Integration */ 00274 *outStream << "Begin checking Gauss integration\n"; 00275 alpha = -0.5; 00276 while(alpha <= 5.0){ 00277 beta = -0.5; 00278 while(beta <= 5.0){ 00279 00280 for(np = NPLOWER; np <= NPUPPER; ++np){ 00281 ipl::zwgj(z,w,np,alpha,beta); 00282 for(n = 2; n < 2*np-1; ++n){ 00283 ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta); 00284 sum = ddot(np,w,1,p,1); 00285 if(std::abs(sum)>TEST_EPS) { 00286 errorFlag = -1000; 00287 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00288 ", np = " << np << ", n = " << n << " integral was " << sum << "\n"; 00289 } 00290 } 00291 } 00292 00293 beta += 0.5; 00294 } 00295 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00296 alpha += 0.5; 00297 } 00298 *outStream << "Finished checking Gauss Integration\n"; 00299 #endif 00300 00301 #if GAUSS_RADAUM_INT 00302 /* Gauss Radau (z = -1) Integration */ 00303 *outStream << "Begin checking Gauss-Radau (z = -1) integration\n"; 00304 alpha = -0.5; 00305 while(alpha <= 5.0){ 00306 beta = -0.5; 00307 while(beta <= 5.0){ 00308 00309 for(np = NPLOWER; np <= NPUPPER; ++np){ 00310 ipl::zwgrjm(z,w,np,alpha,beta); 00311 for(n = 2; n < 2*np-2; ++n){ 00312 ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta); 00313 sum = ddot(np,w,1,p,1); 00314 if(std::abs(sum)>TEST_EPS) { 00315 errorFlag = -1000; 00316 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00317 ", np = " << np << ", n = " << n << " integral was " << sum << "\n"; 00318 } 00319 } 00320 } 00321 00322 beta += 0.5; 00323 } 00324 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00325 alpha += 0.5; 00326 } 00327 *outStream << "Finished checking Gauss-Radau (z = -1) Integration\n"; 00328 #endif 00329 00330 #if GAUSS_RADAUP_INT 00331 /* Gauss Radau (z = +1) Integration */ 00332 *outStream << "Begin checking Gauss-Radau (z = +1) integration\n"; 00333 alpha = -0.5; 00334 while(alpha <= 5.0){ 00335 beta = -0.5; 00336 while(beta <= 5.0){ 00337 00338 for(np = NPLOWER; np <= NPUPPER; ++np){ 00339 ipl::zwgrjp(z,w,np,alpha,beta); 00340 for(n = 2; n < 2*np-2; ++n){ 00341 ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta); 00342 sum = ddot(np,w,1,p,1); 00343 if(std::abs(sum)>TEST_EPS) { 00344 errorFlag = -1000; 00345 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00346 ", np = " << np << ", n = " << n << " integral was " << sum << "\n"; 00347 } 00348 } 00349 } 00350 00351 beta += 0.5; 00352 } 00353 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00354 alpha += 0.5; 00355 } 00356 *outStream << "Finished checking Gauss-Radau (z = +1) Integration\n"; 00357 #endif 00358 00359 #if GAUSS_LOBATTO_INT 00360 /* Gauss Lobatto Integration */ 00361 *outStream << "Begin checking Gauss-Lobatto integration\n"; 00362 alpha = -0.5; 00363 while(alpha <= 5.0){ 00364 beta = -0.5; 00365 while(beta <= 5.0){ 00366 00367 for(np = NPLOWER; np <= NPUPPER; ++np){ 00368 ipl::zwglj(z,w,np,alpha,beta); 00369 for(n = 2; n < 2*np-3; ++n){ 00370 ipl::jacobfd(np,z,p,(double*)0,n,alpha,beta); 00371 sum = ddot(np,w,1,p,1); 00372 if(std::abs(sum)>TEST_EPS) { 00373 errorFlag = -1000; 00374 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00375 ", np = " << np << ", n = " << n << " integral was " << sum << "\n"; 00376 } 00377 } 00378 } 00379 00380 beta += 0.5; 00381 } 00382 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00383 alpha += 0.5; 00384 } 00385 *outStream << "Finished checking Gauss-Lobatto Integration\n"; 00386 #endif 00387 00388 #if GAUSS_DIFF 00389 *outStream << "Begin checking differentiation through Gauss points\n"; 00390 alpha = -0.5; 00391 while(alpha <= 5.0){ 00392 beta = -0.5; 00393 while(beta <= 5.0){ 00394 00395 for(np = NPLOWER; np <= NPUPPER; ++np){ 00396 ipl::zwgj(z,w,np,alpha,beta); 00397 for(n = 2; n < np-1; ++n){ 00398 ipl::Dgj(d,z,np,alpha,beta); 00399 00400 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n); 00401 sum = 0; 00402 for(i = 0; i < np; ++i) 00403 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1)); 00404 sum /= np; 00405 if(std::abs(sum)>TEST_EPS) { 00406 errorFlag = -1000; 00407 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00408 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00409 } 00410 } 00411 } 00412 beta += 0.5; 00413 } 00414 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00415 alpha += 0.5; 00416 } 00417 *outStream << "Finished checking Gauss Jacobi differentiation\n"; 00418 #endif 00419 00420 #if GAUSS_RADAUM_DIFF 00421 *outStream << "Begin checking differentiation through Gauss-Radau (z=-1) points\n"; 00422 alpha = -0.5; 00423 while(alpha <= 5.0){ 00424 beta = -0.5; 00425 while(beta <= 5.0){ 00426 00427 for(np = NPLOWER; np <= NPUPPER; ++np){ 00428 ipl::zwgrjm(z,w,np,alpha,beta); 00429 for(n = 2; n < np-1; ++n){ 00430 ipl::Dgrjm(d,z,np,alpha,beta); 00431 00432 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n); 00433 sum = 0; 00434 for(i = 0; i < np; ++i) 00435 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1)); 00436 sum /= np; 00437 if(std::abs(sum)>TEST_EPS) { 00438 errorFlag = -1000; 00439 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00440 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00441 } 00442 } 00443 } 00444 beta += 0.5; 00445 } 00446 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00447 alpha += 0.5; 00448 } 00449 *outStream << "Finished checking Gauss-Radau (z=-1) differentiation\n"; 00450 #endif 00451 00452 #if GAUSS_RADAUP_DIFF 00453 *outStream << "Begin checking differentiation through Gauss-Radau (z=+1) points\n"; 00454 alpha = -0.5; 00455 while(alpha <= 5.0){ 00456 beta = -0.5; 00457 while(beta <= 5.0){ 00458 00459 for(np = NPLOWER; np <= NPUPPER; ++np){ 00460 ipl::zwgrjp(z,w,np,alpha,beta); 00461 for(n = 2; n < np-1; ++n){ 00462 ipl::Dgrjp(d,z,np,alpha,beta); 00463 00464 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n); 00465 sum = 0; 00466 for(i = 0; i < np; ++i) 00467 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1)); 00468 sum /= np; 00469 if(std::abs(sum)>TEST_EPS) { 00470 errorFlag = -1000; 00471 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00472 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00473 } 00474 } 00475 } 00476 beta += 0.5; 00477 } 00478 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00479 alpha += 0.5; 00480 } 00481 *outStream << "Finished checking Gauss-Radau (z=+1) differentiation\n"; 00482 #endif 00483 00484 #if GAUSS_RADAUP_DIFF 00485 *outStream << "Begin checking differentiation through Gauss-Lobatto (z=+1) points\n"; 00486 alpha = -0.5; 00487 while(alpha <= 5.0){ 00488 beta = -0.5; 00489 while(beta <= 5.0){ 00490 00491 for(np = NPLOWER; np <= NPUPPER; ++np){ 00492 ipl::zwglj(z,w,np,alpha,beta); 00493 for(n = 2; n < np-1; ++n){ 00494 ipl::Dglj(d,z,np,alpha,beta); 00495 00496 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n); 00497 sum = 0; 00498 for(i = 0; i < np; ++i) 00499 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1)); 00500 sum /= np; 00501 if(std::abs(sum)>TEST_EPS) { 00502 errorFlag = -1000; 00503 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00504 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00505 } 00506 } 00507 } 00508 beta += 0.5; 00509 } 00510 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00511 alpha += 0.5; 00512 } 00513 *outStream << "Finished checking Gauss-Lobatto differentiation\n"; 00514 #endif 00515 00516 #if GAUSS_INTERP 00517 *outStream << "Begin checking interpolation through Gauss points\n"; 00518 alpha = -0.5; 00519 while(alpha <= 5.0){ 00520 beta = -0.5; 00521 while(beta <= 5.0){ 00522 00523 for(np = NPLOWER; np <= NPUPPER; ++np){ 00524 ipl::zwgj(z,w,np,alpha,beta); 00525 for(n = 2; n < np-1; ++n){ 00526 for(i = 0; i < np; ++i) { 00527 w[i] = 2.0*i/(double)(np-1)-1.0; 00528 p[i] = std::pow(z[i],n); 00529 } 00530 ipl::Imgj(d,z,w,np,np,alpha,beta); 00531 sum = 0; 00532 for(i = 0; i < np; ++i) 00533 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n)); 00534 sum /= np; 00535 if(std::abs(sum)>TEST_EPS) { 00536 errorFlag = -1000; 00537 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00538 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00539 } 00540 } 00541 } 00542 beta += 0.5; 00543 } 00544 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00545 alpha += 0.5; 00546 } 00547 *outStream << "Finished checking Gauss Jacobi interpolation\n"; 00548 #endif 00549 00550 #if GAUSS_RADAUM_INTERP 00551 *outStream << "Begin checking interpolation through Gauss-Radau (z=-1) points\n"; 00552 alpha = -0.5; 00553 while(alpha <= 5.0){ 00554 beta = -0.5; 00555 while(beta <= 5.0){ 00556 00557 for(np = NPLOWER; np <= NPUPPER; ++np){ 00558 ipl::zwgrjm(z,w,np,alpha,beta); 00559 for(n = 2; n < np-1; ++n){ 00560 for(i = 0; i < np; ++i) { 00561 w[i] = 2.0*i/(double)(np-1)-1.0; 00562 p[i] = std::pow(z[i],n); 00563 } 00564 ipl::Imgrjm(d,z,w,np,np,alpha,beta); 00565 sum = 0; 00566 for(i = 0; i < np; ++i) 00567 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n)); 00568 sum /= np; 00569 if(std::abs(sum)>TEST_EPS) { 00570 errorFlag = -1000; 00571 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00572 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00573 } 00574 } 00575 } 00576 beta += 0.5; 00577 } 00578 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00579 alpha += 0.5; 00580 } 00581 *outStream << "Finished checking Gauss-Radau (z=-1) interpolation\n"; 00582 #endif 00583 00584 #if GAUSS_RADAUP_INTERP 00585 *outStream << "Begin checking interpolation through Gauss-Radau (z=+1) points\n"; 00586 alpha = -0.5; 00587 while(alpha <= 5.0){ 00588 beta = -0.5; 00589 while(beta <= 5.0){ 00590 00591 for(np = NPLOWER; np <= NPUPPER; ++np){ 00592 ipl::zwgrjp(z,w,np,alpha,beta); 00593 for(n = 2; n < np-1; ++n){ 00594 for(i = 0; i < np; ++i) { 00595 w[i] = 2.0*i/(double)(np-1)-1.0; 00596 p[i] = std::pow(z[i],n); 00597 } 00598 ipl::Imgrjp(d,z,w,np,np,alpha,beta); 00599 sum = 0; 00600 for(i = 0; i < np; ++i) 00601 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n)); 00602 sum /= np; 00603 if(std::abs(sum)>TEST_EPS) { 00604 errorFlag = -1000; 00605 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00606 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00607 } 00608 } 00609 } 00610 beta += 0.5; 00611 } 00612 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00613 alpha += 0.5; 00614 } 00615 *outStream << "Finished checking Gauss-Radau (z=+1) interpolation\n"; 00616 #endif 00617 00618 #if GAUSS_LOBATTO_INTERP 00619 *outStream << "Begin checking interpolation through Gauss-Lobatto points\n"; 00620 alpha = -0.5; 00621 while(alpha <= 5.0){ 00622 beta = -0.5; 00623 while(beta <= 5.0){ 00624 00625 for(np = NPLOWER; np <= NPUPPER; ++np){ 00626 ipl::zwglj(z,w,np,alpha,beta); 00627 for(n = 2; n < np-1; ++n){ 00628 for(i = 0; i < np; ++i) { 00629 w[i] = 2.0*i/(double)(np-1)-1.0; 00630 p[i] = std::pow(z[i],n); 00631 } 00632 ipl::Imglj(d,z,w,np,np,alpha,beta); 00633 sum = 0; 00634 for(i = 0; i < np; ++i) 00635 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n)); 00636 sum /= np; 00637 if(std::abs(sum)>TEST_EPS) { 00638 errorFlag = -1000; 00639 *outStream << "ERROR: alpha = " << alpha << ", beta = " << beta << 00640 ", np = " << np << ", n = " << n << " difference " << sum << "\n"; 00641 } 00642 } 00643 } 00644 beta += 0.5; 00645 } 00646 *outStream << "finished checking all beta values for alpha = " << alpha << "\n"; 00647 alpha += 0.5; 00648 } 00649 *outStream << "Finished checking Gauss-Lobatto interpolation\n"; 00650 #endif 00651 00652 free(z); free(w); free(p); free(d); 00653 00654 } // end scope 00655 00656 /******************************************/ 00657 *outStream << "\n"; 00658 } 00659 catch (std::logic_error err) { 00660 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00661 *outStream << err.what() << '\n'; 00662 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00663 errorFlag = -1000; 00664 }; 00665 00666 00667 if (errorFlag != 0) 00668 std::cout << "End Result: TEST FAILED\n"; 00669 else 00670 std::cout << "End Result: TEST PASSED\n"; 00671 00672 // reset format state of std::cout 00673 std::cout.copyfmt(oldFormatState); 00674 00675 return errorFlag; 00676 }
1.7.6.1