Changeset 74e9b5f in sasmodels

Timestamp:

Oct 12, 2018 10:52:48 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:

4de14584

Parents:

b0de252

Message:

autotag functions as device functions for cuda. Refs #1076.

Location:

sasmodels

Files:

• kernel_header.c (modified) (5 diffs)
• kernel_iq.c (modified) (11 diffs)
• kernelcuda.py (modified) (3 diffs)
• model_test.py (modified) (9 diffs)
• models/cylinder.c (modified) (5 diffs)
• models/lib/gauss76.c (modified) (2 diffs)
• models/lib/polevl.c (modified) (2 diffs)
• models/lib/sas_J1.c (modified) (12 diffs)

sasmodels/kernel_header.c

@@ -5 +5 @@
 #elif defined(_OPENMP)
 # define USE_OPENMP
-#elif defined(__CUDACC__)
-# define USE_CUDA
 #endif
 
@@ -17 +15 @@
 
    #define USE_GPU
+   #define pglobal global
+   #define pconstant constant
+
    typedef int int32_t;
-   #define global_par global
-   #define local_par local
-   #define constant_par constant
-   #define global_var global
-   #define local_var local
-   #define constant_var constant
-   #define __device__
 
    #if defined(USE_SINCOS)
@@ -45 +39 @@
 
    #define USE_GPU
-   #define global_par
-   #define local_par
-   #define constant_par const
-   #define global_var
-   #define local_var __shared__
-   #define constant_var __constant__
-
+   #define local __shared__
+   #define pglobal
+   #define constant __constant__
+   #define pconstant const
    #define kernel extern "C" __global__
 
@@ -67 +58 @@
 #else // !USE_OPENCL && !USE_CUDA
 
-   #define global_par
-   #define local_par
-   #define constant_par const
-   #define global_var
-   #define local_var
-   #define constant_var const
-   #define __device__
+   #define local
+   #define pglobal
+   #define constant const
+   #define pconstant const
 
    #ifdef __cplusplus
@@ -193 +181 @@
 #  define M_4PI_3 4.18879020478639
 #endif
-__device__
 inline double square(double x) { return x*x; }
-__device__
 inline double cube(double x) { return x*x*x; }
-__device__
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 

sasmodels/kernel_iq.c

@@ -67 +67 @@
 
 // Return value restricted between low and high
-__device__
 static double clip(double value, double low, double high)
 {
@@ -80 +79 @@
 //     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 weight[6])
 {
@@ -105 +103 @@
 
 // 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=ud.imag
@@ -154 +151 @@
 // 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(
@@ -188 +184 @@
 // Apply the rotation matrix returned from qac_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qc]'
-__device__
 static void
 qac_apply(
@@ -216 +211 @@
 // 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(
@@ -263 +257 @@
 // Apply the rotation matrix returned from qabc_rotation to the point (qx,qy),
 // returning R*[qx,qy]' = [qa,qb,qc]'
-__device__
 static void
 qabc_apply(
@@ -285 +278 @@
     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_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
+    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
     )
@@ -386 +379 @@
   const int n4 = pd_length[4];
   const int p4 = pd_par[4];
-  global_var const double *v4 = pd_value + pd_offset[4];
-  global_var const double *w4 = pd_weight + pd_offset[4];
+  pglobal const double *v4 = pd_value + pd_offset[4];
+  pglobal const double *w4 = pd_weight + pd_offset[4];
   int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start
 
@@ -573 +566 @@
   const int n##_LOOP = details->pd_length[_LOOP]; \
   const int p##_LOOP = details->pd_par[_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]; \
+  pglobal const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  pglobal const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
   int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
 
@@ -598 +591 @@
 // Pointers to the start of the dispersity and weight vectors, if needed.
 #if MAX_PD>0
-  global_var const double *pd_value = values + NUM_VALUES;
-  global_var const double *pd_weight = pd_value + details->num_weights;
+  pglobal const double *pd_value = values + NUM_VALUES;
+  pglobal const double *pd_weight = pd_value + details->num_weights;
 #endif
 

sasmodels/kernelcuda.py

@@ -62 +62 @@
 import logging
 import time
+import re
 
 import numpy as np  # type: ignore
@@ -146 +147 @@
     return dtype in (generate.F32, generate.F64)
 
+
+FUNCTION_PATTERN = re.compile(r"""^
+  (?P<space>\s*)                   # initial space
+  (?P<qualifiers>^(?:\s*\b\w+\b\s*)+) # one or more qualifiers before function
+  (?P<function>\s*\b\w+\b\s*[(])      # function name plus open parens
+  """, re.VERBOSE|re.MULTILINE)
+
+MARKED_PATTERN = re.compile(r"""
+  \b(return|else|kernel|device|__device__)\b
+  """, re.VERBOSE|re.MULTILINE)
+
+def _add_device_tag(match):
+    # type: (None) -> str
+    # Note: should be re.Match, but that isn't a simple type
+    """
+    replace qualifiers with __device__ qualifiers if needed
+    """
+    qualifiers = match.group("qualifiers")
+    if MARKED_PATTERN.search(qualifiers):
+        start, end = match.span()
+        return match.string[start:end]
+    else:
+        function = match.group("function")
+        space = match.group("space")
+        return "".join((space, "__device__ ", qualifiers, function))
+
+def mark_device_functions(source):
+    # type: (str) -> str
+    """
+    Mark all function declarations as __device__ functions (except kernel).
+    """
+    return FUNCTION_PATTERN.sub(_add_device_tag, source)
+
+def show_device_functions(source):
+    # type: (str) -> str
+    """
+    Show all discovered function declarations, but don't change any.
+    """
+    for match in FUNCTION_PATTERN.finditer(source):
+        print(match.group('qualifiers').replace('\n',r'\n'), match.group('function'), '(')
+    return source
+
 def compile_model(source, dtype, fast=False):
     # type: (str, np.dtype, bool) -> SourceModule
@@ -163 +206 @@
     source_list.insert(0, "#define USE_SINCOS\n")
     source = "\n".join(source_list)
-    options = '-use_fast_math' if fast else None
+    #source = show_device_functions(source)
+    source = mark_device_functions(source)
+    #with open('/tmp/kernel.cu', 'w') as fd: fd.write(source)
+    #print(source)
+    #options = ['--verbose', '-E']
+    options = ['--use_fast_math'] if fast else None
     program = SourceModule(source, no_extern_c=True, options=options) # include_dirs=[...]
+
     #print("done with "+program)
     return program

sasmodels/model_test.py

@@ -5 +5 @@
 Usage::
 
-    python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+    python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ...
 
     if model1 is 'all', then all except the remaining models will be tested
@@ -63 +63 @@
 from .modelinfo import expand_pars
 from .kernelcl import use_opencl
+from .kernelcuda import use_cuda
 
 # pylint: disable=unused-import
@@ -80 +81 @@
     Construct the pyunit test suite.
 
-    *loaders* is the list of kernel drivers to use, which is one of
-    *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models,
-    the python driver is always used.
+    *loaders* is the list of kernel drivers to use (dll, opencl or cuda).
+    For python model the python driver is always used.
 
     *models* is the list of models to test, or *["all"]* to test all models.
@@ -135 +135 @@
 
             # test using dll if desired
-            if 'dll' in loaders or not use_opencl():
+            if 'dll' in loaders:
                 test_name = "%s-dll"%model_name
                 test_method_name = "test_%s_dll" % model_info.id
@@ -156 +156 @@
                                      test_method_name,
                                      platform="ocl", dtype=None,
+                                     stash=stash)
+                #print("defining", test_name)
+                suite.addTest(test)
+
+            # test using cuda if desired and available
+            if 'cuda' in loaders and use_cuda():
+                test_name = "%s-cuda"%model_name
+                test_method_name = "test_%s_cuda" % model_info.id
+                # Using dtype=None so that the models that are only
+                # correct for double precision are not tested using
+                # single precision.  The choice is determined by the
+                # presence of *single=False* in the model file.
+                test = ModelTestCase(test_name, model_info,
+                                     test_method_name,
+                                     platform="cuda", dtype=None,
                                      stash=stash)
                 #print("defining", test_name)
@@ -220 +235 @@
 
                 # Check for missing tests.  Only do so for the "dll" tests
-                # to reduce noise from both opencl and dll, and because
+                # to reduce noise from both opencl and cuda, and because
                 # python kernels use platform="dll".
                 if self.platform == "dll":
@@ -368 +383 @@
 
     # Build a test suite containing just the model
-    loaders = ['opencl'] if use_opencl() else ['dll']
+    loader = 'opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'
     models = [model]
     try:
-        suite = make_suite(loaders, models)
+        suite = make_suite([loader], models)
     except Exception:
         import traceback
@@ -434 +449 @@
         loaders = ['opencl']
         models = models[1:]
+    elif models and models[0] == 'cuda':
+        if not use_cuda():
+            print("cuda is not available")
+            return 1
+        loaders = ['cuda']
+        models = models[1:]
     elif models and models[0] == 'dll':
         # TODO: test if compiler is available?
         loaders = ['dll']
         models = models[1:]
-    elif models and models[0] == 'opencl_and_dll':
-        loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
-        models = models[1:]
     else:
-        loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
+        loaders = ['dll']
+        if use_opencl():
+            loaders.append('opencl')
+        if use_cuda():
+            loaders.append('cuda')
     if not models:
         print("""\
 usage:
-  python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ...
+  python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ...
 
 If -v is included on the command line, then use verbose output.
 
-If neither opencl nor dll is specified, then models will be tested with
-both OpenCL and dll; the compute target is ignored for pure python models.
+If no platform is specified, then models will be tested with dll, and
+if available, OpenCL and CUDA; the compute target is ignored for pure python models.
 
 If model1 is 'all', then all except the remaining models will be tested.
@@ -471 +493 @@
     Run "nosetests sasmodels" on the command line to invoke it.
     """
-    loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
+    loaders = ['dll']
+    if use_opencl():
+        loaders.append('opencl')
+    if use_cuda():
+        loaders.append('cuda')
     tests = make_suite(loaders, ['all'])
     def build_test(test):

sasmodels/models/cylinder.c

@@ -1 +1 @@
 #define INVALID(v) (v.radius<0 || v.length<0)
 
-__device__
 static double
 form_volume(double radius, double length)
@@ -8 +7 @@
 }
 
-__device__
 static double
 fq(double qab, double qc, double radius, double length)
@@ -15 +13 @@
 }
 
-__device__
 static double
 orient_avg_1D(double q, double radius, double length)
@@ -36 +33 @@
 }
 
-__device__
 static double
 Iq(double q,
@@ -48 +44 @@
 }
 
-__device__
 static double
 Iqac(double qab, double qc,

sasmodels/models/lib/gauss76.c

@@ -11 +11 @@
 
 // Gaussians
-constant_var double Gauss76Wt[76] = {
+constant double Gauss76Wt[76] = {
         .00126779163408536,             //0
         .00294910295364247,
@@ -90 +90 @@
 };
 
-constant_var double Gauss76Z[76] = {
+constant double Gauss76Z[76] = {
         -.999505948362153,              //0
         -.997397786355355,

sasmodels/models/lib/polevl.c

@@ -51 +51 @@
 */
 
-__device__ static
-double polevl( double x, constant_par double *coef, int N )
+static
+double polevl( double x, pconstant double *coef, int N )
 {
 
@@ -72 +72 @@
  */
 
-__device__ static
-double p1evl( double x, constant_par double *coef, int N )
+static
+double p1evl( double x, pconstant double *coef, int N )
 {
     int i=0;

sasmodels/models/lib/sas_J1.c

@@ -43 +43 @@
 //Cephes double pression function
 
-constant_var double RPJ1[8] = {
+constant double RPJ1[8] = {
     -8.99971225705559398224E8,
     4.52228297998194034323E11,
@@ -53 +53 @@
     0.0 };
 
-constant_var double RQJ1[8] = {
+constant double RQJ1[8] = {
     6.20836478118054335476E2,
     2.56987256757748830383E5,
@@ -64 +64 @@
     };
 
-constant_var double PPJ1[8] = {
+constant double PPJ1[8] = {
     7.62125616208173112003E-4,
     7.31397056940917570436E-2,
@@ -75 +75 @@
 
 
-constant_var double PQJ1[8] = {
+constant double PQJ1[8] = {
     5.71323128072548699714E-4,
     6.88455908754495404082E-2,
@@ -85 +85 @@
     0.0 };
 
-constant_var double QPJ1[8] = {
+constant double QPJ1[8] = {
     5.10862594750176621635E-2,
     4.98213872951233449420E0,
@@ -95 +95 @@
     2.52070205858023719784E1 };
 
-constant_var double QQJ1[8] = {
+constant double QQJ1[8] = {
     7.42373277035675149943E1,
     1.05644886038262816351E3,
@@ -105 +105 @@
     0.0 };
 
-__device__ static
+static
 double cephes_j1(double x)
 {
@@ -155 +155 @@
 #else
 //Single precission version of cephes
-constant_var float JPJ1[8] = {
+constant float JPJ1[8] = {
     -4.878788132172128E-009,
     6.009061827883699E-007,
@@ -166 +166 @@
     };
 
-constant_var float MO1J1[8] = {
+constant float MO1J1[8] = {
     6.913942741265801E-002,
     -2.284801500053359E-001,
@@ -177 +177 @@
     };
 
-constant_var float PH1J1[8] = {
+constant float PH1J1[8] = {
     -4.497014141919556E+001,
     5.073465654089319E+001,
@@ -188 +188 @@
     };
 
-__device__ static
+static
 float cephes_j1f(float xx)
 {
@@ -239 +239 @@
 
 //Finally J1c function that equals 2*J1(x)/x
-__device__ static
+static
 double sas_2J1x_x(double x)
 {