Index: doc/guide/gpu_setup.rst
===================================================================
--- doc/guide/gpu_setup.rst	(revision 63602b15f168fcfb37242dc2e94348a4751243a3)
+++ doc/guide/gpu_setup.rst	(revision 8b31efadb6413d7ae55766f7546c483c7538368e)
@@ -94,4 +94,6 @@
 Device Selection
 ================
+**OpenCL drivers**
+
 If you have multiple GPU devices you can tell the program which device to use.
 By default, the program looks for one GPU and one CPU device from available
@@ -104,8 +106,5 @@
 was used to run the model.
 
-**If you don't want to use OpenCL, you can set** *SAS_OPENCL=None*
-**in your environment settings, and it will only use normal programs.**
-
-If you want to use one of the other devices, you can run the following
+If you want to use a specific driver and devices, you can run the following
 from the python console::
 
@@ -115,5 +114,47 @@
 This will provide a menu of different OpenCL drivers available.
 When one is selected, it will say "set PYOPENCL_CTX=..."
-Use that value as the value of *SAS_OPENCL*.
+Use that value as the value of *SAS_OPENCL=driver:device*.
+
+To use the default OpenCL device (rather than CUDA or None),
+set *SAS_OPENCL=opencl*.
+
+In batch queues, you may need to set *XDG_CACHE_HOME=~/.cache* 
+(Linux only) to a different directory, depending on how the filesystem 
+is configured.  You should also set *SAS_DLL_PATH* for CPU-only modules.
+
+    -DSAS_MODELPATH=path sets directory containing custom models
+    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device
+    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only)
+    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler
+    -DSAS_OPENMP=1 turns on OpenMP for the DLLs
+    -DSAS_DLL_PATH=path sets the path to the compiled modules
+
+
+**CUDA drivers**
+
+If OpenCL drivers are not available on your system, but NVidia CUDA
+drivers are available, then set *SAS_OPENCL=cuda* or
+*SAS_OPENCL=cuda:n* for a particular device number *n*.  If no device
+number is specified, then the CUDA drivers looks for look for
+*CUDA_DEVICE=n* or a file ~/.cuda-device containing n for the device number.
+
+In batch queues, the SLURM command *sbatch --gres=gpu:1 ...* will set
+*CUDA_VISIBLE_DEVICES=n*, which ought to set the correct device
+number for *SAS_OPENCL=cuda*.  If not, then set
+*CUDA_DEVICE=$CUDA_VISIBLE_DEVICES* within the batch script.  You may
+need to set the CUDA cache directory to a folder accessible across the
+cluster with *PYCUDA_CACHE_DIR* (or *PYCUDA_DISABLE_CACHE* to disable
+caching), and you may need to set environment specific compiler flags
+with *PYCUDA_DEFAULT_NVCC_FLAGS*.  You should also set *SAS_DLL_PATH* 
+for CPU-only modules.
+
+**No GPU support**
+
+If you don't want to use OpenCL or CUDA, you can set *SAS_OPENCL=None*
+in your environment settings, and it will only use normal programs.
+
+In batch queues, you may need to set *SAS_DLL_PATH* to a directory
+accessible on the compute node.
+
 
 Device Testing
@@ -154,3 +195,3 @@
 *Document History*
 
-| 2017-09-27 Paul Kienzle
+| 2018-10-15 Paul Kienzle
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 610ef23645ef28d7c1e8c87aab1e03c8cd123a93)
+++ sasmodels/compare.py	(revision 4de145843b2267c4f1e8be16517c48de601151b6)
@@ -41,4 +41,5 @@
 from . import kerneldll
 from . import kernelcl
+from . import kernelcuda
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
@@ -115,5 +116,5 @@
     === environment variables ===
     -DSAS_MODELPATH=path sets directory containing custom models
-    -DSAS_OPENCL=vendor:device|none sets the target OpenCL device
+    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device
     -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only)
     -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler
@@ -725,5 +726,6 @@
         set_integration_size(model_info, ngauss)
 
-    if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl():
+    if (dtype != "default" and not dtype.endswith('!') 
+            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())):
         raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR)
 
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 2dcd6e7a1b6c16cef017142936211c8906750cae)
+++ sasmodels/core.py	(revision b0de25217231a59b280fb18c92467f20fb38c947)
@@ -21,4 +21,5 @@
 from . import mixture
 from . import kernelpy
+from . import kernelcuda
 from . import kernelcl
 from . import kerneldll
@@ -210,4 +211,6 @@
         #print("building dll", numpy_dtype)
         return kerneldll.load_dll(source['dll'], model_info, numpy_dtype)
+    elif platform == "cuda":
+        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast)
     else:
         #print("building ocl", numpy_dtype)
@@ -245,5 +248,5 @@
     # type: (ModelInfo, str, str) -> (np.dtype, bool, str)
     """
-    Interpret dtype string, returning np.dtype and fast flag.
+    Interpret dtype string, returning np.dtype, fast flag and platform.
 
     Possible types include 'half', 'single', 'double' and 'quad'.  If the
@@ -253,7 +256,10 @@
     default for the model and platform.
 
-    Platform preference can be specfied ("ocl" vs "dll"), with the default
-    being OpenCL if it is availabe.  If the dtype name ends with '!' then
-    platform is forced to be DLL rather than OpenCL.
+    Platform preference can be specfied ("ocl", "cuda", "dll"), with the
+    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype
+    name ends with '!' then platform is forced to be DLL rather than GPU.
+    The default platform is set by the environment variable SAS_OPENCL,
+    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA
+    or SAS_OPENCL=none for DLL.
 
     This routine ignores the preferences within the model definition.  This
@@ -268,6 +274,4 @@
     if platform is None:
         platform = "ocl"
-    if not kernelcl.use_opencl() or not model_info.opencl:
-        platform = "dll"
 
     # Check if type indicates dll regardless of which platform is given
@@ -275,4 +279,12 @@
         platform = "dll"
         dtype = dtype[:-1]
+
+    # Make sure model allows opencl/gpu
+    if not model_info.opencl:
+        platform = "dll"
+
+    # Make sure opencl is available, or fallback to cuda then to dll
+    if platform == "ocl" and not kernelcl.use_opencl():
+        platform = "cuda" if kernelcuda.use_cuda() else "dll"
 
     # Convert special type names "half", "fast", and "quad"
@@ -285,18 +297,23 @@
         dtype = "float16"
 
-    # Convert dtype string to numpy dtype.
+    # Convert dtype string to numpy dtype.  Use single precision for GPU
+    # if model allows it, otherwise use double precision.
     if dtype is None or dtype == "default":
-        numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single
+        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda")
                        else generate.F64)
     else:
         numpy_dtype = np.dtype(dtype)
 
-    # Make sure that the type is supported by opencl, otherwise use dll
+    # Make sure that the type is supported by GPU, otherwise use dll
     if platform == "ocl":
         env = kernelcl.environment()
-        if not env.has_type(numpy_dtype):
-            platform = "dll"
-            if dtype is None:
-                numpy_dtype = generate.F64
+    elif platform == "cuda":
+        env = kernelcuda.environment()
+    else:
+        env = None
+    if env is not None and not env.has_type(numpy_dtype):
+        platform = "dll"
+        if dtype is None:
+            numpy_dtype = generate.F64
 
     return numpy_dtype, fast, platform
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/kernel_header.c	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -1,17 +1,28 @@
 #ifdef __OPENCL_VERSION__
 # define USE_OPENCL
+#elif defined(__CUDACC__)
+# define USE_CUDA
 #elif defined(_OPENMP)
 # define USE_OPENMP
 #endif
+
+// Use SAS_DOUBLE to force the use of double even for float kernels
+#define SAS_DOUBLE dou ## ble
 
 // If opencl is not available, then we are compiling a C function
 // Note: if using a C++ compiler, then define kernel as extern "C"
 #ifdef USE_OPENCL
+
+   #define USE_GPU
+   #define pglobal global
+   #define pconstant constant
+
    typedef int int32_t;
-#  if defined(USE_SINCOS)
-#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
-#  else
-#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
-#  endif
+
+   #if defined(USE_SINCOS)
+   #  define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+   #else
+   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+   #endif
    // Intel CPU on Mac gives strange values for erf(); on the verified
    // platforms (intel, nvidia, amd), the cephes erf() is significantly
@@ -24,8 +35,33 @@
    #  define erfcf erfc
    #endif
-#else // !USE_OPENCL
-// Use SAS_DOUBLE to force the use of double even for float kernels
-#  define SAS_DOUBLE dou ## ble
-#  ifdef __cplusplus
+
+#elif defined(USE_CUDA)
+
+   #define USE_GPU
+   #define local __shared__
+   #define pglobal
+   #define constant __constant__
+   #define pconstant const
+   #define kernel extern "C" __global__
+
+   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+   // OpenCL pown(a,b) = C99 pow(a,b), b integer
+   #define powr(a,b) pow(a,b)
+   #define pown(a,b) pow(a,b)
+   //typedef int int32_t;
+   #if defined(USE_SINCOS)
+   #  define SINCOS(angle,svar,cvar) sincos(angle,&svar,&cvar)
+   #else
+   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+   #endif
+
+#else // !USE_OPENCL && !USE_CUDA
+
+   #define local
+   #define pglobal
+   #define constant const
+   #define pconstant const
+
+   #ifdef __cplusplus
       #include <cstdio>
       #include <cmath>
@@ -51,5 +87,5 @@
      #endif
      inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
-#  else // !__cplusplus
+   #else // !__cplusplus
      #include <inttypes.h>  // C99 guarantees that int32_t types is here
      #include <stdio.h>
@@ -76,12 +112,10 @@
      #define kernel
      #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
-#  endif  // !__cplusplus
-#  define global
-#  define local
-#  define constant const
-// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
-// OpenCL pown(a,b) = C99 pow(a,b), b integer
-#  define powr(a,b) pow(a,b)
-#  define pown(a,b) pow(a,b)
+   #endif  // !__cplusplus
+   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+   // OpenCL pown(a,b) = C99 pow(a,b), b integer
+   #define powr(a,b) pow(a,b)
+   #define pown(a,b) pow(a,b)
+
 #endif // !USE_OPENCL
 
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 705307783b90bea1a66cdf1ff38974433eb673cb)
+++ sasmodels/kernel_iq.c	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -278,14 +278,18 @@
     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
+    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
     )
 {
-#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
@@ -340,8 +344,8 @@
   // 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) {
@@ -352,5 +356,5 @@
     }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
-#endif // !USE_OPENCL
+#endif // !USE_GPU
 
 
@@ -375,6 +379,6 @@
   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];
+  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
 
@@ -562,6 +566,6 @@
   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]; \
+  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;
 
@@ -587,6 +591,6 @@
 // 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;
+  pglobal const double *pd_value = values + NUM_VALUES;
+  pglobal const double *pd_weight = pd_value + details->num_weights;
 #endif
 
@@ -648,5 +652,5 @@
       BUILD_ROTATION();
 
-#ifndef USE_OPENCL
+#if !defined(USE_GPU)
       // DLL needs to explicitly loop over the q values.
       #ifdef USE_OPENMP
@@ -654,5 +658,5 @@
       #endif
       for (q_index=0; q_index<nq; q_index++)
-#endif // !USE_OPENCL
+#endif // !USE_GPU
       {
 
@@ -697,9 +701,9 @@
 //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
       }
     }
@@ -725,12 +729,12 @@
 
 // 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 **
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 95f62aa49060757723d07a2113acdde358fa84f5)
+++ sasmodels/kernelcl.py	(revision 95f62aa49060757723d07a2113acdde358fa84f5)
@@ -1,4 +1,6 @@
 """
 GPU driver for C kernels
+
+TODO: docs are out of date
 
 There should be a single GPU environment running on the system.  This
@@ -59,5 +61,5 @@
 
 
-# Attempt to setup opencl. This may fail if the opencl package is not
+# Attempt to setup opencl. This may fail if the pyopencl package is not
 # installed or if it is installed but there are no devices available.
 try:
@@ -132,5 +134,6 @@
 
 def use_opencl():
-    return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none"
+    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower()
+    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda")
 
 ENV = None
@@ -180,27 +183,4 @@
         cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
         queue.device)
-
-def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
-    # type: (np.ndarray, np.dtype, float, int) -> np.ndarray
-    """
-    Stretch an input vector to the correct boundary.
-
-    Performance on the kernels can drop by a factor of two or more if the
-    number of values to compute does not fall on a nice power of two
-    boundary.   The trailing additional vector elements are given a
-    value of *extra*, and so f(*extra*) will be computed for each of
-    them.  The returned array will thus be a subset of the computed array.
-
-    *boundary* should be a power of 2 which is at least 32 for good
-    performance on current platforms (as of Jan 2015).  It should
-    probably be the max of get_warp(kernel,queue) and
-    device.min_data_type_align_size//4.
-    """
-    remainder = vector.size % boundary
-    if remainder != 0:
-        size = vector.size + (boundary - remainder)
-        vector = np.hstack((vector, [extra] * (size - vector.size)))
-    return np.ascontiguousarray(vector, dtype=dtype)
-
 
 def compile_model(context, source, dtype, fast=False):
@@ -342,9 +322,13 @@
     Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment,
     otherwise scans for the most appropriate device using
-    :func:`_get_default_context`
-    """
-    if 'SAS_OPENCL' in os.environ:
-        #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
-        os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"]
+    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which
+    indicates that an OpenCL device should be used without specifying
+    which one (and not a CUDA device, or no GPU).
+    """
+    # Assume we do not get here if SAS_OPENCL is None or CUDA
+    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl')
+    if sas_opencl.lower() != 'opencl':
+        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
+        os.environ["PYOPENCL_CTX"] = sas_opencl
 
     if 'PYOPENCL_CTX' in os.environ:
@@ -644,5 +628,5 @@
                 current_time = time.clock()
                 if current_time - last_nap > 0.5:
-                    time.sleep(0.05)
+                    time.sleep(0.001)
                     last_nap = current_time
         cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for)
Index: sasmodels/kernelcuda.py
===================================================================
--- sasmodels/kernelcuda.py	(revision 8b31efadb6413d7ae55766f7546c483c7538368e)
+++ sasmodels/kernelcuda.py	(revision 8b31efadb6413d7ae55766f7546c483c7538368e)
@@ -0,0 +1,557 @@
+"""
+GPU driver for C kernels (with CUDA)
+
+To select cuda, use SAS_OPENCL=cuda, or SAS_OPENCL=cuda:n for a particular
+device number.  If no device number is specified, then look for CUDA_DEVICE=n
+or a file ~/.cuda-device containing n for the device number.  Otherwise, try
+all available device numbers.
+
+TODO: docs are out of date
+
+There should be a single GPU environment running on the system.  This
+environment is constructed on the first call to :func:`env`, and the
+same environment is returned on each call.
+
+After retrieving the environment, the next step is to create the kernel.
+This is done with a call to :meth:`GpuEnvironment.make_kernel`, which
+returns the type of data used by the kernel.
+
+Next a :class:`GpuData` object should be created with the correct kind
+of data.  This data object can be used by multiple kernels, for example,
+if the target model is a weighted sum of multiple kernels.  The data
+should include any extra evaluation points required to compute the proper
+data smearing.  This need not match the square grid for 2D data if there
+is an index saying which q points are active.
+
+Together the GpuData, the program, and a device form a :class:`GpuKernel`.
+This kernel is used during fitting, receiving new sets of parameters and
+evaluating them.  The output value is stored in an output buffer on the
+devices, where it can be combined with other structure factors and form
+factors and have instrumental resolution effects applied.
+
+In order to use OpenCL for your models, you will need OpenCL drivers for
+your machine.  These should be available from your graphics card vendor.
+Intel provides OpenCL drivers for CPUs as well as their integrated HD
+graphics chipsets.  AMD also provides drivers for Intel CPUs, but as of
+this writing the performance is lacking compared to the Intel drivers.
+NVidia combines drivers for CUDA and OpenCL in one package.  The result
+is a bit messy if you have multiple drivers installed.  You can see which
+drivers are available by starting python and running:
+
+    import pyopencl as cl
+    cl.create_some_context(interactive=True)
+
+Once you have done that, it will show the available drivers which you
+can select.  It will then tell you that you can use these drivers
+automatically by setting the SAS_OPENCL environment variable, which is
+PYOPENCL_CTX equivalent but not conflicting with other pyopnecl programs.
+
+Some graphics cards have multiple devices on the same card.  You cannot
+yet use both of them concurrently to evaluate models, but you can run
+the program twice using a different device for each session.
+
+OpenCL kernels are compiled when needed by the device driver.  Some
+drivers produce compiler output even when there is no error.  You
+can see the output by setting PYOPENCL_COMPILER_OUTPUT=1.  It should be
+harmless, albeit annoying.
+"""
+from __future__ import print_function
+
+import os
+import warnings
+import logging
+import time
+import re
+
+import numpy as np  # type: ignore
+
+
+# Attempt to setup cuda. This may fail if the pycuda package is not
+# installed or if it is installed but there are no devices available.
+try:
+    import pycuda.driver as cuda  # type: ignore
+    from pycuda.compiler import SourceModule
+    from pycuda.tools import make_default_context, clear_context_caches
+    # Ask CUDA for the default context (so that we know that one exists)
+    # then immediately throw it away in case the user doesn't want it.
+    # Note: cribbed from pycuda.autoinit
+    cuda.init()
+    context = make_default_context()
+    context.pop()
+    clear_context_caches()
+    del context
+    HAVE_CUDA = True
+    CUDA_ERROR = ""
+except Exception as exc:
+    HAVE_CUDA = False
+    CUDA_ERROR = str(exc)
+
+from . import generate
+from .kernel import KernelModel, Kernel
+
+# pylint: disable=unused-import
+try:
+    from typing import Tuple, Callable, Any
+    from .modelinfo import ModelInfo
+    from .details import CallDetails
+except ImportError:
+    pass
+# pylint: enable=unused-import
+
+# The max loops number is limited by the amount of local memory available
+# on the device.  You don't want to make this value too big because it will
+# waste resources, nor too small because it may interfere with users trying
+# to do their polydispersity calculations.  A value of 1024 should be much
+# larger than necessary given that cost grows as npts^k where k is the number
+# of polydisperse parameters.
+MAX_LOOPS = 2048
+
+def use_cuda():
+    env = os.environ.get("SAS_OPENCL", "").lower()
+    return HAVE_CUDA and (env == "" or env.startswith("cuda"))
+
+ENV = None
+def reset_environment():
+    """
+    Call to create a new OpenCL context, such as after a change to SAS_OPENCL.
+    """
+    global ENV
+    # Free any previous allocated context.
+    if ENV is not None and ENV.context is not None:
+        ENV.release()
+    ENV = GpuEnvironment() if use_cuda() else None
+
+def environment():
+    # type: () -> "GpuEnvironment"
+    """
+    Returns a singleton :class:`GpuEnvironment`.
+
+    This provides an OpenCL context and one queue per device.
+    """
+    if ENV is None:
+        if not HAVE_CUDA:
+            raise RuntimeError("CUDA startup failed with ***"
+                            + CUDA_ERROR + "***; using C compiler instead")
+        reset_environment()
+        if ENV is None:
+            raise RuntimeError("SAS_OPENCL=None in environment")
+    return ENV
+
+def has_type(dtype):
+    # type: (np.dtype) -> bool
+    """
+    Return true if device supports the requested precision.
+    """
+    # Assume the nvidia card supports 32-bit and 64-bit floats.
+    # TODO: check if pycuda support F16
+    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
+    """
+    Build a model to run on the gpu.
+
+    Returns the compiled program and its type.  The returned type will
+    be float32 even if the desired type is float64 if any of the
+    devices in the context do not support the cl_khr_fp64 extension.
+    """
+    dtype = np.dtype(dtype)
+    if not has_type(dtype):
+        raise RuntimeError("%s not supported for devices"%dtype)
+
+    source_list = [generate.convert_type(source, dtype)]
+
+    source_list.insert(0, "#define USE_SINCOS\n")
+    source = "\n".join(source_list)
+    #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
+
+
+# for now, this returns one device in the context
+# TODO: create a context that contains all devices on all platforms
+class GpuEnvironment(object):
+    """
+    GPU context, with possibly many devices, and one queue per device.
+    """
+    context = None # type: cuda.Context
+    def __init__(self, devnum=None):
+        # type: (int) -> None
+        # Byte boundary for data alignment
+        #self.data_boundary = max(d.min_data_type_align_size
+        #                         for d in self.context.devices)
+        self.compiled = {}
+        env = os.environ.get("SAS_OPENCL", "").lower()
+        if devnum is None and env.startswith("cuda:"):
+            devnum = int(env[5:])
+        # Set the global context to the particular device number if one is
+        # given, otherwise use the default context.  Perhaps this will be set
+        # by an environment variable within autoinit.
+        if devnum is not None:
+            self.context = cuda.Device(devnum).make_context()
+        else:
+            self.context = make_default_context()
+
+    def release(self):
+        if self.context is not None:
+            self.context.pop()
+            self.context = None
+
+    def __del__(self):
+        self.release()
+
+    def has_type(self, dtype):
+        # type: (np.dtype) -> bool
+        """
+        Return True if all devices support a given type.
+        """
+        return has_type(dtype)
+
+    def compile_program(self, name, source, dtype, fast, timestamp):
+        # type: (str, str, np.dtype, bool, float) -> cl.Program
+        """
+        Compile the program for the device in the given context.
+        """
+        # Note: PyOpenCL caches based on md5 hash of source, options and device
+        # so we don't really need to cache things for ourselves.  I'll do so
+        # anyway just to save some data munging time.
+        tag = generate.tag_source(source)
+        key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else ""))
+        # Check timestamp on program
+        program, program_timestamp = self.compiled.get(key, (None, np.inf))
+        if program_timestamp < timestamp:
+            del self.compiled[key]
+        if key not in self.compiled:
+            logging.info("building %s for CUDA", key)
+            program = compile_model(str(source), dtype, fast)
+            self.compiled[key] = (program, timestamp)
+        return program
+
+class GpuModel(KernelModel):
+    """
+    GPU wrapper for a single model.
+
+    *source* and *model_info* are the model source and interface as returned
+    from :func:`generate.make_source` and :func:`generate.make_model_info`.
+
+    *dtype* is the desired model precision.  Any numpy dtype for single
+    or double precision floats will do, such as 'f', 'float32' or 'single'
+    for single and 'd', 'float64' or 'double' for double.  Double precision
+    is an optional extension which may not be available on all devices.
+    Half precision ('float16','half') may be available on some devices.
+    Fast precision ('fast') is a loose version of single precision, indicating
+    that the compiler is allowed to take shortcuts.
+    """
+    info = None # type: ModelInfo
+    source = "" # type: str
+    dtype = None # type: np.dtype
+    fast = False # type: bool
+    program = None # type: SourceModule
+    _kernels = None # type: List[cuda.Function]
+
+    def __init__(self, source, model_info, dtype=generate.F32, fast=False):
+        # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None
+        self.info = model_info
+        self.source = source
+        self.dtype = dtype
+        self.fast = fast
+        self.program = None # delay program creation
+        self._kernels = None
+
+    def __getstate__(self):
+        # type: () -> Tuple[ModelInfo, str, np.dtype, bool]
+        return self.info, self.source, self.dtype, self.fast
+
+    def __setstate__(self, state):
+        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
+        self.info, self.source, self.dtype, self.fast = state
+        self.program = None
+
+    def make_kernel(self, q_vectors):
+        # type: (List[np.ndarray]) -> "GpuKernel"
+        if self.program is None:
+            compile_program = environment().compile_program
+            timestamp = generate.ocl_timestamp(self.info)
+            self.program = compile_program(
+                self.info.name,
+                self.source['opencl'],
+                self.dtype,
+                self.fast,
+                timestamp)
+            variants = ['Iq', 'Iqxy', 'Imagnetic']
+            names = [generate.kernel_name(self.info, k) for k in variants]
+            kernels = [self.program.get_function(k) for k in names]
+            self._kernels = dict((k, v) for k, v in zip(variants, kernels))
+        is_2d = len(q_vectors) == 2
+        if is_2d:
+            kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']]
+        else:
+            kernel = [self._kernels['Iq']]*2
+        return GpuKernel(kernel, self.dtype, self.info, q_vectors)
+
+    def release(self):
+        # type: () -> None
+        """
+        Free the resources associated with the model.
+        """
+        if self.program is not None:
+            self.program = None
+
+    def __del__(self):
+        # type: () -> None
+        self.release()
+
+# TODO: check that we don't need a destructor for buffers which go out of scope
+class GpuInput(object):
+    """
+    Make q data available to the gpu.
+
+    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
+    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
+    to get the best performance on OpenCL, which may involve shifting and
+    stretching the array to better match the memory architecture.  Additional
+    points will be evaluated with *q=1e-3*.
+
+    *dtype* is the data type for the q vectors. The data type should be
+    set to match that of the kernel, which is an attribute of
+    :class:`GpuProgram`.  Note that not all kernels support double
+    precision, so even if the program was created for double precision,
+    the *GpuProgram.dtype* may be single precision.
+
+    Call :meth:`release` when complete.  Even if not called directly, the
+    buffer will be released when the data object is freed.
+    """
+    def __init__(self, q_vectors, dtype=generate.F32):
+        # type: (List[np.ndarray], np.dtype) -> None
+        # TODO: do we ever need double precision q?
+        self.nq = q_vectors[0].size
+        self.dtype = np.dtype(dtype)
+        self.is_2d = (len(q_vectors) == 2)
+        # TODO: stretch input based on get_warp()
+        # not doing it now since warp depends on kernel, which is not known
+        # at this point, so instead using 32, which is good on the set of
+        # architectures tested so far.
+        if self.is_2d:
+            # Note: 16 rather than 15 because result is 1 longer than input.
+            width = ((self.nq+16)//16)*16
+            self.q = np.empty((width, 2), dtype=dtype)
+            self.q[:self.nq, 0] = q_vectors[0]
+            self.q[:self.nq, 1] = q_vectors[1]
+        else:
+            # Note: 32 rather than 31 because result is 1 longer than input.
+            width = ((self.nq+32)//32)*32
+            self.q = np.empty(width, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
+        self.global_size = [self.q.shape[0]]
+        #print("creating inputs of size", self.global_size)
+        self.q_b = cuda.to_device(self.q)
+
+    def release(self):
+        # type: () -> None
+        """
+        Free the memory.
+        """
+        if self.q_b is not None:
+            self.q_b.free()
+            self.q_b = None
+
+    def __del__(self):
+        # type: () -> None
+        self.release()
+
+class GpuKernel(Kernel):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the GpuKernel object to call
+
+    *model_info* is the module information
+
+    *q_vectors* is the q vectors at which the kernel should be evaluated
+
+    *dtype* is the kernel precision
+
+    The resulting call method takes the *pars*, a list of values for
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    vectors for the polydisperse parameters.  *cutoff* determines the
+    integration limits: any points with combined weight less than *cutoff*
+    will not be calculated.
+
+    Call :meth:`release` when done with the kernel instance.
+    """
+    def __init__(self, kernel, dtype, model_info, q_vectors):
+        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
+        q_input = GpuInput(q_vectors, dtype)
+        self.kernel = kernel
+        self.info = model_info
+        self.dtype = dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
+        # plus three for the normalization values
+        self.result = np.empty(q_input.nq+1, dtype)
+
+        # Inputs and outputs for each kernel call
+        # Note: res may be shorter than res_b if global_size != nq
+        self.result_b = cuda.mem_alloc(q_input.global_size[0] * dtype.itemsize)
+        self.q_input = q_input # allocated by GpuInput above
+
+        self._need_release = [self.result_b]
+        self.real = (np.float32 if dtype == generate.F32
+                     else np.float64 if dtype == generate.F64
+                     else np.float16 if dtype == generate.F16
+                     else np.float32)  # will never get here, so use np.float32
+
+    def __call__(self, call_details, values, cutoff, magnetic):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        # Arrange data transfer to card
+        details_b = cuda.to_device(call_details.buffer)
+        values_b = cuda.to_device(values)
+
+        kernel = self.kernel[1 if magnetic else 0]
+        args = [
+            np.uint32(self.q_input.nq), None, None,
+            details_b, values_b, self.q_input.q_b, self.result_b,
+            self.real(cutoff),
+        ]
+        grid = partition(self.q_input.nq)
+        #print("Calling OpenCL")
+        #call_details.show(values)
+        # Call kernel and retrieve results
+        last_nap = time.clock()
+        step = 100000000//self.q_input.nq + 1
+        #step = 1000000000
+        for start in range(0, call_details.num_eval, step):
+            stop = min(start + step, call_details.num_eval)
+            #print("queuing",start,stop)
+            args[1:3] = [np.int32(start), np.int32(stop)]
+            kernel(*args, **grid)
+            if stop < call_details.num_eval:
+                sync()
+                # Allow other processes to run
+                current_time = time.clock()
+                if current_time - last_nap > 0.5:
+                    time.sleep(0.001)
+                    last_nap = current_time
+        sync()
+        cuda.memcpy_dtoh(self.result, self.result_b)
+        #print("result", self.result)
+
+        details_b.free()
+        values_b.free()
+
+        pd_norm = self.result[self.q_input.nq]
+        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
+        background = values[1]
+        #print("scale",scale,values[0],self.result[self.q_input.nq],background)
+        return scale*self.result[:self.q_input.nq] + background
+        # return self.result[:self.q_input.nq]
+
+    def release(self):
+        # type: () -> None
+        """
+        Release resources associated with the kernel.
+        """
+        for p in self._need_release:
+            p.free()
+        self._need_release = []
+
+    def __del__(self):
+        # type: () -> None
+        self.release()
+
+
+def sync():
+    """
+    Overview:
+        Waits for operation in the current context to complete.
+
+    Note: Maybe context.synchronize() is sufficient.
+    """
+    #return # The following works in C++; don't know what pycuda is doing
+    # Create an event with which to synchronize
+    done = cuda.Event()
+
+    # Schedule an event trigger on the GPU.
+    done.record()
+
+    #line added to not hog resources
+    while not done.query():
+        time.sleep(0.01)
+
+    # Block until the GPU executes the kernel.
+    done.synchronize()
+    # Clean up the event; I don't think they can be reused.
+    del done
+
+
+def partition(n):
+    '''
+    Overview:
+        Auto grids the thread blocks to achieve some level of calculation
+    efficiency.
+    '''
+    max_gx, max_gy = 65535, 65535
+    blocksize = 32
+    #max_gx, max_gy = 5, 65536
+    #blocksize = 3
+    block = (blocksize, 1, 1)
+    num_blocks = int((n+blocksize-1)/blocksize)
+    if num_blocks < max_gx:
+        grid = (num_blocks, 1)
+    else:
+        gx = max_gx
+        gy = (num_blocks + max_gx - 1) / max_gx
+        if gy >= max_gy:
+            raise ValueError("vector is too large")
+        grid = (gx, gy)
+    #print("block", block, "grid", grid)
+    #print("waste", block[0]*block[1]*block[2]*grid[0]*grid[1] - n)
+    return dict(block=block, grid=grid)
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 012cd342ba6d27483147e3401f01fb7b93ea589c)
+++ sasmodels/model_test.py	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -5,5 +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,4 +63,5 @@
 from .modelinfo import expand_pars
 from .kernelcl import use_opencl
+from .kernelcuda import use_cuda
 
 # pylint: disable=unused-import
@@ -80,7 +81,6 @@
     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,5 +135,5 @@
 
             # 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,4 +156,19 @@
                                      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,5 +235,5 @@
 
                 # 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,8 +383,8 @@
 
     # 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,22 +449,29 @@
         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,5 +493,9 @@
     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):
Index: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 99b84ec0c1c8127b9cbe4beeecc251208c44071c)
+++ sasmodels/models/lib/gauss76.c	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -11,5 +11,5 @@
 
 // Gaussians
-constant double Gauss76Wt[76]={
+constant double Gauss76Wt[76] = {
 	.00126779163408536,		//0
 	.00294910295364247,
@@ -90,5 +90,5 @@
 };
 
-constant double Gauss76Z[76]={
+constant double Gauss76Z[76] = {
 	-.999505948362153,		//0
 	-.997397786355355,
Index: sasmodels/models/lib/polevl.c
===================================================================
--- sasmodels/models/lib/polevl.c	(revision 447e9aab67fe022330aad0a4c8873028603d4b87)
+++ sasmodels/models/lib/polevl.c	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -51,6 +51,6 @@
 */
 
-double polevl( double x, constant double *coef, int N );
-double polevl( double x, constant double *coef, int N )
+static
+double polevl( double x, pconstant double *coef, int N )
 {
 
@@ -72,6 +72,6 @@
  */
 
-double p1evl( double x, constant double *coef, int N );
-double p1evl( double x, constant double *coef, int N )
+static
+double p1evl( double x, pconstant double *coef, int N )
 {
     int i=0;
Index: sasmodels/models/lib/sas_J1.c
===================================================================
--- sasmodels/models/lib/sas_J1.c	(revision 5181ccc6cbc1277527636856e4db39c6557b26b7)
+++ sasmodels/models/lib/sas_J1.c	(revision 74e9b5fa9e460b3726df690a1411aef195f15b4e)
@@ -42,5 +42,4 @@
 #if FLOAT_SIZE>4
 //Cephes double pression function
-double cephes_j1(double x);
 
 constant double RPJ1[8] = {
@@ -106,4 +105,5 @@
     0.0 };
 
+static
 double cephes_j1(double x)
 {
@@ -155,6 +155,4 @@
 #else
 //Single precission version of cephes
-float cephes_j1f(float x);
-
 constant float JPJ1[8] = {
     -4.878788132172128E-009,
@@ -190,4 +188,5 @@
     };
 
+static
 float cephes_j1f(float xx)
 {
@@ -240,5 +239,5 @@
 
 //Finally J1c function that equals 2*J1(x)/x
-double sas_2J1x_x(double x);
+static
 double sas_2J1x_x(double x)
 {
