Changes in sasmodels/kernelcl.py [da32ec3:3c56da87] in sasmodels

File:

• sasmodels/kernelcl.py (modified) (15 diffs)

sasmodels/kernelcl.py

@@ -30 +30 @@
 try:
     import pyopencl as cl
-except ImportError,exc:
+    # Ask OpenCL for the default context so that we know that one exists
+    cl.create_some_context(interactive=False)
+except Exception, exc:
     warnings.warn(str(exc))
     raise RuntimeError("OpenCL not available")
@@ -37 +39 @@
 
 from . import generate
-from .kernelpy import PyInput, PyKernel
+from .kernelpy import PyModel
 
 F64_DEFS = """\
@@ -61 +63 @@
     """
     source, info = generate.make(kernel_module)
+    if callable(info.get('Iq', None)):
+        return PyModel(info)
     ## for debugging, save source to a .cl file, edit it, and reload as model
     #open(info['name']+'.cl','w').write(source)
@@ -106 +110 @@
     device.min_data_type_align_size//4.
     """
-    remainder = vector.size%boundary
+    remainder = vector.size % boundary
     if remainder != 0:
         size = vector.size + (boundary - remainder)
-        vector = np.hstack((vector, [extra]*(size-vector.size)))
+        vector = np.hstack((vector, [extra] * (size - vector.size)))
     return np.ascontiguousarray(vector, dtype=dtype)
 
@@ -122 +126 @@
     """
     dtype = np.dtype(dtype)
-    if dtype==generate.F64 and not all(has_double(d) for d in context.devices):
+    if dtype == generate.F64 and not all(has_double(d) for d in context.devices):
         raise RuntimeError("Double precision not supported for devices")
 
@@ -131 +135 @@
     if context.devices[0].type == cl.device_type.GPU:
         header += "#define USE_SINCOS\n"
-    program  = cl.Program(context, header+source).build()
+    program = cl.Program(context, header + source).build()
     return program
 
@@ -156 +160 @@
 
         if not self.context:
-            self.context = self._find_context()
+            self.context = _get_default_context()
 
         # Byte boundary for data alignment
@@ -169 +173 @@
         try:
             self.context = cl.create_some_context(interactive=False)
-        except Exception,exc:
+        except Exception, exc:
             warnings.warn(str(exc))
             warnings.warn("pyopencl.create_some_context() failed")
             warnings.warn("the environment variable 'PYOPENCL_CTX' might not be set correctly")
-
-    def _find_context(self):
-        default = None
-        for platform in cl.get_platforms():
-            for device in platform.get_devices():
-                if device.type == cl.device_type.GPU:
-                    return cl.Context([device])
-                if default is None:
-                    default = device
-
-        if not default:
-            raise RuntimeError("OpenCL device not found")
-
-        return cl.Context([default])
 
     def compile_program(self, name, source, dtype):
@@ -199 +189 @@
             del self.compiled[name]
 
+def _get_default_context():
+    default = None
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            if device.type == cl.device_type.GPU:
+                return cl.Context([device])
+            if default is None:
+                default = device
+
+    if not default:
+        raise RuntimeError("OpenCL device not found")
+
+    return cl.Context([default])
+
 
 class GpuModel(object):
@@ -226 +230 @@
         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:
+    def __call__(self, input_value):
+        if self.dtype != input_value.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)
+            compiler = environment().compile_program
+            self.program = compiler(self.info['name'], self.source, self.dtype)
+        kernel_name = generate.kernel_name(self.info, input_value.is_2D)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, input)
+        return GpuKernel(kernel, self.info, input_value)
 
     def release(self):
@@ -254 +253 @@
         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)
+        return GpuInput(q_vectors, dtype=self.dtype)
 
 # TODO: check that we don't need a destructor for buffers which go out of scope
@@ -293 +286 @@
         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)
+            cl.Buffer(env.context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
             for q in self.q_vectors
         ]
@@ -311 +304 @@
     *info* is the module information
 
-    *input* is the DllInput q vectors at which the kernel should be
+    *q_input* is the DllInput q vectors at which the kernel should be
     evaluated.
 
@@ -322 +315 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, input):
-        self.input = input
+    def __init__(self, kernel, info, q_input):
+        self.q_input = q_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]
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_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
@@ -335 +328 @@
         env = environment()
         self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  2*MAX_LOOPS*input.dtype.itemsize)
+                                  2 * MAX_LOOPS * q_input.dtype.itemsize)
                         for _ in env.queues]
         self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                input.global_size[0]*input.dtype.itemsize)
+                                q_input.global_size[0] * q_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")
+    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
+        real = np.float32 if self.q_input.dtype == generate.F32 else np.float64
+
         device_num = 0
+        queuei = environment().queues[device_num]
         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)
+        nq = np.uint32(self.q_input.nq)
+        if pd_pars:
+            cutoff = real(cutoff)
+            loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+            loops = np.hstack(pd_pars) \
+                if pd_pars else np.empty(0, dtype=self.q_input.dtype)
+            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
+            #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")
+
+            loops_bi = self.loops_b[device_num]
+            cl.enqueue_copy(queuei, loops_bi, loops)
+            loops_l = cl.LocalMemory(len(loops.data))
+            #ctx = environment().context
+            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops)
+            dispersed = [loops_bi, loops_l, cutoff] + loops_N
+        else:
+            dispersed = []
+        fixed = [real(p) for p in fixed_pars]
+        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
+        self.kernel(queuei, self.q_input.global_size, None, *args)
         cl.enqueue_copy(queuei, self.res, res_bi)