00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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
00050
00051
00052
00053
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
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 }