Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Shared/IntrepidPolylib/test

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 = 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 (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 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Shared/IntrepidPolylib/test_01.cpp
