Changeset 0db7dbd in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Feb 16, 2018 7:10:04 PM (6 years ago)

Author:

pkienzle

Branches:

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

Children:

47fb816

Parents:

aa90015

Message:

cuda support: allow cylinder model to run under CUDA as well as OpenCL

File:

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

sasmodels/kernel_iq.c

@@ -67 +67 @@
 
 // Return value restricted between low and high
+__device__
 static double clip(double value, double low, double high)
 {
@@ -79 +80 @@
 //     du * (m_sigma_y + 1j*m_sigma_z);
 // weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
+__device__
 static void set_spin_weights(double in_spin, double out_spin, double spins[4])
 {
@@ -92 +94 @@
 
 // Compute the magnetic sld
+__device__
 static double mag_sld(
   const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag
@@ -140 +143 @@
 // jitter angles (dtheta, dphi).  This matrix can be applied to all of the
 // (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qc]'
+__device__
 static void
 qac_rotation(
@@ -173 +177 @@
 // Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qc]'
-static double
+__device__
+static void
 qac_apply(
     QACRotation *rotation,
@@ -200 +205 @@
 // jitter angles (dtheta, dphi, dpsi).  This matrix can be applied to all of the
 // (qx, qy) points in the image to produce R*[qx,qy]' = [qa,qb,qc]'
+__device__
 static void
 qabc_rotation(
@@ -246 +252 @@
 // Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qb,qc]'
-static double
+__device__
+static void
 qabc_apply(
     QABCRotation *rotation,
@@ -267 +274 @@
     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
-    global const ProblemDetails *details,
-    global const double *values,
-    global const double *q, // nq q values, with padding to boundary
-    global double *result,  // nq+1 return values, again with padding
+    global_par const ProblemDetails *details,
+    global_par const double *values,
+    global_par const double *q, // nq q values, with padding to boundary
+    global_par double *result,  // nq+1 return values, again with padding
     const double cutoff     // cutoff in the dispersity weight product
     )
 {
-#ifdef USE_OPENCL
+#if defined(USE_GPU)
   // who we are and what element we are working with
+  #if defined(USE_OPENCL)
   const int q_index = get_global_id(0);
+  #else // USE_CUDA
+  const int q_index = threadIdx.x + blockIdx.x * blockDim.x;
+  #endif
   if (q_index >= nq) return;
 #else
@@ -329 +340 @@
   // seeing one q value (stored in the variable "this_result") while the dll
   // version must loop over all q.
-  #ifdef USE_OPENCL
+  #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]);
-  #else // !USE_OPENCL
+  #else // !USE_GPU
     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
     if (pd_start == 0) {
@@ -341 +352 @@
     }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
-#endif // !USE_OPENCL
+#endif // !USE_GPU
 
 
@@ -364 +375 @@
   const int n4 = pd_length[4];
   const int p4 = pd_par[4];
-  global const double *v4 = pd_value + pd_offset[4];
-  global const double *w4 = pd_weight + pd_offset[4];
+  global_var const double *v4 = pd_value + pd_offset[4];
+  global_var const double *w4 = pd_weight + pd_offset[4];
   int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start
 
@@ -551 +562 @@
   const int n##_LOOP = details->pd_length[_LOOP]; \
   const int p##_LOOP = details->pd_par[_LOOP]; \
-  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
-  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
+  global_var const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  global_var const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
   int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
 
@@ -576 +587 @@
 // Pointers to the start of the dispersity and weight vectors, if needed.
 #if MAX_PD>0
-  global const double *pd_value = values + NUM_VALUES;
-  global const double *pd_weight = pd_value + details->num_weights;
+  global_var const double *pd_value = values + NUM_VALUES;
+  global_var const double *pd_weight = pd_value + details->num_weights;
 #endif
 
@@ -637 +648 @@
       BUILD_ROTATION();
 
-#ifndef USE_OPENCL
+#if !defined(USE_GPU)
       // DLL needs to explicitly loop over the q values.
       #ifdef USE_OPENMP
@@ -643 +654 @@
       #endif
       for (q_index=0; q_index<nq; q_index++)
-#endif // !USE_OPENCL
+#endif // !USE_GPU
       {
 
@@ -684 +695 @@
 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
 
-        #ifdef USE_OPENCL
+        #if defined(USE_GPU)
           this_result += weight * scattering;
-        #else // !USE_OPENCL
+        #else // !USE_GPU
           result[q_index] += weight * scattering;
-        #endif // !USE_OPENCL
+        #endif // !USE_GPU
       }
     }
@@ -712 +723 @@
 
 // Remember the current result and the updated norm.
-#ifdef USE_OPENCL
+#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_OPENCL
+#else // !USE_GPU
   result[nq] = pd_norm;
 //printf("res: %g/%g\n", result[0], pd_norm);
-#endif // !USE_OPENCL
+#endif // !USE_GPU
 
 // ** clear the macros in preparation for the next kernel **