source: sasmodels/sasmodels/kernelcl.py @ f872fd1

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

simplify F32/F64 handling

  • Property mode set to 100644
File size: 24.6 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# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path)
91def quote_path(v):
92    """
93    Quote the path if it is not already quoted.
94
95    If v starts with '-', then assume that it is a -I option or similar
96    and do not quote it.  This is fragile:  -Ipath with space needs to
97    be quoted.
98    """
99    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v
100
101def fix_pyopencl_include():
102    """
103    Monkey patch pyopencl to allow spaces in include file path.
104    """
105    import pyopencl as cl
106    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'):
107        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS]
108
109if HAVE_OPENCL:
110    fix_pyopencl_include()
111
112# The max loops number is limited by the amount of local memory available
113# on the device.  You don't want to make this value too big because it will
114# waste resources, nor too small because it may interfere with users trying
115# to do their polydispersity calculations.  A value of 1024 should be much
116# larger than necessary given that cost grows as npts^k where k is the number
117# of polydisperse parameters.
118MAX_LOOPS = 2048
119
120
121# Pragmas for enable OpenCL features.  Be sure to protect them so that they
122# still compile even if OpenCL is not present.
123_F16_PRAGMA = """\
124#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
125#  pragma OPENCL EXTENSION cl_khr_fp16: enable
126#endif
127"""
128
129_F64_PRAGMA = """\
130#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
131#  pragma OPENCL EXTENSION cl_khr_fp64: enable
132#endif
133"""
134
135def use_opencl():
136    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower()
137    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda")
138
139ENV = None
140def reset_environment():
141    """
142    Call to create a new OpenCL context, such as after a change to SAS_OPENCL.
143    """
144    global ENV
145    ENV = GpuEnvironment() if use_opencl() else None
146
147def environment():
148    # type: () -> "GpuEnvironment"
149    """
150    Returns a singleton :class:`GpuEnvironment`.
151
152    This provides an OpenCL context and one queue per device.
153    """
154    if ENV is None:
155        if not HAVE_OPENCL:
156            raise RuntimeError("OpenCL startup failed with ***"
157                               + OPENCL_ERROR + "***; using C compiler instead")
158        reset_environment()
159        if ENV is None:
160            raise RuntimeError("SAS_OPENCL=None in environment")
161    return ENV
162
163def has_type(device, dtype):
164    # type: (cl.Device, np.dtype) -> bool
165    """
166    Return true if device supports the requested precision.
167    """
168    if dtype == F32:
169        return True
170    elif dtype == F64:
171        return "cl_khr_fp64" in device.extensions
172    else:
173        # Not supporting F16 type since it isn't accurate enough
174        return False
175
176def get_warp(kernel, queue):
177    # type: (cl.Kernel, cl.CommandQueue) -> int
178    """
179    Return the size of an execution batch for *kernel* running on *queue*.
180    """
181    return kernel.get_work_group_info(
182        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
183        queue.device)
184
185def compile_model(context, source, dtype, fast=False):
186    # type: (cl.Context, str, np.dtype, bool) -> cl.Program
187    """
188    Build a model to run on the gpu.
189
190    Returns the compiled program and its type.
191
192    Raises an error if the desired precision is not available.
193    """
194    dtype = np.dtype(dtype)
195    if not all(has_type(d, dtype) for d in context.devices):
196        raise RuntimeError("%s not supported for devices"%dtype)
197
198    source_list = [generate.convert_type(source, dtype)]
199
200    if dtype == generate.F16:
201        source_list.insert(0, _F16_PRAGMA)
202    elif dtype == generate.F64:
203        source_list.insert(0, _F64_PRAGMA)
204
205    # Note: USE_SINCOS makes the intel cpu slower under opencl
206    if context.devices[0].type == cl.device_type.GPU:
207        source_list.insert(0, "#define USE_SINCOS\n")
208    options = (get_fast_inaccurate_build_options(context.devices[0])
209               if fast else [])
210    source = "\n".join(source_list)
211    program = cl.Program(context, source).build(options=options)
212    #print("done with "+program)
213    return program
214
215
216# for now, this returns one device in the context
217# TODO: create a context that contains all devices on all platforms
218class GpuEnvironment(object):
219    """
220    GPU context, with possibly many devices, and one queue per device.
221
222    Because the environment can be reset during a live program (e.g., if the
223    user changes the active GPU device in the GUI), everything associated
224    with the device context must be cached in the environment and recreated
225    if the environment changes.  The *cache* attribute is a simple dictionary
226    which holds keys and references to objects, such as compiled kernels and
227    allocated buffers.  The running program should check in the cache for
228    long lived objects and create them if they are not there.  The program
229    should not hold onto cached objects, but instead only keep them active
230    for the duration of a function call.  When the environment is destroyed
231    then the *release* method for each active cache item is called before
232    the environment is freed.  This means that each cl buffer should be
233    in its own cache entry.
234    """
235    def __init__(self):
236        # type: () -> None
237        # find gpu context
238        context_list = _create_some_context()
239
240        # Find a context for F32 and for F64 (maybe the same one).
241        # F16 isn't good enough.
242        self.context = {}
243        for dtype in (F32, F64):
244            for context in context_list:
245                if has_type(context.devices[0], dtype):
246                    self.context[dtype] = context
247                    break
248            else:
249                self.context[dtype] = None
250
251        # Build a queue for each context
252        self.queue = {}
253        context = self.context[F32]
254        self.queue[F32] = cl.CommandQueue(context, context.devices[0])
255        if self.context[F64] == self.context[F32]:
256            self.queue[F64] = self.queue[F32]
257        else:
258            context = self.context[F64]
259            self.queue[F64] = cl.CommandQueue(context, context.devices[0])
260
261        # Byte boundary for data alignment
262        #self.data_boundary = max(context.devices[0].min_data_type_align_size
263        #                         for context in self.context.values())
264
265        # Cache for compiled programs, and for items in context
266        self.compiled = {}
267        self.cache = {}
268
269    def has_type(self, dtype):
270        # type: (np.dtype) -> bool
271        """
272        Return True if all devices support a given type.
273        """
274        return self.context.get(dtype, None) is not None
275
276    def compile_program(self, name, source, dtype, fast, timestamp):
277        # type: (str, str, np.dtype, bool, float) -> cl.Program
278        """
279        Compile the program for the device in the given context.
280        """
281        # Note: PyOpenCL caches based on md5 hash of source, options and device
282        # so we don't really need to cache things for ourselves.  I'll do so
283        # anyway just to save some data munging time.
284        tag = generate.tag_source(source)
285        key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else ""))
286        # Check timestamp on program
287        program, program_timestamp = self.compiled.get(key, (None, np.inf))
288        if program_timestamp < timestamp:
289            del self.compiled[key]
290        if key not in self.compiled:
291            context = self.context[dtype]
292            logging.info("building %s for OpenCL %s", key,
293                         context.devices[0].name.strip())
294            program = compile_model(self.context[dtype],
295                                    str(source), dtype, fast)
296            self.compiled[key] = (program, timestamp)
297        return program
298
299    def free_buffer(self, key):
300        if key in self.cache:
301            self.cache[key].release()
302            del self.cache[key]
303
304    def __del__(self):
305        for v in self.cache.values():
306            release = getattr(v, 'release', lambda: None)
307            release()
308        self.cache = {}
309
310_CURRENT_ID = 0
311def unique_id():
312    global _CURRENT_ID
313    _CURRENT_ID += 1
314    return _CURRENT_ID
315
316def _create_some_context():
317    # type: () -> cl.Context
318    """
319    Protected call to cl.create_some_context without interactivity.
320
321    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment,
322    otherwise scans for the most appropriate device using
323    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which
324    indicates that an OpenCL device should be used without specifying
325    which one (and not a CUDA device, or no GPU).
326    """
327    # Assume we do not get here if SAS_OPENCL is None or CUDA
328    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl')
329    if sas_opencl.lower() != 'opencl':
330        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
331        os.environ["PYOPENCL_CTX"] = sas_opencl
332
333    if 'PYOPENCL_CTX' in os.environ:
334        try:
335            return [cl.create_some_context(interactive=False)]
336        except Exception as exc:
337            warnings.warn(str(exc))
338            warnings.warn("pyopencl.create_some_context() failed")
339            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly")
340
341    return _get_default_context()
342
343def _get_default_context():
344    # type: () -> List[cl.Context]
345    """
346    Get an OpenCL context, preferring GPU over CPU, and preferring Intel
347    drivers over AMD drivers.
348    """
349    # Note: on mobile devices there is automatic clock scaling if either the
350    # CPU or the GPU is underutilized; probably doesn't affect us, but we if
351    # it did, it would mean that putting a busy loop on the CPU while the GPU
352    # is running may increase throughput.
353    #
354    # Macbook pro, base install:
355    #     {'Apple': [Intel CPU, NVIDIA GPU]}
356    # Macbook pro, base install:
357    #     {'Apple': [Intel CPU, Intel GPU]}
358    # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed
359    #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]}
360    gpu, cpu = None, None
361    for platform in cl.get_platforms():
362        # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it.
363        # If someone has bothered to install the AMD/NVIDIA drivers, prefer
364        # them over the integrated graphics driver that may have been supplied
365        # with the CPU chipset.
366        preferred_cpu = (platform.vendor.startswith('Intel')
367                         or platform.vendor.startswith('Apple'))
368        preferred_gpu = (platform.vendor.startswith('Advanced')
369                         or platform.vendor.startswith('NVIDIA'))
370        for device in platform.get_devices():
371            if device.type == cl.device_type.GPU:
372                # If the existing type is not GPU then it will be CUSTOM
373                # or ACCELERATOR so don't override it.
374                if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU):
375                    gpu = device
376            elif device.type == cl.device_type.CPU:
377                if cpu is None or preferred_cpu:
378                    cpu = device
379            else:
380                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM
381                # Intel Phi for example registers as an accelerator
382                # Since the user installed a custom device on their system
383                # and went through the pain of sorting out OpenCL drivers for
384                # it, lets assume they really do want to use it as their
385                # primary compute device.
386                gpu = device
387
388    # order the devices by gpu then by cpu; when searching for an available
389    # device by data type they will be checked in this order, which means
390    # that if the gpu supports double then the cpu will never be used (though
391    # we may make it possible to explicitly request the cpu at some point).
392    devices = []
393    if gpu is not None:
394        devices.append(gpu)
395    if cpu is not None:
396        devices.append(cpu)
397    return [cl.Context([d]) for d in devices]
398
399
400class GpuModel(KernelModel):
401    """
402    GPU wrapper for a single model.
403
404    *source* and *model_info* are the model source and interface as returned
405    from :func:`generate.make_source` and :func:`generate.make_model_info`.
406
407    *dtype* is the desired model precision.  Any numpy dtype for single
408    or double precision floats will do, such as 'f', 'float32' or 'single'
409    for single and 'd', 'float64' or 'double' for double.  Double precision
410    is an optional extension which may not be available on all devices.
411    Half precision ('float16','half') may be available on some devices.
412    Fast precision ('fast') is a loose version of single precision, indicating
413    that the compiler is allowed to take shortcuts.
414    """
415    def __init__(self, source, model_info, dtype=generate.F32, fast=False):
416        # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None
417        self.info = model_info
418        self.source = source
419        self.dtype = dtype
420        self.fast = fast
421        self.timestamp = generate.ocl_timestamp(self.info)
422        self._cache_key = unique_id()
423
424    def __getstate__(self):
425        # type: () -> Tuple[ModelInfo, str, np.dtype, bool]
426        return self.info, self.source, self.dtype, self.fast
427
428    def __setstate__(self, state):
429        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
430        self.info, self.source, self.dtype, self.fast = state
431
432    def make_kernel(self, q_vectors):
433        # type: (List[np.ndarray]) -> "GpuKernel"
434        return GpuKernel(self, q_vectors)
435
436    @property
437    def Iq(self):
438        return self._fetch_kernel('Iq')
439
440    def fetch_kernel(self, name):
441        # type: (str) -> cl.Kernel
442        """
443        Fetch the kernel from the environment by name, compiling it if it
444        does not already exist.
445        """
446        gpu = environment()
447        key = self._cache_key
448        if key not in gpu.cache:
449            program = gpu.compile_program(
450                self.info.name,
451                self.source['opencl'],
452                self.dtype,
453                self.fast,
454                self.timestamp)
455            variants = ['Iq', 'Iqxy', 'Imagnetic']
456            names = [generate.kernel_name(self.info, k) for k in variants]
457            kernels = [getattr(program, k) for k in names]
458            data = dict((k, v) for k, v in zip(variants, kernels))
459            # keep a handle to program so GC doesn't collect
460            data['program'] = program
461            gpu.cache[key] = data
462        else:
463            data = gpu.cache[key]
464        return data[name]
465
466# TODO: check that we don't need a destructor for buffers which go out of scope
467class GpuInput(object):
468    """
469    Make q data available to the gpu.
470
471    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
472    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
473    to get the best performance on OpenCL, which may involve shifting and
474    stretching the array to better match the memory architecture.  Additional
475    points will be evaluated with *q=1e-3*.
476
477    *dtype* is the data type for the q vectors. The data type should be
478    set to match that of the kernel, which is an attribute of
479    :class:`GpuProgram`.  Note that not all kernels support double
480    precision, so even if the program was created for double precision,
481    the *GpuProgram.dtype* may be single precision.
482
483    Call :meth:`release` when complete.  Even if not called directly, the
484    buffer will be released when the data object is freed.
485    """
486    def __init__(self, q_vectors, dtype=generate.F32):
487        # type: (List[np.ndarray], np.dtype) -> None
488        # TODO: do we ever need double precision q?
489        self.nq = q_vectors[0].size
490        self.dtype = np.dtype(dtype)
491        self.is_2d = (len(q_vectors) == 2)
492        # TODO: stretch input based on get_warp()
493        # not doing it now since warp depends on kernel, which is not known
494        # at this point, so instead using 32, which is good on the set of
495        # architectures tested so far.
496        if self.is_2d:
497            width = ((self.nq+15)//16)*16
498            self.q = np.empty((width, 2), dtype=dtype)
499            self.q[:self.nq, 0] = q_vectors[0]
500            self.q[:self.nq, 1] = q_vectors[1]
501        else:
502            width = ((self.nq+31)//32)*32
503            self.q = np.empty(width, dtype=dtype)
504            self.q[:self.nq] = q_vectors[0]
505        self.global_size = [self.q.shape[0]]
506        self._cache_key = unique_id()
507
508    @property
509    def q_b(self):
510        """Lazy creation of q buffer so it can survive context reset"""
511        env = environment()
512        key = self._cache_key
513        if key not in env.cache:
514            context = env.context[self.dtype]
515            #print("creating inputs of size", self.global_size)
516            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
517                               hostbuf=self.q)
518            env.cache[key] = buffer
519        return env.cache[key]
520
521    def release(self):
522        # type: () -> None
523        """
524        Free the buffer associated with the q value
525        """
526        environment().free_buffer(id(self))
527
528    def __del__(self):
529        # type: () -> None
530        self.release()
531
532class GpuKernel(Kernel):
533    """
534    Callable SAS kernel.
535
536    *model* is the GpuModel object to call
537
538    The following attributes are defined:
539
540    *info* is the module information
541
542    *dtype* is the kernel precision
543
544    *dim* is '1d' or '2d'
545
546    *result* is a vector to contain the results of the call
547
548    The resulting call method takes the *pars*, a list of values for
549    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
550    vectors for the polydisperse parameters.  *cutoff* determines the
551    integration limits: any points with combined weight less than *cutoff*
552    will not be calculated.
553
554    Call :meth:`release` when done with the kernel instance.
555    """
556    def __init__(self, model, q_vectors):
557        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
558        dtype = model.dtype
559        self.q_input = GpuInput(q_vectors, dtype)
560        self._model = model
561        # F16 isn't sufficient, so don't support it
562        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32
563        self._cache_key = unique_id()
564
565        # attributes accessed from the outside
566        self.dim = '2d' if self.q_input.is_2d else '1d'
567        self.info = model.info
568        self.dtype = model.dtype
569
570        # holding place for the returned value
571        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
572        extra_q = 4  # total weight, form volume, shell volume and R_eff
573        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype)
574
575    @property
576    def _result_b(self):
577        """Lazy creation of result buffer so it can survive context reset"""
578        env = environment()
579        key = self._cache_key
580        if key not in env.cache:
581            context = env.context[self.dtype]
582            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize
583            buffer = cl.Buffer(context, mf.READ_WRITE, width)
584            env.cache[key] = buffer
585        return env.cache[key]
586
587    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
588        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
589        env = environment()
590        queue = env.queue[self._model.dtype]
591        context = queue.context
592
593        # Arrange data transfer to/from card
594        q_b = self.q_input.q_b
595        result_b = self._result_b
596        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
597                              hostbuf=call_details.buffer)
598        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
599                             hostbuf=values)
600
601        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy'
602        kernel = self._model.fetch_kernel(name)
603        kernel_args = [
604            np.uint32(self.q_input.nq), None, None,
605            details_b, values_b, q_b, result_b,
606            self._as_dtype(cutoff),
607            np.uint32(effective_radius_type),
608        ]
609        #print("Calling OpenCL")
610        #call_details.show(values)
611        #Call kernel and retrieve results
612        wait_for = None
613        last_nap = time.clock()
614        step = 1000000//self.q_input.nq + 1
615        for start in range(0, call_details.num_eval, step):
616            stop = min(start + step, call_details.num_eval)
617            #print("queuing",start,stop)
618            kernel_args[1:3] = [np.int32(start), np.int32(stop)]
619            wait_for = [kernel(queue, self.q_input.global_size, None,
620                               *kernel_args, wait_for=wait_for)]
621            if stop < call_details.num_eval:
622                # Allow other processes to run
623                wait_for[0].wait()
624                current_time = time.clock()
625                if current_time - last_nap > 0.5:
626                    time.sleep(0.001)
627                    last_nap = current_time
628        cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for)
629        #print("result", self.result)
630
631        # Free buffers
632        for v in (details_b, values_b):
633            if v is not None:
634                v.release()
635
636    def release(self):
637        # type: () -> None
638        """
639        Release resources associated with the kernel.
640        """
641        environment().free_buffer(id(self))
642        self.q_input.release()
643
644    def __del__(self):
645        # type: () -> None
646        self.release()
Note: See TracBrowser for help on using the repository browser.