Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #include <math.h>
00032 #include "SundanceStdMathFunctors.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "Teuchos_Utils.hpp"
00035 
00036 using namespace Sundance;
00037 using namespace Teuchos;
00038 using std::endl;
00039 
00040 PowerFunctor::PowerFunctor(const double& p) 
00041   : UnaryFunctor("pow("+Teuchos::toString(p)+")",
00042      powerDomain(p)), 
00043     p_(p),
00044     powerIsInteger_(p==floor(p))
00045 {}
00046 
00047 RCP<FunctorDomain> PowerFunctor::powerDomain(const double& p)
00048 {
00049   /* There are four cases: 
00050    * (1) p is a positive integer, in which case the domain is all of R.
00051    * (2) p is a negative integer, in which case the domain is R\0.
00052    * (3) p is a positive real \notin Z, in which case the domain is [0,\infty).
00053    * (3) p is a negative real \notin Z, in which case the domain is (0,\infty).
00054    */
00055   bool isInZ = floor(p)==p;
00056   bool isNegative = p < 0.0;
00057 
00058   if (isInZ)
00059     {
00060       if (isNegative) return rcp(new NonzeroDomain());
00061     }
00062   else
00063     {
00064       if (isNegative) return rcp(new StrictlyPositiveDomain());
00065       return rcp(new PositiveDomain());
00066     }
00067   return rcp(new UnboundedDomain());
00068 }
00069 
00070 void PowerFunctor::eval1(const double* const x, 
00071                         int nx, 
00072                         double* f, 
00073                         double* df) const
00074 {
00075   if (p_==2)
00076     {
00077       for (int i=0; i<nx; i++) 
00078         {
00079           df[i] = 2.0*x[i];
00080           f[i] = x[i]*x[i];
00081         }
00082     }
00083   else if (p_==1)
00084     {
00085       for (int i=0; i<nx; i++) 
00086         {
00087           df[i] = 1.0;
00088           f[i] = x[i];
00089         }
00090     }
00091   else if (p_==0)
00092     {
00093       for (int i=0; i<nx; i++) 
00094         {
00095           df[i] = 0.0;
00096           f[i] = 1.0;
00097         }
00098     }
00099   else
00100     {
00101       if (checkResults())
00102         {
00103           for (int i=0; i<nx; i++) 
00104             {
00105         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(1,x[i]), std::runtime_error,
00106          "first deriv of pow(" << x[i] 
00107          << ", " << p_ << ") "
00108          "is undefined");
00109         
00110               double px = ::pow(x[i], p_-1);
00111               df[i] = p_*px;
00112               f[i] = x[i]*px;
00113               //bvbw tried to include math.h, without success
00114 #ifdef REDDISH_PORT_PROBLEM
00115               TEUCHOS_TEST_FOR_EXCEPTION(fpclassify(f[i]) != FP_NORMAL 
00116                                  || fpclassify(df[i]) != FP_NORMAL,
00117                                  std::runtime_error,
00118                                  "Non-normal floating point result detected in "
00119                                  "evaluation of unary functor " << name());
00120 #endif
00121             }
00122         }
00123       else 
00124         {
00125     for (int i=0; i<nx; i++) 
00126       {
00127         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(1,x[i]), std::runtime_error,
00128          "first deriv of pow(" << x[i] 
00129          << ", " << p_ << ") "
00130          "is undefined");
00131         double px = ::pow(x[i], p_-1);
00132         df[i] = p_*px;
00133         f[i] = x[i]*px;
00134       }
00135         }
00136     }
00137 }
00138 
00139 
00140 
00141 
00142 void PowerFunctor::eval3(const double* const x, 
00143                          int nx, 
00144                          double* f, 
00145                          double* df,
00146                          double* d2f_dxx,
00147                          double* d3f_dxxx) const
00148 {
00149   if (p_==3)
00150     {
00151       for (int i=0; i<nx; i++) 
00152         {
00153           d3f_dxxx[i] = 6.0;
00154           d2f_dxx[i] = 6.0*x[i];
00155           df[i] = 3.0*x[i]*x[i];
00156           f[i] = x[i]*x[i]*x[i];
00157         }
00158     }
00159   else if (p_==2)
00160     {
00161       for (int i=0; i<nx; i++) 
00162         {
00163           d3f_dxxx[i] = 0.0;
00164           d2f_dxx[i] = 2.0;
00165           df[i] = 2.0*x[i];
00166           f[i] = x[i]*x[i];
00167         }
00168     }
00169   else if (p_==1)
00170     {
00171        for (int i=0; i<nx; i++) 
00172         {
00173           d3f_dxxx[i] = 0.0;
00174           d2f_dxx[i] = 0.0;
00175           df[i] = 1.0;
00176           f[i] = x[i];
00177         }
00178     }
00179   else if (p_==0)
00180     {
00181       for (int i=0; i<nx; i++) 
00182         {
00183           d3f_dxxx[i] = 0.0;
00184           d2f_dxx[i] = 0.0;
00185           df[i] = 0.0;
00186           f[i] = 1.0;
00187         }
00188     }
00189   else
00190     {
00191       if (checkResults())
00192         {
00193           for (int i=0; i<nx; i++) 
00194             {
00195               double px = ::pow(x[i], p_-3);
00196               d3f_dxxx[i] = p_ * (p_-1) * (p_-2) * px;
00197               d2f_dxx[i] = p_ * (p_-1) * x[i] * px;
00198               df[i] = p_*x[i]*x[i]*px;
00199               f[i] = x[i]*x[i]*x[i]*px;
00200         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(3,x[i]), std::runtime_error,
00201          "third deriv of pow(" << x[i] 
00202          << ", " << p_ << ") "
00203          "is undefined");
00204 
00205 
00206 #ifdef REDDISH_PORT_PROBLEM
00207               TEUCHOS_TEST_FOR_EXCEPTION(fpclassify(f[i]) != FP_NORMAL 
00208                                  || fpclassify(df[i]) != FP_NORMAL,
00209                                  std::runtime_error,
00210                                  "Non-normal floating point result detected in "
00211                                  "evaluation of unary functor " << name());
00212 #endif
00213             }
00214         }
00215       else
00216         {
00217           for (int i=0; i<nx; i++) 
00218             {
00219         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(3,x[i]), std::runtime_error,
00220          "third deriv of pow(" << x[i] 
00221          << ", " << p_ << ") "
00222          "is undefined");
00223 
00224               double px = ::pow(x[i], p_-3);
00225               d3f_dxxx[i] = p_ * (p_-1) * (p_-2) * px;
00226               d2f_dxx[i] = p_ * (p_-1) * x[i] * px;
00227               df[i] = p_*x[i]*x[i]*px;
00228               f[i] = x[i]*x[i]*x[i]*px;
00229             }
00230         }
00231     }
00232 }
00233 
00234 
00235 void PowerFunctor::eval2(const double* const x, 
00236                         int nx, 
00237                         double* f, 
00238                         double* df,
00239                         double* d2f_dxx) const
00240 {
00241   if (p_==2)
00242     {
00243       for (int i=0; i<nx; i++) 
00244         {
00245           d2f_dxx[i] = 2.0;
00246           df[i] = 2.0*x[i];
00247           f[i] = x[i]*x[i];
00248         }
00249     }
00250   else if (p_==1)
00251     {
00252        for (int i=0; i<nx; i++) 
00253         {
00254           d2f_dxx[i] = 0.0;
00255           df[i] = 1.0;
00256           f[i] = x[i];
00257         }
00258     }
00259   else if (p_==0)
00260     {
00261       for (int i=0; i<nx; i++) 
00262         {
00263           d2f_dxx[i] = 0.0;
00264           df[i] = 0.0;
00265           f[i] = 1.0;
00266         }
00267     }
00268   else
00269     {
00270       if (checkResults())
00271         {
00272           for (int i=0; i<nx; i++) 
00273             {
00274         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(2,x[i]), std::runtime_error,
00275          "second deriv of pow(" << x[i] 
00276          << ", " << p_ << ") "
00277          "is undefined");
00278 
00279 
00280               double px = ::pow(x[i], p_-2);
00281               d2f_dxx[i] = p_ * (p_-1) * px;
00282               df[i] = p_*x[i]*px;
00283               f[i] = x[i]*x[i]*px;
00284 #ifdef REDDISH_PORT_PROBLEM
00285               TEUCHOS_TEST_FOR_EXCEPTION(fpclassify(f[i]) != FP_NORMAL 
00286                                  || fpclassify(df[i]) != FP_NORMAL,
00287                                  std::runtime_error,
00288                                  "Non-normal floating point result detected in "
00289                                  "evaluation of unary functor " << name());
00290 #endif
00291             }
00292         }
00293       else
00294         {
00295     for (int i=0; i<nx; i++) 
00296       {
00297         TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(2,x[i]), std::runtime_error,
00298          "second deriv of pow(" << x[i] 
00299          << ", " << p_ << ") "
00300          "is undefined");
00301         
00302         double px = ::pow(x[i], p_-2);
00303         
00304         d2f_dxx[i] = p_ * (p_-1) * px;
00305         df[i] = p_*x[i]*px;
00306         f[i] = x[i]*x[i]*px;
00307       }
00308   }
00309     }
00310 }
00311 
00312 
00313 
00314 void PowerFunctor::eval0(const double* const x, 
00315                         int nx, 
00316                         double* f) const
00317 {
00318   if (checkResults())
00319     {
00320       for (int i=0; i<nx; i++) 
00321         {
00322     TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(0,x[i]), std::runtime_error,
00323            "pow(" << x[i] 
00324            << ", " << p_ << ") "
00325            "is undefined");
00326 
00327 
00328           f[i] = ::pow(x[i], p_);
00329 #ifdef REDDISH_PORT_PROBLEM
00330           TEUCHOS_TEST_FOR_EXCEPTION(fpclassify(f[i]) != FP_NORMAL, 
00331                              std::runtime_error,
00332                              "Non-normal floating point result detected in "
00333                              "evaluation of unary functor " << name());
00334 #endif
00335   }
00336     }
00337   else
00338     {
00339       for (int i=0; i<nx; i++) 
00340   {
00341     TEUCHOS_TEST_FOR_EXCEPTION(!acceptX(0,x[i]), std::runtime_error,
00342            "pow(" << x[i] 
00343            << ", " << p_ << ") "
00344            "is undefined");
00345 
00346     f[i] = ::pow(x[i], p_);
00347   }
00348     }
00349 }
SundanceStdMathFunctors.cpp
