Index: .gitignore
===================================================================
--- .gitignore	(revision 4a82d4dd41f712a03b7f365ffe63e370cf57b792)
+++ .gitignore	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -4,5 +4,4 @@
 /*.csv
 *.pyc
-*.cl
 *.so
 *.obj
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 8d620089dc0fbbdb285c976744abfa4439490599)
+++ sasmodels/compare.py	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -497,4 +497,5 @@
         value = calculator(**pars)
     average_time = toc()*1000./Nevals
+    #print("I(q)",value)
     return value, average_time
 
@@ -604,5 +605,5 @@
     if Nbase > 0 and Ncomp > 0:
         resid = (base_value - comp_value)
-        relerr = resid/comp_value
+        relerr = resid/np.where(comp_value!=0., abs(comp_value), 1.0)
         _print_stats("|%s-%s|"
                      % (base.engine, comp.engine) + (" "*(3+len(comp.engine))),
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision dd7fc126b4eb39f92db8393c7093392becb55218)
+++ sasmodels/core.py	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -142,6 +142,8 @@
             or not HAVE_OPENCL
             or not kernelcl.environment().has_type(numpy_dtype)):
+        #print("building dll", numpy_dtype)
         return kerneldll.load_dll(source, model_info, numpy_dtype)
     else:
+        #print("building ocl", numpy_dtype)
         return kernelcl.GpuModel(source, model_info, numpy_dtype, fast=fast)
 
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
+++ sasmodels/generate.py	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -173,5 +173,5 @@
 
 try:
-    from typing import Tuple, Sequence, Iterator
+    from typing import Tuple, Sequence, Iterator, Dict
     from .modelinfo import ModelInfo
 except ImportError:
@@ -298,4 +298,7 @@
     Return a timestamp for the model corresponding to the most recently
     changed file or dependency.
+
+    Note that this does not look at the time stamps for the OpenCL header
+    information since that need not trigger a recompile of the DLL.
     """
     source_files = (model_sources(model_info)
@@ -304,4 +307,13 @@
     newest = max(getmtime(f) for f in source_files)
     return newest
+
+def model_templates():
+    # type: () -> List[str]
+    # TODO: fails DRY; templates appear two places.
+    # should instead have model_info contain a list of paths
+    # Note: kernel_iq.cl is not on this list because changing it need not
+    # trigger a recompile of the dll.
+    return [joinpath(TEMPLATE_ROOT, filename)
+            for filename in ('kernel_header.c', 'kernel_iq.c')]
 
 def convert_type(source, dtype):
@@ -377,11 +389,4 @@
     return _template_cache[filename][1]
 
-def model_templates():
-    # type: () -> List[str]
-    # TODO: fails DRY; templates are listed in two places.
-    # should instead have model_info contain a list of paths
-    return [joinpath(TEMPLATE_ROOT, filename)
-            for filename in ('kernel_header.c', 'kernel_iq.c')]
-
 
 _FN_TEMPLATE = """\
@@ -390,5 +395,4 @@
     %(body)s
 }
-
 
 """
@@ -467,5 +471,6 @@
     # Load templates and user code
     kernel_header = load_template('kernel_header.c')
-    kernel_code = load_template('kernel_iq.c')
+    dll_code = load_template('kernel_iq.c')
+    ocl_code = load_template('kernel_iq.cl')
     user_code = [open(f).read() for f in model_sources(model_info)]
 
@@ -506,5 +511,5 @@
     if _have_Iqxy(user_code):
         # Call 2D model
-        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters)
+        refs = ["q[2*_i]", "q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters)
         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs))
     else:
@@ -521,22 +526,39 @@
     # TODO: allow mixed python/opencl kernels?
 
-    # define the Iq kernel
-    source.append("#define KERNEL_NAME %s_Iq"%model_info.name)
-    source.append(call_iq)
-    source.append(kernel_code)
-    source.append("#undef CALL_IQ")
-    source.append("#undef KERNEL_NAME")
-
-    # define the Iqxy kernel from the same source with different #defines
-    source.append("#define KERNEL_NAME %s_Iqxy"%model_info.name)
-    source.append(call_iqxy)
-    source.append(kernel_code)
-    source.append("#undef CALL_IQ")
-    source.append("#undef KERNEL_NAME")
-
+    source.append("#if defined(USE_OPENCL)")
+    source.extend(_add_kernels(ocl_code, call_iq, call_iqxy, model_info.name))
+    source.append("#else /* !USE_OPENCL */")
+    source.extend(_add_kernels(dll_code, call_iq, call_iqxy, model_info.name))
+    source.append("#endif /* !USE_OPENCL */")
     return '\n'.join(source)
+
+def _add_kernels(kernel_code, call_iq, call_iqxy, name):
+    # type: (str, str, str, str) -> List[str]
+    source = [
+        # define the Iq kernel
+        "#define KERNEL_NAME %s_Iq"%name,
+        call_iq,
+        kernel_code,
+        "#undef CALL_IQ",
+        "#undef KERNEL_NAME",
+
+        # define the Iqxy kernel from the same source with different #defines
+        "#define KERNEL_NAME %s_Iqxy"%name,
+        call_iqxy,
+        kernel_code,
+        "#undef CALL_IQ",
+        "#undef KERNEL_NAME",
+    ]
+    return source
 
 def load_kernel_module(model_name):
     # type: (str) -> module
+    """
+    Return the kernel module named in *model_name*.
+
+    If the name ends in *.py* then load it as a custom model using
+    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
+    load it from :mod:`sasmodels.models`.
+    """
     if model_name.endswith('.py'):
         kernel_module = load_custom_kernel_module(model_name)
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 6e7ff6d026d856e1eb236a05975a2600ba9ba2c9)
+++ sasmodels/kernel_iq.c	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -54,4 +54,8 @@
   local ParameterBlock local_values;  // current parameter values
   double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
+  double norm;
+
+  // number of active loops
+  const int num_active = details->num_active;
 
   // Fill in the initial variables
@@ -63,4 +67,26 @@
   }
 
+  // Monodisperse computation
+  if (num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+    norm = CALL_VOLUME(local_values);
+
+    const double scale = values[0];
+    const double background = values[1];
+    // result[nq] = norm; // Total volume normalization
+
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      double scattering = CALL_IQ(q, i, local_values);
+      result[i] = (norm>0. ? scale*scattering/norm + background : background);
+    }
+    return;
+  }
+
+#if MAX_PD > 0
   // If it is the first round initialize the result to zero, otherwise
   // assume that the previous result has been passed back.
@@ -75,31 +101,8 @@
       result[i] = 0.0;
     }
-  }
-
-  // Monodisperse computation
-  if (details->num_active == 0) {
-    #ifdef INVALID
-    if (INVALID(local_values)) { return; }
-    #endif
-
-    const double norm = CALL_VOLUME(local_values);
-    const double scale = values[0];
-    const double background = values[1];
-    #ifdef USE_OPENMP
-    #pragma omp parallel for
-    #endif
-    result[nq] = norm; // Total volume normalization
-    for (int i=0; i < nq; i++) {
-      double scattering = CALL_IQ(q, i, local_values);
-      result[i] = (norm>0. ? scale*scattering/norm + background : background);
-    }
-    return;
-  }
-
-#if MAX_PD > 0
-  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
-  // Since we are no longer looping over the entire polydispersity hypercube
-  // for each q, we need to track the normalization values between calls.
-  double norm = 0.0;
+    norm = 0.0;
+  } else {
+    norm = result[nq];
+  }
 
   // need product of weights at every Iq calc, so keep product of
@@ -119,12 +122,15 @@
   pd_index[0] = fast_length;
 
+  // Number of coordinated indices
+  const int num_coord = details->num_coord;
+
   // Loop over the weights then loop over q, accumulating values
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
-    // check if indices need to be updated
+    // check if fast loop needs to be reset
     if (pd_index[0] == fast_length) {
-      //printf("should be here with %d active\n", details->num_active);
+      //printf("should be here with %d active\n", num_active);
 
       // Compute position in polydispersity hypercube
-      for (int k=0; k < details->num_active; k++) {
+      for (int k=0; k < num_active; k++) {
         pd_index[k] = (loop_index/details->pd_stride[k])%details->pd_length[k];
         //printf("pd_index[%d] = %d\n",k,pd_index[k]);
@@ -134,5 +140,5 @@
       partial_weight = 1.0;
       //printf("partial weight %d: ", loop_index);
-      for (int k=1; k < details->num_active; k++) {
+      for (int k=1; k < num_active; k++) {
         double wi = weights[details->pd_offset[k] + pd_index[k]];
         //printf("pd[%d]=par[%d]=%g ", k, details->pd_par[k], wi);
@@ -143,5 +149,5 @@
       // Update parameter offsets in weight vector
       //printf("slow %d: ", loop_index);
-      for (int k=0; k < details->num_coord; k++) {
+      for (int k=0; k < num_coord; k++) {
         int par = details->par_coord[k];
         int coord = details->pd_coord[k];
@@ -160,5 +166,5 @@
         // if theta is not coordinated with fast index, precompute spherical correction
         if (par == details->theta_par && !(details->par_coord[k]&1)) {
-          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
         }
       }
@@ -170,5 +176,5 @@
     double weight = partial_weight*wi;
     //printf("fast %d: ", loop_index);
-    for (int k=0; k < details->num_coord; k++) {
+    for (int k=0; k < num_coord; k++) {
       if (details->pd_coord[k]&1) {
         const int par = details->par_coord[k];
@@ -177,5 +183,5 @@
         // if theta is coordinated with fast index, compute spherical correction each time
         if (par == details->theta_par) {
-          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6);
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
         }
       }
@@ -205,7 +211,4 @@
   }
 
-  // Accumulate norm.
-  result[nq] += norm;
-
   // End of the PD loop we can normalize
   if (pd_stop >= details->total_pd) {
@@ -219,4 +222,7 @@
     }
   }
+
+  // Remember the updated norm.
+  result[nq] = norm;
 #endif // MAX_PD > 0
 }
Index: sasmodels/kernel_iq.cl
===================================================================
--- sasmodels/kernel_iq.cl	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
+++ sasmodels/kernel_iq.cl	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -0,0 +1,243 @@
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
+
+typedef struct {
+#if MAX_PD > 0
+    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
+    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
+    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
+    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
+#endif // MAX_PD > 0
+    int32_t par_offset[NPARS];  // offset of par value blocks in the value & weight vector
+    int32_t par_coord[NPARS];   // ids of the coordination parameters
+    int32_t pd_coord[NPARS];    // polydispersity coordination bitvector
+    int32_t num_active;         // number of non-trivial pd loops
+    int32_t total_pd;           // total number of voxels in hypercube
+    int32_t num_coord;          // number of coordinated parameters
+    int32_t theta_par;          // id of spherical correction variable
+} ProblemDetails;
+
+typedef struct {
+    PARAMETER_TABLE;
+} ParameterBlock;
+#endif
+
+
+kernel
+void KERNEL_NAME(
+    int32_t nq,                 // number of q values
+    const int32_t pd_start,     // where we are in the polydispersity loop
+    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
+    global const ProblemDetails *details,
+    global const double *weights,
+    global const double *values,
+    global const double *q, // nq q values, with padding to boundary
+    global double *result,  // nq+3 return values, again with padding
+    const double cutoff     // cutoff in the polydispersity weight product
+    )
+{
+  double norm;
+
+  // who we are and what element we are working with
+  const int q_index = get_global_id(0);
+
+  // number of active loops
+  const int num_active = details->num_active;
+
+  // Storage for the current parameter values.  These will be updated as we
+  // walk the polydispersity cube.
+  ParameterBlock local_values;  // current parameter values
+  double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
+  // Fill in the initial variables
+  for (int k = 0; k < NPARS; k++) {
+    pvec[k] = values[details->par_offset[k]];
+  }
+
+  // Monodisperse computation
+  if (num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+
+    double scale, background;
+    norm = CALL_VOLUME(local_values);
+    scale = values[0];
+    background = values[1];
+
+    // if (i==0) result[nq] = norm; // Total volume normalization
+
+    if (q_index < nq) {
+      double scattering = CALL_IQ(q, q_index, local_values);
+      result[q_index] = (norm>0. ? scale*scattering/norm + background : background);
+    }
+    return;
+  }
+
+#if MAX_PD > 0
+
+  // If it is the first round initialize the result to zero, otherwise
+  // assume that the previous result has been passed back.
+  // Note: doing this even in the monodisperse case in order to handle the
+  // rare case where the model parameters are invalid and zero is returned.
+  // So slightly increased cost for slightly smaller code size.
+  double this_result;
+
+  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
+  // norm will be shared across all threads.
+
+  // need product of weights at every Iq calc, so keep product of
+  // weights from the outer loops so that weight = partial_weight * fast_weight
+  double partial_weight; // product of weight w4*w3*w2 but not w1
+  double spherical_correction;  // cosine correction for latitude variation
+
+  // Location in the polydispersity hypercube, one index per dimension.
+  int pd_index[MAX_PD];
+
+  // Location of the coordinated parameters in their own sub-cubes.
+  int offset[NPARS];
+
+  // Number of elements in the longest polydispersity loop
+  const int fast_length = details->pd_length[0];
+
+  // Number of coordinated indices
+  const int num_coord = details->num_coord;
+
+  // We could in theory spread this work across different threads, but
+  // lets keep it simple;
+  norm = pd_start == 0 ? 0.0 : result[nq];
+  spherical_correction = 1.0;  // the usual case.
+  // partial_weight = NAN;
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  pd_index[0] = fast_length;
+
+  // Since we are no longer looping over the entire polydispersity hypercube
+  // for each q, we need to track the result and normalization values between
+  // calls.  This means initializing them to 0 at the start and accumulating
+  // them between calls.
+  if (q_index < nq) {
+    this_result = pd_start == 0 ? 0.0 : result[q_index];
+  }
+
+  // Loop over the weights then loop over q, accumulating values
+  for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
+    // check if fast loop needs to be reset
+    if (pd_index[0] == fast_length) {
+      //printf("should be here with %d active\n", num_active);
+
+      // Compute position in polydispersity hypercube
+      for (int k=0; k < num_active; k++) {
+          pd_index[k] = (loop_index/details->pd_stride[k])%details->pd_length[k];
+          //printf("pd_index[%d] = %d\n",k,pd_index[k]);
+      }
+
+      // need to compute the product of the weights.  If the vector were really
+      // long, we could split the work into groups, with each thread taking
+      // every nth weight, but there really is no call for it here.  We could
+      // also do some clever pair-wise multiplication similar to parallel
+      // prefix, but again simpler is probably faster since n is likely small.
+      // Compute partial weights
+      partial_weight = 1.0;
+      //printf("partial weight %d: ", loop_index);
+      for (int k=1; k < num_active; k++) {
+        double wi = weights[details->pd_offset[k] + pd_index[k]];
+        //printf("pd[%d]=par[%d]=%g ", k, details->pd_par[k], wi);
+        partial_weight *= wi;
+      }
+      //printf("\n");
+
+      // Update parameter offsets in weight vector
+      //printf("slow %d: ", loop_index);
+      for (int k=0; k < num_coord; k++) {
+        if (k < num_coord) {
+          int par = details->par_coord[k];
+          int coord = details->pd_coord[k];
+          int this_offset = details->par_offset[par];
+          int block_size = 1;
+          for (int bit=0; coord != 0; bit++) {
+            if (coord&1) {
+                this_offset += block_size * pd_index[bit];
+                block_size *= details->pd_length[bit];
+            }
+            coord >>= 1;
+          }
+          offset[par] = this_offset;
+          pvec[par] = values[this_offset];
+          //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]);
+          // if theta is not coordinated with fast index, precompute spherical correction
+          if (par == details->theta_par && !(details->par_coord[k]&1)) {
+            spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
+          }
+        }
+      }
+      //printf("\n");
+    }
+
+    double weight;
+    const double wi = weights[details->pd_offset[0] + pd_index[0]];
+    weight = partial_weight*wi;
+    pd_index[0]++;
+
+    // Increment fast index
+    //printf("fast %d: ", loop_index);
+    for (int k=0; k < num_coord; k++) {
+      if (k < num_coord) {
+        if (details->pd_coord[k]&1) {
+          const int par = details->par_coord[k];
+          pvec[par] = values[offset[par]++];
+          //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]);
+          // if theta is coordinated with fast index, compute spherical correction each time
+          if (par == details->theta_par) {
+            spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1.e-6);
+          }
+        }
+      }
+    }
+    //printf("\n");
+
+    #ifdef INVALID
+    if (INVALID(local_values)) continue;
+    #endif
+
+    // Accumulate I(q)
+    // Note: weight==0 must always be excluded
+    if (weight > cutoff) {
+      // spherical correction has some nasty effects when theta is +90 or -90
+      // where it becomes zero.  If the entirety of the correction
+      weight *= spherical_correction;
+      norm += weight * CALL_VOLUME(local_values);
+
+      const double scattering = CALL_IQ(q, q_index, local_values);
+      this_result += weight*scattering;
+    }
+  }
+
+  if (q_index < nq) {
+    if (pd_stop >= details->total_pd) {
+      // End of the PD loop we can normalize
+      const double scale = values[0];
+      const double background = values[1];
+      result[q_index] = (norm>0. ? scale*this_result/norm + background : background);
+    } else {
+      // Partial result, so remember it but don't normalize it.
+      result[q_index] = this_result;
+    }
+    // Accumulate norm.
+    if (q_index == 0) result[nq] = norm;
+  }
+
+#endif // MAX_PD > 0
+}
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 8d620089dc0fbbdb285c976744abfa4439490599)
+++ sasmodels/kernelcl.py	(revision f2f67a6b8ebed2920432485a4351f9ac24f3cbf2)
@@ -56,5 +56,5 @@
 
 try:
-    raise NotImplementedError("OpenCL not yet implemented for new kernel template")
+    #raise NotImplementedError("OpenCL not yet implemented for new kernel template")
     import pyopencl as cl  # type: ignore
     # Ask OpenCL for the default context so that we know that one exists
@@ -264,5 +264,5 @@
         key = "%s-%s-%s"%(name, dtype, fast)
         if key not in self.compiled:
-            print("compiling",name)
+            #print("OpenCL compile",name)
             dtype = np.dtype(dtype)
             program = compile_model(self.get_context(dtype),
@@ -373,5 +373,5 @@
         kernel_name = generate.kernel_name(self.info, is_2d)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, q_vectors)
+        return GpuKernel(kernel, self.dtype, self.info, q_vectors)
 
     def release(self):
@@ -443,7 +443,7 @@
         Free the memory.
         """
-        if self.q is not None:
-            self.q.release()
-            self.q = None
+        if self.q_b is not None:
+            self.q_b.release()
+            self.q_b = None
 
     def __del__(self):
@@ -471,36 +471,35 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, model_info, q_vectors):
-        # type: (cl.Kernel, ModelInfo, List[np.ndarray]) -> None
+    def __init__(self, kernel, dtype, model_info, q_vectors):
+        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
         max_pd = model_info.parameters.max_pd
         npars = len(model_info.parameters.kernel_parameters)-2
-        q_input = GpuInput(q_vectors, kernel.dtype)
+        q_input = GpuInput(q_vectors, dtype)
         self.kernel = kernel
         self.info = model_info
-        self.dtype = kernel.dtype
+        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+3, q_input.dtype)
+        self.result = np.empty(q_input.nq+3, dtype)
 
         # Inputs and outputs for each kernel call
         # Note: res may be shorter than res_b if global_size != nq
         env = environment()
-        self.queue = env.get_queue(kernel.dtype)
+        self.queue = env.get_queue(dtype)
 
         # details is int32 data, padded to an 8 integer boundary
         size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8
         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                               q_input.global_size[0] * kernel.dtype.itemsize)
+                               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 self.q_input.dtype == generate.F32
-                     else np.float64 if self.q_input.dtype == generate.F64
-                     else np.float16 if self.q_input.dtype == generate.F16
+        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, weights, values, cutoff):
         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray
-
         context = self.queue.context
         # Arrange data transfer to card
@@ -508,5 +507,5 @@
                               hostbuf=call_details.buffer)
         weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
-                              hostbuf=weights)
+                              hostbuf=weights) if len(weights) else None
         values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
                              hostbuf=values)
@@ -521,5 +520,5 @@
         cl.enqueue_copy(self.queue, self.result, self.result_b)
         for v in (details_b, weights_b, values_b):
-            v.release()
+            if v is not None: v.release()
 
         return self.result[:self.q_input.nq]
