Changeset 03cac08 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Mar 20, 2016 9:44:11 PM (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:

303d8d6

Parents:

d5ac45f

Message:

new generator produces code that compiles

File:

• sasmodels/kernel_iq.c (modified) (11 diffs)

sasmodels/kernel_iq.c

@@ -12 +12 @@
 */
 
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
 
 #define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters
-#define PD_2N 16  // PD_2N is the size of the coordination step table
 
 typedef struct {
@@ -31 +32 @@
 
 typedef struct {
-    PARAMETER_DECL;
+    PARAMETER_TABLE;
 } ParameterBlock;
-
-#define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC
-KERNEL
-void FULL_KERNEL_NAME(
+#endif
+
+
+kernel
+void KERNEL_NAME(
     int nq,                 // number of q values
     global const ProblemDetails *problem,
@@ -42 +44 @@
     global const double *pars,
     global const double *q, // nq q values, with padding to boundary
-    global double *result,  // nq return values, again with padding
+    global double *result,  // nq+3 return values, again with padding
     const double cutoff,    // cutoff in the polydispersity weight product
     const int pd_start,     // where we are in the polydispersity loop
-    const int pd_stop,      // where we are stopping in the polydispersity loop
+    const int pd_stop       // where we are stopping in the polydispersity loop
     )
 {
@@ -52 +54 @@
   // walk the polydispersity cube.
   local ParameterBlock local_pars;  // current parameter values
-  const double *parvec = &local_pars;  // Alias named parameters with a vector
+  double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
 
   local int offset[NPARS-2];
@@ -60 +62 @@
     // Shouldn't need to copy!!
     for (int k=0; k < NPARS; k++) {
-      parvec[k] = pars[k+2];  // skip scale and background
+      pvec[k] = pars[k+2];  // skip scale and background
     }
 
@@ -68 +70 @@
     for (int i=0; i < nq; i++) {
     {
-      const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS);
+      const double scattering = CALL_IQ(q, i, local_pars);
       result[i] += pars[0]*scattering + pars[1];
     }
@@ -99 +101 @@
     norm = result[nq];
     vol = result[nq+1];
-    norm_vol = results[nq+2];
+    norm_vol = result[nq+2];
   }
 
   // Location in the polydispersity hypercube, one index per dimension.
-  local int pd_index[PD_MAX];
+  local int pd_index[MAX_PD];
 
   // Trigger the reset behaviour that happens at the end the fast loop
   // by setting the initial index >= weight vector length.
-  pd_index[0] = pd_length[0];
+  pd_index[0] = problem->pd_length[0];
+
 
   // need product of weights at every Iq calc, so keep product of
@@ -120 +123 @@
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
     // check if indices need to be updated
-    if (pd_index[0] >= pd_length[0]) {
+    if (pd_index[0] >= problem->pd_length[0]) {
 
       // RESET INDICES
-      pd_index[0] = loop_index%pd_length[0];
+      pd_index[0] = loop_index%problem->pd_length[0];
       partial_weight = 1.0;
       partial_volweight = 1.0;
       for (int k=1; k < MAX_PD; k++) {
-        pd_index[k] = (loop_index%pd_length[k])/pd_stride[k];
-        const double wi = weights[pd_offset[0]+pd_index[0]];
+        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
+        const double wi = weights[problem->pd_offset[0]+pd_index[0]];
         partial_weight *= wi;
-        if (pd_isvol[k]) partial_volweight *= wi;
+        if (problem->pd_isvol[k]) partial_volweight *= wi;
       }
       for (int k=0; k < NPARS; k++) {
-        int coord = par_coord[k];
-        int this_offset = par_offset[k];
+        int coord = problem->par_coord[k];
+        int this_offset = problem->par_offset[k];
         int block_size = 1;
         for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
           if (coord&1) {
               this_offset += block_size * pd_index[bit];
-              block_size *= pd_length[bit];
+              block_size *= problem->pd_length[bit];
           }
           coord /= 2;
         }
         offset[k] = this_offset;
-        parvec[k] = pars[this_offset];
-      }
-      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
-      if (theta_var >= 0) {
-        spherical_correction = fabs(cos(M_PI_180*parvec[theta_var]));
-      }
-      if (!fast_theta) weight *= spherical_correction;
+        pvec[k] = pars[this_offset];
+      }
+      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
+      if (problem->theta_var >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
+      if (!problem->fast_theta) {
+        weight *= spherical_correction;
+      }
 
     } else {
@@ -156 +161 @@
       // INCREMENT INDICES
       pd_index[0] += 1;
-      const double wi = weights[pd_offset[0]+pd_index[0]];
+      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
       weight = partial_weight*wi;
-      if (pd_isvol[0]) vol_weight *= wi;
+      if (problem->pd_isvol[0]) vol_weight *= wi;
       for (int k=0; k < problem->fast_coord_count; k++) {
-        parvec[ fast_coord_index[k]]
-            = pars[offset[fast_coord_index[k]] + pd_index[0]];
-      }
-      if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var]));
-
+        pvec[problem->fast_coord_index[k]]
+            = pars[offset[problem->fast_coord_index[k]]++];
+      }
+      if (problem->fast_theta) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
     }
 
@@ -180 +186 @@
       #endif
       for (int i=0; i < nq; i++) {
-      {
-        const double scattering = CALL_IQ(q, nq, i, local_pars);
+        const double scattering = CALL_IQ(q, i, local_pars);
         result[i] += weight*scattering;
       }
+    }
   }
   //Makes a normalization avialable for the next round
@@ -191 +197 @@
 
   //End of the PD loop we can normalize
-  if (pd_stop == pd_stride[MAX_PD-1]) {}
+  if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
     #ifdef USE_OPENMP
     #pragma omp parallel for