Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 3221de0ed9b3480163246ceb0503432126bd5134)
+++ sasmodels/core.py	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -13,4 +13,7 @@
 from glob import glob
 import re
+
+# Set "SAS_OPENCL=cuda" in the environment to use the CUDA rather than OpenCL
+USE_CUDA = os.environ.get("SAS_OPENCL", "") == "cuda"
 
 import numpy as np # type: ignore
@@ -21,5 +24,8 @@
 from . import mixture
 from . import kernelpy
-from . import kernelcl
+if USE_CUDA:
+    from . import kernelcuda
+else:
+    from . import kernelcl
 from . import kerneldll
 from . import custom
@@ -210,4 +216,7 @@
         #print("building dll", numpy_dtype)
         return kerneldll.load_dll(source['dll'], model_info, numpy_dtype)
+    elif USE_CUDA:
+        #print("building cuda", numpy_dtype)
+        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast)
     else:
         #print("building ocl", numpy_dtype)
@@ -268,6 +277,12 @@
     if platform is None:
         platform = "ocl"
-    if not kernelcl.use_opencl() or not model_info.opencl:
+    if not model_info.opencl:
         platform = "dll"
+    elif USE_CUDA:
+        if not kernelcuda.use_cuda():
+            platform = "dll"
+    else:
+        if not kernelcl.use_opencl():
+            platform = "dll"
 
     # Check if type indicates dll regardless of which platform is given
@@ -294,5 +309,8 @@
     # Make sure that the type is supported by opencl, otherwise use dll
     if platform == "ocl":
-        env = kernelcl.environment()
+        if USE_CUDA:
+            env = kernelcuda.environment()
+        else:
+            env = kernelcl.environment()
         if not env.has_type(numpy_dtype):
             platform = "dll"
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/kernel_header.c	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -1,4 +1,6 @@
 #ifdef __OPENCL_VERSION__
 # define USE_OPENCL
+#elif defined(__CUDACC__)
+# define USE_CUDA
 #elif defined(_OPENMP)
 # define USE_OPENMP
@@ -8,10 +10,20 @@
 // Note: if using a C++ compiler, then define kernel as extern "C"
 #ifdef USE_OPENCL
+
+   #define USE_GPU
    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
+   #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)
+   #  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 +36,39 @@
    #  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 global_par
+   #define local_par
+   #define constant_par const
+   #define global_var
+   #define local_var __shared__
+   #define constant_var __constant__
+
+   #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 global_par
+   #define local_par
+   #define constant_par const
+   #define global_var
+   #define local_var
+   #define constant_var const
+   #define __device__
+
+   #ifdef __cplusplus
       #include <cstdio>
       #include <cmath>
@@ -51,5 +94,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,13 +119,14 @@
      #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
+
+// Use SAS_DOUBLE to force the use of double even for float kernels
+#define SAS_DOUBLE dou ## ble
 
 #if defined(NEED_EXPM1)
@@ -147,6 +191,9 @@
 #  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; }
 
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision aadec17f5a4db1c04d3dbefd23a3ac36b2dd5b96)
+++ sasmodels/kernel_iq.c	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -67,4 +67,5 @@
 
 // Return value restricted between low and high
+__device__
 static double clip(double value, double low, double high)
 {
@@ -79,4 +80,5 @@
 //     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,4 +94,5 @@
 
 // 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,4 +143,5 @@
 // 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,5 +177,6 @@
 // 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,4 +205,5 @@
 // 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,5 +252,6 @@
 // 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,14 +274,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
+    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,8 +340,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) {
@@ -341,5 +352,5 @@
     }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
-#endif // !USE_OPENCL
+#endif // !USE_GPU
 
 
@@ -364,6 +375,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];
+  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,6 +562,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]; \
+  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,6 +587,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;
+  global_var const double *pd_value = values + NUM_VALUES;
+  global_var const double *pd_weight = pd_value + details->num_weights;
 #endif
 
@@ -637,5 +648,5 @@
       BUILD_ROTATION();
 
-#ifndef USE_OPENCL
+#if !defined(USE_GPU)
       // DLL needs to explicitly loop over the q values.
       #ifdef USE_OPENMP
@@ -643,5 +654,5 @@
       #endif
       for (q_index=0; q_index<nq; q_index++)
-#endif // !USE_OPENCL
+#endif // !USE_GPU
       {
 
@@ -684,9 +695,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
       }
     }
@@ -712,12 +723,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/kernelcuda.py
===================================================================
--- sasmodels/kernelcuda.py	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
+++ sasmodels/kernelcuda.py	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -0,0 +1,491 @@
+"""
+GPU driver for C kernels
+
+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 numpy as np  # type: ignore
+
+
+# Attempt to setup opencl. This may fail if the opencl package is not
+# installed or if it is installed but there are no devices available.
+try:
+    import pycuda.autoinit
+    import pycuda.driver as cuda  # type: ignore
+    from pycuda.compiler import SourceModule
+    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
+
+
+# Pragmas for enable OpenCL features.  Be sure to protect them so that they
+# still compile even if OpenCL is not present.
+_F16_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
+#  pragma OPENCL EXTENSION cl_khr_fp16: enable
+#endif
+"""
+
+_F64_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
+#  pragma OPENCL EXTENSION cl_khr_fp64: enable
+#endif
+"""
+
+def use_cuda():
+    return HAVE_CUDA
+
+ENV = None
+def reset_environment():
+    """
+    Call to create a new OpenCL context, such as after a change to SAS_OPENCL.
+    """
+    global ENV
+    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("OpenCL 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 _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(source, dtype, fast=False):
+    # type: (str, np.dtype, bool) -> cl.Program
+    """
+    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.
+    """
+    source_list = [generate.convert_type(source, dtype)]
+
+    if dtype == generate.F16:
+        source_list.insert(0, _F16_PRAGMA)
+    elif dtype == generate.F64:
+        source_list.insert(0, _F64_PRAGMA)
+
+    source_list.insert(0, "#define USE_SINCOS\n")
+    source = "\n".join(source_list)
+    program = SourceModule(source) # no_extern_c=True, include_dirs=[...]
+    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.
+    """
+    def __init__(self, devnum=0):
+        # type: () -> None
+        # Byte boundary for data alignment
+        #self.data_boundary = max(d.min_data_type_align_size
+        #                         for d in self.context.devices)
+        self.compiled = {}
+        #self.device = cuda.Device(devnum)
+        #self.context = self.device.make_context()
+
+    def has_type(self, dtype):
+        # type: (np.dtype) -> bool
+        """
+        Return True if all devices support a given type.
+        """
+        return True
+
+    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.
+    """
+    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.q_input]
+        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 = 1000000//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.05)
+                    last_nap = current_time
+        sync()
+        details_b.free()
+        values_b.free()
+        cuda.memcpy_dtoh(self.result, self.result_b)
+        #print("result", self.result)
+
+        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.
+        """
+        if self.result_b is not None:
+            self.result_b.free()
+            self.result_b = None
+
+    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/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/models/cylinder.c	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -1,4 +1,5 @@
 #define INVALID(v) (v.radius<0 || v.length<0)
 
+__device__
 static double
 form_volume(double radius, double length)
@@ -7,4 +8,5 @@
 }
 
+__device__
 static double
 fq(double qab, double qc, double radius, double length)
@@ -13,4 +15,5 @@
 }
 
+__device__
 static double
 orient_avg_1D(double q, double radius, double length)
@@ -33,4 +36,5 @@
 }
 
+__device__
 static double
 Iq(double q,
@@ -44,4 +48,5 @@
 }
 
+__device__
 static double
 Iqac(double qab, double qc,
Index: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 99b84ec0c1c8127b9cbe4beeecc251208c44071c)
+++ sasmodels/models/lib/gauss76.c	(revision 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -11,5 +11,5 @@
 
 // Gaussians
-constant double Gauss76Wt[76]={
+constant_var double Gauss76Wt[76] = {
 	.00126779163408536,		//0
 	.00294910295364247,
@@ -90,5 +90,5 @@
 };
 
-constant double Gauss76Z[76]={
+constant_var 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 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -51,6 +51,6 @@
 */
 
-double polevl( double x, constant double *coef, int N );
-double polevl( double x, constant double *coef, int N )
+__device__ static
+double polevl( double x, constant_par 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 )
+__device__ static
+double p1evl( double x, constant_par 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 0db7dbd2e3d48bc33cf18b4a6955d3db0416c80a)
@@ -42,7 +42,6 @@
 #if FLOAT_SIZE>4
 //Cephes double pression function
-double cephes_j1(double x);
-
-constant double RPJ1[8] = {
+
+constant_var double RPJ1[8] = {
     -8.99971225705559398224E8,
     4.52228297998194034323E11,
@@ -54,5 +53,5 @@
     0.0 };
 
-constant double RQJ1[8] = {
+constant_var double RQJ1[8] = {
     6.20836478118054335476E2,
     2.56987256757748830383E5,
@@ -65,5 +64,5 @@
     };
 
-constant double PPJ1[8] = {
+constant_var double PPJ1[8] = {
     7.62125616208173112003E-4,
     7.31397056940917570436E-2,
@@ -76,5 +75,5 @@
 
 
-constant double PQJ1[8] = {
+constant_var double PQJ1[8] = {
     5.71323128072548699714E-4,
     6.88455908754495404082E-2,
@@ -86,5 +85,5 @@
     0.0 };
 
-constant double QPJ1[8] = {
+constant_var double QPJ1[8] = {
     5.10862594750176621635E-2,
     4.98213872951233449420E0,
@@ -96,5 +95,5 @@
     2.52070205858023719784E1 };
 
-constant double QQJ1[8] = {
+constant_var double QQJ1[8] = {
     7.42373277035675149943E1,
     1.05644886038262816351E3,
@@ -106,4 +105,5 @@
     0.0 };
 
+__device__ static
 double cephes_j1(double x)
 {
@@ -155,7 +155,5 @@
 #else
 //Single precission version of cephes
-float cephes_j1f(float x);
-
-constant float JPJ1[8] = {
+constant_var float JPJ1[8] = {
     -4.878788132172128E-009,
     6.009061827883699E-007,
@@ -168,5 +166,5 @@
     };
 
-constant float MO1J1[8] = {
+constant_var float MO1J1[8] = {
     6.913942741265801E-002,
     -2.284801500053359E-001,
@@ -179,5 +177,5 @@
     };
 
-constant float PH1J1[8] = {
+constant_var float PH1J1[8] = {
     -4.497014141919556E+001,
     5.073465654089319E+001,
@@ -190,4 +188,5 @@
     };
 
+__device__ 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);
+__device__ static
 double sas_2J1x_x(double x)
 {
