Changeset c036ddb in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Aug 7, 2018 8:45:45 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:

7e923c2

Parents:

7b0abf8

Message:

refactor so Iq is not needed if Fq is defined

File:

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

sasmodels/kernel_iq.c

@@ -29 +29 @@
 //  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
@@ -85 +87 @@
   out_spin = clip(out_spin, 0.0, 1.0);
   // Previous version of this function took the square root of the weights,
-  // under the assumption that 
+  // under the assumption that
   //
   //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
@@ -272 +274 @@
 
 // ==================== KERNEL CODE ========================
+#define COMPUTE_F1_F2 defined(CALL_FQ)
 kernel
 void KERNEL_NAME(
@@ -338 +341 @@
   // 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
-  // version must loop over all q
-  #if BETA
-  double *beta_result = &(result[nq+1]); // = result + nq + 1
-  double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-  #endif
+  // version must loop over all q.
   #ifdef USE_OPENCL
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
-    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
-    #if BETA
-      double this_beta_result = (pd_start == 0 ? 0.0 : result[nq + q_index];
+    #if COMPUTE_F1_F2
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]);
+      double this_F1 = (pd_start == 0 ? 0.0 : result[q_index+nq+1]);
+    #else
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
+    #endif
   #else // !USE_OPENCL
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
+    #if COMPUTE_F1_F2
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+    #else
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    #endif
     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 BETA
-      for (int q_index=0; q_index < nq; q_index++) beta_result[q_index] = 0.0;
+      #if COMPUTE_F1_F2
+          for (int q_index=0; q_index < 2*nq+2; 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]);
 #endif // !USE_OPENCL
-
-
-
 
 
@@ -454 +461 @@
 
 
-#if defined(CALL_IQ)
-  // unoriented 1D
+#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
+  // unoriented 1D returning <F> and <F^2>
   double qk;
-  #if BETA == 0
-    #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)
-
-  // unoriented 1D Beta
-  #elif BETA == 1
-    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_IQ(qk,F1,F2,local_values.table)
-  #endif
+  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)
 
 #elif defined(CALL_IQ_A)
@@ -504 +519 @@
   #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
@@ -517 +533 @@
 // 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
@@ -668 +685 @@
     if (weight > cutoff) {
       pd_norm += weight * CALL_VOLUME(local_values.table);
+      #if COMPUTE_F1_F2
+      weight_norm += weight;
+      #endif
       BUILD_ROTATION();
-#if BETA
-    if (weight > cutoff) {
-      weight_norm += weight;}
-#endif
+
 #ifndef USE_OPENCL
       // DLL needs to explicitly loop over the q values.
@@ -718 +735 @@
           }
         #else  // !MAGNETIC
-          #if defined(CALL_IQ) && BETA == 1
-            CALL_KERNEL();
-            const double scatteringF1 = F1;
-            const double scatteringF2 = F2;
+          #if COMPUTE_F1_F2
+            CALL_KERNEL(); // sets F1 and F2 by reference
           #else
             const double scattering = CALL_KERNEL();
@@ -729 +744 @@
 
         #ifdef USE_OPENCL
-          #if defined(CALL_IQ)&& BETA == 1
-             this_result += weight * scatteringF2;
-             this_beta_result += weight * scatteringF1;
-            #else
-              this_result += weight * scattering;
+          #if COMPUTE_F1_F2
+            this_F1 += weight * F1;
+            this_F2 += weight * F2;
+          #else
+            this_result += weight * scattering;
           #endif
         #else // !USE_OPENCL
-          #if defined(CALL_IQ)&& BETA == 1
-            result[q_index] += weight * scatteringF2;
-            beta_result[q_index] += weight * scatteringF1;
-            #endif
-            #else
+          #if COMPUTE_F1_F2
+            result[q_index] += weight * F2;
+            result[q_index+nq+1] += weight * F1;
+          #else
             result[q_index] += weight * scattering;
           #endif
@@ -767 +781 @@
 // Remember the current result and the updated norm.
 #ifdef USE_OPENCL
-  result[q_index] = this_result;
-  if (q_index == 0) result[nq] = pd_norm;
-  #if BETA
-  beta_result[q_index] = this_beta_result;
+  #if COMPUTE_F1_F2
+    result[q_index] = this_F2;
+    result[q_index+nq+1] = this_F1;
+    if (q_index == 0) {
+      result[nq] = pd_norm;
+      result[2*nq+1] = weight_norm;
+    }
+  #else
+    result[q_index] = this_result;
+    if (q_index == 0) result[nq] = pd_norm;
   #endif
-  if (q_index == 0) result[nq] = pd_norm;
 
 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
 #else // !USE_OPENCL
-  result[nq] = pd_norm;
-  #if BETA
-  result[2*nq+1] = weight_norm;
+  #if COMPUTE_F1_F2
+    result[nq] = pd_norm;
+    result[2*nq+1] = weight_norm;
+  #else
+    result[nq] = pd_norm;
   #endif
 //printf("res: %g/%g\n", result[0], pd_norm);
@@ -784 +805 @@
 
 // ** clear the macros in preparation for the next kernel **
+#undef COMPUTE_F1_F2
 #undef PD_INIT
 #undef PD_OPEN