Changeset b3796fa in sasmodels

Timestamp:

Aug 4, 2016 10:11:51 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

38ce0ab

Parents:

4fd2c63

Message:

Allow use of cephes erf() even if erf() is available in the math library

Location:

sasmodels

Files:

• kernel_header.c (modified) (2 diffs)
• models/lib/sas_erf.c (modified) (12 diffs)

sasmodels/kernel_header.c

@@ -15 +15 @@
 #  endif
    // OpenCL only has type generic math
+   #define expf exp
    #define erff erf
    #define erfcf erfc
+   // Intel CPU on Mac gives strange values for erf(), so maybe don't use it
+   //#define NEED_ERF
 #else // !USE_OPENCL
 // Use SAS_DOUBLE to force the use of double even for float kernels
@@ -123 +126 @@
 #  define M_E 2.718281828459045091
 #endif
+#ifndef M_SQRT1_2
+#  define M_SQRT1_2 0.70710678118654746
+#endif
 
 // Non-standard function library

sasmodels/models/lib/sas_erf.c

@@ -89 +89 @@
 
 #if FLOAT_SIZE>4  // DOUBLE_PRECISION
-double erf(double x);
-double erfc(double a);
+
+double cephes_erf(double x);
+double cephes_erfc(double a);
 
 constant double PD[] = {
@@ -152 +153 @@
 };
 
-double erfc(double a)
+double cephes_erfc(double a)
 {
     double MAXLOG = 88.72283905206835;
@@ -158 +159 @@
 
 
-    /*if (a < 0.0)
-        x = -a;
-    else
-        x = a;*/
-
     x = fabs(a);
 
-
     if (x < 1.0) {
-        //The line bellow is a troublemaker for GPU, so sas_erf function
-        //is explicit here for the case < 1.0
-        //return (1.0 - sas_erf(a));
+        // The line below causes problems on the GPU, so inline
+        // the erf function instead and z < 1.0.
+        //return (1.0 - cephes_erf(a));
         z = x * x;
         y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
@@ -211 +206 @@
 
 
-double erf(double x)
+double cephes_erf(double x)
 {
     double y, z;
 
     if (fabs(x) > 1.0)
-        return (1.0 - erfc(x));
+        return (1.0 - cephes_erfc(x));
 
     z = x * x;
-    #if FLOAT_SIZE>4
-        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
-    #else
-        y = x * polevl( z, TF, 6 );
-    #endif
+    y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
 
     return y;
@@ -230 +221 @@
 #else // SINGLE PRECISION
 
-double erff(double x);
-double erfcf(double a);
+float cephes_erff(float x);
+float cephes_erfcf(float a);
 
 /* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */
-constant double PF[] = {
+constant float PF[] = {
     2.326819970068386E-002,
     -1.387039388740657E-001,
@@ -247 +238 @@
 
 /* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */
-constant double RF[] = {
+constant float RF[] = {
     -1.047766399936249E+001,
     1.297719955372516E+001,
@@ -259 +250 @@
 
 /* erf(x) = x P(x^2), 0 < x < 1 */
- constant double TF[] = {
+ constant float TF[] = {
     7.853861353153693E-005,
     -8.010193625184903E-004,
@@ -270 +261 @@
 
 
-float erfcf(float a)
+float cephes_erfcf(float a)
 {
     float MAXLOG = 88.72283905206835;
@@ -327 +318 @@
 
 
-float erff(float x)
+float cephes_erff(float x)
 {
     float y, z;
@@ -333 +324 @@
     // TODO: tinycc does not support fabsf
     if (fabs(x) > 1.0)
-        return (1.0 - erfcf(x));
+        return (1.0 - cephes_erfcf(x));
 
     z = x * x;
@@ -342 +333 @@
 
 #endif // SINGLE_PRECISION
-#endif // NEED_ERF
 
 #if FLOAT_SIZE>4
+//static double sas_erf(double x) { return erf(x); }
+//static double sas_erfc(double x) { return erfc(x); }
+#define sas_erf cephes_erf
+#define sas_erfc cephes_erfc
+#else
+#define sas_erf cephes_erff
+#define sas_erfc cephes_erfcf
+#endif
+
+#else // !NEED_ERF
+
+#if FLOAT_SIZE>4
+//static double sas_erf(double x) { return erf(x); }
+//static double sas_erfc(double x) { return erfc(x); }
 #define sas_erf erf
 #define sas_erfc erfc
@@ -351 +355 @@
 #define sas_erfc erfcf
 #endif
+#endif // !NEED_ERF