source: sasmodels/sasmodels/kernelcl.py @ 880a2ed

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 880a2ed was 880a2ed, checked in by wojciech, 8 years ago

Handling for SAS_OPENCL enviromental variable added

  • Property mode set to 100644
File size: 22.9 KB
Line 
1"""
2GPU driver for C kernels
3
4There should be a single GPU environment running on the system.  This
5environment is constructed on the first call to :func:`env`, and the
6same environment is returned on each call.
7
8After retrieving the environment, the next step is to create the kernel.
9This is done with a call to :meth:`GpuEnvironment.make_kernel`, which
10returns the type of data used by the kernel.
11
12Next a :class:`GpuData` object should be created with the correct kind
13of data.  This data object can be used by multiple kernels, for example,
14if the target model is a weighted sum of multiple kernels.  The data
15should include any extra evaluation points required to compute the proper
16data smearing.  This need not match the square grid for 2D data if there
17is an index saying which q points are active.
18
19Together the GpuData, the program, and a device form a :class:`GpuKernel`.
20This kernel is used during fitting, receiving new sets of parameters and
21evaluating them.  The output value is stored in an output buffer on the
22devices, where it can be combined with other structure factors and form
23factors and have instrumental resolution effects applied.
24
25In order to use OpenCL for your models, you will need OpenCL drivers for
26your machine.  These should be available from your graphics card vendor.
27Intel provides OpenCL drivers for CPUs as well as their integrated HD
28graphics chipsets.  AMD also provides drivers for Intel CPUs, but as of
29this writing the performance is lacking compared to the Intel drivers.
30NVidia combines drivers for CUDA and OpenCL in one package.  The result
31is a bit messy if you have multiple drivers installed.  You can see which
32drivers are available by starting python and running:
33
34    import pyopencl as cl
35    cl.create_some_context(interactive=True)
36
37Once you have done that, it will show the available drivers which you
38can select.  It will then tell you that you can use these drivers
39automatically by setting the SAS_OPENCL environment variable, which is
40PYOPENCL_CTX equivalent but not conflicting with other pyopnecl programs.
41
42Some graphics cards have multiple devices on the same card.  You cannot
43yet use both of them concurrently to evaluate models, but you can run
44the program twice using a different device for each session.
45
46OpenCL kernels are compiled when needed by the device driver.  Some
47drivers produce compiler output even when there is no error.  You
48can see the output by setting PYOPENCL_COMPILER_OUTPUT=1.  It should be
49harmless, albeit annoying.
50"""
51from __future__ import print_function
52
53import os
54import warnings
55import logging
56
57import numpy as np  # type: ignore
58
59try:
60    #raise NotImplementedError("OpenCL not yet implemented for new kernel template")
61    import pyopencl as cl  # type: ignore
62    # Ask OpenCL for the default context so that we know that one exists
63    cl.create_some_context(interactive=False)
64except Exception as exc:
65    warnings.warn("OpenCL startup failed with ***"
66                  + str(exc) + "***; using C compiler instead")
67    raise RuntimeError("OpenCL not available")
68
69from pyopencl import mem_flags as mf
70from pyopencl.characterize import get_fast_inaccurate_build_options
71
72from . import generate
73from .kernel import KernelModel, Kernel
74
75try:
76    from typing import Tuple, Callable, Any
77    from .modelinfo import ModelInfo
78    from .details import CallDetails
79except ImportError:
80    pass
81
82# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path)
83def quote_path(v):
84    """
85    Quote the path if it is not already quoted.
86
87    If v starts with '-', then assume that it is a -I option or similar
88    and do not quote it.  This is fragile:  -Ipath with space needs to
89    be quoted.
90    """
91    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v
92
93def fix_pyopencl_include():
94    """
95    Monkey patch pyopencl to allow spaces in include file path.
96    """
97    import pyopencl as cl
98    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'):
99        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS]
100
101fix_pyopencl_include()
102
103
104# The max loops number is limited by the amount of local memory available
105# on the device.  You don't want to make this value too big because it will
106# waste resources, nor too small because it may interfere with users trying
107# to do their polydispersity calculations.  A value of 1024 should be much
108# larger than necessary given that cost grows as npts^k where k is the number
109# of polydisperse parameters.
110MAX_LOOPS = 2048
111
112
113# Pragmas for enable OpenCL features.  Be sure to protect them so that they
114# still compile even if OpenCL is not present.
115_F16_PRAGMA = """\
116#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
117#  pragma OPENCL EXTENSION cl_khr_fp16: enable
118#endif
119"""
120
121_F64_PRAGMA = """\
122#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
123#  pragma OPENCL EXTENSION cl_khr_fp64: enable
124#endif
125"""
126
127
128ENV = None
129def environment():
130    # type: () -> "GpuEnvironment"
131    """
132    Returns a singleton :class:`GpuEnvironment`.
133
134    This provides an OpenCL context and one queue per device.
135    """
136    global ENV
137    if ENV is None:
138        ENV = GpuEnvironment()
139    return ENV
140
141def has_type(device, dtype):
142    # type: (cl.Device, np.dtype) -> bool
143    """
144    Return true if device supports the requested precision.
145    """
146    if dtype == generate.F32:
147        return True
148    elif dtype == generate.F64:
149        return "cl_khr_fp64" in device.extensions
150    elif dtype == generate.F16:
151        return "cl_khr_fp16" in device.extensions
152    else:
153        return False
154
155def get_warp(kernel, queue):
156    # type: (cl.Kernel, cl.CommandQueue) -> int
157    """
158    Return the size of an execution batch for *kernel* running on *queue*.
159    """
160    return kernel.get_work_group_info(
161        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
162        queue.device)
163
164def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
165    # type: (np.ndarray, np.dtype, float, int) -> np.ndarray
166    """
167    Stretch an input vector to the correct boundary.
168
169    Performance on the kernels can drop by a factor of two or more if the
170    number of values to compute does not fall on a nice power of two
171    boundary.   The trailing additional vector elements are given a
172    value of *extra*, and so f(*extra*) will be computed for each of
173    them.  The returned array will thus be a subset of the computed array.
174
175    *boundary* should be a power of 2 which is at least 32 for good
176    performance on current platforms (as of Jan 2015).  It should
177    probably be the max of get_warp(kernel,queue) and
178    device.min_data_type_align_size//4.
179    """
180    remainder = vector.size % boundary
181    if remainder != 0:
182        size = vector.size + (boundary - remainder)
183        vector = np.hstack((vector, [extra] * (size - vector.size)))
184    return np.ascontiguousarray(vector, dtype=dtype)
185
186
187def compile_model(context, source, dtype, fast=False):
188    # type: (cl.Context, str, np.dtype, bool) -> cl.Program
189    """
190    Build a model to run on the gpu.
191
192    Returns the compiled program and its type.  The returned type will
193    be float32 even if the desired type is float64 if any of the
194    devices in the context do not support the cl_khr_fp64 extension.
195    """
196    dtype = np.dtype(dtype)
197    if not all(has_type(d, dtype) for d in context.devices):
198        raise RuntimeError("%s not supported for devices"%dtype)
199
200    source_list = [generate.convert_type(source, dtype)]
201
202    if dtype == generate.F16:
203        source_list.insert(0, _F16_PRAGMA)
204    elif dtype == generate.F64:
205        source_list.insert(0, _F64_PRAGMA)
206
207    # Note: USE_SINCOS makes the intel cpu slower under opencl
208    if context.devices[0].type == cl.device_type.GPU:
209        source_list.insert(0, "#define USE_SINCOS\n")
210    options = (get_fast_inaccurate_build_options(context.devices[0])
211               if fast else [])
212    source = "\n".join(source_list)
213    program = cl.Program(context, source).build(options=options)
214    #print("done with "+program)
215    return program
216
217
218# for now, this returns one device in the context
219# TODO: create a context that contains all devices on all platforms
220class GpuEnvironment(object):
221    """
222    GPU context, with possibly many devices, and one queue per device.
223    """
224    def __init__(self):
225        # type: () -> None
226        # find gpu context
227        #self.context = cl.create_some_context()
228
229        self.context = None
230        if 'SAS_OPENCL' in os.environ:
231            #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
232            os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"]
233            self._create_some_context()
234
235        if not self.context:
236            self.context = _get_default_context()
237
238        # Byte boundary for data alignment
239        #self.data_boundary = max(d.min_data_type_align_size
240        #                         for d in self.context.devices)
241        self.queues = [cl.CommandQueue(context, context.devices[0])
242                       for context in self.context]
243        self.compiled = {}
244
245    def has_type(self, dtype):
246        # type: (np.dtype) -> bool
247        """
248        Return True if all devices support a given type.
249        """
250        return any(has_type(d, dtype)
251                   for context in self.context
252                   for d in context.devices)
253
254    def get_queue(self, dtype):
255        # type: (np.dtype) -> cl.CommandQueue
256        """
257        Return a command queue for the kernels of type dtype.
258        """
259        for context, queue in zip(self.context, self.queues):
260            if all(has_type(d, dtype) for d in context.devices):
261                return queue
262
263    def get_context(self, dtype):
264        # type: (np.dtype) -> cl.Context
265        """
266        Return a OpenCL context for the kernels of type dtype.
267        """
268        for context in self.context:
269            if all(has_type(d, dtype) for d in context.devices):
270                return context
271
272    def _create_some_context(self):
273        # type: () -> cl.Context
274        """
275        Protected call to cl.create_some_context without interactivity.  Use
276        this if SAS_OPENCL is set in the environment.  Sets the *context*
277        attribute.
278        """
279        try:
280            self.context = [cl.create_some_context(interactive=False)]
281        except Exception as exc:
282            warnings.warn(str(exc))
283            warnings.warn("pyopencl.create_some_context() failed")
284            warnings.warn("the environment variable 'SAS_OPENCL' might not be set correctly")
285
286    def compile_program(self, name, source, dtype, fast, timestamp):
287        # type: (str, str, np.dtype, bool, float) -> cl.Program
288        """
289        Compile the program for the device in the given context.
290        """
291        # Note: PyOpenCL caches based on md5 hash of source, options and device
292        # so we don't really need to cache things for ourselves.  I'll do so
293        # anyway just to save some data munging time.
294        key = "%s-%s%s"%(name, dtype, ("-fast" if fast else ""))
295        # Check timestamp on program
296        program, program_timestamp = self.compiled.get(key, (None, np.inf))
297        if program_timestamp < timestamp:
298            del self.compiled[key]
299        if key not in self.compiled:
300            context = self.get_context(dtype)
301            logging.info("building %s for OpenCL %s", key,
302                         context.devices[0].name.strip())
303            program = compile_model(self.get_context(dtype),
304                                    str(source), dtype, fast)
305            self.compiled[key] = (program, timestamp)
306        return program
307
308def _get_default_context():
309    # type: () -> List[cl.Context]
310    """
311    Get an OpenCL context, preferring GPU over CPU, and preferring Intel
312    drivers over AMD drivers.
313    """
314    # Note: on mobile devices there is automatic clock scaling if either the
315    # CPU or the GPU is underutilized; probably doesn't affect us, but we if
316    # it did, it would mean that putting a busy loop on the CPU while the GPU
317    # is running may increase throughput.
318    #
319    # Macbook pro, base install:
320    #     {'Apple': [Intel CPU, NVIDIA GPU]}
321    # Macbook pro, base install:
322    #     {'Apple': [Intel CPU, Intel GPU]}
323    # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed
324    #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]}
325    gpu, cpu = None, None
326    for platform in cl.get_platforms():
327        # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it.
328        # If someone has bothered to install the AMD/NVIDIA drivers, prefer
329        # them over the integrated graphics driver that may have been supplied
330        # with the CPU chipset.
331        preferred_cpu = (platform.vendor.startswith('Intel')
332                         or platform.vendor.startswith('Apple'))
333        preferred_gpu = (platform.vendor.startswith('Advanced')
334                         or platform.vendor.startswith('NVIDIA'))
335        for device in platform.get_devices():
336            if device.type == cl.device_type.GPU:
337                # If the existing type is not GPU then it will be CUSTOM
338                # or ACCELERATOR so don't override it.
339                if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU):
340                    gpu = device
341            elif device.type == cl.device_type.CPU:
342                if cpu is None or preferred_cpu:
343                    cpu = device
344            else:
345                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM
346                # Intel Phi for example registers as an accelerator
347                # Since the user installed a custom device on their system
348                # and went through the pain of sorting out OpenCL drivers for
349                # it, lets assume they really do want to use it as their
350                # primary compute device.
351                gpu = device
352
353    # order the devices by gpu then by cpu; when searching for an available
354    # device by data type they will be checked in this order, which means
355    # that if the gpu supports double then the cpu will never be used (though
356    # we may make it possible to explicitly request the cpu at some point).
357    devices = []
358    if gpu is not None:
359        devices.append(gpu)
360    if cpu is not None:
361        devices.append(cpu)
362    return [cl.Context([d]) for d in devices]
363
364
365class GpuModel(KernelModel):
366    """
367    GPU wrapper for a single model.
368
369    *source* and *model_info* are the model source and interface as returned
370    from :func:`generate.make_source` and :func:`generate.make_model_info`.
371
372    *dtype* is the desired model precision.  Any numpy dtype for single
373    or double precision floats will do, such as 'f', 'float32' or 'single'
374    for single and 'd', 'float64' or 'double' for double.  Double precision
375    is an optional extension which may not be available on all devices.
376    Half precision ('float16','half') may be available on some devices.
377    Fast precision ('fast') is a loose version of single precision, indicating
378    that the compiler is allowed to take shortcuts.
379    """
380    def __init__(self, source, model_info, dtype=generate.F32, fast=False):
381        # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None
382        self.info = model_info
383        self.source = source
384        self.dtype = dtype
385        self.fast = fast
386        self.program = None # delay program creation
387        self._kernels = None
388
389    def __getstate__(self):
390        # type: () -> Tuple[ModelInfo, str, np.dtype, bool]
391        return self.info, self.source, self.dtype, self.fast
392
393    def __setstate__(self, state):
394        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
395        self.info, self.source, self.dtype, self.fast = state
396        self.program = None
397
398    def make_kernel(self, q_vectors):
399        # type: (List[np.ndarray]) -> "GpuKernel"
400        if self.program is None:
401            compile_program = environment().compile_program
402            timestamp = generate.ocl_timestamp(self.info)
403            self.program = compile_program(
404                self.info.name,
405                self.source['opencl'],
406                self.dtype,
407                self.fast,
408                timestamp)
409            variants = ['Iq', 'Iqxy', 'Imagnetic']
410            names = [generate.kernel_name(self.info, k) for k in variants]
411            kernels = [getattr(self.program, k) for k in names]
412            self._kernels = dict((k, v) for k, v in zip(variants, kernels))
413        is_2d = len(q_vectors) == 2
414        if is_2d:
415            kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']]
416        else:
417            kernel = [self._kernels['Iq']]*2
418        return GpuKernel(kernel, self.dtype, self.info, q_vectors)
419
420    def release(self):
421        # type: () -> None
422        """
423        Free the resources associated with the model.
424        """
425        if self.program is not None:
426            self.program = None
427
428    def __del__(self):
429        # type: () -> None
430        self.release()
431
432# TODO: check that we don't need a destructor for buffers which go out of scope
433class GpuInput(object):
434    """
435    Make q data available to the gpu.
436
437    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
438    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
439    to get the best performance on OpenCL, which may involve shifting and
440    stretching the array to better match the memory architecture.  Additional
441    points will be evaluated with *q=1e-3*.
442
443    *dtype* is the data type for the q vectors. The data type should be
444    set to match that of the kernel, which is an attribute of
445    :class:`GpuProgram`.  Note that not all kernels support double
446    precision, so even if the program was created for double precision,
447    the *GpuProgram.dtype* may be single precision.
448
449    Call :meth:`release` when complete.  Even if not called directly, the
450    buffer will be released when the data object is freed.
451    """
452    def __init__(self, q_vectors, dtype=generate.F32):
453        # type: (List[np.ndarray], np.dtype) -> None
454        # TODO: do we ever need double precision q?
455        env = environment()
456        self.nq = q_vectors[0].size
457        self.dtype = np.dtype(dtype)
458        self.is_2d = (len(q_vectors) == 2)
459        # TODO: stretch input based on get_warp()
460        # not doing it now since warp depends on kernel, which is not known
461        # at this point, so instead using 32, which is good on the set of
462        # architectures tested so far.
463        if self.is_2d:
464            # Note: 17 rather than 15 because results is 2 elements
465            # longer than input.
466            width = ((self.nq+17)//16)*16
467            self.q = np.empty((width, 2), dtype=dtype)
468            self.q[:self.nq, 0] = q_vectors[0]
469            self.q[:self.nq, 1] = q_vectors[1]
470        else:
471            # Note: 33 rather than 31 because results is 2 elements
472            # longer than input.
473            width = ((self.nq+33)//32)*32
474            self.q = np.empty(width, dtype=dtype)
475            self.q[:self.nq] = q_vectors[0]
476        self.global_size = [self.q.shape[0]]
477        context = env.get_context(self.dtype)
478        #print("creating inputs of size", self.global_size)
479        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
480                             hostbuf=self.q)
481
482    def release(self):
483        # type: () -> None
484        """
485        Free the memory.
486        """
487        if self.q_b is not None:
488            self.q_b.release()
489            self.q_b = None
490
491    def __del__(self):
492        # type: () -> None
493        self.release()
494
495class GpuKernel(Kernel):
496    """
497    Callable SAS kernel.
498
499    *kernel* is the GpuKernel object to call
500
501    *model_info* is the module information
502
503    *q_vectors* is the q vectors at which the kernel should be evaluated
504
505    *dtype* is the kernel precision
506
507    The resulting call method takes the *pars*, a list of values for
508    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
509    vectors for the polydisperse parameters.  *cutoff* determines the
510    integration limits: any points with combined weight less than *cutoff*
511    will not be calculated.
512
513    Call :meth:`release` when done with the kernel instance.
514    """
515    def __init__(self, kernel, dtype, model_info, q_vectors):
516        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
517        q_input = GpuInput(q_vectors, dtype)
518        self.kernel = kernel
519        self.info = model_info
520        self.dtype = dtype
521        self.dim = '2d' if q_input.is_2d else '1d'
522        # plus three for the normalization values
523        self.result = np.empty(q_input.nq+3, dtype)
524
525        # Inputs and outputs for each kernel call
526        # Note: res may be shorter than res_b if global_size != nq
527        env = environment()
528        self.queue = env.get_queue(dtype)
529
530        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
531                                  q_input.global_size[0] * dtype.itemsize)
532        self.q_input = q_input # allocated by GpuInput above
533
534        self._need_release = [self.result_b, self.q_input]
535        self.real = (np.float32 if dtype == generate.F32
536                     else np.float64 if dtype == generate.F64
537                     else np.float16 if dtype == generate.F16
538                     else np.float32)  # will never get here, so use np.float32
539
540    def __call__(self, call_details, values, cutoff, magnetic):
541        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
542        context = self.queue.context
543        # Arrange data transfer to card
544        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
545                              hostbuf=call_details.buffer)
546        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
547                             hostbuf=values)
548
549        kernel = self.kernel[1 if magnetic else 0]
550        args = [
551            np.uint32(self.q_input.nq), None, None,
552            details_b, values_b, self.q_input.q_b, self.result_b,
553            self.real(cutoff),
554        ]
555        #print("Calling OpenCL")
556        #call_details.show(values)
557        # Call kernel and retrieve results
558        last_call = None
559        step = 100
560        for start in range(0, call_details.num_eval, step):
561            stop = min(start + step, call_details.num_eval)
562            #print("queuing",start,stop)
563            args[1:3] = [np.int32(start), np.int32(stop)]
564            last_call = [kernel(self.queue, self.q_input.global_size,
565                                None, *args, wait_for=last_call)]
566        cl.enqueue_copy(self.queue, self.result, self.result_b)
567        #print("result", self.result)
568
569        # Free buffers
570        for v in (details_b, values_b):
571            if v is not None: v.release()
572
573        pd_norm = self.result[self.q_input.nq]
574        scale = values[0]/(pd_norm if pd_norm!=0.0 else 1.0)
575        background = values[1]
576        #print("scale",scale,values[0],self.result[self.q_input.nq],background)
577        return scale*self.result[:self.q_input.nq] + background
578        # return self.result[:self.q_input.nq]
579
580    def release(self):
581        # type: () -> None
582        """
583        Release resources associated with the kernel.
584        """
585        for v in self._need_release:
586            v.release()
587        self._need_release = []
588
589    def __del__(self):
590        # type: () -> None
591        self.release()
Note: See TracBrowser for help on using the repository browser.