Changeset f2f67a6 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Apr 15, 2016 7:26:24 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:

ae2b6b5, 38a9b07, eb97b11

Parents:

0ff62d4

Message:

reenable opencl; works on cpu but not gpu

File:

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

sasmodels/kernel_iq.c

@@ -54 +54 @@
   local ParameterBlock local_values;  // current parameter values
   double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
+  double norm;
+
+  // number of active loops
+  const int num_active = details->num_active;
 
   // Fill in the initial variables
@@ -63 +67 @@
   }
 
+  // Monodisperse computation
+  if (num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+    norm = CALL_VOLUME(local_values);
+
+    const double scale = values[0];
+    const double background = values[1];
+    // result[nq] = norm; // Total volume normalization
+
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      double scattering = CALL_IQ(q, i, local_values);
+      result[i] = (norm>0. ? scale*scattering/norm + background : background);
+    }
+    return;
+  }
+
+#if MAX_PD > 0
   // If it is the first round initialize the result to zero, otherwise
   // assume that the previous result has been passed back.
@@ -75 +101 @@
       result[i] = 0.0;
     }
-  }
-
-  // Monodisperse computation
-  if (details->num_active == 0) {
-    #ifdef INVALID
-    if (INVALID(local_values)) { return; }
-    #endif
-
-    const double norm = CALL_VOLUME(local_values);
-    const double scale = values[0];
-    const double background = values[1];
-    #ifdef USE_OPENMP
-    #pragma omp parallel for
-    #endif
-    result[nq] = norm; // Total volume normalization
-    for (int i=0; i < nq; i++) {
-      double scattering = CALL_IQ(q, i, local_values);
-      result[i] = (norm>0. ? scale*scattering/norm + background : background);
-    }
-    return;
-  }
-
-#if MAX_PD > 0
-  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
-  // Since we are no longer looping over the entire polydispersity hypercube
-  // for each q, we need to track the normalization values between calls.
-  double norm = 0.0;
+    norm = 0.0;
+  } else {
+    norm = result[nq];
+  }
 
   // need product of weights at every Iq calc, so keep product of
@@ -119 +122 @@
   pd_index[0] = fast_length;
 
+  // Number of coordinated indices
+  const int num_coord = details->num_coord;
+
   // Loop over the weights then loop over q, accumulating values
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
-    // check if indices need to be updated
+    // check if fast loop needs to be reset
     if (pd_index[0] == fast_length) {
-      //printf("should be here with %d active\n", details->num_active);
+      //printf("should be here with %d active\n", num_active);
 
       // Compute position in polydispersity hypercube
-      for (int k=0; k < details->num_active; k++) {
+      for (int k=0; k < num_active; k++) {
         pd_index[k] = (loop_index/details->pd_stride[k])%details->pd_length[k];
         //printf("pd_index[%d] = %d\n",k,pd_index[k]);
@@ -134 +140 @@
       partial_weight = 1.0;
       //printf("partial weight %d: ", loop_index);
-      for (int k=1; k < details->num_active; k++) {
+      for (int k=1; k < num_active; k++) {
         double wi = weights[details->pd_offset[k] + pd_index[k]];
         //printf("pd[%d]=par[%d]=%g ", k, details->pd_par[k], wi);
@@ -143 +149 @@
       // Update parameter offsets in weight vector
       //printf("slow %d: ", loop_index);
-      for (int k=0; k < details->num_coord; k++) {
+      for (int k=0; k < num_coord; k++) {
         int par = details->par_coord[k];
         int coord = details->pd_coord[k];
@@ -160 +166 @@
         // if theta is not coordinated with fast index, precompute spherical correction
         if (par == details->theta_par && !(details->par_coord[k]&1)) {
-          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
         }
       }
@@ -170 +176 @@
     double weight = partial_weight*wi;
     //printf("fast %d: ", loop_index);
-    for (int k=0; k < details->num_coord; k++) {
+    for (int k=0; k < num_coord; k++) {
       if (details->pd_coord[k]&1) {
         const int par = details->par_coord[k];
@@ -177 +183 @@
         // if theta is coordinated with fast index, compute spherical correction each time
         if (par == details->theta_par) {
-          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
         }
       }
@@ -205 +211 @@
   }
 
-  // Accumulate norm.
-  result[nq] += norm;
-
   // End of the PD loop we can normalize
   if (pd_stop >= details->total_pd) {
@@ -219 +222 @@
     }
   }
+
+  // Remember the updated norm.
+  result[nq] = norm;
 #endif // MAX_PD > 0
 }