source: sasmodels/sasmodels/src @ ba32cdd

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ba32cdd was 9c490bc, checked in by wojciech, 8 years ago

Some intermediate changes
  • Property mode set to 100644
File size: 10.4 KB
Line 
1// GENERATED CODE --- DO NOT EDIT ---
2// Code is produced by sasmodels.gen from sasmodels/models/MODEL.c
3
4#ifdef __OPENCL_VERSION__
5# define USE_OPENCL
6#endif
7
8#define USE_KAHAN_SUMMATION 0
9
10// If opencl is not available, then we are compiling a C function
11// Note: if using a C++ compiler, then define kernel as extern "C"
12#ifndef USE_OPENCL
13#  ifdef __cplusplus
14      #include <cstdio>
15      #include <cmath>
16      using namespace std;
17      #if defined(_MSC_VER)
18         #include <limits>
19         #include <float.h>
20         #define kernel extern "C" __declspec( dllexport )
21         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
22         inline double fmin(double x, double y) { return x>y ? y : x; }
23         inline double fmax(double x, double y) { return x<y ? y : x; }
24         inline double isnan(double x) { return _isnan(x); }
25         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
26         static double cephes_expm1(double x) {
27            // Adapted from the cephes math library.
28            // Copyright 1984 - 1992 by Stephen L. Moshier
29            if (x != x || x == 0.0) {
30               return x; // NaN and +/- 0
31            } else if (x < -0.5 || x > 0.5) {
32               return exp(x) - 1.0;
33            } else {
34               const double xsq = x*x;
35               const double p = (((
36                  +1.2617719307481059087798E-4)*xsq
37                  +3.0299440770744196129956E-2)*xsq
38                  +9.9999999999999999991025E-1);
39               const double q = ((((
40                  +3.0019850513866445504159E-6)*xsq
41                  +2.5244834034968410419224E-3)*xsq
42                  +2.2726554820815502876593E-1)*xsq
43                  +2.0000000000000000000897E0);
44               double r = x * p;
45               r =  r / (q - r);
46               return r+r;
47             }
48         }
49         #define expm1 cephes_expm1
50     #else
51         #define kernel extern "C"
52     #endif
53     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
54#  else
55     #include <stdio.h>
56     #include <tgmath.h> // C99 type-generic math, so sin(float) => sinf
57     // MSVC doesn't support C99, so no need for dllexport on C99 branch
58     #define kernel
59     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
60#  endif
61#  define global
62#  define local
63#  define constant const
64// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
65// OpenCL pown(a,b) = C99 pow(a,b), b integer
66#  define powr(a,b) pow(a,b)
67#  define pown(a,b) pow(a,b)
68#else
69#  if defined(USE_SINCOS)
70#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
71#  else
72#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
73#  endif
74#endif
75
76// Standard mathematical constants:
77//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
78//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
79// OpenCL defines M_constant_F for float constants, and nothing if double
80// is not enabled on the card, which is why these constants may be missing
81#ifndef M_PI
82#  define M_PI 3.141592653589793
83#endif
84#ifndef M_PI_2
85#  define M_PI_2 1.570796326794897
86#endif
87#ifndef M_PI_4
88#  define M_PI_4 0.7853981633974483
89#endif
90#ifndef M_E
91#  define M_E 2.718281828459045091
92#endif
93
94// Non-standard function library
95// pi/180, used for converting between degrees and radians
96// 4/3 pi for computing sphere volumes
97// square and cube for computing squares and cubes
98#ifndef M_PI_180
99#  define M_PI_180 0.017453292519943295
100#endif
101#ifndef M_4PI_3
102#  define M_4PI_3 4.18879020478639
103#endif
104//inline double square(double x) { return pow(x,2.0); }
105//inline double square(double x) { return pown(x,2); }
106inline double square(double x) { return x*x; }
107inline double cube(double x) { return x*x*x; }
108inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
109
110
111#define VOLUME_PARAMETER_DECLARATIONS void
112#define IQ_KERNEL_NAME bessel_Iq
113#define IQ_PARAMETERS ignored
114#define IQ_FIXED_PARAMETER_DECLARATIONS const double scale, \
115    const double background, \
116    const double ignored
117#define IQ_PARAMETER_DECLARATIONS double ignored
118#define IQXY_KERNEL_NAME bessel_Iqxy
119#define IQXY_PARAMETERS ignored
120#define IQXY_FIXED_PARAMETER_DECLARATIONS const double scale, \
121    const double background, \
122    const double ignored
123#define IQXY_PARAMETER_DECLARATIONS double ignored
124
125/*
126The wrapper for gamma function from OpenCL and standard libraries
127The OpenCL gamma function fails misserably on values lower than 1.0
128while works fine on larger values.
129We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
130to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1)
131*/
132
133
134inline double sas_gamma( double x) { return tgamma(x+1)/x; }
135
136double form_volume(VOLUME_PARAMETER_DECLARATIONS);
137double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
138   
139   
140}
141
142
143double Iq(double q, IQ_PARAMETER_DECLARATIONS);
144double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
145   
146    return sas_gamma(q);
147   
148}
149
150
151double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
152double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
153   
154    // never called since no orientation or magnetic parameters.
155    //return -1.0;
156   
157}
158
159
160/*
161    ##########################################################
162    #                                                        #
163    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
164    #   !!                                              !!   #
165    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
166    #   !!                                              !!   #
167    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
168    #                                                        #
169    ##########################################################
170*/
171
172#ifdef IQ_KERNEL_NAME
173kernel void IQ_KERNEL_NAME(
174    global const double *q,
175    global double *result,
176    const int Nq,
177#ifdef IQ_OPEN_LOOPS
178  #ifdef USE_OPENCL
179    global double *loops_g,
180  #endif
181    local double *loops,
182    const double cutoff,
183    IQ_DISPERSION_LENGTH_DECLARATIONS,
184#endif
185    IQ_FIXED_PARAMETER_DECLARATIONS
186    )
187{
188#ifdef USE_OPENCL
189  #ifdef IQ_OPEN_LOOPS
190  // copy loops info to local memory
191  event_t e = async_work_group_copy(loops, loops_g, (IQ_DISPERSION_LENGTH_SUM)*2, 0);
192  wait_group_events(1, &e);
193  #endif
194
195  int i = get_global_id(0);
196  if (i < Nq)
197#else
198  #pragma omp parallel for
199  for (int i=0; i < Nq; i++)
200#endif
201  {
202    const double qi = q[i];
203#ifdef IQ_OPEN_LOOPS
204    double ret=0.0, norm=0.0;
205  #ifdef VOLUME_PARAMETERS
206    double vol=0.0, norm_vol=0.0;
207  #endif
208    IQ_OPEN_LOOPS
209    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
210    //  const double radius = loops[2*(radius_i)];
211    //  const double radius_w = loops[2*(radius_i)+1];
212
213    const double weight = IQ_WEIGHT_PRODUCT;
214    if (weight > cutoff) {
215      const double scattering = Iq(qi, IQ_PARAMETERS);
216      // allow kernels to exclude invalid regions by returning NaN
217      if (!isnan(scattering)) {
218        ret += weight*scattering;
219        norm += weight;
220      #ifdef VOLUME_PARAMETERS
221        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
222        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
223        norm_vol += vol_weight;
224      #endif
225      }
226    //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
227    }
228    IQ_CLOSE_LOOPS
229  #ifdef VOLUME_PARAMETERS
230    if (vol*norm_vol != 0.0) {
231      ret *= norm_vol/vol;
232    }
233  #endif
234    result[i] = scale*ret/norm+background;
235#else
236    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background;
237#endif
238  }
239}
240#endif
241
242
243#ifdef IQXY_KERNEL_NAME
244kernel void IQXY_KERNEL_NAME(
245    global const double *qx,
246    global const double *qy,
247    global double *result,
248    const int Nq,
249#ifdef IQXY_OPEN_LOOPS
250  #ifdef USE_OPENCL
251    global double *loops_g,
252  #endif
253    local double *loops,
254    const double cutoff,
255    IQXY_DISPERSION_LENGTH_DECLARATIONS,
256#endif
257    IQXY_FIXED_PARAMETER_DECLARATIONS
258    )
259{
260#ifdef USE_OPENCL
261  #ifdef IQXY_OPEN_LOOPS
262  // copy loops info to local memory
263  event_t e = async_work_group_copy(loops, loops_g, (IQXY_DISPERSION_LENGTH_SUM)*2, 0);
264  wait_group_events(1, &e);
265  #endif
266
267  int i = get_global_id(0);
268  if (i < Nq)
269#else
270  #pragma omp parallel for
271  for (int i=0; i < Nq; i++)
272#endif
273  {
274    const double qxi = qx[i];
275    const double qyi = qy[i];
276    #if USE_KAHAN_SUMMATION
277    double accumulated_error = 0.0;
278    #endif
279#ifdef IQXY_OPEN_LOOPS
280    double ret=0.0, norm=0.0;
281    #ifdef VOLUME_PARAMETERS
282    double vol=0.0, norm_vol=0.0;
283    #endif
284    IQXY_OPEN_LOOPS
285    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
286    //  const double radius = loops[2*(radius_i)];
287    //  const double radius_w = loops[2*(radius_i)+1];
288
289    const double weight = IQXY_WEIGHT_PRODUCT;
290    if (weight > cutoff) {
291
292      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
293      if (!isnan(scattering)) { // if scattering is bad, exclude it from sum
294      //if (scattering >= 0.0) { // scattering cannot be negative
295        // TODO: use correct angle for spherical correction
296        // Definition of theta and phi are probably reversed relative to the
297        // equation which gave rise to this correction, leading to an
298        // attenuation of the pattern as theta moves through pi/2.  Either
299        // reverse the meanings of phi and theta in the forms, or use phi
300        // rather than theta in this correction.  Current code uses cos(theta)
301        // so that values match those of sasview.
302      #if defined(IQXY_HAS_THETA) // && 0
303        const double spherical_correction
304          = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0);
305        const double next = spherical_correction * weight * scattering;
306      #else
307        const double next = weight * scattering;
308      #endif
309      #if USE_KAHAN_SUMMATION
310        const double y = next - accumulated_error;
311        const double t = ret + y;
312        accumulated_error = (t - ret) - y;
313        ret = t;
314      #else
315        ret += next;
316      #endif
317        norm += weight;
318      #ifdef VOLUME_PARAMETERS
319        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
320        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
321        norm_vol += vol_weight;
322      #endif
323      }
324      //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
325    }
326    IQXY_CLOSE_LOOPS
327  #ifdef VOLUME_PARAMETERS
328    if (vol*norm_vol != 0.0) {
329      ret *= norm_vol/vol;
330    }
331  #endif
332    result[i] = scale*ret/norm+background;
333#else
334    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background;
335#endif
336  }
337}
338#endif
339
Note: See TracBrowser for help on using the repository browser.