diff --git a/ADOL-C/boost-test/traceOperatorScalar.cpp b/ADOL-C/boost-test/traceOperatorScalar.cpp index f58e4ff5..0136bfda 100644 --- a/ADOL-C/boost-test/traceOperatorScalar.cpp +++ b/ADOL-C/boost-test/traceOperatorScalar.cpp @@ -4514,6 +4514,99 @@ BOOST_AUTO_TEST_CASE(ErfOperator_FOS_Reverse) myfree1(z); } +BOOST_AUTO_TEST_CASE(ErfcOperator_ZOS_Forward) +{ + double a = 7.1, aout; + adouble ad; + + trace_on(1); + ad <<= a; + + ad = erfc(ad); + + ad >>= aout; + trace_off(); + + a = std::erfc(a); + + BOOST_TEST(aout == a, tt::tolerance(tol)); + + double *x = myalloc1(1); + double *y = myalloc1(1); + + *x = 7.1; + + zos_forward(1, 1, 1, 0, x, y); + + BOOST_TEST(*y == a, tt::tolerance(tol)); + + myfree1(x); + myfree1(y); +} + +BOOST_AUTO_TEST_CASE(ErfcOperator_FOS_Forward) +{ + double a = 7.1, aout; + adouble ad; + + trace_on(1); + ad <<= a; + + ad = erfc(ad); + + ad >>= aout; + trace_off(); + + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a * a); + a = std::erfc(a); + + double *x = myalloc1(1); + double *xd = myalloc1(1); + double *y = myalloc1(1); + double *yd = myalloc1(1); + + *x = 7.1; + *xd = 1.; + + fos_forward(1, 1, 1, 0, x, xd, y, yd); + + BOOST_TEST(*y == a, tt::tolerance(tol)); + BOOST_TEST(*yd == aDerivative, tt::tolerance(tol)); + + myfree1(x); + myfree1(xd); + myfree1(y); + myfree1(yd); +} + +BOOST_AUTO_TEST_CASE(ErfcOperator_FOS_Reverse) +{ + double a = 7.1, aout; + adouble ad; + + trace_on(1, 1); + ad <<= a; + + ad = erfc(ad); + + ad >>= aout; + trace_off(); + + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a * a); + + double *u = myalloc1(1); + double *z = myalloc1(1); + + *u = 1.; + + fos_reverse(1, 1, 1, u, z); + + BOOST_TEST(*z == aDerivative, tt::tolerance(tol)); + + myfree1(u); + myfree1(z); +} + BOOST_AUTO_TEST_CASE(EqOperator_ZOS_Forward) { double a = 10.01, aout; diff --git a/ADOL-C/boost-test/traceOperatorVector.cpp b/ADOL-C/boost-test/traceOperatorVector.cpp index 3867ba0a..0fb6c594 100644 --- a/ADOL-C/boost-test/traceOperatorVector.cpp +++ b/ADOL-C/boost-test/traceOperatorVector.cpp @@ -3394,6 +3394,75 @@ BOOST_AUTO_TEST_CASE(ErfOperator_FOV_Reverse) myfree2(z); } +BOOST_AUTO_TEST_CASE(ErfcOperator_FOV_Forward) +{ + double a = 7.1, aout; + adouble ad; + + trace_on(1); + ad <<= a; + + ad = erfc(ad); + + ad >>= aout; + trace_off(); + + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a*a); + a = std::erfc(a); + + double *x = myalloc1(1); + double **xd = myalloc2(1, 2); + double *y = myalloc1(1); + double **yd = myalloc2(1, 2); + + x[0] = 7.1; + + for (int i = 0; i < 2; i++) { + xd[0][i] = std::pow(3., i*2.); + } + + fov_forward(1, 1, 1, 2, x, xd, y, yd); + + BOOST_TEST(*y == a, tt::tolerance(tol)); + BOOST_TEST(yd[0][0] == aDerivative, tt::tolerance(tol)); + BOOST_TEST(yd[0][1] == aDerivative * std::pow(3., 2.), tt::tolerance(tol)); + + myfree1(x); + myfree2(xd); + myfree1(y); + myfree2(yd); +} + +BOOST_AUTO_TEST_CASE(ErfcOperator_FOV_Reverse) +{ + double a = 7.1, aout; + adouble ad; + + trace_on(1, 1); + ad <<= a; + + ad = erfc(ad); + + ad >>= aout; + trace_off(); + + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a*a); + + double **u = myalloc2(2, 1); + double **z = myalloc2(2, 1); + + u[0][0] = 1.; + u[1][0] = -1.1; + + fov_reverse(1, 1, 1, 2, u, z); + + BOOST_TEST(z[0][0] == aDerivative, tt::tolerance(tol)); + BOOST_TEST(z[1][0] == -1.1*aDerivative, tt::tolerance(tol)); + + myfree2(u); + myfree2(z); +} + BOOST_AUTO_TEST_CASE(EqPlusOperator_FOV_Forward) { double a = 5.132, aout; diff --git a/ADOL-C/boost-test/tracelessCompositeTests.cpp b/ADOL-C/boost-test/tracelessCompositeTests.cpp index 1ad29cdc..47172f70 100644 --- a/ADOL-C/boost-test/tracelessCompositeTests.cpp +++ b/ADOL-C/boost-test/tracelessCompositeTests.cpp @@ -282,6 +282,32 @@ BOOST_AUTO_TEST_CASE(CompositeErfFabs_Traceless) BOOST_TEST(ax1.getADValue(1) == x2Derivative, tt::tolerance(tol)); } +BOOST_AUTO_TEST_CASE(CompositeErfcFabs_Traceless) +{ + double x1 = 4.56, x2 = 6.45; + adouble ax1 = x1, ax2 = x2; + + ax1.setADValue(0, 1.); + ax2.setADValue(1, 1.); + + ax1 = adtl::erfc(adtl::fabs(ax1 - 2.)*adtl::sinh(ax2)); + + double x1Derivative = -2. / std::sqrt(std::acos(-1.)) + * std::exp(-std::pow(std::fabs(x1 - 2.) + * std::sinh(x2), 2)) + * std::sinh(x2); + double x2Derivative = -2. / std::sqrt(std::acos(-1.)) + * std::exp(-std::pow(std::fabs(x1 - 2.) + * std::sinh(x2), 2)) + * std::fabs(x1 - 2.) * std::cosh(x2); + + x1 = std::erfc(std::fabs(x1 - 2.)*std::sinh(x2)); + + BOOST_TEST(ax1.getValue() == x1, tt::tolerance(tol)); + BOOST_TEST(ax1.getADValue(0) == x1Derivative, tt::tolerance(tol)); + BOOST_TEST(ax1.getADValue(1) == x2Derivative, tt::tolerance(tol)); +} + BOOST_AUTO_TEST_CASE(ExpTrigSqrt_Traceless) { double x1 = 1.2, x2 = 2.1; diff --git a/ADOL-C/boost-test/tracelessOperatorScalar.cpp b/ADOL-C/boost-test/tracelessOperatorScalar.cpp index b543f6e3..00737902 100644 --- a/ADOL-C/boost-test/tracelessOperatorScalar.cpp +++ b/ADOL-C/boost-test/tracelessOperatorScalar.cpp @@ -1436,6 +1436,36 @@ BOOST_AUTO_TEST_CASE(ErfOperatorDerivative) BOOST_TEST(ad.getADValue(0) == aDerivative, tt::tolerance(tol)); } +/* The complementary error function erfc(a) is defined as + * 1.0 - erf(a). + */ +BOOST_AUTO_TEST_CASE(ErfcOperatorPrimal) +{ + double a = 7.1; + adouble ad = a; + + a = std::erfc(a); + ad = adtl::erfc(ad); + + BOOST_TEST(ad.getValue() == a, tt::tolerance(tol)); +} + +BOOST_AUTO_TEST_CASE(ErfcOperatorDerivative) +{ + double a = 7.1; + adouble ad = a; + + /* The inverse cosine is used for pi so that no additional math library + * import is necessary, see also the implementation in adtl.h. + */ + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a * a); + + ad.setADValue(0, 1.); + ad = adtl::erfc(ad); + + BOOST_TEST(ad.getADValue(0) == aDerivative, tt::tolerance(tol)); +} + /* Test the primitive non-temporary operations =, +=, -=, *=, /=. */ BOOST_AUTO_TEST_CASE(EqOperatorPrimal_Derivative) diff --git a/ADOL-C/boost-test/tracelessOperatorVector.cpp b/ADOL-C/boost-test/tracelessOperatorVector.cpp index 24d4d7fd..e060098b 100644 --- a/ADOL-C/boost-test/tracelessOperatorVector.cpp +++ b/ADOL-C/boost-test/tracelessOperatorVector.cpp @@ -868,6 +868,21 @@ BOOST_AUTO_TEST_CASE(ErfOperatorDerivativeVectorMode) BOOST_TEST(ad.getADValue(1) == aDerivative*(-2.5), tt::tolerance(tol)); } +BOOST_AUTO_TEST_CASE(ErfcOperatorDerivativeVectorMode) +{ + double a = 7.1; + adouble ad = a; + + double aDerivative = -2. / std::sqrt(std::acos(-1.)) * std::exp(-a * a); + + ad.setADValue(0, 1.); + ad.setADValue(1, -2.5); + ad = adtl::erfc(ad); + + BOOST_TEST(ad.getADValue(0) == aDerivative, tt::tolerance(tol)); + BOOST_TEST(ad.getADValue(1) == aDerivative*(-2.5), tt::tolerance(tol)); +} + BOOST_AUTO_TEST_CASE(EqOperatorDerivativeVectorMode) { double a = 10.01; diff --git a/ADOL-C/c_interface/ADOLC_TB_interface.cpp b/ADOL-C/c_interface/ADOLC_TB_interface.cpp index 1e7036d4..6dfcea24 100644 --- a/ADOL-C/c_interface/ADOLC_TB_interface.cpp +++ b/ADOL-C/c_interface/ADOLC_TB_interface.cpp @@ -261,6 +261,10 @@ extern "C" { return new adouble(erf(*static_cast(a))); } + TBAdoubleHandle tb_erfc(TBAdoubleHandle a) + { + return new adouble(erfc(*static_cast(a))); + } } /* diff --git a/ADOL-C/c_interface/ADOLC_TB_interface.h b/ADOL-C/c_interface/ADOLC_TB_interface.h index 438d42e4..bce1d0ef 100644 --- a/ADOL-C/c_interface/ADOLC_TB_interface.h +++ b/ADOL-C/c_interface/ADOLC_TB_interface.h @@ -115,6 +115,7 @@ extern "C" TBAdoubleHandle tb_floor(TBAdoubleHandle a); TBAdoubleHandle tb_ldexp(TBAdoubleHandle a, int n); TBAdoubleHandle tb_erf(TBAdoubleHandle a); + TBAdoubleHandle tb_erfc(TBAdoubleHandle a); #ifdef __cplusplus } #endif @@ -139,4 +140,4 @@ extern "C" } #endif -#endif // ADOLC_TB_INTERFACE_H \ No newline at end of file +#endif // ADOLC_TB_INTERFACE_H diff --git a/ADOL-C/c_interface/ADOLC_TL_interface.cpp b/ADOL-C/c_interface/ADOLC_TL_interface.cpp index de85c659..05beab39 100644 --- a/ADOL-C/c_interface/ADOLC_TL_interface.cpp +++ b/ADOL-C/c_interface/ADOLC_TL_interface.cpp @@ -292,4 +292,8 @@ extern "C" { return new adtl::adouble(erf(*static_cast(a))); } + TLAdoubleHandle tl_erfc(TLAdoubleHandle a) + { + return new adtl::adouble(erfc(*static_cast(a))); + } } diff --git a/ADOL-C/c_interface/ADOLC_TL_interface.h b/ADOL-C/c_interface/ADOLC_TL_interface.h index 5558463e..e329644b 100644 --- a/ADOL-C/c_interface/ADOLC_TL_interface.h +++ b/ADOL-C/c_interface/ADOLC_TL_interface.h @@ -122,8 +122,9 @@ extern "C" TLAdoubleHandle tl_floor(TLAdoubleHandle a); TLAdoubleHandle tl_ldexp(TLAdoubleHandle a, const int n); TLAdoubleHandle tl_erf(TLAdoubleHandle a); + TLAdoubleHandle tl_erfc(TLAdoubleHandle a); #ifdef __cplusplus } #endif -#endif // ADOLC_TL_INTERFACE_H \ No newline at end of file +#endif // ADOLC_TL_INTERFACE_H diff --git a/ADOL-C/include/adolc/adoublecuda.h b/ADOL-C/include/adolc/adoublecuda.h index 78dc15b8..32acc231 100644 --- a/ADOL-C/include/adolc/adoublecuda.h +++ b/ADOL-C/include/adolc/adoublecuda.h @@ -152,6 +152,7 @@ class adouble { CUDADEV friend adouble ldexp (const double v, const adouble &a); CUDADEV friend double frexp (const adouble &a, int* v); CUDADEV friend adouble erf (const adouble &a); + CUDADEV friend adouble erfc (const adouble &a); /******************* nontemporary results ***************************/ @@ -794,6 +795,16 @@ CUDADEV adouble erf (const adouble &a) { return tmp; } +CUDADEV adouble erfc (const adouble &a) { + adouble tmp; + tmp.val=ADOLC_MATH_NSP_ERF::erfc(a.val); + double tmp2 = -2.0 / + ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) * + ADOLC_MATH_NSP_ERF::exp(-a.val*a.val); + FOR_I_EQ_0_LT_NUMDIR + tmp.ADVAL_I=tmp2*a.ADVAL_I; + return tmp; +} /******************* nontemporary results *********************************/ CUDADEV void adouble::operator = (const double v) { diff --git a/ADOL-C/include/adolc/adtl_hov.h b/ADOL-C/include/adolc/adtl_hov.h index 451e3e6d..3a9e2986 100644 --- a/ADOL-C/include/adolc/adtl_hov.h +++ b/ADOL-C/include/adolc/adtl_hov.h @@ -1938,6 +1938,26 @@ inline adouble erf (const adouble &a) { return tmp; } +inline adouble erfc (const adouble &a) { + adouble tmp; + if (unlikely(!adouble::_do_val() && adouble::_do_adval())) { + fprintf(DIAG_OUT, "ADOL-C error: Traceless: Incorrect mode, call setMode(enum Mode mode)\n"); + throw logic_error("incorrect function call, errorcode=1"); + } + if (do_val()) + tmp.val=ADOLC_MATH_NSP_ERF::erfc(a.val); + if (likely(adouble::_do_adval() && adouble::_do_val())) { + double tmp2 = -2.0 / + ADOLC_MATH_NSP_ERF::sqrt(ADOLC_MATH_NSP::acos(-1.0)) * + ADOLC_MATH_NSP_ERF::exp(-a.val*a.val); + FOR_I_EQ_0_LT_NUMDIR + tmp.ADVAL_I=tmp2*a.ADVAL_I; + } + if (do_indo()) + tmp.add_to_pattern( a.get_pattern() ) ; + return tmp; +} + inline void condassign( adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2 ) { if (no_do_val()) { diff --git a/ADOL-C/include/adolc/adtl_indo.h b/ADOL-C/include/adolc/adtl_indo.h index 349ec7cd..f9444268 100644 --- a/ADOL-C/include/adolc/adtl_indo.h +++ b/ADOL-C/include/adolc/adtl_indo.h @@ -129,6 +129,7 @@ class adouble { inline friend adouble ldexp (const double v, const adouble &a); inline friend double frexp (const adouble &a, int* v); inline friend adouble erf (const adouble &a); + inline friend adouble erfc (const adouble &a); inline friend void condassign( adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2 ); @@ -690,6 +691,13 @@ inline adouble erf (const adouble &a) { return tmp; } +inline adouble erfc (const adouble &a) { + adouble tmp; + tmp.val=ADOLC_MATH_NSP_ERF::erfc(a.val); + tmp.add_to_pattern( a.get_pattern() ) ; + return tmp; +} + inline void condassign( adouble &res, const adouble &cond, const adouble &arg1, const adouble &arg2 ) { if (cond.getValue() > 0) diff --git a/ADOL-C/include/adolc/param.h b/ADOL-C/include/adolc/param.h index 4ff1faed..162a81f7 100644 --- a/ADOL-C/include/adolc/param.h +++ b/ADOL-C/include/adolc/param.h @@ -185,6 +185,7 @@ inline adub asinh ( const pdouble& p) { return asinh(adub(p)); } inline adub acosh ( const pdouble& p) { return acosh(adub(p)); } inline adub atanh ( const pdouble& p) { return atanh(adub(p)); } inline adub erf ( const pdouble& p) { return erf(adub(p)); } +inline adub erfc ( const pdouble& p) { return erfc(adub(p)); } inline adub fabs ( const pdouble& p) { return fabs(adub(p)); } inline adub ceil ( const pdouble& p) { return ceil(adub(p)); } inline adub floor ( const pdouble& p) { return floor(adub(p)); } diff --git a/ADOL-C/src/tapedoc/tapedoc.c b/ADOL-C/src/tapedoc/tapedoc.c index 26300c82..c2abd59f 100644 --- a/ADOL-C/src/tapedoc/tapedoc.c +++ b/ADOL-C/src/tapedoc/tapedoc.c @@ -968,6 +968,25 @@ void tape_doc(short tnum, /* tape id */ filewrite(operation,"erf op",3,loc_a,val_a,0,cst_d); break; + /*--------------------------------------------------------------------------*/ + case erfc_op: /* erfc_op */ + arg1 = get_locint_f(); + arg2 = get_locint_f(); + res = get_locint_f(); + loc_a[0]=arg1; + loc_a[1]=arg2; + loc_a[2]=res; +#ifdef ADOLC_TAPE_DOC_VALUES + val_a[0]=dp_T0[arg1]; + dp_T0[res] = erfc(dp_T0[arg1]); + ADOLC_OPENMP_RESTORE_THREAD_NUMBER; + val_a[1]=dp_T0[arg2]; + val_a[2]=dp_T0[res]; +#endif + filewrite(operation,"erfc op",3,loc_a,val_a,0,cst_d); + break; + + /*--------------------------------------------------------------------------*/ case log_op: /* log_op */ arg = get_locint_f(); diff --git a/ADOL-C/swig/adubswigfuncs.hpp b/ADOL-C/swig/adubswigfuncs.hpp index 8036d3fc..8424e8bf 100644 --- a/ADOL-C/swig/adubswigfuncs.hpp +++ b/ADOL-C/swig/adubswigfuncs.hpp @@ -108,6 +108,7 @@ adub* asinh ( const badouble& ); adub* acosh ( const badouble& ); adub* atanh ( const badouble& ); adub* erf ( const badouble& ); +adub* erfc ( const badouble& ); adub* fabs ( const badouble& ); adub* ceil ( const badouble& ); @@ -147,6 +148,7 @@ adub* asinh ( const pdouble& ); adub* acosh ( const pdouble& ); adub* atanh ( const pdouble& ); adub* erf ( const pdouble& ); +adub* erfc ( const pdouble& ); adub* fabs ( const pdouble& ); adub* ceil ( const pdouble& );