Changeset 149eb53 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Oct 25, 2018 1:43:27 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

6d5601c

Parents:

07646b6 (diff), 508475a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'cuda-test' into beta_approx

File:

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

sasmodels/kernel_iq.c

@@ -26 +26 @@
 //  MAGNETIC_PARS : a comma-separated list of indices to the sld
 //      parameters in the parameter table.
-//  CALL_VOLUME(table) : call the form volume function
+//  CALL_VOLUME(form, shell, table) : assign form and shell values
+//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
 //  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
+//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs.
+//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data.
 //  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
 //  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
@@ -36 +39 @@
 //  PROJECTION : equirectangular=1, sinusoidal=2
 //      see explore/jitter.py for definitions.
+
 
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
@@ -270 +274 @@
 #endif // _QABC_SECTION
 
-
 // ==================== KERNEL CODE ========================
-
 kernel
 void KERNEL_NAME(
-    int32_t nq,                 // number of q values
-    const int32_t pd_start,     // where we are in the dispersity loop
-    const int32_t pd_stop,      // where we are stopping in the dispersity loop
+    int32_t nq,                   // number of q values
+    const int32_t pd_start,       // where we are in the dispersity loop
+    const int32_t pd_stop,        // where we are stopping in the dispersity loop
     pglobal const ProblemDetails *details,
-    pglobal const double *values,
-    pglobal const double *q, // nq q values, with padding to boundary
-    pglobal double *result,  // nq+1 return values, again with padding
-    const double cutoff     // cutoff in the dispersity weight product
+    pglobal const double *values, // parameter values and distributions
+    pglobal const double *q,      // nq q values, with padding to boundary
+    pglobal double *result,       // nq+1 return values, again with padding
+    const double cutoff,          // cutoff in the dispersity weight product
+    int32_t effective_radius_type // which effective radius to compute
     )
 {
@@ -342 +345 @@
   //
   // The code differs slightly between opencl and dll since opencl is only
-  // seeing one q value (stored in the variable "this_result") while the dll
+  // seeing one q value (stored in the variable "this_F2") while the dll
   // version must loop over all q.
+  #if defined(CALL_FQ)
+    double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
+    double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+    double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+    double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]);
+  #else
+    double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]);
+    double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]);
+    double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]);
+  #endif
   #if defined(USE_GPU)
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
-    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
+    #if defined(CALL_FQ)
+      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]);
+      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]);
+    #else
+      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]);
+    #endif
   #else // !USE_GPU
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
     if (pd_start == 0) {
       #ifdef USE_OPENMP
       #pragma omp parallel for
       #endif
-      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+      #if defined(CALL_FQ)
+          // 2*nq for F^2,F pairs
+          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0;
+      #else
+          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+      #endif
     }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
@@ -415 +437 @@
           FETCH_Q         // set qx,qy from the q input vector
           APPLY_ROTATION  // convert qx,qy to qa,qb,qc
-          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
+          CALL_KERNEL     // F2 = Iqxy(qa, qb, qc, p1, p2, ...)
 
       ++step;  // increment counter representing position in dispersity mesh
@@ -446 +468 @@
 // inner loop and defines the macros that use them.
 
-#if defined(CALL_IQ)
-  // unoriented 1D
+
+#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
+  // unoriented 1D returning <F> and <F^2>
+  double qk;
+  double F1, F2;
+  #define FETCH_Q() do { qk = q[q_index]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table)
+
+#elif defined(CALL_FQ_A)
+  // unoriented 2D return <F> and <F^2>
+  double qx, qy;
+  double F1, F2;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table)
+
+#elif defined(CALL_IQ)
+  // unoriented 1D return <F^2>
   double qk;
   #define FETCH_Q() do { qk = q[q_index]; } while (0)
   #define BUILD_ROTATION() do {} while(0)
   #define APPLY_ROTATION() do {} while(0)
-  #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
+  #define CALL_KERNEL() CALL_IQ(qk,local_values.table)
 
 #elif defined(CALL_IQ_A)
@@ -486 +527 @@
   #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
+
 #elif defined(CALL_IQ_XY)
   // direct call to qx,qy calculator
@@ -499 +541 @@
 // logic in the IQ_AC and IQ_ABC branches.  This will become more important
 // if we implement more projections, or more complicated projections.
-#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
+#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A)
+  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
 #elif defined(CALL_IQ_XY) // pass orientation to the model
@@ -649 +692 @@
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
-      pd_norm += weight * CALL_VOLUME(local_values.table);
+      double form, shell;
+      CALL_VOLUME(form, shell, local_values.table);
+      weight_norm += weight;
+      weighted_form += weight * form;
+      weighted_shell += weight * shell;
+      if (effective_radius_type != 0) {
+        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
+      }
       BUILD_ROTATION();
 
@@ -667 +717 @@
         #if defined(MAGNETIC) && NUM_MAGNETIC > 0
           // Compute the scattering from the magnetic cross sections.
-          double scattering = 0.0;
+          double F2 = 0.0;
           const double qsq = qx*qx + qy*qy;
           if (qsq > 1.e-16) {
@@ -692 +742 @@
 //  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
                 }
-                scattering += xs_weight * CALL_KERNEL();
+                F2 += xs_weight * CALL_KERNEL();
               }
             }
           }
         #else  // !MAGNETIC
-          const double scattering = CALL_KERNEL();
+          #if defined(CALL_FQ)
+            CALL_KERNEL(); // sets F1 and F2 by reference
+          #else
+            const double F2 = CALL_KERNEL();
+          #endif
         #endif // !MAGNETIC
-//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
+//printf("q_index:%d %g %g %g %g\n", q_index, F2, weight0);
 
         #if defined(USE_GPU)
-          this_result += weight * scattering;
-        #else // !USE_GPU
-          result[q_index] += weight * scattering;
-        #endif // !USE_GPU
+          #if defined(CALL_FQ)
+            this_F2 += weight * F2;
+            this_F1 += weight * F1;
+          #else
+            this_F2 += weight * F2;
+          #endif
+        #else // !USE_OPENCL
+          #if defined(CALL_FQ)
+            result[2*q_index+0] += weight * F2;
+            result[2*q_index+1] += weight * F1;
+          #else
+            result[q_index] += weight * F2;
+          #endif
+        #endif // !USE_OPENCL
       }
     }
   }
-
 // close nested loops
 ++step;
@@ -728 +791 @@
 #endif
 
-// Remember the current result and the updated norm.
+// Remember the results and the updated norm.
 #if defined(USE_GPU)
-  result[q_index] = this_result;
-  if (q_index == 0) result[nq] = pd_norm;
-//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
-#else // !USE_GPU
-  result[nq] = pd_norm;
-//printf("res: %g/%g\n", result[0], pd_norm);
-#endif // !USE_GPU
+  #if defined(CALL_FQ)
+  result[2*q_index+0] = this_F2;
+  result[2*q_index+1] = this_F1;
+  #else
+  result[q_index] = this_F2;
+  #endif
+  if (q_index == 0)
+#endif
+  {
+#if defined(CALL_FQ)
+    result[2*nq] = weight_norm;
+    result[2*nq+1] = weighted_form;
+    result[2*nq+2] = weighted_shell;
+    result[2*nq+3] = weighted_radius;
+#else
+    result[nq] = weight_norm;
+    result[nq+1] = weighted_form;
+    result[nq+2] = weighted_shell;
+    result[nq+3] = weighted_radius;
+#endif
+  }
 
 // ** clear the macros in preparation for the next kernel **