Changeset 380e8c9 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Mar 24, 2016 10:56:44 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:

151f3bc

Parents:

60eab2a

Message:

progress on new polydispersity loop, but still broken

File:

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

sasmodels/kernel_iq.c

@@ -54 +54 @@
   double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
 
-#if MAX_PD > 0
-  if (problem->pd_length[0] == 1) {
-#endif // MAX_PD > 0
+  // Monodisperse computation
+  if (pd_stop == 1) {
     // Shouldn't need to copy!!
     for (int k=0; k < NPARS; k++) {
@@ -72 +71 @@
     }
     return;
+  }
+
 #if MAX_PD > 0
-  }
-
-  // polydispersity loop index positions
-  local int offset[NPARS];  // NPARS excludes scale/background
-
-  printf("Entering polydispersity\n");
+
+printf("Entering polydispersity\n");
   // Since we are no longer looping over the entire polydispersity hypercube
   // for each q, we need to track the normalization values for each q in a
@@ -107 +104 @@
   // Location in the polydispersity hypercube, one index per dimension.
   local int pd_index[MAX_PD];
+
+  // polydispersity loop index positions
+  local int offset[NPARS];  // NPARS excludes scale/background
 
   // Trigger the reset behaviour that happens at the end the fast loop
@@ -132 +132 @@
       for (int k=1; k < MAX_PD; k++) {
         pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
-        const double wi = weights[problem->pd_offset[0]+pd_index[0]];
+        const double wi = weights[problem->pd_offset[k]+pd_index[k]];
         partial_weight *= wi;
         if (problem->pd_isvol[k]) partial_volweight *= wi;
       }
+      printf("slow %d: ", loop_index);
       for (int k=0; k < NPARS; k++) {
         int coord = problem->par_coord[k];
@@ -149 +150 @@
         offset[k] = this_offset;
         pvec[k] = values[this_offset];
-      }
+        printf("p[%d]=v[%d]=%g ", k, offset[k], pvec[k]);
+      }
+      printf("\n");
       weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
       if (problem->theta_par >= 0) {
@@ -157 +160 @@
         weight *= spherical_correction;
       }
+      pd_index[0] += 1;
 
     } else {
 
       // INCREMENT INDICES
-      pd_index[0] += 1;
       const double wi = weights[problem->pd_offset[0]+pd_index[0]];
       weight = partial_weight*wi;
       if (problem->pd_isvol[0]) vol_weight *= wi;
+      printf("fast %d: ", loop_index);
       for (int k=0; k < problem->fast_coord_count; k++) {
-         printf("fast loop %d coord %d idx %d ?%d offset %d %g\n",
-         k,problem->fast_coord_pars[k],pd_index[0],
-        problem->fast_coord_pars[k],
-            offset[problem->fast_coord_pars[k]],
-            values[offset[problem->fast_coord_pars[k]]]
-          );
-        pvec[problem->fast_coord_pars[k]]
-            = values[offset[problem->fast_coord_pars[k]]++];
-
-      }
-      if (problem->theta_par ==problem->pd_par[0]) {
+        const int pindex = problem->fast_coord_pars[k];
+        pvec[pindex] = values[++offset[pindex]];
+        printf("p[%d]=v[%d]=%g ", pindex, offset[pindex], pvec[pindex]);
+      }
+      printf("\n");
+      if (problem->theta_par == problem->pd_par[0]) {
         weight *= fabs(cos(M_PI_180*pvec[problem->theta_par]));
       }
+      pd_index[0] += 1;
     }
     #ifdef INVALID
@@ -201 +201 @@
   }
 
-  //Makes a normalization available for the next round
+  // Make normalization available for the next round
   result[nq] = norm;
   result[nq+1] = vol;
   result[nq+2] = norm_vol;
 
-  //End of the PD loop we can normalize
+  // End of the PD loop we can normalize
   if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
     #ifdef USE_OPENMP