source: sasmodels/sasmodels/kernel_iq.cl @ 56547a8

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 56547a8 was 56547a8, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

remove build error for Intel HD 4600 on windows

  • Property mode set to 100644
File size: 11.0 KB
Line 
1
2/*
3    ##########################################################
4    #                                                        #
5    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
6    #   !!                                              !!   #
7    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
8    #   !!                                              !!   #
9    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
10    #                                                        #
11    ##########################################################
12*/
13
14#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
15#define _PAR_BLOCK_
16
17typedef struct {
18#if MAX_PD > 0
19    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
20    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
21    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
22    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
23#endif // MAX_PD > 0
24    int32_t pd_prod;            // total number of voxels in hypercube
25    int32_t pd_sum;             // total length of the weights vector
26    int32_t num_active;         // number of non-trivial pd loops
27    int32_t theta_par;          // id of spherical correction variable
28} ProblemDetails;
29
30typedef struct {
31    PARAMETER_TABLE
32} ParameterBlock;
33#endif // _PAR_BLOCK_
34
35
36#if defined(MAGNETIC) && NUM_MAGNETIC>0
37
38// Return value restricted between low and high
39static double clip(double value, double low, double high)
40{
41  return (value < low ? low : (value > high ? high : value));
42}
43
44// Compute spin cross sections given in_spin and out_spin
45// To convert spin cross sections to sld b:
46//     uu * (sld - m_sigma_x);
47//     dd * (sld + m_sigma_x);
48//     ud * (m_sigma_y + 1j*m_sigma_z);
49//     du * (m_sigma_y - 1j*m_sigma_z);
50static void set_spins(double in_spin, double out_spin, double spins[4])
51{
52  in_spin = clip(in_spin, 0.0, 1.0);
53  out_spin = clip(out_spin, 0.0, 1.0);
54  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
55  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
56  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
57  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
58}
59
60static double mag_sld(double qx, double qy, double p,
61                       double mx, double my, double sld)
62{
63    const double perp = qy*mx - qx*my;
64    return sld + perp*p;
65}
66
67#endif // MAGNETIC
68
69kernel
70void KERNEL_NAME(
71    int32_t nq,                 // number of q values
72    const int32_t pd_start,     // where we are in the polydispersity loop
73    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
74    global const ProblemDetails *details,
75    global const double *values,
76    global const double *q, // nq q values, with padding to boundary
77    global double *result,  // nq+1 return values, again with padding
78    const double cutoff     // cutoff in the polydispersity weight product
79    )
80{
81
82  // who we are and what element we are working with
83  const int q_index = get_global_id(0);
84  if (q_index >= nq) return;
85
86  // Storage for the current parameter values.  These will be updated as we
87  // walk the polydispersity cube.  local_values will be aliased to pvec.
88  ParameterBlock local_values;
89  double *pvec = (double *)&local_values;
90
91  // Fill in the initial variables
92  for (int i=0; i < NUM_PARS; i++) {
93    pvec[i] = values[2+i];
94//if (q_index==0) printf("p%d = %g\n",i, pvec[i]);
95  }
96
97#if defined(MAGNETIC) && NUM_MAGNETIC>0
98  // Location of the sld parameters in the parameter pvec.
99  // These parameters are updated with the effective sld due to magnetism.
100  #if NUM_MAGNETIC > 3
101  const int32_t slds[] = { MAGNETIC_PARS };
102  #endif
103
104  // TODO: could precompute these outside of the kernel.
105  // Interpret polarization cross section.
106  //     up_frac_i = values[NUM_PARS+2];
107  //     up_frac_f = values[NUM_PARS+3];
108  //     up_angle = values[NUM_PARS+4];
109  double spins[4];
110  double cos_mspin, sin_mspin;
111  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
112  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
113#endif // MAGNETIC
114
115  double pd_norm, this_result;
116  if (pd_start == 0) {
117    pd_norm = this_result = 0.0;
118  } else {
119    pd_norm = result[nq];
120    this_result = result[q_index];
121  }
122//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result);
123
124#if MAX_PD>0
125  global const double *pd_value = values + NUM_VALUES + 2;
126  global const double *pd_weight = pd_value + details->pd_sum;
127#endif
128
129  // Jump into the middle of the polydispersity loop
130#if MAX_PD>4
131  int n4=details->pd_length[4];
132  int i4=(pd_start/details->pd_stride[4])%n4;
133  const int p4=details->pd_par[4];
134  global const double *v4 = pd_value + details->pd_offset[4];
135  global const double *w4 = pd_weight + details->pd_offset[4];
136#endif
137#if MAX_PD>3
138  int n3=details->pd_length[3];
139  int i3=(pd_start/details->pd_stride[3])%n3;
140  const int p3=details->pd_par[3];
141  global const double *v3 = pd_value + details->pd_offset[3];
142  global const double *w3 = pd_weight + details->pd_offset[3];
143//if (q_index==0) printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES);
144#endif
145#if MAX_PD>2
146  int n2=details->pd_length[2];
147  int i2=(pd_start/details->pd_stride[2])%n2;
148  const int p2=details->pd_par[2];
149  global const double *v2 = pd_value + details->pd_offset[2];
150  global const double *w2 = pd_weight + details->pd_offset[2];
151#endif
152#if MAX_PD>1
153  int n1=details->pd_length[1];
154  int i1=(pd_start/details->pd_stride[1])%n1;
155  const int p1=details->pd_par[1];
156  global const double *v1 = pd_value + details->pd_offset[1];
157  global const double *w1 = pd_weight + details->pd_offset[1];
158#endif
159#if MAX_PD>0
160  int n0=details->pd_length[0];
161  int i0=(pd_start/details->pd_stride[0])%n0;
162  const int p0=details->pd_par[0];
163  global const double *v0 = pd_value + details->pd_offset[0];
164  global const double *w0 = pd_weight + details->pd_offset[0];
165#endif
166
167
168  double spherical_correction=1.0;
169  const int theta_par = details->theta_par;
170#if MAX_PD>0
171  const int fast_theta = (theta_par == p0);
172  const int slow_theta = (theta_par >= 0 && !fast_theta);
173#else
174  const int slow_theta = (theta_par >= 0);
175#endif
176
177  int step = pd_start;
178
179
180#if MAX_PD>4
181  const double weight5 = 1.0;
182  while (i4 < n4) {
183    pvec[p4] = v4[i4];
184    double weight4 = w4[i4] * weight5;
185//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4);
186#elif MAX_PD>3
187    const double weight4 = 1.0;
188#endif
189#if MAX_PD>3
190  while (i3 < n3) {
191    pvec[p3] = v3[i3];
192    double weight3 = w3[i3] * weight4;
193//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3);
194#elif MAX_PD>2
195    const double weight3 = 1.0;
196#endif
197#if MAX_PD>2
198  while (i2 < n2) {
199    pvec[p2] = v2[i2];
200    double weight2 = w2[i2] * weight3;
201//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2);
202#elif MAX_PD>1
203    const double weight2 = 1.0;
204#endif
205#if MAX_PD>1
206  while (i1 < n1) {
207    pvec[p1] = v1[i1];
208    double weight1 = w1[i1] * weight2;
209//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1);
210#elif MAX_PD>0
211    const double weight1 = 1.0;
212#endif
213    if (slow_theta) { // Theta is not in inner loop
214      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
215    }
216#if MAX_PD>0
217  while(i0 < n0) {
218    pvec[p0] = v0[i0];
219    double weight0 = w0[i0] * weight1;
220//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0);
221    if (fast_theta) { // Theta is in inner loop
222      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6);
223    }
224#else
225    const double weight0 = 1.0;
226#endif
227
228//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); }
229//if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
230
231    #ifdef INVALID
232    if (!INVALID(local_values))
233    #endif
234    {
235      // Accumulate I(q)
236      // Note: weight==0 must always be excluded
237      if (weight0 > cutoff) {
238        // spherical correction is set at a minimum of 1e-6, otherwise there
239        // would be problems looking at models with theta=90.
240        const double weight = weight0 * spherical_correction;
241        pd_norm += weight * CALL_VOLUME(local_values);
242
243#if defined(MAGNETIC) && NUM_MAGNETIC > 0
244        const double qx = q[2*q_index];
245        const double qy = q[2*q_index+1];
246        const double qsq = qx*qx + qy*qy;
247
248        // Constant across orientation, polydispersity for given qx, qy
249        double scattering = 0.0;
250        // TODO: what is the magnetic scattering at q=0
251        if (qsq > 1.e-16) {
252          double p[4];  // spin_i, spin_f
253          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
254          p[3] = -p[0];
255          p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
256
257          for (int index=0; index<4; index++) {
258            const double xs = spins[index];
259            if (xs > 1.e-8) {
260              const int spin_flip = (index==1) || (index==2);
261              const double pk = p[index];
262              for (int axis=0; axis<=spin_flip; axis++) {
263                #define M1 NUM_PARS+5
264                #define M2 NUM_PARS+8
265                #define M3 NUM_PARS+13
266                #define SLD(_M_offset, _sld_offset) \
267                    pvec[_sld_offset] = xs * (axis \
268                    ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
269                    : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
270                              (spin_flip ? 0.0 : values[_sld_offset+2])))
271                #if NUM_MAGNETIC==1
272                    SLD(M1, MAGNETIC_PAR1);
273                #elif NUM_MAGNETIC==2
274                    SLD(M1, MAGNETIC_PAR1);
275                    SLD(M2, MAGNETIC_PAR2);
276                #elif NUM_MAGNETIC==3
277                    SLD(M1, MAGNETIC_PAR1);
278                    SLD(M2, MAGNETIC_PAR2);
279                    SLD(M3, MAGNETIC_PAR3);
280                #else
281                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
282                    SLD(M1+3*sk, slds[sk]);
283                }
284                #endif
285                scattering += CALL_IQ(q, q_index, local_values);
286              }
287            }
288          }
289        }
290#else  // !MAGNETIC
291        const double scattering = CALL_IQ(q, q_index, local_values);
292#endif // !MAGNETIC
293        this_result += weight * scattering;
294      }
295    }
296    ++step;
297#if MAX_PD>0
298    if (step >= pd_stop) break;
299    ++i0;
300  }
301  i0 = 0;
302#endif
303#if MAX_PD>1
304    if (step >= pd_stop) break;
305    ++i1;
306  }
307  i1 = 0;
308#endif
309#if MAX_PD>2
310    if (step >= pd_stop) break;
311    ++i2;
312  }
313  i2 = 0;
314#endif
315#if MAX_PD>3
316    if (step >= pd_stop) break;
317    ++i3;
318  }
319  i3 = 0;
320#endif
321#if MAX_PD>4
322    if (step >= pd_stop) break;
323    ++i4;
324  }
325  i4 = 0;
326#endif
327
328//if (q_index==0) printf("res: %g/%g\n", this_result, pd_norm);
329  // Remember the current result and the updated norm.
330  result[q_index] = this_result;
331  if (q_index == 0) result[nq] = pd_norm;
332}
Note: See TracBrowser for help on using the repository browser.