source: sasmodels/sasmodels/kernel_iq.c @ d5ac45f

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

refactoring generate/kernel_template in process…

  • Property mode set to 100644
File size: 6.8 KB
Line 
1
2/*
3    ##########################################################
4    #                                                        #
5    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
6    #   !!                                              !!   #
7    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
8    #   !!                                              !!   #
9    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
10    #                                                        #
11    ##########################################################
12*/
13
14
15#define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters
16#define PD_2N 16  // PD_2N is the size of the coordination step table
17
18typedef struct {
19    int pd_par[MAX_PD];     // index of the nth polydispersity variable
20    int pd_length[MAX_PD];  // length of the nth polydispersity weight vector
21    int pd_offset[MAX_PD];  // offset of pd weights in the par & weight vector
22    int pd_stride[MAX_PD];  // stride to move to the next index at this level
23    int pd_isvol[MAX_PD];   // True if parameter is a volume weighting parameter
24    int par_offset[NPARS];  // offset of par values in the par & weight vector
25    int par_coord[NPARS];   // polydispersity coordination bitvector
26    int fast_coord_index[NPARS]; // index of the fast coordination parameters
27    int fast_coord_count;   // number of parameters coordinated with pd 1
28    int theta_var;          // id of spherical correction variable
29    int fast_theta;         // true if spherical correction depends on pd 1
30} ProblemDetails;
31
32typedef struct {
33    PARAMETER_DECL;
34} ParameterBlock;
35
36#define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC
37KERNEL
38void FULL_KERNEL_NAME(
39    int nq,                 // number of q values
40    global const ProblemDetails *problem,
41    global const double *weights,
42    global const double *pars,
43    global const double *q, // nq q values, with padding to boundary
44    global double *result,  // nq return values, again with padding
45    const double cutoff,    // cutoff in the polydispersity weight product
46    const int pd_start,     // where we are in the polydispersity loop
47    const int pd_stop,      // where we are stopping in the polydispersity loop
48    )
49{
50
51  // Storage for the current parameter values.  These will be updated as we
52  // walk the polydispersity cube.
53  local ParameterBlock local_pars;  // current parameter values
54  const double *parvec = &local_pars;  // Alias named parameters with a vector
55
56  local int offset[NPARS-2];
57
58#if defined(USE_SHORTCUT_OPTIMIZATION)
59  if (pd_length[0] == 1) {
60    // Shouldn't need to copy!!
61    for (int k=0; k < NPARS; k++) {
62      parvec[k] = pars[k+2];  // skip scale and background
63    }
64
65    #ifdef USE_OPENMP
66    #pragma omp parallel for
67    #endif
68    for (int i=0; i < nq; i++) {
69    {
70      const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS);
71      result[i] += pars[0]*scattering + pars[1];
72    }
73    return;
74  }
75#endif
76
77
78  // Since we are no longer looping over the entire polydispersity hypercube
79  // for each q, we need to track the normalization values for each q in a
80  // separate work vector.
81  double norm;   // contains sum over weights
82  double vol; // contains sum over volume
83  double norm_vol; // contains weights over volume
84
85  // Initialize the results to zero
86  if (pd_start == 0) {
87    norm_vol = 0.0;
88    norm = 0.0;
89    vol = 0.0;
90
91    #ifdef USE_OPENMP
92    #pragma omp parallel for
93    #endif
94    for (int i=0; i < nq; i++) {
95      result[i] = 0.0;
96    }
97  } else {
98    //Pulling values from previous segment
99    norm = result[nq];
100    vol = result[nq+1];
101    norm_vol = results[nq+2];
102  }
103
104  // Location in the polydispersity hypercube, one index per dimension.
105  local int pd_index[PD_MAX];
106
107  // Trigger the reset behaviour that happens at the end the fast loop
108  // by setting the initial index >= weight vector length.
109  pd_index[0] = pd_length[0];
110
111  // need product of weights at every Iq calc, so keep product of
112  // weights from the outer loops so that weight = partial_weight * fast_weight
113  double partial_weight = NAN; // product of weight w4*w3*w2 but not w1
114  double partial_volweight = NAN;
115  double weight = 1.0;        // set to 1 in case there are no weights
116  double vol_weight = 1.0;    // set to 1 in case there are no vol weights
117  double spherical_correction = 1.0;  // correction for latitude variation
118
119  // Loop over the weights then loop over q, accumulating values
120  for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
121    // check if indices need to be updated
122    if (pd_index[0] >= pd_length[0]) {
123
124      // RESET INDICES
125      pd_index[0] = loop_index%pd_length[0];
126      partial_weight = 1.0;
127      partial_volweight = 1.0;
128      for (int k=1; k < MAX_PD; k++) {
129        pd_index[k] = (loop_index%pd_length[k])/pd_stride[k];
130        const double wi = weights[pd_offset[0]+pd_index[0]];
131        partial_weight *= wi;
132        if (pd_isvol[k]) partial_volweight *= wi;
133      }
134      for (int k=0; k < NPARS; k++) {
135        int coord = par_coord[k];
136        int this_offset = par_offset[k];
137        int block_size = 1;
138        for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
139          if (coord&1) {
140              this_offset += block_size * pd_index[bit];
141              block_size *= pd_length[bit];
142          }
143          coord /= 2;
144        }
145        offset[k] = this_offset;
146        parvec[k] = pars[this_offset];
147      }
148      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
149      if (theta_var >= 0) {
150        spherical_correction = fabs(cos(M_PI_180*parvec[theta_var]));
151      }
152      if (!fast_theta) weight *= spherical_correction;
153
154    } else {
155
156      // INCREMENT INDICES
157      pd_index[0] += 1;
158      const double wi = weights[pd_offset[0]+pd_index[0]];
159      weight = partial_weight*wi;
160      if (pd_isvol[0]) vol_weight *= wi;
161      for (int k=0; k < problem->fast_coord_count; k++) {
162        parvec[ fast_coord_index[k]]
163            = pars[offset[fast_coord_index[k]] + pd_index[0]];
164      }
165      if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var]));
166
167    }
168
169    #ifdef INVALID
170    if (INVALID(local_pars)) continue;
171    #endif
172
173    if (weight > cutoff) {
174      norm += weight;
175      vol += vol_weight * CALL_VOLUME(local_pars);
176      norm_vol += vol_weight;
177
178      #ifdef USE_OPENMP
179      #pragma omp parallel for
180      #endif
181      for (int i=0; i < nq; i++) {
182      {
183        const double scattering = CALL_IQ(q, nq, i, local_pars);
184        result[i] += weight*scattering;
185      }
186  }
187  //Makes a normalization avialable for the next round
188  result[nq] = norm;
189  result[nq+1] = vol;
190  result[nq+2] = norm_vol;
191
192  //End of the PD loop we can normalize
193  if (pd_stop == pd_stride[MAX_PD-1]) {}
194    #ifdef USE_OPENMP
195    #pragma omp parallel for
196    #endif
197    for (int i=0; i < nq; i++) {
198      if (vol*norm_vol != 0.0) {
199        result[i] *= norm_vol/vol;
200      }
201      result[i] = pars[0]*result[i]/norm + pars[1];
202    }
203  }
204}
Note: See TracBrowser for help on using the repository browser.