boost tests for erfc

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;

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;

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;

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)

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;

ADOL-C/c_interface/ADOLC_TB_interface.cpp

@@ -261,6 +261,10 @@ extern "C"
     {
         return new adouble(erf(*static_cast<adouble *>(a)));
     }
+    TBAdoubleHandle tb_erfc(TBAdoubleHandle a)
+    {
+        return new adouble(erfc(*static_cast<adouble *>(a)));
+    }
 }
 
 /*

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
+#endif // ADOLC_TB_INTERFACE_H

ADOL-C/c_interface/ADOLC_TL_interface.cpp

@@ -292,4 +292,8 @@ extern "C"
     {
         return new adtl::adouble(erf(*static_cast<adtl::adouble *>(a)));
     }
+    TLAdoubleHandle tl_erfc(TLAdoubleHandle a)
+    {
+        return new adtl::adouble(erfc(*static_cast<adtl::adouble *>(a)));
+    }
 }

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
+#endif // ADOLC_TL_INTERFACE_H

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) {

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()) {

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) 

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)); }