Index: compare.py
===================================================================
--- compare.py	(revision 87c722ebaa7a1cf12d7bf15ea0c81e18122e3d22)
+++ compare.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -8,5 +8,5 @@
 
 from sasmodels.bumps_model import BumpsModel, plot_data, tic
-from sasmodels import gpu, dll
+from sasmodels import kernelcl, kerneldll
 from sasmodels.convert import revert_model
 
@@ -42,5 +42,5 @@
     sasmodels = __import__('sasmodels.models.'+modelname)
     module = getattr(sasmodels.models, modelname, None)
-    kernel = gpu.load_model(module, dtype=dtype)
+    kernel = kernelcl.load_model(module, dtype=dtype)
     return kernel
 
@@ -48,5 +48,5 @@
     sasmodels = __import__('sasmodels.models.'+modelname)
     module = getattr(sasmodels.models, modelname, None)
-    kernel = dll.load_model(module, dtype=dtype)
+    kernel = kerneldll.load_model(module, dtype=dtype)
     return kernel
 
Index: compare_many.py
===================================================================
--- compare_many.py	(revision 34756fdd2ca589745b9e48738e90455e9bd78d86)
+++ compare_many.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -5,5 +5,5 @@
 import numpy as np
 
-from sasmodels.gpu import environment
+from sasmodels.kernelcl import environment
 from compare import (MODELS, randomize_model, suppress_pd, eval_sasview,
                      eval_opencl, eval_ctypes, make_data)
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 87c722ebaa7a1cf12d7bf15ea0c81e18122e3d22)
+++ sasmodels/bumps_model.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -14,10 +14,10 @@
 
 try:
-    from .gpu import load_model as _loader
+    from .kernelcl import load_model as _loader
 except ImportError,exc:
     import warnings
     warnings.warn(str(exc))
     warnings.warn("OpenCL not available --- using ctypes instead")
-    from .dll import load_model as _loader
+    from .kerneldll import load_model as _loader
 
 def load_model(modelname, dtype='single'):
Index: smodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
+++ 	(revision )
@@ -1,191 +1,0 @@
-"""
-C types wrapper for sasview models.
-"""
-import sys
-import os
-import tempfile
-import ctypes as ct
-from ctypes import c_void_p, c_int, c_double
-
-import numpy as np
-
-from . import generate
-from .pykernel import PyInput, PyKernel
-
-from .generate import F32, F64
-# Compiler platform details
-if sys.platform == 'darwin':
-    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
-    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
-elif os.name == 'nt':
-    # make sure vcvarsall.bat is called first in order to set compiler, headers, lib paths, etc.
-    ##COMPILER = r'"C:\Program Files (x86)\Common Files\Microsoft\Visual C++ for Python\9.0\VC\Bin\cl.exe"'
-    # Can't find VCOMP90.DLL (don't know why), so remove openmp support from windows compiler build
-    #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-    COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-    #COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
-else:
-    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
-
-DLL_PATH = tempfile.gettempdir()
-
-
-def dll_path(info):
-    """
-    Path to the compiled model defined by *info*.
-    """
-    from os.path import join as joinpath, split as splitpath, splitext
-    basename = splitext(splitpath(info['filename'])[1])[0]
-    return joinpath(DLL_PATH, basename+'.so')
-
-
-def load_model(kernel_module, dtype=None):
-    """
-    Load the compiled model defined by *kernel_module*.
-
-    Recompile if any files are newer than the model file.
-
-    *dtype* is ignored.  Compiled files are always double.
-
-    The DLL is not loaded until the kernel is called so models an
-    be defined without using too many resources.
-    """
-    import tempfile
-
-    source, info = generate.make(kernel_module)
-    source_files = generate.sources(info) + [info['filename']]
-    newest = max(os.path.getmtime(f) for f in source_files)
-    dllpath = dll_path(info)
-    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
-        # Replace with a proper temp file
-        fid, filename = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
-        os.fdopen(fid,"w").write(source)
-        command = COMPILE%{"source":filename, "output":dllpath}
-        print "Compile command:",command
-        status = os.system(command)
-        if status != 0:
-            print "compile failed.  File is in %r"%filename
-        else:
-            ## uncomment the following to keep the generated c file
-            #os.unlink(filename); print "saving compiled file in %r"%filename
-            pass
-    return DllModel(dllpath, info)
-
-
-IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double]
-IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int, c_void_p, c_double]
-
-class DllModel(object):
-    """
-    ctypes wrapper for a single model.
-
-    *source* and *info* are the model source and interface as returned
-    from :func:`gen.make`.
-
-    *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.
-
-    Call :meth:`release` when done with the kernel.
-    """
-    def __init__(self, dllpath, info):
-        self.info = info
-        self.dllpath = dllpath
-        self.dll = None
-
-    def _load_dll(self):
-        Nfixed1d = len(self.info['partype']['fixed-1d'])
-        Nfixed2d = len(self.info['partype']['fixed-2d'])
-        Npd1d = len(self.info['partype']['pd-1d'])
-        Npd2d = len(self.info['partype']['pd-2d'])
-
-        #print "dll",self.dllpath
-        self.dll = ct.CDLL(self.dllpath)
-
-        self.Iq = self.dll[generate.kernel_name(self.info, False)]
-        self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed1d + [c_int]*Npd1d
-
-        self.Iqxy = self.dll[generate.kernel_name(self.info, True)]
-        self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed2d + [c_int]*Npd2d
-
-    def __getstate__(self):
-        return {'info': self.info, 'dllpath': self.dllpath, 'dll': None}
-
-    def __setstate__(self, state):
-        self.__dict__ = state
-
-    def __call__(self, input):
-        # Support pure python kernel call
-        if input.is_2D and callable(self.info['Iqxy']):
-            return PyKernel(self.info['Iqxy'], self.info, input)
-        elif not input.is_2D and callable(self.info['Iq']):
-            return PyKernel(self.info['Iq'], self.info, input)
-
-        if self.dll is None: self._load_dll()
-        kernel = self.Iqxy if input.is_2D else self.Iq
-        return DllKernel(kernel, self.info, input)
-
-    def make_input(self, q_vectors):
-        """
-        Make q input vectors available to the model.
-
-        Note that each model needs its own q vector even if the case of
-        mixture models because some models may be OpenCL, some may be
-        ctypes and some may be pure python.
-        """
-        return PyInput(q_vectors, dtype=F64)
-
-    def release(self):
-        pass # TODO: should release the dll
-
-
-class DllKernel(object):
-    """
-    Callable SAS kernel.
-
-    *kernel* is the c function to call.
-
-    *info* is the module information
-
-    *input* is the DllInput q vectors at which the kernel should be
-    evaluated.
-
-    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, info, input):
-        self.info = info
-        self.input = input
-        self.kernel = kernel
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
-        self.fixed_pars = info['partype']['fixed-'+dim]
-        self.pd_pars = info['partype']['pd-'+dim]
-
-        # In dll kernel, but not in opencl kernel
-        self.p_res = self.res.ctypes.data
-
-    def __call__(self, pars, pd_pars, cutoff):
-        real = np.float32 if self.input.dtype == F32 else np.float64
-        fixed = [real(p) for p in pars]
-        cutoff = real(cutoff)
-        loops = np.hstack(pd_pars)
-        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
-        loops_N = [np.uint32(len(p[0])) for p in pd_pars]
-
-        nq = c_int(self.input.nq)
-        p_loops = loops.ctypes.data
-        args = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loops_N
-        #print pars
-        self.kernel(*args)
-
-        return self.res
-
-    def release(self):
-        pass
Index: smodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
+++ 	(revision )
@@ -1,340 +1,0 @@
-"""
-GPU support through OpenCL
-
-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.
-"""
-import numpy as np
-
-import pyopencl as cl
-from pyopencl import mem_flags as mf
-
-from . import generate
-from .pykernel import PyInput, PyKernel
-
-F64_DEFS = """\
-#ifdef cl_khr_fp64
-#  pragma OPENCL EXTENSION cl_khr_fp64: enable
-#endif
-"""
-
-# 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 load_model(kernel_module, dtype="single"):
-    """
-    Load the OpenCL model defined by *kernel_module*.
-
-    Access to the OpenCL device is delayed until the kernel is called
-    so models can be defined without using too many resources.
-    """
-    source, info = generate.make(kernel_module)
-    ## for debugging, save source to a .cl file, edit it, and reload as model
-    #open(info['name']+'.cl','w').write(source)
-    #source = open(info['name']+'.cl','r').read()
-    return GpuModel(source, info, dtype)
-
-ENV = None
-def environment():
-    """
-    Returns a singleton :class:`GpuEnvironment`.
-
-    This provides an OpenCL context and one queue per device.
-    """
-    global ENV
-    if ENV is None:
-        ENV = GpuEnvironment()
-    return ENV
-
-def has_double(device):
-    """
-    Return true if device supports double precision.
-    """
-    return "cl_khr_fp64" in device.extensions
-
-def get_warp(kernel, queue):
-    """
-    Return the size of an execution batch for *kernel* running on *queue*.
-    """
-    return kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
-                                      queue.device)
-
-def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
-    """
-    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):
-    """
-    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 dtype==generate.F64 and not all(has_double(d) for d in context.devices):
-        raise RuntimeError("Double precision not supported for devices")
-
-    header = F64_DEFS if dtype == generate.F64 else ""
-    if dtype == generate.F32:
-        source = generate.use_single(source)
-    # Note: USE_SINCOS makes the intel cpu slower under opencl
-    if context.devices[0].type == cl.device_type.GPU:
-        header += "#define USE_SINCOS\n"
-    program  = cl.Program(context, header+source).build()
-    return program
-
-
-def make_result(self, size):
-    self.res = np.empty(size, dtype=self.dtype)
-    self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes)
-    return self.res, self.res_b
-
-
-# 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):
-        self.context = cl.create_some_context()
-        self.queues = [cl.CommandQueue(self.context, d)
-                       for d in self.context.devices]
-        # Byte boundary for data alignment
-        #self.data_boundary = max(d.min_data_type_align_size
-        #                         for d in self.context.devices)
-        self.has_double = all(has_double(d) for d in self.context.devices)
-        self.compiled = {}
-
-    def compile_program(self, name, source, dtype):
-        if name not in self.compiled:
-            #print "compiling",name
-            self.compiled[name] = compile_model(self.context, source, dtype)
-        return self.compiled[name]
-
-    def release_program(self, name):
-        if name in self.compiled:
-            self.compiled[name].release()
-            del self.compiled[name]
-
-class GpuModel(object):
-    """
-    GPU wrapper for a single model.
-
-    *source* and *info* are the model source and interface as returned
-    from :func:`gen.make`.
-
-    *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.
-    """
-    def __init__(self, source, info, dtype=generate.F32):
-        self.info = info
-        self.source = source
-        self.dtype = dtype
-        self.program = None # delay program creation
-
-    def __getstate__(self):
-        state = self.__dict__.copy()
-        state['program'] = None
-        return state
-
-    def __setstate__(self, state):
-        self.__dict__ = state.copy()
-
-    def __call__(self, input):
-        # Support pure python kernel call
-        if input.is_2D and callable(self.info['Iqxy']):
-            return PyKernel(self.info['Iqxy'], self.info, input)
-        elif not input.is_2D and callable(self.info['Iq']):
-            return PyKernel(self.info['Iq'], self.info, input)
-
-        if self.dtype != input.dtype:
-            raise TypeError("data and kernel have different types")
-        if self.program is None:
-            self.program = environment().compile_program(self.info['name'],self.source, self.dtype)
-        kernel_name = generate.kernel_name(self.info, input.is_2D)
-        kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, input)
-
-    def release(self):
-        if self.program is not None:
-            environment().release_program(self.info['name'])
-            self.program = None
-
-    def make_input(self, q_vectors):
-        """
-        Make q input vectors available to the model.
-
-        Note that each model needs its own q vector even if the case of
-        mixture models because some models may be OpenCL, some may be
-        ctypes and some may be pure python.
-        """
-        # Support pure python kernel call
-        if len(q_vectors) == 1 and callable(self.info['Iq']):
-            return PyInput(q_vectors, dtype=self.dtype)
-        elif callable(self.info['Iqxy']):
-            return PyInput(q_vectors, dtype=self.dtype)
-        else:
-            return GpuInput(q_vectors, dtype=self.dtype)
-
-# 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):
-        env = environment()
-        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.
-        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
-        self.q_buffers = [
-            cl.Buffer(env.context,  mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
-            for q in self.q_vectors
-        ]
-        self.global_size = [self.q_vectors[0].size]
-
-    def release(self):
-        for b in self.q_buffers:
-            b.release()
-        self.q_buffers = []
-
-class GpuKernel(object):
-    """
-    Callable SAS kernel.
-
-    *kernel* is the GpuKernel object to call.
-
-    *info* is the module information
-
-    *input* is the DllInput q vectors at which the kernel should be
-    evaluated.
-
-    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, info, input):
-        self.input = input
-        self.kernel = kernel
-        self.info = info
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
-        self.fixed_pars = info['partype']['fixed-'+dim]
-        self.pd_pars = info['partype']['pd-'+dim]
-
-        # Inputs and outputs for each kernel call
-        # Note: res may be shorter than res_b if global_size != nq
-        env = environment()
-        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  2*MAX_LOOPS*input.dtype.itemsize)
-                        for _ in env.queues]
-        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                input.global_size[0]*input.dtype.itemsize)
-                      for _ in env.queues]
-
-
-    def __call__(self, pars, pd_pars, cutoff=1e-5):
-        real = np.float32 if self.input.dtype == generate.F32 else np.float64
-        fixed = [real(p) for p in pars]
-        cutoff = real(cutoff)
-        loops = np.hstack(pd_pars)
-        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
-        Nloops = [np.uint32(len(p[0])) for p in pd_pars]
-        #print "loops",Nloops, loops
-
-        #import sys; print >>sys.stderr,"opencl eval",pars
-        #print "opencl eval",pars
-        if len(loops) > 2*MAX_LOOPS:
-            raise ValueError("too many polydispersity points")
-        device_num = 0
-        res_bi = self.res_b[device_num]
-        queuei = environment().queues[device_num]
-        loops_bi = self.loops_b[device_num]
-        loops_l = cl.LocalMemory(len(loops.data))
-        cl.enqueue_copy(queuei, loops_bi, loops)
-        #ctx = environment().context
-        #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
-        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + Nloops
-        self.kernel(queuei, self.input.global_size, None, *args)
-        cl.enqueue_copy(queuei, self.res, res_bi)
-
-        return self.res
-
-    def release(self):
-        for b in self.loops_b:
-            b.release()
-        self.loops_b = []
-        for b in self.res_b:
-            b.release()
-        self.res_b = []
-
-    def __del__(self):
-        self.release()
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
+++ sasmodels/kernelcl.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -0,0 +1,340 @@
+"""
+GPU support through OpenCL
+
+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.
+"""
+import numpy as np
+
+import pyopencl as cl
+from pyopencl import mem_flags as mf
+
+from . import generate
+from .kernelpy import PyInput, PyKernel
+
+F64_DEFS = """\
+#ifdef cl_khr_fp64
+#  pragma OPENCL EXTENSION cl_khr_fp64: enable
+#endif
+"""
+
+# 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 load_model(kernel_module, dtype="single"):
+    """
+    Load the OpenCL model defined by *kernel_module*.
+
+    Access to the OpenCL device is delayed until the kernel is called
+    so models can be defined without using too many resources.
+    """
+    source, info = generate.make(kernel_module)
+    ## for debugging, save source to a .cl file, edit it, and reload as model
+    #open(info['name']+'.cl','w').write(source)
+    #source = open(info['name']+'.cl','r').read()
+    return GpuModel(source, info, dtype)
+
+ENV = None
+def environment():
+    """
+    Returns a singleton :class:`GpuEnvironment`.
+
+    This provides an OpenCL context and one queue per device.
+    """
+    global ENV
+    if ENV is None:
+        ENV = GpuEnvironment()
+    return ENV
+
+def has_double(device):
+    """
+    Return true if device supports double precision.
+    """
+    return "cl_khr_fp64" in device.extensions
+
+def get_warp(kernel, queue):
+    """
+    Return the size of an execution batch for *kernel* running on *queue*.
+    """
+    return kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
+                                      queue.device)
+
+def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
+    """
+    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):
+    """
+    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 dtype==generate.F64 and not all(has_double(d) for d in context.devices):
+        raise RuntimeError("Double precision not supported for devices")
+
+    header = F64_DEFS if dtype == generate.F64 else ""
+    if dtype == generate.F32:
+        source = generate.use_single(source)
+    # Note: USE_SINCOS makes the intel cpu slower under opencl
+    if context.devices[0].type == cl.device_type.GPU:
+        header += "#define USE_SINCOS\n"
+    program  = cl.Program(context, header+source).build()
+    return program
+
+
+def make_result(self, size):
+    self.res = np.empty(size, dtype=self.dtype)
+    self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes)
+    return self.res, self.res_b
+
+
+# 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):
+        self.context = cl.create_some_context()
+        self.queues = [cl.CommandQueue(self.context, d)
+                       for d in self.context.devices]
+        # Byte boundary for data alignment
+        #self.data_boundary = max(d.min_data_type_align_size
+        #                         for d in self.context.devices)
+        self.has_double = all(has_double(d) for d in self.context.devices)
+        self.compiled = {}
+
+    def compile_program(self, name, source, dtype):
+        if name not in self.compiled:
+            #print "compiling",name
+            self.compiled[name] = compile_model(self.context, source, dtype)
+        return self.compiled[name]
+
+    def release_program(self, name):
+        if name in self.compiled:
+            self.compiled[name].release()
+            del self.compiled[name]
+
+class GpuModel(object):
+    """
+    GPU wrapper for a single model.
+
+    *source* and *info* are the model source and interface as returned
+    from :func:`gen.make`.
+
+    *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.
+    """
+    def __init__(self, source, info, dtype=generate.F32):
+        self.info = info
+        self.source = source
+        self.dtype = dtype
+        self.program = None # delay program creation
+
+    def __getstate__(self):
+        state = self.__dict__.copy()
+        state['program'] = None
+        return state
+
+    def __setstate__(self, state):
+        self.__dict__ = state.copy()
+
+    def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, input)
+
+        if self.dtype != input.dtype:
+            raise TypeError("data and kernel have different types")
+        if self.program is None:
+            self.program = environment().compile_program(self.info['name'],self.source, self.dtype)
+        kernel_name = generate.kernel_name(self.info, input.is_2D)
+        kernel = getattr(self.program, kernel_name)
+        return GpuKernel(kernel, self.info, input)
+
+    def release(self):
+        if self.program is not None:
+            environment().release_program(self.info['name'])
+            self.program = None
+
+    def make_input(self, q_vectors):
+        """
+        Make q input vectors available to the model.
+
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
+        """
+        # Support pure python kernel call
+        if len(q_vectors) == 1 and callable(self.info['Iq']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        elif callable(self.info['Iqxy']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        else:
+            return GpuInput(q_vectors, dtype=self.dtype)
+
+# 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):
+        env = environment()
+        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.
+        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
+        self.q_buffers = [
+            cl.Buffer(env.context,  mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
+            for q in self.q_vectors
+        ]
+        self.global_size = [self.q_vectors[0].size]
+
+    def release(self):
+        for b in self.q_buffers:
+            b.release()
+        self.q_buffers = []
+
+class GpuKernel(object):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the GpuKernel object to call.
+
+    *info* is the module information
+
+    *input* is the DllInput q vectors at which the kernel should be
+    evaluated.
+
+    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, info, input):
+        self.input = input
+        self.kernel = kernel
+        self.info = info
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        self.fixed_pars = info['partype']['fixed-'+dim]
+        self.pd_pars = info['partype']['pd-'+dim]
+
+        # Inputs and outputs for each kernel call
+        # Note: res may be shorter than res_b if global_size != nq
+        env = environment()
+        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
+                                  2*MAX_LOOPS*input.dtype.itemsize)
+                        for _ in env.queues]
+        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
+                                input.global_size[0]*input.dtype.itemsize)
+                      for _ in env.queues]
+
+
+    def __call__(self, pars, pd_pars, cutoff=1e-5):
+        real = np.float32 if self.input.dtype == generate.F32 else np.float64
+        fixed = [real(p) for p in pars]
+        cutoff = real(cutoff)
+        loops = np.hstack(pd_pars)
+        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+        Nloops = [np.uint32(len(p[0])) for p in pd_pars]
+        #print "loops",Nloops, loops
+
+        #import sys; print >>sys.stderr,"opencl eval",pars
+        #print "opencl eval",pars
+        if len(loops) > 2*MAX_LOOPS:
+            raise ValueError("too many polydispersity points")
+        device_num = 0
+        res_bi = self.res_b[device_num]
+        queuei = environment().queues[device_num]
+        loops_bi = self.loops_b[device_num]
+        loops_l = cl.LocalMemory(len(loops.data))
+        cl.enqueue_copy(queuei, loops_bi, loops)
+        #ctx = environment().context
+        #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
+        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + Nloops
+        self.kernel(queuei, self.input.global_size, None, *args)
+        cl.enqueue_copy(queuei, self.res, res_bi)
+
+        return self.res
+
+    def release(self):
+        for b in self.loops_b:
+            b.release()
+        self.loops_b = []
+        for b in self.res_b:
+            b.release()
+        self.res_b = []
+
+    def __del__(self):
+        self.release()
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
+++ sasmodels/kerneldll.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -0,0 +1,191 @@
+"""
+C types wrapper for sasview models.
+"""
+import sys
+import os
+import tempfile
+import ctypes as ct
+from ctypes import c_void_p, c_int, c_double
+
+import numpy as np
+
+from . import generate
+from .kernelpy import PyInput, PyKernel
+
+from .generate import F32, F64
+# Compiler platform details
+if sys.platform == 'darwin':
+    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
+    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
+elif os.name == 'nt':
+    # make sure vcvarsall.bat is called first in order to set compiler, headers, lib paths, etc.
+    ##COMPILER = r'"C:\Program Files (x86)\Common Files\Microsoft\Visual C++ for Python\9.0\VC\Bin\cl.exe"'
+    # Can't find VCOMP90.DLL (don't know why), so remove openmp support from windows compiler build
+    #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
+    COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
+    #COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
+else:
+    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
+
+DLL_PATH = tempfile.gettempdir()
+
+
+def dll_path(info):
+    """
+    Path to the compiled model defined by *info*.
+    """
+    from os.path import join as joinpath, split as splitpath, splitext
+    basename = splitext(splitpath(info['filename'])[1])[0]
+    return joinpath(DLL_PATH, basename+'.so')
+
+
+def load_model(kernel_module, dtype=None):
+    """
+    Load the compiled model defined by *kernel_module*.
+
+    Recompile if any files are newer than the model file.
+
+    *dtype* is ignored.  Compiled files are always double.
+
+    The DLL is not loaded until the kernel is called so models an
+    be defined without using too many resources.
+    """
+    import tempfile
+
+    source, info = generate.make(kernel_module)
+    source_files = generate.sources(info) + [info['filename']]
+    newest = max(os.path.getmtime(f) for f in source_files)
+    dllpath = dll_path(info)
+    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
+        # Replace with a proper temp file
+        fid, filename = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
+        os.fdopen(fid,"w").write(source)
+        command = COMPILE%{"source":filename, "output":dllpath}
+        print "Compile command:",command
+        status = os.system(command)
+        if status != 0:
+            print "compile failed.  File is in %r"%filename
+        else:
+            ## uncomment the following to keep the generated c file
+            #os.unlink(filename); print "saving compiled file in %r"%filename
+            pass
+    return DllModel(dllpath, info)
+
+
+IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double]
+IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int, c_void_p, c_double]
+
+class DllModel(object):
+    """
+    ctypes wrapper for a single model.
+
+    *source* and *info* are the model source and interface as returned
+    from :func:`gen.make`.
+
+    *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.
+
+    Call :meth:`release` when done with the kernel.
+    """
+    def __init__(self, dllpath, info):
+        self.info = info
+        self.dllpath = dllpath
+        self.dll = None
+
+    def _load_dll(self):
+        Nfixed1d = len(self.info['partype']['fixed-1d'])
+        Nfixed2d = len(self.info['partype']['fixed-2d'])
+        Npd1d = len(self.info['partype']['pd-1d'])
+        Npd2d = len(self.info['partype']['pd-2d'])
+
+        #print "dll",self.dllpath
+        self.dll = ct.CDLL(self.dllpath)
+
+        self.Iq = self.dll[generate.kernel_name(self.info, False)]
+        self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed1d + [c_int]*Npd1d
+
+        self.Iqxy = self.dll[generate.kernel_name(self.info, True)]
+        self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed2d + [c_int]*Npd2d
+
+    def __getstate__(self):
+        return {'info': self.info, 'dllpath': self.dllpath, 'dll': None}
+
+    def __setstate__(self, state):
+        self.__dict__ = state
+
+    def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, input)
+
+        if self.dll is None: self._load_dll()
+        kernel = self.Iqxy if input.is_2D else self.Iq
+        return DllKernel(kernel, self.info, input)
+
+    def make_input(self, q_vectors):
+        """
+        Make q input vectors available to the model.
+
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
+        """
+        return PyInput(q_vectors, dtype=F64)
+
+    def release(self):
+        pass # TODO: should release the dll
+
+
+class DllKernel(object):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the c function to call.
+
+    *info* is the module information
+
+    *input* is the DllInput q vectors at which the kernel should be
+    evaluated.
+
+    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, info, input):
+        self.info = info
+        self.input = input
+        self.kernel = kernel
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        self.fixed_pars = info['partype']['fixed-'+dim]
+        self.pd_pars = info['partype']['pd-'+dim]
+
+        # In dll kernel, but not in opencl kernel
+        self.p_res = self.res.ctypes.data
+
+    def __call__(self, pars, pd_pars, cutoff):
+        real = np.float32 if self.input.dtype == F32 else np.float64
+        fixed = [real(p) for p in pars]
+        cutoff = real(cutoff)
+        loops = np.hstack(pd_pars)
+        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+        loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+
+        nq = c_int(self.input.nq)
+        p_loops = loops.ctypes.data
+        args = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loops_N
+        #print pars
+        self.kernel(*args)
+
+        return self.res
+
+    def release(self):
+        pass
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
+++ sasmodels/kernelpy.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -0,0 +1,209 @@
+import numpy as np
+from numpy import pi, sin, cos, sqrt
+
+from .generate import F32, F64
+
+class PyInput(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):
+        self.nq = q_vectors[0].size
+        self.dtype = dtype
+        self.is_2D = (len(q_vectors) == 2)
+        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
+        self.q_pointers = [q.ctypes.data for q in q_vectors]
+
+    def release(self):
+        self.q_vectors = []
+
+class PyKernel(object):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the DllKernel object to call.
+
+    *info* is the module information
+
+    *input* is the DllInput q vectors at which the kernel should be
+    evaluated.
+
+    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, info, input):
+        self.info = info
+        self.input = input
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        # Loop over q unless user promises that the kernel is vectorized by
+        # taggining it with vectorized=True
+        if not getattr(kernel, 'vectorized', False):
+            if dim == '2d':
+                def vector_kernel(qx, qy, *args):
+                    return np.array([kernel(qxi,qyi,*args) for qxi,qyi in zip(qx,qy)])
+            else:
+                def vector_kernel(q, *args):
+                    return np.array([kernel(qi,*args) for qi in q])
+            self.kernel = vector_kernel
+        else:
+            self.kernel = kernel
+        fixed_pars = info['partype']['fixed-'+dim]
+        pd_pars = info['partype']['pd-'+dim]
+        vol_pars = info['partype']['volume']
+
+        # First two fixed pars are scale and background
+        pars = [p[0] for p in info['parameters'][2:]]
+        offset = len(self.input.q_vectors)
+        self.args = self.input.q_vectors + [None]*len(pars)
+        self.fixed_index = np.array([pars.index(p)+offset for p in fixed_pars[2:] ])
+        self.pd_index = np.array([pars.index(p)+offset for p in pd_pars])
+        self.vol_index = np.array([pars.index(p)+offset for p in vol_pars])
+        try: self.theta_index = pars.index('theta')+offset
+        except ValueError: self.theta_index = -1
+
+        # Caller needs fixed_pars and pd_pars
+        self.fixed_pars = fixed_pars
+        self.pd_pars = pd_pars
+
+    def __call__(self, fixed, pd, cutoff):
+        #print "fixed",fixed
+        #print "pd", pd
+        args = self.args[:]  # grab a copy of the args
+        form, form_volume = self.kernel, self.info['form_volume']
+        # First two fixed
+        scale, background = fixed[:2]
+        for index,value in zip(self.fixed_index, fixed[2:]):
+            args[index] = float(value)
+        res = _loops(form, form_volume, cutoff, scale, background,  args,
+                     pd, self.pd_index, self.vol_index, self.theta_index)
+
+        return res
+
+    def release(self):
+        self.input = None
+
+def _loops(form, form_volume, cutoff, scale, background,
+           args, pd, pd_index, vol_index, theta_index):
+    """
+    *form* is the name of the form function, which should be vectorized over
+    q, but otherwise have an interface like the opencl kernels, with the
+    q parameters first followed by the individual parameters in the order
+    given in model.parameters (see :mod:`sasmodels.generate`).
+
+    *form_volume* calculates the volume of the shape.  *vol_index* gives
+    the list of volume parameters
+
+    *cutoff* ignores the corners of the dispersion hypercube
+
+    *scale*, *background* multiplies the resulting form and adds an offset
+
+    *args* is the prepopulated set of arguments to the form function, starting
+    with the q vectors, and including slots for all the parameters.  The
+    values for the parameters will be substituted with values from the
+    dispersion functions.
+
+    *pd* is the list of dispersion parameters
+
+    *pd_index* are the indices of the dispersion parameters in *args*
+
+    *vol_index* are the indices of the volume parameters in *args*
+
+    *theta_index* is the index of the theta parameter for the sasview
+    spherical correction, or -1 if there is no angular dispersion
+    """
+
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH GENERATE.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+
+    # weight vector, to be populated by polydispersity loops
+    weight = np.empty(len(pd), 'd')
+
+    # identify which pd parameters are volume parameters
+    vol_weight_index = np.array([(index in vol_index) for index in pd_index])
+
+    # Sort parameters in decreasing order of pd length
+    Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
+    order = np.argsort(Npd)[::-1]
+    stride = np.cumprod(Npd[order])
+    pd = [pd[index] for index in order]
+    pd_index = pd_index[order]
+    vol_weight_index = vol_weight_index[order]
+
+    fast_value = pd[0][0]
+    fast_weight = pd[0][1]
+
+    ret = np.zeros_like(args[0])
+    norm = np.zeros_like(ret)
+    vol = np.zeros_like(ret)
+    vol_norm = np.zeros_like(ret)
+    for k in range(stride[-1]):
+        # update polydispersity parameter values
+        fast_index = k%stride[0]
+        if fast_index == 0:  # bottom loop complete ... check all other loops
+            for i,index, in enumerate(k%stride):
+                args[pd_index[i]] = pd[i][0][index]
+                weight[i] = pd[i][1][index]
+        else:
+            args[pd_index[0]] = fast_value[fast_index]
+            weight[0] = fast_weight[fast_index]
+        # This computes the weight, and if it is sufficient, calls the scattering
+        # function and adds it to the total.  If there is a volume normalization,
+        # it will also be added here.
+        # Note: make sure this is consistent with the code in PY_LOOP_BODY!!
+        # Note: can precompute w1*w2*...*wn
+        w = np.prod(weight)
+        if w > cutoff:
+            I = form(*args)
+            positive = (I>=0.0)
+
+            # Note: can precompute spherical correction if theta_index is not the fast index
+            # Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
+            #spherical_correction = abs(sin(pi*args[theta_index])) if theta_index>=0 else 1.0
+            #spherical_correction = abs(cos(pi*args[theta_index]))*pi/2 if theta_index>=0 else 1.0
+            spherical_correction = 1.0
+            ret += w*I*spherical_correction*positive
+            norm += w*positive
+
+            # Volume normalization.
+            # If there are "volume" polydispersity parameters, then these will be used
+            # to call the form_volume function from the user supplied kernel, and accumulate
+            # a normalized weight.
+            # Note: can precompute volume norm if the fast index is not a volume index
+            if form_volume:
+                vol_args = [args[index] for index in vol_index]
+                vol_weight = np.prod(weight[vol_weight_index])
+                vol += vol_weight*form_volume(*vol_args)*positive
+                vol_norm += vol_weight*positive
+
+    positive = (vol*vol_norm != 0.0)
+    ret[positive] *= vol_norm[positive] / vol[positive]
+    result = scale*ret/norm+background
+    return result
Index: smodels/pykernel.py
===================================================================
--- sasmodels/pykernel.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
+++ 	(revision )
@@ -1,209 +1,0 @@
-import numpy as np
-from numpy import pi, sin, cos, sqrt
-
-from .generate import F32, F64
-
-class PyInput(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):
-        self.nq = q_vectors[0].size
-        self.dtype = dtype
-        self.is_2D = (len(q_vectors) == 2)
-        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
-        self.q_pointers = [q.ctypes.data for q in q_vectors]
-
-    def release(self):
-        self.q_vectors = []
-
-class PyKernel(object):
-    """
-    Callable SAS kernel.
-
-    *kernel* is the DllKernel object to call.
-
-    *info* is the module information
-
-    *input* is the DllInput q vectors at which the kernel should be
-    evaluated.
-
-    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, info, input):
-        self.info = info
-        self.input = input
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
-        # Loop over q unless user promises that the kernel is vectorized by
-        # taggining it with vectorized=True
-        if not getattr(kernel, 'vectorized', False):
-            if dim == '2d':
-                def vector_kernel(qx, qy, *args):
-                    return np.array([kernel(qxi,qyi,*args) for qxi,qyi in zip(qx,qy)])
-            else:
-                def vector_kernel(q, *args):
-                    return np.array([kernel(qi,*args) for qi in q])
-            self.kernel = vector_kernel
-        else:
-            self.kernel = kernel
-        fixed_pars = info['partype']['fixed-'+dim]
-        pd_pars = info['partype']['pd-'+dim]
-        vol_pars = info['partype']['volume']
-
-        # First two fixed pars are scale and background
-        pars = [p[0] for p in info['parameters'][2:]]
-        offset = len(self.input.q_vectors)
-        self.args = self.input.q_vectors + [None]*len(pars)
-        self.fixed_index = np.array([pars.index(p)+offset for p in fixed_pars[2:] ])
-        self.pd_index = np.array([pars.index(p)+offset for p in pd_pars])
-        self.vol_index = np.array([pars.index(p)+offset for p in vol_pars])
-        try: self.theta_index = pars.index('theta')+offset
-        except ValueError: self.theta_index = -1
-
-        # Caller needs fixed_pars and pd_pars
-        self.fixed_pars = fixed_pars
-        self.pd_pars = pd_pars
-
-    def __call__(self, fixed, pd, cutoff):
-        #print "fixed",fixed
-        #print "pd", pd
-        args = self.args[:]  # grab a copy of the args
-        form, form_volume = self.kernel, self.info['form_volume']
-        # First two fixed
-        scale, background = fixed[:2]
-        for index,value in zip(self.fixed_index, fixed[2:]):
-            args[index] = float(value)
-        res = _loops(form, form_volume, cutoff, scale, background,  args,
-                     pd, self.pd_index, self.vol_index, self.theta_index)
-
-        return res
-
-    def release(self):
-        self.input = None
-
-def _loops(form, form_volume, cutoff, scale, background,
-           args, pd, pd_index, vol_index, theta_index):
-    """
-    *form* is the name of the form function, which should be vectorized over
-    q, but otherwise have an interface like the opencl kernels, with the
-    q parameters first followed by the individual parameters in the order
-    given in model.parameters (see :mod:`sasmodels.generate`).
-
-    *form_volume* calculates the volume of the shape.  *vol_index* gives
-    the list of volume parameters
-
-    *cutoff* ignores the corners of the dispersion hypercube
-
-    *scale*, *background* multiplies the resulting form and adds an offset
-
-    *args* is the prepopulated set of arguments to the form function, starting
-    with the q vectors, and including slots for all the parameters.  The
-    values for the parameters will be substituted with values from the
-    dispersion functions.
-
-    *pd* is the list of dispersion parameters
-
-    *pd_index* are the indices of the dispersion parameters in *args*
-
-    *vol_index* are the indices of the volume parameters in *args*
-
-    *theta_index* is the index of the theta parameter for the sasview
-    spherical correction, or -1 if there is no angular dispersion
-    """
-
-    ##########################################################
-    #                                                        #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #   !!                                              !!   #
-    #   !!  KEEP THIS CODE CONSISTENT WITH GENERATE.PY  !!   #
-    #   !!                                              !!   #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #                                                        #
-    ##########################################################
-
-    # weight vector, to be populated by polydispersity loops
-    weight = np.empty(len(pd), 'd')
-
-    # identify which pd parameters are volume parameters
-    vol_weight_index = np.array([(index in vol_index) for index in pd_index])
-
-    # Sort parameters in decreasing order of pd length
-    Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
-    order = np.argsort(Npd)[::-1]
-    stride = np.cumprod(Npd[order])
-    pd = [pd[index] for index in order]
-    pd_index = pd_index[order]
-    vol_weight_index = vol_weight_index[order]
-
-    fast_value = pd[0][0]
-    fast_weight = pd[0][1]
-
-    ret = np.zeros_like(args[0])
-    norm = np.zeros_like(ret)
-    vol = np.zeros_like(ret)
-    vol_norm = np.zeros_like(ret)
-    for k in range(stride[-1]):
-        # update polydispersity parameter values
-        fast_index = k%stride[0]
-        if fast_index == 0:  # bottom loop complete ... check all other loops
-            for i,index, in enumerate(k%stride):
-                args[pd_index[i]] = pd[i][0][index]
-                weight[i] = pd[i][1][index]
-        else:
-            args[pd_index[0]] = fast_value[fast_index]
-            weight[0] = fast_weight[fast_index]
-        # This computes the weight, and if it is sufficient, calls the scattering
-        # function and adds it to the total.  If there is a volume normalization,
-        # it will also be added here.
-        # Note: make sure this is consistent with the code in PY_LOOP_BODY!!
-        # Note: can precompute w1*w2*...*wn
-        w = np.prod(weight)
-        if w > cutoff:
-            I = form(*args)
-            positive = (I>=0.0)
-
-            # Note: can precompute spherical correction if theta_index is not the fast index
-            # Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
-            #spherical_correction = abs(sin(pi*args[theta_index])) if theta_index>=0 else 1.0
-            #spherical_correction = abs(cos(pi*args[theta_index]))*pi/2 if theta_index>=0 else 1.0
-            spherical_correction = 1.0
-            ret += w*I*spherical_correction*positive
-            norm += w*positive
-
-            # Volume normalization.
-            # If there are "volume" polydispersity parameters, then these will be used
-            # to call the form_volume function from the user supplied kernel, and accumulate
-            # a normalized weight.
-            # Note: can precompute volume norm if the fast index is not a volume index
-            if form_volume:
-                vol_args = [args[index] for index in vol_index]
-                vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight*form_volume(*vol_args)*positive
-                vol_norm += vol_weight*positive
-
-    positive = (vol*vol_norm != 0.0)
-    ret[positive] *= vol_norm[positive] / vol[positive]
-    result = scale*ret/norm+background
-    return result
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 87c722ebaa7a1cf12d7bf15ea0c81e18122e3d22)
+++ sasmodels/sasview_model.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
@@ -25,9 +25,9 @@
 
 try:
-    from .gpu import load_model
+    from .kernelcl import load_model
 except ImportError,exc:
     warnings.warn(str(exc))
     warnings.warn("OpenCL not available --- using ctypes instead")
-    from .dll import load_model
+    from .kerneldll import load_model
 
 
