source: sasmodels/sasmodels/kernelcl.py @ 3199b17

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3199b17 was 3199b17, checked in by Paul Kienzle <pkienzle@…>, 4 months ago

PEP 8 changes and improved consistency between OpenCL/CUDA/DLL code

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